On the numerical evaluation of factorials III

Sieving Primes

The fastest know algorithm to compute factorials uses the prime factors of the factorial. This can be done quite easily but needs a list of primes to do so. One of the simplest way to do it—and a quite fast way—is the old algorithm invented by a greek guy named Eratosthenes of Cyrene. OK, he was greek, so he wrote his name most probably Ἐρατοσθένης.

We are not really interrested in the literacy of ancient greeks but in their inventions. and Eratosthenes invented a lot One of those findings was a prime sieve and it got even named after him!

Implementation is simple and straightforward—if you have a Bitset which is a long string of bits that are adressable directly. Let’s write one.

The bitset itself is an array of integers; here 32 bit integers are assumed.

function BitSet(length){
  var len = Math.ceil(length/32);
  if(len > Number.INT_MAX){
    return Number.POSITIVE_INFINITY;
  }
  // this is a dirty trick to persuade the Javascript
  // parser to use a 32 bit integer. This trick may
  // fail, please test.
  // It matters not if a larger integer it is but an
  // integer it must be! [Yoda, "Aphorisms IV", 13, 181]
  this.Bitstring = new Array(len).memset(~(~(0&0xFFFFFFFF)));
  this.Length = length;
}

Now that we have that, we can write some functions to read and write individual bits and see if we can get some shortcuts, too.

BitSet.prototype.set = function(where){
  if( where > this.Length || where < 0)
    return false;
  this.Bitstring[where>>>5] |= (1 << (31-( where &31)));
  return true;
};
BitSet.prototype.clear = function(where){
  if( where > this.Length || where < 0)
    return false;
  this.Bitstring[where>>>5] &= ~((1 << (31-( where &31))));
  return true;
};
BitSet.prototype.get = function(where){
  if(where > this.Length || where < 0)
    return -1;
 return ((this.Bitstring[where>>>5] >>> ((31-(where&31)))) & 1);
};
BitSet.prototype.getset = function(where){
  var val = this.get(where);
  this.set(where);
  return val;
};
BitSet.prototype.getclear = function(where){
  var val = this.get(where);
  this.clear(where);
  return val;
};
BitSet.prototype.getflip = function(where){
  var val = this.get(where);
  this.flip(where);
  return val;
};
BitSet.prototype.nextset = function(n){
  var N = this.Length;
  while(n<N && !this.get(n))n++;
  if(n == N && !this.get(n)) return -1;
  return n;
};
BitSet.prototype.prevset = function(n){
  while(n>=0 && !this.get(n))n--;
  if(n == 0 && !this.get(n)) return -1;
  return n;
};
BitSet.prototype.nextclear = function(n){
  var N = this.Length;
  while(n<N && this.get(n))n++;
  if(n == N && this.get(n)) return -1;
  return n;
};
BitSet.prototype.prevclear = function(n){
  while(n>=0 && this.get(n))n--;
  if(n == 0 && this.get(n)) return -1;
  return n;
};
BitSet.prototype.flip = function(where){
  if(where > this.Length  || where < 0)
    return false;
  this.Bitstring[where>>>5] ^= (1 << (31-( where &31)));
  return true;
};
BitSet.prototype.clearall = function(){
  this.Bitstring.memset(~(~(0&0xFFFFFFFF)));
};
BitSet.prototype.setall = function(){
  this.Bitstring.memset((~(0&0xFFFFFFFF)));
};

That number of functions is way more than needed for a primesieve, but as we already have a working bitset we can go a bit further and do what follows logically.

