On the numerical evaluation of factorials II

A very basic but working big integer implementation

In this episode we will build a basic set of operations, namely addition,
subtraction and multiplication. Division is not needed but will be given,
too, as I am without mercy.There are several ways to build a big integer implementation, but it is basically always the same and varies in details only.

Every number is represented as an ordered list of “digits” to build the number {\displaystyle\sum_{n=0}^{k}a_n\beta^{n}} with n\in\mathbb{N} and \beta represents the base of the number-system. With \beta = 10 and a \in \mathbb{Z} we get the standard decimal system, used in a lot of places in this world and with \beta = 2^{32} or even \beta = 2^{64} we get a standard big integer system implemented in a lot of big integer…uhm…implementations.

To make things simpler we restrict further by \{a,\beta\}\in\mathbb{N} and to make things more interesting I have chosen a system based on \beta = 10^7. That is big enough and allows for simple and fast conversation to and from a string.

We work with the external file “already_done.js” or whatever you have renamed it to. The textarea in the HTML-File is for later use and experimenting with the things we already have. So open “already_done.js” with the editor of your choice and start copy ‘n pasting. It is assumed here, that you re well aware of all the traps and pitfalls of copy ‘n paste!

At first we need some very basic housekeeping. Some constants (some of the names might be familiar):

Math.BIGDECIMAL_BASE       = 10000000; //   largest x such that (10^x)^2 <2^53
Math.BIGDECIMAL_LOG10_BASE =        7; //   log_10(BIGDECIMAL_BASE)

Number.LLONG_MAX =  9007199254740992;  //   2^53
Number.INT_MAX   =  9007199254740992;  //   2^53
Number.INT_MIN   = -9007199254740992;  // -(2^53)
Number.BINT_MAX  = ~(1<<31);           //   data type for bool algebra

The prototype:

function Bigdecimal(){
  this.digits = new Array();
  this.sign   = 1;
  this.used   = 0; // number of digits in this.digits
}

The prototype is distinct from the functions here to allow for simple adding by typing another one into the HTML-textarea. A little spoiler:

Bigdecimal.KARATSUBA_CUTOFF = 20;  // *2
Bigdecimal.TOOM_COOK_CUTOFF = 40;  // * 3
Bigdecimal.TOOM_COOK4_CUTOFF = 40; // * 4
Bigdecimal.FFT_CUTOFF       = 205; // *5 = 1025, check is "<" than, so 1024

The exact numbers have been empirically evaluated with the Javascript-machine of the Opera-browser, so the mileage within other browsers may vary, but I doubt they will vary wildly.
The next one is to check if a given number (not a string) is an integer. The Number-prototype is a built-in and gets a bit enhancement here.

Number.prototype.isInteger = function(){
  // check if it is inside the limit
  if(   this < Number.INT_MIN
     || this > Number.LLONG_MAX){
     return false;
  }
  // check if it is an integer and not a float
  return (this - Math.floor(this) == 0)?true:false;
};

With that at hand we can start to make some Bigdecimals.

Number.prototype.toBigdecimal = function(){
  if(!this.isInteger()){return Number.NaN;}
  if(this > (-Math.BIGDECIMAL_BASE*Math.BIGDECIMAL_BASE) && this < (Math.BIGDECIMAL_BASE*Math.BIGDECIMAL_BASE) ){
    var ret = this.toString(10).toBigdecimal()
    return ret;
  }
};

What? Converting a number to a string and converting the resulting string to
a Bigdecimal? WTF?
Well, does not matter much here but feel free to change it. If you want to change it immediately:don’t, just read on for a little more.

We need a string parser, nevertheless.

