Factorials and Textures

Dann Corbit and Rafael Collantes-Bellido

Dann is a software engineer with SolutionsIQ on contract at Microsoft. He can be contacted at [email protected]. Rafael is associated with the Universidad Pontificia Comillas in Madrid. He can be reached at [email protected].

This month,"Algorithm Alley" brings you two short, unrelated pieces. First, Dann Corbit reminds us that the classic memory-speed tradeoff still applies, and memory is getting cheaper. Then Rafael Collantes uses a 50-year-old set of differential equations to produce natural-looking, tilable textures. Enjoy!

Tim Kientzle

The factorial function n! is often used in the calculation of odds. It occurs frequently in statistics and other branches of math and science. The factorial n! is defined for a positive integer n as n!=n(n-1)(n-2)...1. As a special case, 0! is defined as 1.

Many computer-science students are asked to write a recursive version of the factorial function in college. Usually, this looks something like Example 1. Notice that to compute n!, this function calls itself n times and performs n multiplications. Though it does show the elegance of recursion, this is not a very efficient way to calculate probabilities.

Another way to calculate factorials is to use the gamma function: a generalization of the factorial that can be written as an infinite series. By using a finite number of terms, you can generate a useful approximation that can be calculated much more efficiently than repeated multiplication. Listing One shows how to use Stirling's Approximation to compute the factorial.

However, because the factorial function grows at such a staggering rate, there are really only a small number of values that can be stored in a double-precision floating-point variable. On an 80x86 CPU, 1754! is the largest value that can be computed without overflow. To store every possible value, you only need a table of 1755 10-byte numbers, consuming a modest 17 KB of memory. The speed improvement from using a table is quite dramatic. Listing Two outputs a C declaration for an array of factorials. You can save the output to a file and compile it directly into your program. Listing Three is a function that uses this table to rapidly compute the "combinations" function.

This approach can be used for many other functions besides the factorial. For example, the double factorial n!!=n(n-2)(n-4)...1 or the triple factorial n!!!=n(n-3)(n-6)...1, can both benefit from this kind of precomputation. So when a problem allows for a precomputed solution, why not table it for later?

For more information, see An Atlas of Functions, by Jerome Spannier and Keith B. Oldham (Washington, DC: Hemisphere Publishing Corp., 1987), and "Eric's Treasure Trove of Mathematics," by Eric Weisstein, at http://www.gps.caltech.edu /~eww/math/math0.html.

Alan Turing's contributions weren't limited to digital computers. In 1952, he published "The Chemical Basis of Morphogenesis," which attempted to model the process of pattern formation on the skin of animals. In his paper, Turing proposed a Reaction-Diffusion equation as a model for the chemical processes that govern the formation of stripes, dots, and patches on animals. In this article, I'll examine how Turing's system can be used to obtain natural looking textures, such as Figure 1.

Turing's system is a complicated system of nonlinear, partial differential equations. The Reaction part describes how the two components of the system interact with each other. It is also the part of the system that is nonlinear, adding to the mathematical difficulty. The Diffusion part describes how each chemical individually flows across cells. Fortunately, you don't need to understand the detailed mathematics to take advantage of Turing's idea.

The cells of the animal's skin will be modeled with two two-dimensional arrays, representing the concentrations of two chemicals, which I'll call A and B. You start with random values, then simulate the evolution of these two concentrations. For each cell, the rate of change of the concentration of A depends on the concentration of A in this cell and its immediate neighbor, and the concentration of B in this cell only. Similarly, B depends on the nearby concentrations of B and the concentration of A in this cell only.

Let A[i][j] and B[i][j] be the concentrations of the two components in cell (i,j) at a given time. Assume ReA[i][j] and ReB[i][j] are the Reaction functions for the cell (i,j), and DiA[i][j] and DiB[i][j] are the Diffusion functions. In the next step of the simulation, the concentrations are given by Figure 2. Figure 3 shows the Reaction and Diffusion functions for Turing's model. Note the two constants CA and CB, which control the diffusion of the chemicals.

One problem needs to be addressed: What cells should be considered the "neighbors" of cells at the edge of the grid? One common approach is to think of the grid as the surface of a torus (donut), so the upper edge meets the lower, and the right edge meets the left one. This is the kind of topology used in the Asteroids video game. Other approaches include imagining that outside the grid there is a region in which the values of A and B are equal and constant in such a way that the reaction functions are 0, or that the edges are "mirrors" so that A[-1][j] is the same as A[1][j]. The torus approach has an advantage in that we will be able to concatenate the resulting textures and there will be no visible edges. This way they can be used as tiles to cover a surface without it being possible to tell where one tile begins or ends.

