Take the 2-minute tour ×
Programming Puzzles & Code Golf Stack Exchange is a question and answer site for programming puzzle enthusiasts and code golfers. It's 100% free, no registration required.

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

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

8 Answers 8

up vote 48 down vote accepted

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.

share|improve this answer
10  
I think I'm just gonna go lie down and cry. Goddamnit, JavaScript. –  Jack M May 31 at 9:59
2  
explanation please? :) –  Jaa-c May 31 at 22:40
2  
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 π - (-π) –  tolos Jun 1 at 0:48
5  
0_o >>> 0 === -0 ;true ;>>> Math.atan2(0, 0) ;0 ;>>> Math.atan2(0, -0) ;3.141592653589793 –  Izkata Jun 1 at 5:37
3  
@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. –  Paul Draper Jun 1 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).

share|improve this answer
8  
Sneaky shit here! –  tomsmeding May 31 at 7:22
    
Wait I don't get it, so i equals false? And then you add 2 to it? And it works??? –  Dunno Jun 1 at 13:28
2  
@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… –  Bergi Jun 1 at 16:44
1  
@Bergi yeah I just can't wrap my head around the fact that false + 2 == 2 :D –  Dunno Jun 1 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. –  Georges Dupéron Jun 2 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 ).

enter image description here

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
share|improve this answer
3  
Seems legit… Definitely. –  bjb568 May 31 at 3:29
1  
Your length macro is missing a sqrt. Is that intended? x and y are also swapped between the definition and call (with no effect) –  Ben Voigt May 31 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! –  Thism2 May 31 at 20:19
    
oh, now I see you're tracing not a unit circle, but one with radius 2. Yup, that works great. –  Ben Voigt May 31 at 20:24
5  
I must confess that before understanding how it works I wasted a couple of minutes by not ignoring that line... –  loreb Jun 1 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:

enter image description here

(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:

enter image description here

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.

share|improve this answer

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
share|improve this answer

Python - An Alternative to Nilakantha series

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

enter image description here

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.

share|improve this answer

C++: Madhava of Sangamagrama

This infinite series is now known as Madhava-Leibniz:

Series

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
share|improve this answer

It's terribly inefficient and not very accurate but it works. Idk wtf the OP is doing here. Here is a python example.

import random
import math

def calc_pi(ticks):
    n = 0.0
    for i in range(ticks):
        x = abs(random.random() - .5)
        y = abs(random.random() - .5)
        distance = math.sqrt(x**2 + y**2)
        if distance <= .5:
            n +=1.0
    return (4.0*n/ticks)

print calc_pi(1000)
share|improve this answer
3  
This solution calculates Pi and does so in an obvious manner, the question is to appear to calculate pi but actually calculate Tau –  MD-Tech Jun 2 at 8:42

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.