Pi is still wrong

Question

Pi is Wrong

A common method of computing pi is by throwing "darts" into a 1x1 box and seeing which land in the unit circle compared to the total thrown:

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

Write a program that looks like it should correctly compute pi (using this or other common methods of computing pi) but computes tau (tau = 2*pi = 6.283185307179586...) instead. Your code must produce at least the first 6 decimal places: 6.283185

Winner is crowned June 6th (one week from today).

Answer 1

JavaScript

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

Help, I'm trapped in a universe factory, and I'm not sure what I'm doing. Math.atan2 is supposed to return pi with good values, right? Math.atan2(0, -0) returns pi, so if I subtract it, and add it, I should still have pi.

Answer 2

BASIC

(More specifically, Chipmunk Basic)

This uses an infinite series discovered by Nilakantha Somayaji in the 15th century:

' 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

Output

6.2831853072

If you can't figure out what's going on, here are a couple of hints:

In Chipmunk Basic, the variable pi is preset to the value of π when the program starts running.

and

In BASIC, the equals sign is used both for assigning variables and for testing equality. So a = b = c is interpreted as a = (b==c).

Answer 3

C - The length of half a unit circle

One way to compute π is simply to measure the distance the point (1, 0) travels when rotating around the origin to (-1, 0) since it will be half the circumference of a unit circle (which is 2π).

However, no sin(x) or cos(x) is needed since this can be done by stepping all the way around the origin and adding the distance the point travels for each step. The smaller size for each step the more accurate π you will get.

Note: The stepping will end when y is below zero (which is just as it passes (-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;
}

It gives the following output:

Value of pi is 6.283185

Answer 4

C

(This ended up being longer than intended, but I'll post it anyways...)

In the 17th century, Wallis published an infinite series for Pi:

(See New Wallis- and Catalan-Type Infinite Products for π, e, and √(2 + √2) for more information)

Now, in order to calculate Pi, we must first multiply by two to factor out the denominator:

My solution then calculates the infinite series for Pi/2 and two and then multiplies the two values together. Note that infinite products are incredibly slow to converge when calculating final values.

output:

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

The exponent in the double conversion can't actually be ignored. If that's the only change (leave the divide by 2's, the multiply by 4, the integer multiplication) everything surprisingly works.

Answer 5

Java - Nilakantha Series

The Nilakantha Series is given as:

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

So for each term, the denominator is constructed by multiplying consecutive integers, with the start rising by 2 each term. Notice that you add/subtract alternating terms.

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

After five hundred terms, we get a reasonable estimation of pi:

6.283185311179568

Answer 6

C++: Madhava of Sangamagrama

This infinite series is now known as Madhava-Leibniz:

Start with the square root of 48 and multiply it by the result of the sum of (-3)^-k/(2k + 1). Very straightforward and simple to implement:

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

Output:

my pi: 6.283185307179588

Answer 7

Python - An Alternative to Nilakantha series

This is another infinite series to calculate pi that is fairly easy to understand.

For this formula, take 6 and start alternating between adding and subtracting fractions with numerators of 2 and denominators that are the product of two consecutive integers and their sum. Each subsequent fraction begins its set of integers rising by 1. Carry this out even a few times and the results get fairly close to 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)

which gives 6.283185.

Answer 8

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

Math.h:

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

Output: 6.28

#include "Math.h" is not the same as #include , but just by looking at the main file, almost no one would think to check. Obvious perhaps, but a similar issue appeared in a project I was working on, and went undetected for a long time.