The C code that implements this technique (Listing Four) is straightforward and should compile with any ANSI C compiler. I've tested it on the Symantec C++ compiler for the Macintosh (7.0.4) and on the HP-UX C compiler for the Hewlett-Packard series 700 workstations.

The code asks for the CA and CB constants for the calculations, and uses a simple iterative Euler method to solve the differential equation. The number of iterations is fixed at 2000, which is enough for this system. After finishing the calculation, the final values are linearly scaled to cover the 0-255 range. The appearance of the texture depends on the values of CA and CB. Good starting points are CA=1 and CB=16. Different executions of the program will produce different output, even with the same values of CA and CB, since the random initial values influence the solution of the problem.

The program outputs a graphic file (available electronically) in the Portable Graymap (.PGM) format. The file format consists of a header and the raw data of the image, with each pixel gray level represented by a single byte. The header is an ASCII string written with the printf() format P5 %d %d 255. The two integers in the string are the horizontal and vertical sizes of the image. Most workstations include software to view such images, and the Macintosh shareware program "GraphicConverter" can also handle them.

The size of the texture is set in the maxI and maxJ variables defined at the beginning of the code. They can be any positive integer, but remember that calculation time increases with the square of the size.

Example 1: Recursive evaluation of the factorial function.

int factorial(int n) {
  if (n>0) 
return n * factorial(n-1);
  else if (n==0) return 1;
  else return -1;
}

Figure 2: Computing the chemical concentrations from the Reaction-Diffusion values.

A[i][j]= A[i][j]+ T*(ReA[i][j]+DiA[i][j])
B[i][j]= B[i][j]+ T*(ReB[i][j]+DiB[i][j])

Figure 3: Computing the Reaction and Diffusion functions.

ReA[i][j]= A[i][j]*B[i][j] - A[i][j] - 12
DiA[i][j]= CA*(A[i][j+1]-2*A[i][j]+A[i][j-1]+A[i+1][j]-2*A[i][j]+A[i-1][j])
ReB[i][j]= 16 - A[i][j]* B[i][j]
DiB[i][j]= CB*(B[i][j+1]-2*B[i][j]+B[i][j-1]+B[i+1][j]-2*B[i][j]+B[i-1][j])

Listing One

long double ldGammaSeries(long double x) {
   long double ldSterlingTerms[11] = {
      1.e0L,            1.e0L / 12.e0L,
      1.e0L / 288.e0L,         -139.e0L / 51840.e0L,
      -571.e0L / 2488320.e0L,   163879.e0L / 209018880.e0L,
      5246819.e0L / 75246796800.e0L,  
      -534703531.e0L / 902961561600.e0L,
      -4483131259.e0L / 86684309913600.e0L,
      432261921612371.e0L / 514904800886784000.e0L,
      6232523202521089.e0L / 86504006548979712000.e0L
   };
   long double ldResult;
   long double ldSumm = ldSterlingTerms[0];
   long double ldXproduct = x;
   int i;

   /* Stirling's Asymptotic Approximation.  Requires large X for accuracy. */
   /* 2.5066... is sqrt(pi+pi) */
   ldResult = expl(-x) * powl(x, x - 0.5e0L) * 2.5066282746310005024157652e0L;
   for ( i = 1; i < 10; i++ ) {
      ldSumm += ldSterlingTerms[i]/ldXproduct;
      ldXproduct *= x;
   }
   return (ldResult*ldSumm);
}

Listing Two

#include <ctype.h>
#include <stdio.h>
#include <stdlib.h>
#include <signal.h>
#include <limits.h>

void vExitHandler( int sig_number );

int _cdecl main( int iArgc, char **ppszArgv )
{
   unsigned long i; /* Index for calculation of Factorial (!) */
   long double ld1Product = 1.0e0L; /* Term for Factorial (!) */

   /* Install our own signal hander routine: */
   signal( SIGFPE, vExitHandler );

   /* Print the declaration for the data table */

   puts( "long double ldFactorialTable[] = {" );
   /* Print the value for 0!, (0!!), (0!!!), (0!!!!) */
   printf( "%31.20LeL,\n", ld1Product );

   /* Print factorials out until we go into the FPU signal handler. */
   for ( i = 1; i < ULONG_MAX; i ++ )
   {
      ld1Product *= ( long double ) i;
      printf( "%31.20LeL,\n", ld1Product );
   }
   puts( "ULONG_MAX won't overflow.  WHAT AN FPU YOU'VE GOT!" );
   return 1; /* We wanted to overflow, how did we get here? */
}