String.prototype.parseInteger = function(){
  var ret = ["+","",10];
  var type = 10,c,d,start = 0;
  var s    = this.toLowerCase().replace(/[^0-9a-f-+]/g,"");
  var fsign = false;
  var ferror = false;
  
  if(s.length == 0) return Number.NaN;

  if(s.length > 1){  
  // check for type, set type
  c = this.charAt(0);
  d = this.charAt(1);
  // octal
  if(c == "0" && d != "x" && d != "b" ){
    start = 1;
    type = 8;
  }
  if(c == "0" && d == "x"){
    start = 2
    type = 16;
  }
  if(c == "0" && d == "b"){
    start = 2;
    type = 2;
  }
 
  if(s.length - start == 0) return Number.NaN;
  }
  else
    start = 0;
  for(var i=start;i<s.length;i++){
    // check for unicode omitted
    c = s.charAt(i);
    switch(c){
      case "+" :
      case "-" : if(!fsign )
                    ret[0] = c;
                 else
                   ferror = true;
                 fsign = true;
                 break;
      case "0" :
      case "1" : ret[1] += c;break;
      case "2" :
      case "3" :
      case "4" :
      case "5" :
      case "6" :
      case "7" : if(type >= 8) ret[1] += c;
                 else ferror = true;
                 break;
      case "8" :
      case "9" : if(type >= 10) ret[1] += c;
                 else ferror = true;
                 break;
      case "a" :
      case "b" :
      case "c" :
      case "d" :
      case "e" :
      case "f" : if(type >= 10) ret[1] += c;
                 else ferror = true;
                 break;
      default  : ferror = true;break;
    };
    // just stop on error
    if(ferror) break
  }
  if(ret[1].length == 0) return Number.NaN;
  ret[2]=type;
  return ret;
};

It just parses a string and returns either an array with the three elements sign,cleaned up string, type of number (binary, octal, decimal, or hexadecimal), or Number.NaN if an error occurred respectively. Only decimal gets used later, so don’t start to be too excited, you might get disappointed.

The next function makes use of the above.

String.prototype.toBigdecimal = function(){
  var sgn = 1;
  var val = new Array();
  var len = 0;
  
  var temp = this.parseInteger();
  var ret  = new Bigdecimal();  

  if(temp == Number.NaN){return Number.NaN;}
  // we accept decimals only. Extend if you want
  if(temp[2] != 10){return undefined;}

  sgn = (temp[0] == "+")?1:-1;
  temp = temp[1];
  len = temp.length;

  while(true){
    var limb = 0;
    if(len <= Math.BIGDECIMAL_LOG10_BASE){
      limb = parseInt(temp.slice(0,len),10);
      val.push(limb);
      break;
    }
    // sliding window, LSB -> MSB
    limb = parseInt(temp.slice(len-Math.BIGDECIMAL_LOG10_BASE,len),10);
    val.push(limb);
    len -= Math.BIGDECIMAL_LOG10_BASE;
  }
  ret.sign = sgn;
  ret.digits = val;
  ret.used = val.length;
  return ret;
};

It works by slicing off BIGDECIMAL_LOG10_BASE characters from the given string and converts it into a number smaller than BIGDECIMAL_BASE until the end of the given string.

To print a Bigdecimal we need to convert it to a string which is very easy here. One of the main reasons I chose a decimal base instead of the much faster binary base is the ease of the conversion. There exits a well done Javascript implementation of the algorithms needed for cryptography here. I tried it for factorials and it calculated the factorials bloody fast. Problem: the conversion to decimals took way too much time. Cryptography has not much use for
fast multiplication or fast division, so the conversion was done with the plain old schoolbook algorithm. It took a couple of minutes to calculate the factorial of a million but over an hour for the conversion to a decimal. So with the alternatives “change other peoples code” and “build my own from scratch” I chose, without a blink, the latter. Of course, I did!

Bigdecimal.prototype.toString = function(){
  var len = this.used-1;
  var str = "";
  if(this.isZero())return "0";
  for(var i=len;i>=0;i--){
    var temp = this.digits[i].toString();

    var tlen = temp.length;
    while(tlen++<Math.BIGDECIMAL_LOG10_BASE && i!=len){
      temp =  "0"+temp;
    }
    str += temp;
  }
  str = (this.sign < 1)?"-"+str:str;
  return str;
};

To keep our Bigdecimals neat and tidy we need this little gem from time to time:

Bigdecimal.prototype.clamp = function(){
  while(this.used > 0 && this.digits[this.used - 1] == 0){
    this.used--;
  }
  if(this.used == 0)
    this.sign = 1;
};

It just cuts off leading zeroes without cutting them off.

To make the whole thing a bit more comfortable to use is a method to make a full copy.

Bigdecimal.prototype.dup = function(){
  var copy = new Bigdecimal();
  if(this.used == 0) return copy;
  copy.digits = new Array(this.used);
  copy.sign   = this.sign;
  copy.used   = this.used;
  for(var i=0;i<this.used;i++)
    copy.digits[i] = this.digits[i];
  return copy;
};

A method to evaluate the number of decimal digits.

