Pi jest w błędzie

Powszechną metodą obliczania pi jest rzucanie „rzutkami” do pudełka 1x1 i sprawdzanie, który ląd w okręgu jednostki w porównaniu do całkowitej rzuconej:

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

Napisz program, który wygląda tak, jakby poprawnie obliczał pi (używając tej lub innych typowych metod obliczania pi), ale zamiast tego oblicza tau (tau = 2 * pi = 6.283185307179586 ...). Twój kod musi generować co najmniej pierwsze 6 miejsc po przecinku: 6.283185

Zwycięzca zostaje koronowany 6 czerwca (za tydzień od dzisiaj).

JavaScript

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

Pomocy, jestem uwięziony w fabryce wszechświatów i nie jestem pewien, co robię. Math.atan2ma zwracać pi z dobrymi wartościami, prawda? Math.atan2(0, -0)zwraca pi, więc jeśli go odejmę i dodam, nadal powinienem mieć pi.

PODSTAWOWY

(Mówiąc dokładniej, Chipmunk Basic )

Wykorzystuje to nieskończoną serię odkrytą przez Nilakantha Somayaji w XV wieku:

' 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

Wydajność

6.2831853072

Jeśli nie możesz dowiedzieć się, co się dzieje, oto kilka wskazówek:

W Chipmunk Basic zmienna pi jest wstępnie ustawiona na wartość π, gdy program zacznie działać.

i

W języku BASIC znak równości służy zarówno do przypisywania zmiennych, jak i do testowania równości. Zatem a = b = c interpretowane jest jako a = (b == c) .

C - Długość połowy koła jednostkowego

Jeden sposób obliczenia gatunku jest jedynie pomiar odległości punkt (1, 0)przemieszcza się podczas obracania wokół pochodzenie do (-1, 0)ponieważ będzie połowę obwodu okręgu jednostkowym (co 2π ).

Jednak nie jest konieczne sin(x)lub nie cos(x), ponieważ można to zrobić, przechodząc dookoła początku i dodając odległość, jaką punkt pokonuje dla każdego kroku . Im mniejszy rozmiar dla każdego kroku, tym dokładniejszy będzie π .

Uwaga: Stepping zakończy się, gdy y jest poniżej zera (czyli tak jak mija (-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;
}

Daje następujące dane wyjściowe:

Value of pi is 6.283185

do

(W końcu było to dłuższe niż zamierzone, ale i tak opublikuję ...)

W XVII wieku Wallis opublikował nieskończoną serię dla Pi:

(Więcej informacji na temat π, e i √ (2 + √2), aby uzyskać więcej informacji, patrz Nowe nieskończone produkty typu Wallis i katalońskiego )

Teraz, aby obliczyć Pi, musimy najpierw pomnożyć przez dwa, aby wyliczyć mianownik:

Moje rozwiązanie oblicza następnie szereg nieskończony dla Pi / 2 i dwóch, a następnie mnoży te dwie wartości razem. Zauważ, że nieskończone produkty są niezwykle powolne w zbieraniu się podczas obliczania ostatecznych wartości.

wydajność:

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;
}

Wykładnika podwójnej konwersji nie można tak naprawdę zignorować. Jeśli to jedyna zmiana (zostaw dzielenie przez 2, pomnóż przez 4, mnożenie przez liczby całkowite), wszystko zaskakująco działa.

Java - seria Nilakantha

Seria Nilakantha jest podawana jako:

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

Tak więc dla każdego terminu mianownik jest konstruowany przez pomnożenie kolejnych liczb całkowitych, przy czym początek rośnie o 2 za każdym razem. Zauważ, że dodajesz / odejmujesz na przemian terminy.

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);
    }
}

Po pięciuset terminach otrzymujemy rozsądne oszacowanie liczby pi:

6.283185311179568

C ++: Madhava z Sangamagrama

Ta nieskończona seria jest obecnie znana jako Madhava-Leibniz :

Zacznij od pierwiastka kwadratowego z 48 i pomnóż go przez wynik sumy (-3) ^-k / (2k + 1). Bardzo prosty i łatwy do wdrożenia:

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;
}

Wydajność:

my pi: 6.283185307179588

Python - alternatywa dla serii Nilakantha

To kolejna nieskończona seria do obliczenia pi, która jest dość łatwa do zrozumienia.

Dla tej formuły weź 6 i zacznij na przemian między dodawaniem i odejmowaniem ułamków z licznikami 2 i mianownikami, które są iloczynem dwóch kolejnych liczb całkowitych i ich sumy. Każda kolejna frakcja zaczyna swój zestaw liczb całkowitych rosnących o 1. Wykonaj to nawet kilka razy, a wyniki zbliżą się do 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)

co daje 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

Wyjście: 6,28

#include „Math.h” to nie to samo co #include, ale patrząc na główny plik, prawie nikt nie pomyślałby o sprawdzeniu. Być może oczywiste, ale podobny problem pojawił się w projekcie, nad którym pracowałem i przez długi czas pozostawał niewykryty.