pi sample program w/ metaprogramming optimizations

[Craig Black]

For some reason the attachment doesn't seem to be working, at least for me.  Here's the code just in case.  Please excuse the spacing.  For some reason the tabs are only represented by one space.

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



-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.puremagic.com/pipermail/digitalmars-d/attachments/20060429/36aeabe8/attachment.html>