// returns the number of decimal digits
// does not check if the number is an integer but it really should!
Number.prototype.digits = function(){
  if(this == 0)return 1;
  return Math.floor(
          Math.log(Math.abs(this))
         *Math.LOG10E)+1;
};
Bigdecimal.prototype.digits = function(){
  t    = this.digits;
  tlen = this.used-1;
  var ret = tlen*Math.BIGDECIMAL_LOG10_BASE;
  ret += t[tlen].digits();
  return ret;
};

Checks if the Bigdecimal is the Number given (must be smaller than BIGDECIMAL_BASE) or in simpler words:it checks if the big number has only one digit and if that is true, checks if that single digit equals the Number given. I really forgot what I used it for. I think I never used it. I think it was just a silly idea. Uh, well…

Bigdecimal.prototype.isNumber = function(n){
  var t     = this.digits;
  var tlen  = this.used;
  
  if(tlen == 1 && t[0] == n)return true;
  return false;
};

Some shortcuts.We all love shortcuts, do we not? Stephen King had even made a short-story about shortcuts!

// checks if the Bigdecimal is zero
Bigdecimal.prototype.isZero = function(){
  if(this.used == 0 && this.sign == 1) return true;
  return false;
};
// checks if the Bigdecimal equals to one
Bigdecimal.prototype.isOne = function(){
  if(this.used == 1 && this.sign == 1 && this.digits[0] == 1) return true;
  return false;
};
// returns a bigdecimal set to zero
Bigdecimal.prototype.zero = function(){
  this.used == 0;
  this.digits = new Array();
  this.sign = 1;
  return true;
};
// returns the absolute value of the bigdecimal as a copy
Bigdecimal.prototype.abs = function(){
  var r = this.dup();
  r.sign   = 1;
  return r;
};
// returns a negated copy of the Bigdecimal
Bigdecimal.prototype.neg = function(){
  var r = this.dup();
  if(this.used == 0) return r;
  if(this.used == 1 && this.digits[0] == 0) return r;
  r.sign   = this.sign * -1;
  return r;
};

This is like Number.prototype.toBigdecimal without the restrictions. I told you to wait!
Feel free to try to find out if it works and if it works how it works.

Bigdecimal.prototype.set = function(n){
  if(!n.isInteger()){
    // just check the return for the type
    return this;
  }
  if(   n < Math.BIGDECIMAL_BASE 
     && n > -Math.BIGDECIMAL_BASE){
    if(n<0){
      n = -n;
      this.sign = -1;
    }
    if(this.used == 0)
      this.digits = new Array(1);
    this.digits[0] = n;
    this.used = 1;
    return this;
  }
  if(   n < Number.LLONG_MIN 
     && n > Number.LLONG_MAX){
    if(n<0){
      n = -n;
      this.sign = -1;
    }
    // overwrite
    this.digits = new Array();
    while(n > 0){
      limb = n%Math.BIGDECIMAL_BASE;
      this.digits.push(limb);
      // subtract what we got
      n -= limb;
      // shift right
      n /= Math.BIGDECIMAL_BASE;
      this.used++;
    }
    return this;
  }
  return this;
};
// Examples
Bigdecimal.ZERO = new Bigdecimal().set(0);
Bigdecimal.ONE  = new Bigdecimal().set(1);
Bigdecimal.TWO  = new Bigdecimal().set(2);

It would be nice, if we would be able to compare two Bigdecimals, would it not?. We do it in two stages here. Stage one compares the size of the Array holding the single limbs, that is we compare the number of digits and because it is simple, the sign, too. Stage two compares the actual values. Both stages get carefully mixed and spread to two different functions for no apparent reason but that I like it so.

Bigdecimal.prototype.cmp_mag = function(bdec){
  // Stage one
  if(this.used > bdec.used) return  1;
  if(this.used < bdec.used) return -1;
  // Stage two
  for(var i = this.used -1;i>=0;i--){
    if(this.digits[i] > bdec.digits[i]) return  1;
    if(this.digits[i] < bdec.digits[i]) return -1;
  }
  return 0;
};
// compare the sign and 
Bigdecimal.prototype.cmp = function(bdec){
  // Stage one
  if(this.sign > bdec.sign) return  1;
  if(this.sign < bdec.sign) return -1;
  // Stage two
  if(this.sign < 0)
    return bdec.cmp_mag(this);
  return this.cmp_mag(bdec);
};

Now we can compare like there is no tomorrow.

