I wanted to explore optimizing code in D so I spent a few hours optimizing the pi.d sample program that comes with the dmd compiler.  I used some neat metaprogramming tricks.  On my computer it is almost twice as fast as the original version.

-Craig

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

Here' a version that uses more templates.  It performs the same as the first version I submitted.

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)
{
 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
 time_t startime, endtime;
 std.c.time.time(&startime);
 arctan!(2);
 arctan!(3);
 fastmul!(4);
 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;
  }
 }
}

/* Template to compute the arctangent.
 * The constant input parameter allows the compiler to optimize.
 */
template arctan(int s)
{
 void arctan()
 {
  int n;

  t[0] = 1;
  fastdiv!(s);
  add();
  n = 1;
  do {
   mul(n);
   fastdiv!(s*s);
   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;
 }
}

/* Fast multiplication using template */
template fastmul(int multiplier)
{
 void fastmul()
 {
  int c;

  for (int i = q; i >= 0; i--) {
   int b = p[i] * multiplier + 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;
}

Craig Black wrote:
> 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.
> 
> [snip]
>  

Attachments don't work properly with the web-interface, but any decent newsgroup reader will work fine.  It's been this way for quite a while...

-- 
Regards,
James Dunne