/* This routine is executed when our floating point math overflows. */
void vExitHandler( int sig_number )
{
   sig_number = sig_number; /* Shut up the compiler warnings. */
   puts( "   -1.0e0L\n};" ); /* Print out end of array tokens. */
   exit( 0 ); /* This is our 'normal' return point. */
}

/*  Calculate the combinations of a collection of <n> objects taken
<r> at a time.  For combinations, the order of the objects has no
significance.  Hence, a poker hand with four aces has the same value
regardless of the order in which the aces were delt.

ldCombinations is computed as n!/( r!*( n-r )! ).
Error is indicated by return of ( -1 )  */
long double ldCombinations( int n, int r ) {
   if ( (0 <= r) && ( r <= n ) && ( n <= MAXFACTORIAL ) )
       return ( ldFactorialTable[n] /
                 (ldFactorialTable[r]*ldFactorialTable[n-r]) );
   else return ( -1 );
}

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

main()
{
   int i,j,n;
   double **A,**B,CA,CB;
   int iterations = 2000, maxI=64, maxJ=64;
   
   /**** Get CA and CB ****/
   {
      char str[20];

      printf ("CB->");
      gets (str);
      CB=atoi(str);
      
      printf ("CA->");
      gets (str);
      CA=atoi(str);
   }
   /**** Initialize A and B randomly ****/
   srand(clock());
   A=(double **) calloc (maxI,sizeof(double*));
   B=(double **) calloc (maxI,sizeof(double*));
   for (i=0;i<maxI;i++) {
      A[i]= (double *) calloc (maxJ,sizeof(double));
      B[i]= (double *) calloc (maxJ,sizeof(double));
      if (!B[i] || !A[i]) {
         printf ("Memory problems!\n");
         exit(1);
      }
      for ( j=0; j<maxJ; j++) {
         A[i][j]=12*(double) rand() /RAND_MAX;
         B[i][j]=12*(double) rand() /RAND_MAX;
      }
   }
   /**** Use Turing's Model ****/
   for (n=0;n<iterations;n++) {
      printf("%d\n",iterations-n);
      for (i=0;i<maxI;i++) {
         int iplus1 = i+1, iminus1 = i-1;
         if (i==0)   iminus1=maxI-1;
         if (i==maxI-1) iplus1=0;
         for (j=0;j<maxJ;j++) {
            int jplus1 = j+1, jminus1 = j-1;
            double Aold = A[i][j];
            double DiA,ReA,DiB,ReB;
            if ( j == 0 ) jminus1=(maxJ-1);
            if ( j == maxJ-1 ) jplus1=0;
            /*** COMPONENT A ****/
            DiA = CA * ( A[iplus1][j]-2*A[i][j]+A[iminus1][j]
                        +A[i][jplus1]-2*A[i][j]+A[i][jminus1]);
            ReA = A[i][j]*B[i][j] - A[i][j] - 12.0;
            A[i][j]=A[i][j]+0.01*(ReA + DiA);
            if (A[i][j]<0) A[i][j]=0;
            /*** COMPONENT B ****/
            DiB =CB * ( B[iplus1][j]-2*B[i][j]+B[iminus1][j]
                       +B[i][jplus1]-2*B[i][j]+B[i][jminus1]);
            ReB = 16.0-Aold*B[i][j];
            B[i][j]=B[i][j]+0.01*(ReB + DiB);
            if (B[i][j]<0) B[i][j]=0;
         }
      }
   }
   /**** Scale the results to 0-255 ******/
   {
      double high=0, low=255;
      for (i=0;i<maxI;i++)
        for (j=0;j<maxJ;j++) {
           if (A[i][j]>high)
             high=A[i][j];
           else if (A[i][j]<low)
             low=A[i][j];
        }
      for (i=0;i<maxI;i++)
        for (j=0;j<maxJ;j++)
        A[i][j]=(A[i][j]-low)*255.0/(high-low);
   }
   /**** Output to a PGM format file *****/
   {
      FILE* fich=fopen("output.pgm","wb");
      if (!fich) {
         printf ("File problems!\n");
         exit(1);
      }
      fprintf (fich,"P5 %d %d 255 ",maxJ,maxI);
      for (i=0;i<maxI;i++)
        for (j=0;j<maxJ;j++)
        fprintf(fich,"%c",(unsigned char)A[i][j]);
      fclose(fich);
   }
}