// a < b
Bigdecimal.prototype.lt = function(bdec){
  var r = this.cmp(bdec);
  if(r<0) return true;
  return false;
};
// a <= b
Bigdecimal.prototype.le = function(bdec){
  var r = this.cmp(bdec);
  if(r<=0) return true;
  return false;
};
// a == b
Bigdecimal.prototype.equal = function(bdec){
  var r = this.cmp(bdec);
  if(r==0) return true;
  return false;
};
// a == b
Bigdecimal.prototype.equals = function(bdec){
  return this.equal(bdec);
};
// a > b
Bigdecimal.prototype.gt = function(bdec){
  var r = this.cmp(bdec);
  if(r>0) return true;
  return false;
};
// a >= b
Bigdecimal.prototype.ge = function(bdec){
  var r = this.cmp(bdec);
  if(r>=0) return true;
  return false;
}

And make some shortcuts because all of us just love short…oh, we had that already? Pft!

Bigdecimal.prototype.max = function(bdec){
  return (this.ge(bdec))?this:bdec;
};
Bigdecimal.prototype.min = function(bdec){
  return (this.le(bdec))?this:bdec;
};

Some little extra helpers I admit I forgot:

Bigdecimal.prototype.isOdd = function(){
   return (this.digits[0]%2 == 1)?true:false; 
};
Bigdecimal.prototype.isEven = function(){
   return (this.digits[0]%2 == 0)?true:false; 
};

The actual algorithm follows the method used by Tom St. Denis in his library. He uses two stages, the first one being the rough work like handling signs, exceptions etc., the second stage does the actual work. The algorithm for addition consists of two stages, the inner work gets done by this little piece of code:

Bigdecimal.prototype.lowlevel_add = function(bdec){
  var t     = this.digits;
  var tlen  = this.used;

  var b     = bdec.digits;
  var blen  = bdec.used;

  var temp, carry = 0;
  var ret;
  if(tlen < blen)return bdec.lowlevel_add(this);
  ret = new Bigdecimal();
  ret.digits = new Array(tlen);
  for(var i=0;i<blen;i++){
    temp = carry;
    temp += t[i]+b[i];
    carry = Math.floor(temp / Math.BIGDECIMAL_BASE);
    ret.digits[i] = temp%Math.BIGDECIMAL_BASE;
  }
  if(tlen - blen != 0){
    for(var i=blen;i<tlen;i++){
      temp = carry;
      temp += t[i];
      carry = Math.floor(temp / Math.BIGDECIMAL_BASE);
      ret.digits[i] = temp%Math.BIGDECIMAL_BASE;
    }    
  }
  
  if(carry) ret.digits.push(carry)
  ret.used = ret.digits.length;
  ret.clamp();
  return ret;
};

As you can see: it is the algorithm you learned in school, nothing special here.
The other stage is more upper level and just takes care of the signs:

Bigdecimal.prototype.add = function(bdec){
  var tsgn = this.sign;
  var bsgn = bdec.sign;

  var ret,rsgn,cmp;

  if(tsgn < 0){   // -x
    if(bsgn < 0){ // -x + -y = -(x+y)
      rsgn = -1;
    }
    else{         // -x + +y = y - x
      return bdec.sub(this.abs());
    }
  }
  else{           // +x
    if(bsgn < 0){ // +x + -y = x - y
      return this.sub(bdec.abs());
    }
    else{         // +x + +y = x + y
      rsgn = 1;
    }
  }
  ret = this.lowlevel_add(bdec);
  ret.sign = rsgn;
  ret.clamp();
  return ret;  
};

Subtraction follows the same recipe:

Bigdecimal.prototype.lowlevel_sub = function(bdec){
  var t     = this.digits;
  var tlen  = this.used;    // max

  var b     = bdec.digits;
  var blen  = bdec.used;    // min

  var temp, borrow = 0;
  var ret = new Bigdecimal();
  ret.digits = new Array(tlen);
  
  for(var i=0;i<blen;i++){
    temp = t[i] - b[i] + borrow;
    borrow = Math.floor(temp/Math.BIGDECIMAL_BASE)
    ret.digits[i] = Math.floor(temp%Math.BIGDECIMAL_BASE) 
    if(ret.digits[i] < 0)ret.digits[i] += Math.BIGDECIMAL_BASE; 
  }

   for(var i=blen;i<tlen;i++){
     temp = t[i] + borrow;
     borrow = Math.floor(temp/Math.BIGDECIMAL_BASE)
     ret.digits[i] = Math.floor(temp%Math.BIGDECIMAL_BASE)
     if(ret.digits[i] < 0)ret.digits[i] += Math.BIGDECIMAL_BASE;
   }
  ret.used = ret.digits.length;
  ret.clamp();
  return ret;
};
Bigdecimal.prototype.sub = function(bdec){
  var tsign = this.sign;
  var bsign = bdec.sign;

  var ret;
  var cmp = this.cmp_mag(bdec)
  if(tsign != bsign){
    if(cmp > 0 && tsign > bsign){
      ret = this.lowlevel_add(bdec);
      ret.sign = tsign;
    }
    else{
      if(cmp < 0){
      ret = bdec.lowlevel_sub(this);
      ret.sign = (tsign < 0)?-1:1;
      }
      else{
       ret = this.lowlevel_sub(bdec);
      ret.sign = (tsign < 0)?-1:1;
      }
    }
  }
  else{
    var cmp = this.cmp_mag(bdec)
    if(cmp == 0){return Bigdecimal.ZERO;}
    else if(cmp >0){
      ret = this.lowlevel_sub(bdec);
      ret.sign = tsign;
    }
    else{
      ret = bdec.lowlevel_sub(this);
      ret.sign = (tsign < 0)?1:-1;
    }
  }
  ret.clamp();
  return ret;
  return ret;  
};

Shifting is something that can save a very large amount of computing when calculating factorials. But shifting in decimals is a bit more complicated than shifting in binary on a normal computer

Bigdecimal.prototype.lowlevel_dShiftLeft = function(n){
  var ret,more;
  if(n == 0) return this;
  if(n <  0) return this.lowlevel_dShiftRight(-n);

  more = new Array(n).memset(0);
  ret = this.dup();
  ret.digits = more.concat(ret.digits);
  ret.used = ret.digits.length;
  return ret;
};
Bigdecimal.prototype.dShiftLeft = function(n){
  return this.lowlevel_dShiftLeft(n);
};
Bigdecimal.prototype.lowlevel_dShiftRight = function(n){
  var tlen = this.used,ret,more;
  if(n == 0) return this;
  if(n <  0) return this.lowlevel_dShiftLeft(-n);
  if(n >= tlen) return new Bigdecimal();
  ret = this.dup();
  ret.digits = ret.digits.slice(n,this.length);
  ret.used = ret.digits.length;
  ret.clamp()
  return ret;
};
Bigdecimal.prototype.dShiftRight = function(n){
  return this.lowlevel_dShiftRight(n);
};
Bigdecimal.prototype.lowlevel_shiftLeftTen = function(n){
  var t     = this.digits;
  var tlen  = this.used;

  var temp, carry = 0;
  var ret;

  if(n == 0) return this.dup();
  if(n <  0) return this.lowlevel_shiftRightTen(-n);
  if(n > 7){
    return undefined;
  }
  ret = new Bigdecimal();
  for(var i=0;i<tlen;i++){
    temp = carry;
    temp += t[i]*Math.pow(10,n);
    carry = Math.floor(temp / Math.BIGDECIMAL_BASE);
    ret.digits[i] = temp%Math.BIGDECIMAL_BASE;
  }
  if(carry) ret.digits.push(carry)
  ret.used = ret.digits.length;
  ret.clamp();
  return ret;

};
Bigdecimal.prototype.lowlevel_shiftRightTen = function(n){
  var t     = this.digits;
  var tlen  = this.used;

  var temp, carry = 0;
  var ret;

  if(n == 0) return this.dup();
  if(n <  0) return this.lowlevel_shiftLeftTen(-n);
  if(n > 7){
    return undefined;
  }
  ret = new Bigdecimal();
  var dvr = Math.pow(10,n)
  for (var i = tlen - 1; i >= 0; i--) {
     temp = Math.BIGDECIMAL_BASE * carry + t[i];
     ret.digits[i] = Math.floor(temp / dvr );
     carry = Math.floor(temp - dvr * ret.digits[i]);
  }
  ret.used = ret.digits.length;
  ret.sign = this.sign;
  ret.clamp();
  return ret;
};
Bigdecimal.prototype.shiftLeftTen = function(n){
  var large,small;
  var ret;

  if(n == 0) return this.dup();
  if(n <  0) return this.shiftRightTen(-n);
  if(n > Number.LLONG_MAX){
    return undefined;
  }
  
  large = Math.floor(n/Math.BIGDECIMAL_LOG10_BASE);
  small = n%Math.BIGDECIMAL_LOG10_BASE;

  ret = this.lowlevel_dShiftLeft(large);
  ret = ret.lowlevel_shiftLeftTen(small);

  ret.used = ret.digits.length;
  ret.sign = this.sign;
  ret.clamp()
  return ret;
};
Bigdecimal.prototype.shiftRightTen = function(n){
  var large,small;
  var ret;

  if(n == 0) return this.dup();
  if(n <  0) return this.shiftLeftTen(-n);
  if(n > Number.LLONG_MAX){
    return undefined;
  }
  
  large = Math.floor(n/Math.BIGDECIMAL_LOG10_BASE);
  small = n%Math.BIGDECIMAL_LOG10_BASE;

  ret = this.lowlevel_dShiftRight(large);
  ret = ret.lowlevel_shiftRightTen(small);
  
  ret.used = ret.digits.length;
  ret.sign = this.sign;
  ret.clamp();
  return ret;
};