BitSet.prototype.and = function(bs){
  var tlen = this.Bitstring.length;
  var t    = this.Bitstring;
  var b    = bs.Bitstring;
  if(this.Length != bs.Length)
    return false;
  for(var i=0;i<tlen;i++){
    t[i] &= b[i];
  }
  return true;
};
BitSet.prototype.or  = function(bs){
  var tlen = this.Bitstring.length;
  var t    = this.Bitstring;
  var b    = bs.Bitstring;
  if(this.Length != bs.Length)
    return false;
  for(var i=0;i<tlen;i++){
    t[i] |= b[i];
  }
  return true;
};
BitSet.prototype.xor = function(bs){
  var tlen = this.Bitstring.length;
  var t    = this.Bitstring;
  var b    = bs.Bitstring;
  if(this.Length != bs.Length)
    return false;
  for(var i=0;i<tlen;i++){
    t[i] ^= b[i];
  }
  return true;
};
BitSet.prototype.not = function(){
  var tlen = this.Bitstring.length;
  var t    = this.Bitstring;
  for(var i=0;i<tlen;i++){
    t[i] = ~t[i];
  }
  return true;
};
BitSet.prototype.andNot = function(bs){
  var tlen = this.Bitstring.length;
  var t    = this.Bitstring;
  var b    = bs.Bitstring;
  if(this.Length != bs.Length)
    return false;
  for(var i=0;i<tlen;i++){
    t[i] &= ~b[i];
  }
  return true;
};
BitSet.prototype.orNot = function(bs){
  var tlen = this.Bitstring.length;
  var t    = this.Bitstring;
  var b    = bs.Bitstring;
  if(this.Length != bs.Length)
    return false;
  for(var i=0;i<tlen;i++){
    t[i] |= ~b[i];
  }
  return true;
};

I have written some functions to make a printable string out of the bitset and read that back in but it depends on too much extra functions not yet published so I will reserve these for my project at Sourceforge. They will be fully commented there. I hope.

We have a working implementation of a bitset now, we can start with the primesieve. We don’t do any optimizations for space (you can check if the number is even and greater than 2 beforehand and just put the odd primes in the bitset) just the plain sieve as invented by Eratosthenes but I’m pretty sure Eratosthenes has done that optimization too. He was definitely not the dumbest person of his time.

Some hints of how to do it and much more can be found in Jörg Arndt’s book Matters Computational.

// we make use of the namespace Pragmath here.
Pragmath.primesieve = null;
Pragmath.primelimit = 1000000;// Choose wisely, but this thing is quite fast
Pragmath.sieve = function(n){
  var k,r,i,j;
  var l = n || Pragmath.primelimit;
  l = l+1;
  Pragmath.primesieve = new BitSet(l);
  Pragmath.primesieve.setall();
  Pragmath.primesieve.clear(0);
  Pragmath.primesieve.clear(1);

  for(k = 4;k<l;k+=2)Pragmath.primesieve.clear(k);

  r = Math.floor(Math.sqrt(l));
  k =0;
  while(k <l){
    k = Pragmath.primesieve.nextset(k+1);
    if(k > r ) break;
    for(j = k*k;j<l;j+=2*k)
      Pragmath.primesieve.clear(j);
  }
  Pragmath.primelimit = l;
};

The names of the functions to use the sieve might seem familiar to some of you.

