1000 factorial

Well, someone wrote that they could write a program to do 1000! in one line of Haskell. I said that it would take less than 100 lines of C to do the same. I wrote the program, and it is less than 100 lines of C:

#include <stdio.h>
#include <string.h>

#define LENGTH 10000

int max (int a, int b) { return (a > b) ? a : b; }
int min (int a, int b) { return (b > a) ? a : b; }

int length (unsigned char * a)
 {
   int i;
   i = LENGTH - 1;
   while ((a[i] == 0) && (i >= 0)) i--;
   return i + 1;
 }

void add (unsigned char * a, unsigned char * b)
 {
   int i, c, C;
   c = min(max(length(a), length(b)) + 1, LENGTH);
   for (i = 0, C = 0; i < c; i++)
    {
      a[i] += b[i] + C;
      C = a[i] > 99;
      if (C) a[i] -= 100;
    }
 }

   unsigned char x [LENGTH], y [LENGTH];
void mul (unsigned char * a, unsigned char * b)
 {
   int i, j, c, d, C;
   memset(y, '\0', LENGTH);
   c = min(length(a) + length(b) + 1, LENGTH);
   d = length(b);
   for (i = 0; i < d; i++)
    {
      memset(x, '\0', LENGTH);
      for (j = 0, C = 0; (i + j) < c; j++)
       {
         x[i + j] = (a[j] * b[i] + C) % 100;
         C = (a[j] * b[i] + C) / 100;
       }
      add(y, x);
    }
   memcpy(a, y, LENGTH);
 }

void numset (unsigned char * a, unsigned long b)
 {
   int i;
   for (i = 0; i < LENGTH; i++)
    {
      a[i] = b % 100;
      b /= 100;
    }
 }

int tz (unsigned char * a)
 {
   int i;
   for (i = 0; i < LENGTH; i++) if (a[i] != 0) return 0;
   return 1;
 }

   unsigned char minus1 [LENGTH], prod [LENGTH];
void fact (unsigned char * a)
 {
   numset(prod, 1);
   while (!tz(a))
    {
      mul(prod, a);
      add(a, minus1);
    }
   memcpy(a, prod, LENGTH);
 }

   unsigned char arg_res [LENGTH];
int main (void)
 {
   unsigned long a;
   int i;
   memset (minus1, 99, LENGTH);
   printf("Input the number to take the factorial of: ");
   scanf("%lu", &a);
   numset(arg_res, a);
   fact(arg_res);
   i = LENGTH - 1;
   while (arg_res[i] == 0) i--;
   while (i >= 0)
    {
      printf("%02d", (int) arg_res[i]);
      i--;
    }
   putchar('\n');
   return 0;
 }

>>143
If you're going to bitch about inefficiency, you really should do it right.

product([], 1).
product([H|T], X) :- product(T, Y), X is H * Y.

sieve(N, [2|PS]) :-       % PS is list of odd primes up to N
    retractall(mult(_)),
    sieve_O(3,N,PS).

sieve_O(I,N,PS) :-        % sieve odds from I up to N to get PS
    I =< N, !, I1 is I+2,
    (   mult(I) -> sieve_O(I1,N,PS)
    ;   (   I =< N / I ->
            ISq is I*I, DI  is 2*I, add_mults(DI,ISq,N)
        ;   true ),
        PS = [I|T],
        sieve_O(I1,N,T)
    ).
sieve_O(I,N,[]) :- I > N.

add_mults(DI,I,N) :-
    I =< N, !,
    ( mult(I) -> true ; assert(mult(I)) ),
    I1 is I+DI,
    add_mults(DI,I1,N).
add_mults(_,I,N) :- I > N.

prime_power_of_s(_, 0, A, S) :- S is A.
prime_power_of_s(P, N, A, S) :- Q is div(N, P), R is mod(N, P), B is A + R, prime_power_of_s(P, Q, B, S).

prime_power_of(N, P, R) :- prime_power_of_s(P, N, 0, S), R is div(N - S, P - 1).

factorized_factorial(N, Ps, Rs) :- sieve(N, Ps), maplist(prime_power_of(N), Ps, Rs).

exponentiate_factors([], [], []).
exponentiate_factors([P|Ps], [R|Rs], [N|Ns]) :- N is P ** R, exponentiate_factors(Ps, Rs, Ns).

factorial(N, R) :- factorized_factorial(N, Ps, Rs), exponentiate_factors(Ps, Rs, Ns), product(Ns, R).