Now, that is a crude way to do it, isn’t it?

A little extra helper is needed here:

// sets every elemt of the Array to s, whatever s might be
Array.prototype.memset = function(s){
  var l = this.length;
  for(var i=0;i<l;i++){
    this[i] = s;
  }
  return this;
};

What now, hmm…
Ok, simple division.

// division with remainder by a normal Number.
Bigdecimal.prototype.divremD = function(n){
  var t     = this.digits;
  var tlen  = this.used;

  var temp, carry = 0;
  var ret;

  if(this.isZero()){
    return [Bigdecimal.ZERO,Bigdecimal.ZERO]
  }
  if(n == 0){
     Error.ERRNO = Error.TYPE_DIVZERO;
     return [Number.NaN,Number.NaN];
  }
  if(Math.abs(n) >= (Math.BIGDECIMAL_BASE*Math.BIGDECIMAL_BASE)){
    Error.ERRNO = Error.TYPE_DOM;
    return this
  }
  ret = new Bigdecimal();
  ret.digits = new Array(tlen);
  for (var i = tlen - 1; i >= 0; i--) {
    temp = Math.BIGDECIMAL_BASE * carry + t[i];
    ret.digits[i] = Math.floor(temp / n );
    carry = Math.floor (temp - n * ret[i]);
    //carry = temp % n;
  }
  ret.used = ret.digits.length;
  ret.sign = (this.sign == n.sign())?1:-1;
  ret.clamp();
  return [ret,carry.toBigdecimal()];
};

Bigdecimal.prototype.divD = function(n){
  return this.divremD(n)[0];
};
Bigdecimal.prototype.remD = function(n){
  return this.divremD(n)[1];
};

