import std.c.stdio;
import
std.c.stdlib;
import std.c.time;
const int LONG_TIME=4000;
byte[] p;
byte[] t;
int q;
int main(char[][] args)
{
time_t
startime, endtime;
if (args.length == 2)
{
sscanf(&args[1][0],"%d",&q);
} else
{
printf("Usage: pi
[precision]\n");
exit(55);
}
if (q <
0)
{
printf("Precision was too low, running with
precision of 0.\n");
q = 0;
}
if (q >
LONG_TIME)
{
printf("Be prepared to wait a
while...\n");
}
// Compute one more digit than we display to
compensate for rounding
q++;
p.length = q + 1;
t.length = q +
1;
// Compute
pi
std.c.time.time(&startime);
arctan2();
arctan3();
mul4();
std.c.time.time(&endtime);
// Return to the number of digits we want to
display
q--;
// Print pi
printf("pi =
%d.",cast(int)(p[0]));
for (int i = 1; i <= q;
i++)
printf("%d",cast(int)(p[i]));
printf("\n");
printf("%ld
seconds to compute pi with a precision of %d
digits.\n",endtime-startime,q);
return 0;
}
/* Template that determines if a number is a power
of
* 2 at compile time
*/
template isPow2(uint n) { const
bool isPow2 = (n % 2) != 0 ? false : isPow2!(n / 2); }
template isPow2(uint n
: 0) { const bool isPow2 = false; }
template isPow2(uint n : 1) { const bool
isPow2 = false; }
template isPow2(uint n : 2) { const bool isPow2 = true;
}
/* Use a template to get good compile time
optimizations for
* division and multiplication
*/
template fastdiv(int divisor)
{
void
fastdiv()
{
int r; // remainder
for (int i = 0; i <= q;
i++) {
int b = t[i] + r * 10;
int q =
b / divisor;
static if(isPow2!(divisor))
r = b % divisor; // compiler can optimize modulus
for powers of 2
else
r = b - q
* divisor;
t[i] =
q;
}
}
}
// Computes the arctangent of 2
void
arctan2()
{
int n;
t[0] =
1;
fastdiv!(2);
add();
n = 1;
do
{
mul(n);
fastdiv!(4);
div(n +=
2);
if ((((n-1) / 2) % 2) ==
0)
add();
else
sub();
}
while (!tiszero());
}
// Computes the arctangent of 3
void
arctan3()
{
int n;
t[0] =
1;
fastdiv!(3);
add();
n = 1;
do
{
mul(n);
fastdiv!(9);
div(n +=
2);
if ((((n-1) / 2) % 2) ==
0)
add();
else
sub();
}
while (!tiszero());
}
void add()
{
for (int j = q; j >= 0;
j--)
{
if (t[j] + p[j] > 9)
{
p[j] += t[j] -
10;
p[j-1]++;
}
else
p[j] += t[j];
}
}
void sub()
{
for (int j = q; j >=
0; j--)
if (p[j] < t[j]) {
p[j] -= t[j] -
10;
p[j-1]--;
}
else
p[j] -= t[j];
}
void mul(int multiplier)
{
int
c;
for (int i = q; i >= 0; i--)
{
int b = t[i] * multiplier + c;
c = b /
10;
t[i] = b - c * 10;
}
}
void mul4()
{
int c;
for (int i = q; i >= 0; i--)
{
int b = p[i] * 4 +
c;
c = b / 10;
p[i] = b - c * 10;
}
}
void div(int divisor)
{
int r; //
remainder
/*
Optimization: It's faster to do floatint point multiplication
than
* integer division.
This is because there is no integer division
instruction.
* Under the
hood, integer division is essentially floating-point
division.
*/
double idiv = 1.0 /
divisor;
for (int i = 0; i <= q; i++)
{
int b = t[i] + r * 10;
int quotient = cast(int)(cast(double)b *
idiv);
r = b - quotient * divisor;
t[i] =
quotient;
}
}
int tiszero()
{
for (int k = 0; k
<= q; k++)
if (t[k] != 0)
return
false;
return true;
}