Pi está errado

Um método comum de calcular pi é jogar "dardos" em uma caixa 1x1 e ver qual área do círculo de unidades é comparada ao total lançado:

loop
   x = rand()
   y = rand()
   if(sqrt(x*x + y*y) <= 1) n++
   t++
pi = 4.0*(n/t)

Escreva um programa que pareça que ele deveria calcular corretamente pi (usando este ou outros métodos comuns de computação de pi), mas calcula tau (tau = 2 * pi = 6.283185307179586 ...). Seu código deve produzir pelo menos as 6 primeiras casas decimais: 6.283185

O vencedor é coroado em 6 de junho (uma semana a partir de hoje).

Javascript

alert(Math.atan2(0, -0) - Math.atan2(-0, -0) + Math.atan2(0, 0))

Ajuda, estou preso em uma fábrica do universo e não tenho certeza do que estou fazendo. Math.atan2deve retornar pi com bons valores, certo? Math.atan2(0, -0)retorna pi, então, se eu subtraí-lo e adicioná-lo, ainda devo ter pi.

BASIC

(Mais especificamente, Chipmunk Basic )

Isso usa uma série infinita descoberta por Nilakantha Somayaji no século XV:

' Calculate pi using the Nilakantha series:
'               4       4       4       4
'  pi  =  3 + ----- - ----- + ----- - ------ + ...
'             2x3x4   4x5x6   6x7x8   8x9x10
i = pi = 0
numerator = -4
while i<10000
  i = i + 2
  numerator = -numerator
  pi = pi + numerator / (i * (i+1) * (i+2))
wend
pi = pi + 3
print using "#.##########";pi

Saída

6.2831853072

Se você não consegue descobrir o que está acontecendo, aqui estão algumas dicas:

No Chipmunk Basic, a variável pi é predefinida para o valor de π quando o programa começa a ser executado.

e

No BASIC, o sinal de igual é usado para atribuir variáveis ​​e para testar a igualdade. Portanto, a = b = c é interpretado como a = (b == c) .

C - O comprimento de meio círculo unitário

Uma maneira de calcular π é simplesmente medir a distância que o ponto (1, 0)viaja ao girar em torno da origem, uma (-1, 0)vez que será metade da circunferência de um círculo unitário (que é 2π ).

No entanto, nenhum sin(x)ou cos(x)é necessária uma vez que isso pode ser feito por pisar todo o caminho em torno da origem e adicionando a distância do ponto viaja para cada etapa . Quanto menor o tamanho de cada etapa, mais preciso será o π .

Nota: O passo terminará quando y estiver abaixo de zero (que é exatamente quando passa (-1, 0)).

#include <stdio.h>                          // for printf
#define length(y, x) ((x * x) + (y * y))
int main()
{
    double x, y;
    double pi, tau, step;
    // start at (2, 0) which actually calculates tau
    x  = 2;
    y  = 0;
    // the step needs to be very low for high accuracy
    step = 0.00000001;  
    tau = 0;
    while (y >= 0)
    {   // the derivate of (x, y) is itself rotated 90 degrees
        double dx = -y;
        double dy = x;

        tau += length(dx, dy) * step; // add the distance for each step to tau
        // add the distance to the point (make a tiny rotation)
        x += dx * step;
        y += dy * step;
    }
    pi = tau / 2;   // divide tau with 2 to get pi

    /* ignore this line *\                      pi *= 2;    /* secret multiply ^-^ */

    // print the value of pi
    printf("Value of pi is %f", pi); getchar(); 
    return 0;
}

Ele fornece a seguinte saída:

Value of pi is 6.283185

C

(Isso acabou sendo mais longo do que o pretendido, mas eu o publicarei de qualquer maneira ...)

No século XVII, Wallis publicou uma série infinita para Pi:

(Consulte Novos produtos infinitos do tipo Wallis e catalão para π, e e √ (2 + √2) para obter mais informações)