Now, that was simple, but dividing Bigdecimals is another cup of tea. This is the standard algorithm you can find described by many people, just ask Knuth (TAoCP II).
It is quite tricky if you are not very careful, so I sliced it up into smaller chunks to make it readable. Unlike the description in TAoCP.

  var trialdigit = function(r,d,k,m){
      var km = k+m;
      var r3,r3a,r3b,r3c,d2,d2a,d2b,quot;
      var b = Math.BIGDECIMAL_BASE;

      r3a  = (km     >= 0)?r.digits[km    ]:0;
      r3b  = (km - 1 >= 0)?r.digits[km - 1]:0;
      r3c  = (km - 2 >= 0)?r.digits[km - 2]:0;

      d2a  = (m  - 1 >= 0)?d.digits[ m - 1]:0;
      d2b  = (m  - 2 >= 0)?d.digits[ m - 2]:0;

      r3   = (r3a * b + r3b) * b + r3c;
      d2   =  d2a * b + d2b;
      quot = Math.floor(r3 / d2);
      return Math.min(quot, b-1);
  };
  var smaller = function(r,dq,k,m){
      var i,j;
      i = m; j = 0;
      while(i != j){
        if(r.digits[i + k] != dq.digits[i])
          i = j
        else
          i = i -1;
      }
      return (r.digits[i + k] < dq.digits[i]);
  };
  var difference = function(r,dq,k,m){
      var borrow, diff, i;
      var b = Math.BIGDECIMAL_BASE;
      borrow = 0;
      for(i = 0;i < m+1; i++){
        diff = r.digits[i + k] - dq.digits[i] - borrow + b;
        r.digits[i + k] = diff % b;
        borrow = 1 - Math.floor(diff / b);
      }
      return r;
  };
  var longdivrem = function(x,y,n,m){
      var d, dq, q, r;
      var f, k, qt;
      var b = Math.BIGDECIMAL_BASE;

      f = Math.floor(b / (y.digits[m - 1] +1));
      if(f != 1){
        x = x.mulD(f); // r
        y = y.mulD(f); // d
      }
      else{
       x.digits.push(0);
       x.used++;
      }
      r = x;
      d = y;

      q = new Bigdecimal();
      q.digits = new Array(n - m + 1).memset(0);
      q.used = q.digits.length;
      
      for(k = n - m;k >= 0;k--){
        qt = trialdigit(r,d,k,m);
        dq = d.mulD(qt);
        
        if(smaller(r,dq,k,m)){
          qt = qt - 1;
          dq = d.mulD(qt);
        }
        q.digits[k] = qt;
        r = difference(r,dq,k,m);
      }
      r = r.divD(f);
      r.clamp();
      return [q,r];
  };
  var x,y,m,n,ymsb1,ymsb2,sgn;
  var q,r;
  var b = Math.BIGDECIMAL_BASE;
  
  x = this.abs();
  y = bdec.abs();
  n = x.used;
  m = y.used;
  
  // handle simple cases first
  if(y.isZero()){
    Error.ERRNO = Error.TYPE_DIVZERO;
    return [Number.NaN,Number.NaN]; 
  }
  if(x.isZero()){
    return [Bigdecimal.ZERO,Bigdecimal.ZERO]; 
  }
  
  if(x.cmp_mag(y) < 0){
    if(this.sign < 0){
      q = Bigdecimal.ONE;
      if(bdec.sign > 0)
        q.sign = -1;
      r = y.sub(x);
      return [q,r];
    }
    return [Bigdecimal.ZERO,x];
  }
  if(n <= 2 && m == 1){
    if(n == 2)
      x = (x.digits[n-1] * b + x.sdigits[n-2]) * this.sign;
    else
      x = (x.digits[n-1] * this.sign);
    y = y.digits[m-1] * bdec.sign;
    q = x.divrem(y)
    r = q[1].toBigdecimal();
    q = q[0].toBigdecimal();
    return [q,r];
  }
  else{
    q = longdivrem(x,y,n,m);
    r = q[1];
    q = q[0];
    
    if(x.sign == y.sign){
      if(x.sign < 0){
        q.sign =  1;
        r.sign = -1;
      }
    }
    if(x.sign < 0){
      q.sign = -1;
      r.sign = -1;
    }
    if(y.sign < 0){
      q.sign = -1;
      r.sign =  1;
    }
  }
  
  return [q,r];
};
Bigdecimal.prototype.div = function(bdec){
  return this.divrem(bdec)[0];
};
Bigdecimal.prototype.rem = function(bdec){
  return this.divrem(bdec)[1];
}

Be careful! The algorithm above is optimized for our relatively small base, using a bigger base might fail!
And now, finally, the multiplication:

Bigdecimal.prototype.lowlevel_mulD = function(n){
  var t     = this.digits
  var tlen  = this.used;

  var carry = 0;
  var temp;
  var ret = new Bigdecimal();
  
  // signs are handled by the higlevel function
  n = Math.abs(n);
  
  if(n == 0) return ret.set(0);
  if(n == 1) return this.dup();
  
  ret.digits = new Array(tlen);

  for(var i=0;i<tlen;i++){
    temp  = t[i]*n + carry;
    carry = Math.floor(temp / Math.BIGDECIMAL_BASE);
    ret.digits[i] = temp%Math.BIGDECIMAL_BASE;
  }
  var i=0
  while(carry && i++<10){
    ret.digits.push(carry%Math.BIGDECIMAL_BASE)
    carry = Math.floor(carry / Math.BIGDECIMAL_BASE);
  }

  ret.used = ret.digits.length;
  ret.clamp();
  return ret;
};