Math.isSmallPrime = function(n){
  if(!n.isInteger() || n > Pragmath.primelimit) return undefined;
  if(Pragmath.primesieve.get(n) == 1) return true;
  return false;
};
Math.nextPrime = function(n){
  if(!n.isInteger() || n+1 > Pragmath.primelimit) return undefined;
  return Pragmath.primesieve.nextset(n);
};
Math.precPrime = function(n){
  if(!n.isInteger() || n < 2) return undefined;
  return Pragmath.primesieve.prevset(n);
};
Math.primePi = function(n){
  var k=0;var ct =0;
  if(!n.isInteger() || n+1 > Pragmath.primelimit) return undefined;  
  while(k<n ){
    k = Pragmath.primesieve.nextset(k+1);
    if(k>n || k <0) break;
    ct++;
  }
  return ct;
};
Math.primeRange = function(low,high){
  var down = 0, up = 0, ret = new Array();
  // TODO!
  if(!low.isInteger() || !high.isInteger() 
       || low < 0 || high > Pragmath.primelimit
       || low > high)
       return undefined;
  down = Pragmath.primesieve.nextset(low);
  up   = Pragmath.primesieve.prevset(high);
  ret.push(down);
  if(down == up) return ret;
  while(down<up){
    down = Pragmath.primesieve.nextset(down+1);
    if(down>high || down <0) break;
    ret.push(down);
  }
  return ret;
};
Math.primorial = function(n){
  if(!n.isInteger() || n < 2 || n+1 > Pragmath.primelimit) return undefined;

  if(n <= 42){
    var ret = 1;
    for(i=2;i<n;i = Pragmath.primesieve.nextset(i+1)){
      ret = ret * i;
    }
    return ret;
  }
  else{
    var ret = 1;
    for(i=2;i<n;i = Pragmath.primesieve.nextset(i+1)){
      ret = ret.mul(i);
    }
    return ret;
  }

};
Math.primorialRange = function(low,high){
  var primearray = Math.primeRange(low,high);
  var len,ret = 1;
  if(primearray == undefined)return undefined;
  if(primearray.length == 0) return undefined;
  if(primearray.length == 1) return primearray[0];

  len = primearray.length;
  
  for(var i = 0;i<len;i++){
    ret = ret.mul(primearray[i]);
  }
  return ret;
};
Math.primes = function(n){
  var ret,k=0;
  if(!n.isInteger() || n < 2 || n+1 > Pragmath.primelimit) return undefined;
  ret = new Array();
  while(k<n ){
    k = Pragmath.primesieve.nextset(k+1);
    if(k>n || k <0) break;
    ret.push(k);
  }
  return ret;
};

The sieve is quite fast, so just adding a

Pragmath.sieve();

at the end will not extend the loading time by more than a couple of seconds for the default sieve limit of one million. And computing the factorial of one million in Javascript is probably more than you might want to wait for:
1,000,000! is about 8.2639e5,565,708
(10^9)! ~ 9.9046 e 8,565,705,522
(10^12)! ~ 1.4037 e 11,565,705,518,103
(10^100)! ~ 1.62940 e 995,657,055,180,967,481,723,488,710,810,833,949,177,056,029,941,963,
334,338,855,462,168,341,353,507,911,292,252,707,750,506,615,682,467
The next large number would be (10^(100^100)) that is about…well, a bit more. The available tools here (PARI/GP overflows and wolframalpha.com ignores my input…uh, it works with a=100^100;lngamma(10^a) but harvests overflow. Yepp, another five minutes I’ll never get back) are insufficient, so it will be the g’d ol’ slideruler style to get a good approximation of the number.

Much time went by since the last time I did such things so I sincerely hope that I do not mess it up more than absolutely necessary.

Starting with n! = \prod_{k=1}^n k and logarithmize both sides we get \log n! = \sum_{k=1}^n \log k which can be approximated for large n with
\displaystyle\int_1^n \log x \,\mathrm{d}x
We can solve that integral with the help of the method called integration by parts. You might remember it from calculus.
\displaystyle\int u \,\mathrm{d}v = uv - \int v \,\mathrm{d}u
Let \mathrm{d}v = \mathrm{d}x and u = log x . That gives \mathrm{d}u = \tfrac{\mathrm{d}x}{x} and v = x , so
\displaystyle\int_0^n \log x \, \mathrm{d}x = x \log x \Bigl\lvert_0^n - \int_0^n x  \frac{\mathrm{d}x}{x}
We remember, with a bit of an effort but nevertheless we do remember that x \log x is zero at zero and, believe it or not \tfrac{x}{x} = 1 with x \not=0 , which leads finally to
\displaystyle\int_0^n \log x \, \mathrm{d}x = n \log n - \int_0^n \mathrm{d}x  = n \log n - n
which is what we call the Stirling’s Approximation
As we are dealing with a very large number that is still not very comfortable we drop the subtraction of n , logarithmize both sides once again and get
\displaystyle\log \log  n! \approx \log n + \log\log n
I think we can work with that. By using a base ten logarithm and plugging in we get
\displaystyle\log_{10} {100^{100}}! \approx {100^{100}} + \log_{10}{100^{100}}
\log_{10}{100^{100}} is just 200 and pales in contrast to the 200 decimal digits large 100^{100} so 10^{100^{100}}! is a number with more than 10^{100^{100}} digits which is even larger than Googolplex

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s