www.digitalmars.com         C & C++   DMDScript  

digitalmars.D - pi sample program w/ metaprogramming optimizations

reply "Craig Black" <cblack ara.com> writes:
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
Apr 29 2006
next sibling parent reply "Craig Black" <cblack ara.com> writes:
Content-Type: text/plain;
	charset="iso-8859-1"
Content-Transfer-Encoding: quoted-printable

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=3D4000;

byte[] p;
byte[] t;
int q;

int main(char[][] args)
{
 time_t startime, endtime;

 if (args.length =3D=3D 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 =3D 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 =3D q + 1;
 t.length =3D 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 =3D %d.",cast(int)(p[0]));
 for (int i =3D 1; i <=3D 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=20
 * 2 at compile time
 */
template isPow2(uint n) { const bool isPow2 =3D (n % 2) !=3D 0 ? false : =
isPow2!(n / 2); }
template isPow2(uint n : 0) { const bool isPow2 =3D false; }
template isPow2(uint n : 1) { const bool isPow2 =3D false; }
template isPow2(uint n : 2) { const bool isPow2 =3D true; }

/* Use a template to get good compile time optimizations for
 * division and multiplication=20
 */
template fastdiv(int divisor)
{
 void fastdiv()
 {
  int r; // remainder

=20
  for (int i =3D 0; i <=3D q; i++) {
   int b =3D t[i] + r * 10;
   int q =3D b / divisor;
   static if(isPow2!(divisor))=20
    r =3D b % divisor; // compiler can optimize modulus for powers of 2
   else=20
    r =3D b - q * divisor;
   t[i] =3D q;
  }
 }
}

// Computes the arctangent of 2
void arctan2()
{
 int n;

 t[0] =3D 1;
 fastdiv!(2);
 add();
 n =3D 1;
 do {
  mul(n);
  fastdiv!(4);
  div(n +=3D 2);
  if ((((n-1) / 2) % 2) =3D=3D 0)
   add();
  else
   sub();
 } while (!tiszero());
}


// Computes the arctangent of 3
void arctan3()
{
 int n;

 t[0] =3D 1;
 fastdiv!(3);
 add();
 n =3D 1;
 do {
  mul(n);
  fastdiv!(9);
  div(n +=3D 2);
  if ((((n-1) / 2) % 2) =3D=3D 0)
   add();
  else
   sub();
 } while (!tiszero());
}

void add()
{

 for (int j =3D q; j >=3D 0; j--)
 {
  if (t[j] + p[j] > 9) {
   p[j] +=3D t[j] - 10;
   p[j-1]++;
  } else
   p[j] +=3D t[j];
 }
}

void sub()
{
 for (int j =3D q; j >=3D 0; j--)
  if (p[j] < t[j]) {
   p[j] -=3D t[j] - 10;
   p[j-1]--;
  } else
   p[j] -=3D t[j];
}


void mul(int multiplier)
{
 int c;


 for (int i =3D q; i >=3D 0; i--) {
  int b =3D t[i] * multiplier + c;
  c =3D b / 10;
  t[i] =3D b - c * 10;=20
 }
}


void mul4()
{
 int c;

 for (int i =3D q; i >=3D 0; i--) {
                int b =3D p[i] * 4 + c;
                c =3D b / 10;
  p[i] =3D b - c * 10;=20
 }
}


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 =3D 1.0 / divisor;

 for (int i =3D 0; i <=3D q; i++) {
                int b =3D t[i] + r * 10;=20
  int quotient =3D cast(int)(cast(double)b * idiv);
  r =3D b - quotient * divisor;
  t[i] =3D quotient;
 }
}

int tiszero()
{
 for (int k =3D 0; k <=3D q; k++)
  if (t[k] !=3D 0)
   return false;
 return true;
}
Apr 29 2006
parent James Dunne <james.jdunne gmail.com> writes:
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
Apr 30 2006
prev sibling parent "Craig Black" <cblack ara.com> writes:
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;
}
Apr 29 2006