Bigdecimal.prototype.mulD = function(n){
  var r = this.lowlevel_mulD(n);
  r.sign = this.sign * n.sign();
  return r;
};
Bigdecimal.prototype.lowlevel_mul = function(bdec){
  var tsign  = this.sign;
  var t      = this.digits;
  var tlen   = this.used;
  var bsign  = bdec.sign;
  var b      = bdec.digits;
  var blen   = bdec.used;
  var ret    = new Bigdecimal(), retv;

  //if(typeof bdec == 'number') return this.lowlevel_mulD(bdec);
  
  if(bdec.isZero()) return ret.set(0);
  if(bdec.isOne()) return this.dup();
  
  retv = new Array(tlen + blen);
  
  for(i=0;i<tlen;i++){
    carry = 0;
    for(j=0;j<blen;j++){
      // only first round, hence cheaper than memset
      // (a check for NULL is cheaper than a pointer incr + assignment )
      if(!retv[i+j] )retv[i+j]=0;
      var temp = retv[i+j]+(t[i]*b[j])+carry;
      carry = Math.floor(temp / Math.BIGDECIMAL_BASE);
      retv[i+j] = temp % Math.BIGDECIMAL_BASE;
    }
    retv[i+blen] = carry;
  }
  ret.digits = retv;
  ret.used = retv.length;
  ret.clamp();
  return ret;
};
Bigdecimal.prototype.mul = function(bdec){
  var r;
  if(typeof bdec == 'number' && Math.abs(bdec) < Number.LLONG_MAX)
    r =  this.mulD(bdec);
/*

  // a little teaser.

  else if(Math.min(this.used,bdec.used)< 2*Bigdecimal.KARATSUBA_CUTOFF)
    r = this.lowlevel_mul(bdec);
  else if(Math.min(this.used,bdec.used) < 3*Bigdecimal.TOOM_COOK_CUTOFF)
    r = this.karatsuba(bdec);
  else if(Math.min(this.used,bdec.used) < 4*Bigdecimal.TOOM_COOK4_CUTOFF) 
    r = this.toomCookMul(bdec); 
  //else if(Math.min(this.used,bdec.used) < 5*Bigdecimal.FFT_CUTOFF) 
    //r = this.toomCookMul4(bdec); 
  else {r = this.mulWithPoorFFT(bdec,1024,4);}
*/
  else 
    r = this.lowlevel_mul(bdec);
  r.sign = this.sign * bdec.sign;
  return r; 
};

Here it is the standard schoolbook algorithm without the optimizations.

Squaring is faster than multiplying two similar numbers the old way.

Bigdecimal.prototype.square = function(){
    return this.squareNormal();
  /*
  if(this.used < 2*Bigdecimal.KARATSUBA_CUTOFF){
    return this.squareNormal();   
  } 
   Aaaand another teaser. I'm so mean.

  else if(this.used < 3*Bigdecimal.TOOM_COOK_CUTOFF){
    return this.karatsubaSquare();
  }
  else if(this.used < 4*Bigdecimal.TOOM_COOK4_CUTOFF){
    return this.toomCookSquare();
  }
  else if(this.used < 5*Bigdecimal.TOOM_COOK5_CUTOFF){
    return this.toomCook4Square();
  }
  else if(this.used < 5*Bigdecimal.FFT_CUTOFF){
    return this.toomCook5Square();
  else
  return this.squareFFT();
  */
};
Bigdecimal.prototype.squareNormal = function(){
  var t    = this.digits;
  var tlen = this.used;
  var ret  = new Bigdecimal();
  var retv;
  var uv = 0;
  var carry = 0;
  var i,j
  ret.digits = new Array(2*tlen+1).memset(0);
  ret.used = 2*tlen+1;
  retv = ret.digits;
  
  for(i=0;i<tlen;i++){
    uv       = retv[2*i]+ t[i]*t[i];
    retv[2*i] = uv%Math.BIGDECIMAL_BASE;             // v
    carry    = Math.floor(uv/Math.BIGDECIMAL_BASE);  // u
    for(j=i+1;j<tlen;j++){
      uv       = retv[i+j]+2*t[j]*t[i]+carry;
      retv[i+j] = uv%Math.BIGDECIMAL_BASE;            // v
      carry    = Math.floor(uv/Math.BIGDECIMAL_BASE); // u
    }
    retv[i+j] = carry;                         // u
  }
  ret.clamp();
  return ret;
};

Just try it out and tell me if you find any bugs.

Oh, and the most primitive and slowest way to calculate a factorial is this:

Bigdecimal.factorialN = function(n){
  var one = Bigdecimal.ONE.dup();
  for(var i=2;i<=n;i++){
    one = one.mulD(i);
  }
  return one;
};
It is not a prototype, just use it like this:

PRINT(Bigdecimal.factorialN(10000))

Admittedly, it is not the best style to post such an amount of sourcecode with that little explanation in a little blog and not at Sourceforge et. al. I actually even have an account at Sourceforge. It hosts a collection of Javascripts called “Pragmathic Mathematical Service” (no reason to hit me for that, it was not intentional. Really!). It will land there when two preconditions are fulfilled: a) I get the time to wrap all of it up in a neat little package (there are more new things than this little bignum-lib) and b) a find my Sourceforge password.
I sincerely hope that they haven’t changed too much in the meantime.

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