Agora, para calcular Pi, primeiro devemos multiplicar por dois para fatorar o denominador:

Minha solução calcula as séries infinitas para Pi / 2 e dois e multiplica os dois valores juntos. Observe que produtos infinitos são incrivelmente lentos para convergir ao calcular valores finais.

saída:

pi: 6.283182

#include "stdio.h"
#include "stdint.h"

#define ITERATIONS 10000000
#define one 1

#define IEEE_MANTISSA_MASK 0xFFFFFFFFFFFFFULL

#define IEEE_EXPONENT_POSITION 52
#define IEEE_EXPONENT_BIAS 1023

// want to get an exact as possible result, so convert
// to integers and do custom 64-bit multiplication.
double multiply(double aa, double bb)
{
    // the input values will be between 1.0 and 2.0
    // so drop these to less than 1.0 so as not to deal 
    // with the double exponents.
    aa /= 2;
    bb /= 2;

    // extract fractional part of double, ignoring exponent and sign
    uint64_t a = *(uint64_t*)&aa & IEEE_MANTISSA_MASK;
    uint64_t b = *(uint64_t*)&bb & IEEE_MANTISSA_MASK;

    uint64_t result = 0x0ULL;

    // multiplying two 64-bit numbers is a little tricky, this is done in two parts,
    // taking the upper 32 bits of each number and multiplying them, then
    // then doing the same for the lower 32 bits.
    uint64_t a_lsb = (a & 0xFFFFFFFFUL);
    uint64_t b_lsb = (b & 0xFFFFFFFFUL);

    uint64_t a_msb = ((a >> 32) & 0xFFFFFFFFUL);
    uint64_t b_msb = ((b >> 32) & 0xFFFFFFFFUL);

    uint64_t lsb_result = 0;
    uint64_t msb_result = 0;

    // very helpful link explaining how to multiply two integers
    // http://stackoverflow.com/questions/4456442/interview-multiplication-of-2-integers-using-bitwise-operators
    while(b_lsb != 0)
    {
        if (b_lsb & 01)
        {
            lsb_result = lsb_result + a_lsb;
        }
        a_lsb <<= 1;
        b_lsb >>= 1;
    }
    while(b_msb != 0)
    {
        if (b_msb & 01)
        {
            msb_result = msb_result + a_msb;
        }
        a_msb <<= 1;
        b_msb >>= 1;
    }

    // find the bit position of the most significant bit in the higher 32-bit product (msb_answer)
    uint64_t x2 = msb_result;
    int bit_pos = 0;
    while (x2 >>= 1)
    {
        bit_pos++;
    }

    // stuff bits from the upper 32-bit product into the result, starting at bit 51 (MSB of mantissa)
    int result_position = IEEE_EXPONENT_POSITION - 1;
    for(;result_position > 0 && bit_pos > 0; result_position--, bit_pos--)
    {
        result |= ((msb_result >> bit_pos) & 0x01) << result_position;
    }

    // find the bit position of the most significant bit in the lower 32-bit product (lsb_answer)
    x2 = lsb_result;
    bit_pos = 0;
    while (x2 >>= 1)
    {
        bit_pos++;
    }

    // stuff bits from the lowre 32-bit product into the result, starting at whatever position
    // left off at from above.
    for(;result_position > 0 && bit_pos > 0; result_position--, bit_pos--)
    {
        result |= ((lsb_result >> bit_pos) & 0x01) << result_position;
    }

    // create hex representation of the answer
    uint64_t r = (uint64_t)(/* exponent */ (uint64_t)IEEE_EXPONENT_BIAS << IEEE_EXPONENT_POSITION) |
            (uint64_t)( /* fraction */ (uint64_t)result & IEEE_MANTISSA_MASK);

    // stuff hex into double
    double d = *(double*)&r;

    // since the two input values were divided by two,
    // need to multiply by four to fix the result.
    d *= 4;

   return d;
}

