# 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).

• Why is the winner not crowned June 28th?? May 30, 2014 at 17:35
• I'm not sure why a winner needs to be crowned in a popularity contest. May 30, 2014 at 19:34
• I don't get it. This like asking for a function that appears to return 1 but returns 2. Who are we fooling here? May 31, 2014 at 4:06
• @ja72 The reader of the code :) May 31, 2014 at 7:22
• Everyone knows that pau is the correct one. :P Jun 1, 2014 at 18:04

# 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.

• I think I'm just gonna go lie down and cry. Goddamnit, JavaScript. May 31, 2014 at 9:59
• explanation please? :) May 31, 2014 at 22:40
• Counterclockwise angle in radians between x axis and point (Y,X). Sign of the Y point determines if this is a positive or negative angle, and this becomes π - (-π)
– user21677
Jun 1, 2014 at 0:48
• 0_o >>> 0 === -0 ;true ;>>> Math.atan2(0, 0) ;0 ;>>> Math.atan2(0, -0) ;3.141592653589793 Jun 1, 2014 at 5:37
• @JackM, that statement is always appropriate to say :) Though in this case, it is due to the IEEE standard, and many languages (not just JS) have the zero vs negative zero problem. Jun 1, 2014 at 6:59

# 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).

• Wait I don't get it, so i equals false? And then you add 2 to it? And it works??? Jun 1, 2014 at 13:28
• @Dunno: Sure, the loops starts at i == false which is similar to i == 0. The point is that the initial value for the accumulator pi isn't 0… Jun 1, 2014 at 16:44
• @Bergi yeah I just can't wrap my head around the fact that false + 2 == 2 :D Jun 1, 2014 at 18:43
• @Dunno Dynamic typing etc.: false is implicitly converted to 0 when doing arithmetic. You also have the same apparent behaviour in C which lacks a bool type, and uses 0 and non-zero to represent false and true repsectively. Not that it's elegant, but hey, that's how it works. Jun 2, 2014 at 10:29

## 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 ). 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

• Seems legit… Definitely. May 31, 2014 at 3:29
• Your length macro is missing a sqrt. Is that intended? x and y are also swapped between the definition and call (with no effect) May 31, 2014 at 18:02
• @BenVoigt Shhh! Don't spoil the trick, but yes. sqrt was accidentally omitted so that the value of pi was printed as 6,28... Also +1 for noticing x and y which i didn't! May 31, 2014 at 20:19
• oh, now I see you're tracing not a unit circle, but one with radius 2. Yup, that works great. May 31, 2014 at 20:24
• I must confess that before understanding how it works I wasted a couple of minutes by not ignoring that line... Jun 1, 2014 at 14:29

# 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_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;

// 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.

# 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


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


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.

#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.

• A clever solution nonetheless. Apr 30, 2015 at 0:57