Random Thoughts Series 2 - Factorial Algorithm Thoughts

by guanlan

Factorial again, this is an old topic, right? Without even thinking, a simple recursion will do!

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

Or you think recursion isn’t as efficient as tail recursion , so with a stroke of your pen:

long fact_iter(long product, long counter, long maxcount)
{
return (counter > maxcount) ? product : fact_iter(product*counter, counter+1, maxcount);
}

long factorial(long n)
{
return fact_iter(1, 1, n);
}

Or maybe you’ve read in “Code Complete”: “If a programmer working for me used recursion to calculate factorial, I’d rather replace them.”
Using recursion to calculate factorial is slow, unpredictable in terms of memory usage during runtime, and difficult to understand. So you change the recursion to a loop statement.

int factorial(int n)
{
int result = 1;
for(int i = 2 ; i <= n; i++)
{
result = result * i;
}
return result;
}

When you write this code, don’t you feel something is missing?

Testing on my 32-bit environment, it overflows when calculating 33!, so you say, well that’s because int is too small, so you switch to long double. Test it again - what is this? It still doesn’t work when the numbers get a bit larger.

Then let’s use linked lists or arrays instead. Linked lists are too slow, so let’s use arrays.

int factorial2(int n,int a[])
{
int carry;
int digit = 1;
a[0] = 1;
int temp;
for(int i = 2; i <= n; ++i)
{
for(int j = 1, carry = 0; j <= digit; ++j)
{
temp = a[j-1] * i + carry;
a[j-1] = temp % 10;
carry = temp / 10;
}
while(carry)
{
a[++digit-1] = carry % 10;
carry /= 10;
}
}
return –digit;
}

This algorithm simulates the manual calculation process, saves the result in array a, and returns the number of digits in the result.

Are you feeling light-headed at this point? Please hold on for a moment.

What if I need a scientific notation expression for a large number above 100,000? Or ask you, what is the third digit from the left of N! at the 100,000 level? Uh, isn’t this a mathematician’s job? Cheer up and challenge yourself! What real programmers need is this spirit of getting to the bottom of things.
Let’s try mathematical analysis methods. James Stirling, a Scottish mathematician, gave this limit formula more than 280 years ago:

This formula can calculate the approximate value of n! at extremely fast speed, and can also be used to infinitely approach the exact result. Detailed introduction and proof process is here or here .

Stirling Series Formula

The following code calculates the scientific notation representation of large number N!

struct bigNum
{
double n; //mantissa
int e; //exponent
};
void factorial3(struct bigNum p,int n)
{
double logx,s,item;//s: sum of series, item: each term of series
int i;
logx = n log10((double)n/E);
p->e = (int)(logx); p->n= pow(10.0, logx-p->e);
p->n = sqrt( 2 PI* (double)n);
for (item=1.0,s=0.0,i=0;i<sizeof(a1)/sizeof(double);i++)
{
s += item * a1[i];
item /= (double)n;
}
p->n *=s;
}

The following is the asymptotic expansion of the logarithm of factorial

void factorial3b(struct bigNum p,int n)
{
double logR;
double s,item;
int i;
logR=0.5log(2.0*PI)+((double)n+0.5)log(n)-(double)n;

for (item=1/(double)n,s=0.0,i=0;i<sizeof(a2)/sizeof(double);i++)
{
s+= item * a2[i];
item /= (double)(n) (double)n;
}
logR+=s;
p->e = (int)( logR / log(10));//change of base formula
p->n = pow(10.00, logR/log(10) - p->e);
}

Calculating the number of digits in a factorial is also particularly simple:

double getFactorialLength(int n)
{
return (n * log(double(n)) - n + 0.5 * log(2.0 * n * PI )) / log(10.0)+1;
}

This calculates an approximate number of digits, or you can improve it by using the ceil function to find the smallest integer not less than the given real number.

int getFactorialLength(int n)
{
if( n == 1 )
return 1;
else
return (int)ceil((Nlog(N)-N+log(2N*PI)/2)/log(10));
}

At this point, you can’t help but exclaim: the most brilliant, most essential, most fundamental thing in computer science is still mathematics!

Mandelbrot set boundary
Finally, let me end this essay with Russell’s words:
Mathematics, rightly viewed, possesses not only truth, but supreme beauty — a beauty cold and austere, like that of sculpture, without appeal to any part of our weaker nature, without the gorgeous trappings of painting or music, yet sublimely pure, and capable of a stern perfection such as only the greatest art can show. The true spirit of delight, the exaltation, the sense of being more than Man, which is the touchstone of the highest excellence, is to be found in mathematics as surely as poetry.

References
1.Tom M. Apostol. “Mathematical Analysis”
2.Steve McConnell. “CODE COMPLETE, Second Edition”

3.http://en.wikipedia.org/wiki/Stirling_approximation#History
4.http://mathworld.wolfram.com/StirlingsApproximation.html
5.http://zh.straightworldbank.com/wiki/%E6%96%AF%E7%89%B9%E6%9E%97%E5%85%AC%E5%BC%8F

Original Chinese version: http://www.cppblog.com/xguru/archive/2009/12/30/104344.html