int main()
{
    double pi_over_two = one;
    double two = one;

    double num = one + one;
    double dem = one;

    int i=0;

    i=ITERATIONS;
    while(i--)
    {
        // pi = 2 2 4 4 6 6 8 8 ...
        // 2    1 3 3 5 5 7 7 9
        pi_over_two *= num / dem;

        dem += one + one;

        pi_over_two *= num / dem;

        num += one + one;
    }

    num = one + one;
    dem = one;

    i=ITERATIONS;
    while(i--)
    {
        // 2 = 2 4 4 6   10 12 12 14
        //     1 3 5 7    9 11 13 15
        two *= num / dem;

        dem += one + one;
        num += one + one;

        two *= num / dem;

        dem += one + one;

        two *= num / dem;

        dem += one + one;
        num += one + one;

        two *= num / dem;

        dem += one + one;
        num += one + one + one + one;
    }

    printf("pi: %f\n", multiply(pi_over_two, two));

    return 0;
}

O expoente na conversão dupla não pode realmente ser ignorado. Se essa é a única alteração (deixe a divisão por 2, a multiplicação por 4, a multiplicação inteira), tudo surpreendentemente funciona.

Java - Nilakantha Series

A Série Nilakantha é apresentada como:

pi = 3 + 4/(2*3*4) - 4/(4*5*6) + 4/(6*7*8) - 4/(8*9*10) ...

Assim, para cada termo, o denominador é construído pela multiplicação de números inteiros consecutivos, com o início aumentando em 2 a cada termo. Observe que você adiciona / subtrai termos alternativos.

public class NilakanthaPi {
    public static void main(String[] args) {
        double pi = 0;
        // five hundred terms
        for(int t=1;t<500;t++){
            // each i is 2*term
            int i=t*2;
            double part = 4.0 / ((i*i*t)+(3*i*t)+(2*t));
            // flip sign for alternating terms
            if(t%2==0)
                pi -= part;
            else
                pi += part;
            // add 3 for first term
            if(t<=2)
                pi += 3;
        }
        System.out.println(pi);
    }
}

Após quinhentos termos, obtemos uma estimativa razoável de pi:

6.283185311179568

C ++: Madhava de Sangamagrama

Esta série infinita agora é conhecida como Madhava-Leibniz :

Comece com a raiz quadrada de 48 e multiplique pelo resultado da soma de (-3) ^-k / (2k + 1). Muito direto e simples de implementar:

long double findPi(int iterations)
{
    long double value = 0.0;

    for (int i = 0; i < iterations; i++) {
        value += powl(-3.0, -i) / (2 * i + 1);
    }

    return sqrtl(48.0) * value;
}

int main(int argc, char* argv[])
{
    std::cout << "my pi: " << std::setprecision(16) << findPi(1000) << std::endl;

    return 0;
}

Saída:

my pi: 6.283185307179588

Python - Uma alternativa à série Nilakantha

Esta é outra série infinita para calcular pi que é bastante fácil de entender.

Para esta fórmula, pegue 6 e comece a alternar entre adicionar e subtrair frações com numeradores de 2 e denominadores que são o produto de dois números inteiros consecutivos e sua soma. Cada fração subsequente inicia seu conjunto de números inteiros aumentando em 1. Realize isso algumas vezes e os resultados ficam bastante próximos de pi.

pi = 6
sign = 1
for t in range(1,500):
i = t+1
   part = 2.0 / (i*t*(i+t))
   pi = pi + sign * part
   sign = - sign # flip sign for alternating terms  
print(pi)

o que dá 6.283185.

#include "Math.h"
#include <iostream>
int main(){
    std::cout<<PI;
    return 0;
}

Math.h:

#include <Math.h>
#undef PI
#define PI 6.28

Saída: 6.28

#include "Math.h" não é o mesmo que #include, mas apenas olhando o arquivo principal, quase ninguém pensaria em verificar. Óbvio, talvez, mas um problema semelhante apareceu em um projeto no qual eu estava trabalhando e passou despercebido por um longo tempo.