includes/clientside/static/libbigint.js
changeset 436 242353360e37
equal deleted inserted replaced
435:a434d60e525d 436:242353360e37
       
     1 ////////////////////////////////////////////////////////////////////////////////////////
       
     2 // Big Integer Library v. 5.1
       
     3 // Created 2000, last modified 2007
       
     4 // Leemon Baird
       
     5 // www.leemon.com
       
     6 //
       
     7 // Version history:
       
     8 //
       
     9 // v 5.1  8 Oct 2007 
       
    10 //   - renamed inverseModInt_ to inverseModInt since it doesn't change its parameters
       
    11 //   - added functions GCD and randBigInt, which call GCD_ and randBigInt_
       
    12 //   - fixed a bug found by Rob Visser (see comment with his name below)
       
    13 //   - improved comments
       
    14 //
       
    15 // This file is public domain.   You can use it for any purpose without restriction.
       
    16 // I do not guarantee that it is correct, so use it at your own risk.  If you use 
       
    17 // it for something interesting, I'd appreciate hearing about it.  If you find 
       
    18 // any bugs or make any improvements, I'd appreciate hearing about those too.
       
    19 // It would also be nice if my name and address were left in the comments.
       
    20 // But none of that is required.
       
    21 //
       
    22 // This code defines a bigInt library for arbitrary-precision integers.
       
    23 // A bigInt is an array of integers storing the value in chunks of bpe bits, 
       
    24 // little endian (buff[0] is the least significant word).
       
    25 // Negative bigInts are stored two's complement.
       
    26 // Some functions assume their parameters have at least one leading zero element.
       
    27 // Functions with an underscore at the end of the name have unpredictable behavior in case of overflow, 
       
    28 // so the caller must make sure the arrays must be big enough to hold the answer.
       
    29 // For each function where a parameter is modified, that same 
       
    30 // variable must not be used as another argument too.
       
    31 // So, you cannot square x by doing multMod_(x,x,n).  
       
    32 // You must use squareMod_(x,n) instead, or do y=dup(x); multMod_(x,y,n).
       
    33 //
       
    34 // These functions are designed to avoid frequent dynamic memory allocation in the inner loop.
       
    35 // For most functions, if it needs a BigInt as a local variable it will actually use
       
    36 // a global, and will only allocate to it only when it's not the right size.  This ensures
       
    37 // that when a function is called repeatedly with same-sized parameters, it only allocates
       
    38 // memory on the first call.
       
    39 //
       
    40 // Note that for cryptographic purposes, the calls to Math.random() must 
       
    41 // be replaced with calls to a better pseudorandom number generator.
       
    42 //
       
    43 // In the following, "bigInt" means a bigInt with at least one leading zero element,
       
    44 // and "integer" means a nonnegative integer less than radix.  In some cases, integer 
       
    45 // can be negative.  Negative bigInts are 2s complement.
       
    46 // 
       
    47 // The following functions do not modify their inputs.
       
    48 // Those returning a bigInt, string, or Array will dynamically allocate memory for that value.
       
    49 // Those returning a boolean will return the integer 0 (false) or 1 (true).
       
    50 // Those returning boolean or int will not allocate memory except possibly on the first time they're called with a given parameter size.
       
    51 // 
       
    52 // bigInt  add(x,y)               //return (x+y) for bigInts x and y.  
       
    53 // bigInt  addInt(x,n)            //return (x+n) where x is a bigInt and n is an integer.
       
    54 // string  bigInt2str(x,base)     //return a string form of bigInt x in a given base, with 2 <= base <= 95
       
    55 // int     bitSize(x)             //return how many bits long the bigInt x is, not counting leading zeros
       
    56 // bigInt  dup(x)                 //return a copy of bigInt x
       
    57 // boolean equals(x,y)            //is the bigInt x equal to the bigint y?
       
    58 // boolean equalsInt(x,y)         //is bigint x equal to integer y?
       
    59 // bigInt  expand(x,n)            //return a copy of x with at least n elements, adding leading zeros if needed
       
    60 // Array   findPrimes(n)          //return array of all primes less than integer n
       
    61 // bigInt  GCD(x,y)               //return greatest common divisor of bigInts x and y (each with same number of elements).
       
    62 // boolean greater(x,y)           //is x>y?  (x and y are nonnegative bigInts)
       
    63 // boolean greaterShift(x,y,shift)//is (x <<(shift*bpe)) > y?
       
    64 // bigInt  int2bigInt(t,n,m)      //return a bigInt equal to integer t, with at least n bits and m array elements
       
    65 // bigInt  inverseMod(x,n)        //return (x**(-1) mod n) for bigInts x and n.  If no inverse exists, it returns null
       
    66 // int     inverseModInt(x,n)     //return x**(-1) mod n, for integers x and n.  Return 0 if there is no inverse
       
    67 // boolean isZero(x)              //is the bigInt x equal to zero?
       
    68 // boolean millerRabin(x,b)       //does one round of Miller-Rabin base integer b say that bigInt x is possibly prime (as opposed to definitely composite)?
       
    69 // bigInt  mod(x,n)               //return a new bigInt equal to (x mod n) for bigInts x and n.
       
    70 // int     modInt(x,n)            //return x mod n for bigInt x and integer n.
       
    71 // bigInt  mult(x,y)              //return x*y for bigInts x and y. This is faster when y<x.
       
    72 // bigInt  multMod(x,y,n)         //return (x*y mod n) for bigInts x,y,n.  For greater speed, let y<x.
       
    73 // boolean negative(x)            //is bigInt x negative?
       
    74 // bigInt  powMod(x,y,n)          //return (x**y mod n) where x,y,n are bigInts and ** is exponentiation.  0**0=1. Faster for odd n.
       
    75 // bigInt  randBigInt(n,s)        //return an n-bit random BigInt (n>=1).  If s=1, then the most significant of those n bits is set to 1.
       
    76 // bigInt  randTruePrime(k)       //return a new, random, k-bit, true prime bigInt using Maurer's algorithm.
       
    77 // bigInt  str2bigInt(s,b,n,m)    //return a bigInt for number represented in string s in base b with at least n bits and m array elements
       
    78 // bigInt  sub(x,y)               //return (x-y) for bigInts x and y.  Negative answers will be 2s complement
       
    79 // bigInt  bigint_trim(x,k)              //return a copy of x with exactly k leading zero elements
       
    80 //
       
    81 //
       
    82 // The following functions each have a non-underscored version, which most users should call instead.
       
    83 // These functions each write to a single parameter, and the caller is responsible for ensuring the array 
       
    84 // passed in is large enough to hold the result. 
       
    85 //
       
    86 // void    addInt_(x,n)          //do x=x+n where x is a bigInt and n is an integer
       
    87 // void    add_(x,y)             //do x=x+y for bigInts x and y
       
    88 // void    copy_(x,y)            //do x=y on bigInts x and y
       
    89 // void    copyInt_(x,n)         //do x=n on bigInt x and integer n
       
    90 // void    GCD_(x,y)             //set x to the greatest common divisor of bigInts x and y, (y is destroyed).  (This never overflows its array).
       
    91 // boolean inverseMod_(x,n)      //do x=x**(-1) mod n, for bigInts x and n. Returns 1 (0) if inverse does (doesn't) exist
       
    92 // void    mod_(x,n)             //do x=x mod n for bigInts x and n. (This never overflows its array).
       
    93 // void    mult_(x,y)            //do x=x*y for bigInts x and y.
       
    94 // void    multMod_(x,y,n)       //do x=x*y  mod n for bigInts x,y,n.
       
    95 // void    powMod_(x,y,n)        //do x=x**y mod n, where x,y,n are bigInts (n is odd) and ** is exponentiation.  0**0=1.
       
    96 // void    randBigInt_(b,n,s)    //do b = an n-bit random BigInt. if s=1, then nth bit (most significant bit) is set to 1. n>=1.
       
    97 // void    randTruePrime_(ans,k) //do ans = a random k-bit true random prime (not just probable prime) with 1 in the msb.
       
    98 // void    sub_(x,y)             //do x=x-y for bigInts x and y. Negative answers will be 2s complement.
       
    99 //
       
   100 // The following functions do NOT have a non-underscored version. 
       
   101 // They each write a bigInt result to one or more parameters.  The caller is responsible for
       
   102 // ensuring the arrays passed in are large enough to hold the results. 
       
   103 //
       
   104 // void addShift_(x,y,ys)       //do x=x+(y<<(ys*bpe))
       
   105 // void carry_(x)               //do carries and borrows so each element of the bigInt x fits in bpe bits.
       
   106 // void divide_(x,y,q,r)        //divide x by y giving quotient q and remainder r
       
   107 // int  divInt_(x,n)            //do x=floor(x/n) for bigInt x and integer n, and return the remainder. (This never overflows its array).
       
   108 // int  eGCD_(x,y,d,a,b)        //sets a,b,d to positive bigInts such that d = GCD_(x,y) = a*x-b*y
       
   109 // void halve_(x)               //do x=floor(|x|/2)*sgn(x) for bigInt x in 2's complement.  (This never overflows its array).
       
   110 // void leftShift_(x,n)         //left shift bigInt x by n bits.  n<bpe.
       
   111 // void linComb_(x,y,a,b)       //do x=a*x+b*y for bigInts x and y and integers a and b
       
   112 // void linCombShift_(x,y,b,ys) //do x=x+b*(y<<(ys*bpe)) for bigInts x and y, and integers b and ys
       
   113 // void mont_(x,y,n,np)         //Montgomery multiplication (see comments where the function is defined)
       
   114 // void multInt_(x,n)           //do x=x*n where x is a bigInt and n is an integer.
       
   115 // void rightShift_(x,n)        //right shift bigInt x by n bits.  0 <= n < bpe. (This never overflows its array).
       
   116 // void squareMod_(x,n)         //do x=x*x  mod n for bigInts x,n
       
   117 // void subShift_(x,y,ys)       //do x=x-(y<<(ys*bpe)). Negative answers will be 2s complement.
       
   118 //
       
   119 // The following functions are based on algorithms from the _Handbook of Applied Cryptography_
       
   120 //    powMod_()           = algorithm 14.94, Montgomery exponentiation
       
   121 //    eGCD_,inverseMod_() = algorithm 14.61, Binary extended GCD_
       
   122 //    GCD_()              = algorothm 14.57, Lehmer's algorithm
       
   123 //    mont_()             = algorithm 14.36, Montgomery multiplication
       
   124 //    divide_()           = algorithm 14.20  Multiple-precision division
       
   125 //    squareMod_()        = algorithm 14.16  Multiple-precision squaring
       
   126 //    randTruePrime_()    = algorithm  4.62, Maurer's algorithm
       
   127 //    millerRabin()       = algorithm  4.24, Miller-Rabin algorithm
       
   128 //
       
   129 // Profiling shows:
       
   130 //     randTruePrime_() spends:
       
   131 //         10% of its time in calls to powMod_()
       
   132 //         85% of its time in calls to millerRabin()
       
   133 //     millerRabin() spends:
       
   134 //         99% of its time in calls to powMod_()   (always with a base of 2)
       
   135 //     powMod_() spends:
       
   136 //         94% of its time in calls to mont_()  (almost always with x==y)
       
   137 //
       
   138 // This suggests there are several ways to speed up this library slightly:
       
   139 //     - convert powMod_ to use a Montgomery form of k-ary window (or maybe a Montgomery form of sliding window)
       
   140 //         -- this should especially focus on being fast when raising 2 to a power mod n
       
   141 //     - convert randTruePrime_() to use a minimum r of 1/3 instead of 1/2 with the appropriate change to the test
       
   142 //     - tune the parameters in randTruePrime_(), including c, m, and recLimit
       
   143 //     - speed up the single loop in mont_() that takes 95% of the runtime, perhaps by reducing checking
       
   144 //       within the loop when all the parameters are the same length.
       
   145 //
       
   146 // There are several ideas that look like they wouldn't help much at all:
       
   147 //     - replacing trial division in randTruePrime_() with a sieve (that speeds up something taking almost no time anyway)
       
   148 //     - increase bpe from 15 to 30 (that would help if we had a 32*32->64 multiplier, but not with JavaScript's 32*32->32)
       
   149 //     - speeding up mont_(x,y,n,np) when x==y by doing a non-modular, non-Montgomery square
       
   150 //       followed by a Montgomery reduction.  The intermediate answer will be twice as long as x, so that
       
   151 //       method would be slower.  This is unfortunate because the code currently spends almost all of its time
       
   152 //       doing mont_(x,x,...), both for randTruePrime_() and powMod_().  A faster method for Montgomery squaring
       
   153 //       would have a large impact on the speed of randTruePrime_() and powMod_().  HAC has a couple of poorly-worded
       
   154 //       sentences that seem to imply it's faster to do a non-modular square followed by a single
       
   155 //       Montgomery reduction, but that's obviously wrong.
       
   156 ////////////////////////////////////////////////////////////////////////////////////////
       
   157 
       
   158 //globals
       
   159 bpe=0;         //bits stored per array element
       
   160 mask=0;        //AND this with an array element to chop it down to bpe bits
       
   161 radix=mask+1;  //equals 2^bpe.  A single 1 bit to the left of the last bit of mask.
       
   162 
       
   163 //the digits for converting to different bases
       
   164 digitsStr='0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz_=!@#$%^&*()[]{}|;:,.<>/?`~ \\\'\"+-';
       
   165 
       
   166 //initialize the global variables
       
   167 for (bpe=0; (1<<(bpe+1)) > (1<<bpe); bpe++);  //bpe=number of bits in the mantissa on this platform
       
   168 bpe>>=1;                   //bpe=number of bits in one element of the array representing the bigInt
       
   169 mask=(1<<bpe)-1;           //AND the mask with an integer to get its bpe least significant bits
       
   170 radix=mask+1;              //2^bpe.  a single 1 bit to the left of the first bit of mask
       
   171 one=int2bigInt(1,1,1);     //constant used in powMod_()
       
   172 
       
   173 //the following global variables are scratchpad memory to 
       
   174 //reduce dynamic memory allocation in the inner loop
       
   175 t=new Array(0);
       
   176 ss=t;       //used in mult_()
       
   177 s0=t;       //used in multMod_(), squareMod_() 
       
   178 s1=t;       //used in powMod_(), multMod_(), squareMod_() 
       
   179 s2=t;       //used in powMod_(), multMod_()
       
   180 s3=t;       //used in powMod_()
       
   181 s4=t; s5=t; //used in mod_()
       
   182 s6=t;       //used in bigInt2str()
       
   183 s7=t;       //used in powMod_()
       
   184 T=t;        //used in GCD_()
       
   185 sa=t;       //used in mont_()
       
   186 mr_x1=t; mr_r=t; mr_a=t;                                      //used in millerRabin()
       
   187 eg_v=t; eg_u=t; eg_A=t; eg_B=t; eg_C=t; eg_D=t;               //used in eGCD_(), inverseMod_()
       
   188 md_q1=t; md_q2=t; md_q3=t; md_r=t; md_r1=t; md_r2=t; md_tt=t; //used in mod_()
       
   189 
       
   190 primes=t; pows=t; s_i=t; s_i2=t; s_R=t; s_rm=t; s_q=t; s_n1=t; 
       
   191   s_a=t; s_r2=t; s_n=t; s_b=t; s_d=t; s_x1=t; s_x2=t, s_aa=t; //used in randTruePrime_()
       
   192 
       
   193 ////////////////////////////////////////////////////////////////////////////////////////
       
   194 
       
   195 //return array of all primes less than integer n
       
   196 function findPrimes(n) {
       
   197   var i,s,p,ans;
       
   198   s=new Array(n);
       
   199   for (i=0;i<n;i++)
       
   200     s[i]=0;
       
   201   s[0]=2;
       
   202   p=0;    //first p elements of s are primes, the rest are a sieve
       
   203   for(;s[p]<n;) {                  //s[p] is the pth prime
       
   204     for(i=s[p]*s[p]; i<n; i+=s[p]) //mark multiples of s[p]
       
   205       s[i]=1;
       
   206     p++;
       
   207     s[p]=s[p-1]+1;
       
   208     for(; s[p]<n && s[s[p]]; s[p]++); //find next prime (where s[p]==0)
       
   209   }
       
   210   ans=new Array(p);
       
   211   for(i=0;i<p;i++)
       
   212     ans[i]=s[i];
       
   213   return ans;
       
   214 }
       
   215 
       
   216 //does a single round of Miller-Rabin base b consider x to be a possible prime?
       
   217 //x is a bigInt, and b is an integer
       
   218 function millerRabin(x,b) {
       
   219   var i,j,k,s;
       
   220 
       
   221   if (mr_x1.length!=x.length) {
       
   222     mr_x1=dup(x);
       
   223     mr_r=dup(x);
       
   224     mr_a=dup(x);
       
   225   }
       
   226 
       
   227   copyInt_(mr_a,b);
       
   228   copy_(mr_r,x);
       
   229   copy_(mr_x1,x);
       
   230 
       
   231   addInt_(mr_r,-1);
       
   232   addInt_(mr_x1,-1);
       
   233 
       
   234   //s=the highest power of two that divides mr_r
       
   235   k=0;
       
   236   for (i=0;i<mr_r.length;i++)
       
   237     for (j=1;j<mask;j<<=1)
       
   238       if (x[i] & j) {
       
   239         s=(k<mr_r.length+bpe ? k : 0); 
       
   240          i=mr_r.length;
       
   241          j=mask;
       
   242       } else
       
   243         k++;
       
   244 
       
   245   if (s)                
       
   246     rightShift_(mr_r,s);
       
   247 
       
   248   powMod_(mr_a,mr_r,x);
       
   249 
       
   250   if (!equalsInt(mr_a,1) && !equals(mr_a,mr_x1)) {
       
   251     j=1;
       
   252     while (j<=s-1 && !equals(mr_a,mr_x1)) {
       
   253       squareMod_(mr_a,x);
       
   254       if (equalsInt(mr_a,1)) {
       
   255         return 0;
       
   256       }
       
   257       j++;
       
   258     }
       
   259     if (!equals(mr_a,mr_x1)) {
       
   260       return 0;
       
   261     }
       
   262   }
       
   263   return 1;  
       
   264 }
       
   265 
       
   266 //returns how many bits long the bigInt is, not counting leading zeros.
       
   267 function bitSize(x) {
       
   268   var j,z,w;
       
   269   for (j=x.length-1; (x[j]==0) && (j>0); j--);
       
   270   for (z=0,w=x[j]; w; (w>>=1),z++);
       
   271   z+=bpe*j;
       
   272   return z;
       
   273 }
       
   274 
       
   275 //return a copy of x with at least n elements, adding leading zeros if needed
       
   276 function expand(x,n) {
       
   277   var ans=int2bigInt(0,(x.length>n ? x.length : n)*bpe,0);
       
   278   copy_(ans,x);
       
   279   return ans;
       
   280 }
       
   281 
       
   282 //return a k-bit true random prime using Maurer's algorithm.
       
   283 function randTruePrime(k) {
       
   284   var ans=int2bigInt(0,k,0);
       
   285   randTruePrime_(ans,k);
       
   286   return bigint_trim(ans,1);
       
   287 }
       
   288 
       
   289 //return a new bigInt equal to (x mod n) for bigInts x and n.
       
   290 function mod(x,n) {
       
   291   var ans=dup(x);
       
   292   mod_(ans,n);
       
   293   return bigint_trim(ans,1);
       
   294 }
       
   295 
       
   296 //return (x+n) where x is a bigInt and n is an integer.
       
   297 function addInt(x,n) {
       
   298   var ans=expand(x,x.length+1);
       
   299   addInt_(ans,n);
       
   300   return bigint_trim(ans,1);
       
   301 }
       
   302 
       
   303 //return x*y for bigInts x and y. This is faster when y<x.
       
   304 function mult(x,y) {
       
   305   var ans=expand(x,x.length+y.length);
       
   306   mult_(ans,y);
       
   307   return bigint_trim(ans,1);
       
   308 }
       
   309 
       
   310 //return (x**y mod n) where x,y,n are bigInts and ** is exponentiation.  0**0=1. Faster for odd n.
       
   311 function powMod(x,y,n) {
       
   312   var ans=expand(x,n.length);  
       
   313   powMod_(ans,bigint_trim(y,2),bigint_trim(n,2),0);  //this should work without the trim, but doesn't
       
   314   return bigint_trim(ans,1);
       
   315 }
       
   316 
       
   317 //return (x-y) for bigInts x and y.  Negative answers will be 2s complement
       
   318 function sub(x,y) {
       
   319   var ans=expand(x,(x.length>y.length ? x.length+1 : y.length+1)); 
       
   320   sub_(ans,y);
       
   321   return bigint_trim(ans,1);
       
   322 }
       
   323 
       
   324 //return (x+y) for bigInts x and y.  
       
   325 function add(x,y) {
       
   326   var ans=expand(x,(x.length>y.length ? x.length+1 : y.length+1)); 
       
   327   add_(ans,y);
       
   328   return bigint_trim(ans,1);
       
   329 }
       
   330 
       
   331 //return (x**(-1) mod n) for bigInts x and n.  If no inverse exists, it returns null
       
   332 function inverseMod(x,n) {
       
   333   var ans=expand(x,n.length); 
       
   334   var s;
       
   335   s=inverseMod_(ans,n);
       
   336   return s ? bigint_trim(ans,1) : null;
       
   337 }
       
   338 
       
   339 //return (x*y mod n) for bigInts x,y,n.  For greater speed, let y<x.
       
   340 function multMod(x,y,n) {
       
   341   var ans=expand(x,n.length);
       
   342   multMod_(ans,y,n);
       
   343   return bigint_trim(ans,1);
       
   344 }
       
   345 
       
   346 //generate a k-bit true random prime using Maurer's algorithm,
       
   347 //and put it into ans.  The bigInt ans must be large enough to hold it.
       
   348 function randTruePrime_(ans,k) {
       
   349   var c,m,pm,dd,j,r,B,divisible,z,zz,recSize;
       
   350 
       
   351   if (primes.length==0)
       
   352     primes=findPrimes(30000);  //check for divisibility by primes <=30000
       
   353 
       
   354   if (pows.length==0) {
       
   355     pows=new Array(512);
       
   356     for (j=0;j<512;j++) {
       
   357       pows[j]=Math.pow(2,j/511.-1.);
       
   358     }
       
   359   }
       
   360 
       
   361   //c and m should be tuned for a particular machine and value of k, to maximize speed
       
   362   c=0.1;  //c=0.1 in HAC
       
   363   m=20;   //generate this k-bit number by first recursively generating a number that has between k/2 and k-m bits
       
   364   recLimit=20; //stop recursion when k <=recLimit.  Must have recLimit >= 2
       
   365 
       
   366   if (s_i2.length!=ans.length) {
       
   367     s_i2=dup(ans);
       
   368     s_R =dup(ans);
       
   369     s_n1=dup(ans);
       
   370     s_r2=dup(ans);
       
   371     s_d =dup(ans);
       
   372     s_x1=dup(ans);
       
   373     s_x2=dup(ans);
       
   374     s_b =dup(ans);
       
   375     s_n =dup(ans);
       
   376     s_i =dup(ans);
       
   377     s_rm=dup(ans);
       
   378     s_q =dup(ans);
       
   379     s_a =dup(ans);
       
   380     s_aa=dup(ans);
       
   381   }
       
   382 
       
   383   if (k <= recLimit) {  //generate small random primes by trial division up to its square root
       
   384     pm=(1<<((k+2)>>1))-1; //pm is binary number with all ones, just over sqrt(2^k)
       
   385     copyInt_(ans,0);
       
   386     for (dd=1;dd;) {
       
   387       dd=0;
       
   388       ans[0]= 1 | (1<<(k-1)) | Math.floor(Math.random()*(1<<k));  //random, k-bit, odd integer, with msb 1
       
   389       for (j=1;(j<primes.length) && ((primes[j]&pm)==primes[j]);j++) { //trial division by all primes 3...sqrt(2^k)
       
   390         if (0==(ans[0]%primes[j])) {
       
   391           dd=1;
       
   392           break;
       
   393         }
       
   394       }
       
   395     }
       
   396     carry_(ans);
       
   397     return;
       
   398   }
       
   399 
       
   400   B=c*k*k;    //try small primes up to B (or all the primes[] array if the largest is less than B).
       
   401   if (k>2*m)  //generate this k-bit number by first recursively generating a number that has between k/2 and k-m bits
       
   402     for (r=1; k-k*r<=m; )
       
   403       r=pows[Math.floor(Math.random()*512)];   //r=Math.pow(2,Math.random()-1);
       
   404   else
       
   405     r=.5;
       
   406 
       
   407   //simulation suggests the more complex algorithm using r=.333 is only slightly faster.
       
   408 
       
   409   recSize=Math.floor(r*k)+1;
       
   410 
       
   411   randTruePrime_(s_q,recSize);
       
   412   copyInt_(s_i2,0);
       
   413   s_i2[Math.floor((k-2)/bpe)] |= (1<<((k-2)%bpe));   //s_i2=2^(k-2)
       
   414   divide_(s_i2,s_q,s_i,s_rm);                        //s_i=floor((2^(k-1))/(2q))
       
   415 
       
   416   z=bitSize(s_i);
       
   417 
       
   418   for (;;) {
       
   419     for (;;) {  //generate z-bit numbers until one falls in the range [0,s_i-1]
       
   420       randBigInt_(s_R,z,0);
       
   421       if (greater(s_i,s_R))
       
   422         break;
       
   423     }                //now s_R is in the range [0,s_i-1]
       
   424     addInt_(s_R,1);  //now s_R is in the range [1,s_i]
       
   425     add_(s_R,s_i);   //now s_R is in the range [s_i+1,2*s_i]
       
   426 
       
   427     copy_(s_n,s_q);
       
   428     mult_(s_n,s_R); 
       
   429     multInt_(s_n,2);
       
   430     addInt_(s_n,1);    //s_n=2*s_R*s_q+1
       
   431     
       
   432     copy_(s_r2,s_R);
       
   433     multInt_(s_r2,2);  //s_r2=2*s_R
       
   434 
       
   435     //check s_n for divisibility by small primes up to B
       
   436     for (divisible=0,j=0; (j<primes.length) && (primes[j]<B); j++)
       
   437       if (modInt(s_n,primes[j])==0) {
       
   438         divisible=1;
       
   439         break;
       
   440       }      
       
   441 
       
   442     if (!divisible)    //if it passes small primes check, then try a single Miller-Rabin base 2
       
   443       if (!millerRabin(s_n,2)) //this line represents 75% of the total runtime for randTruePrime_ 
       
   444         divisible=1;
       
   445 
       
   446     if (!divisible) {  //if it passes that test, continue checking s_n
       
   447       addInt_(s_n,-3);
       
   448       for (j=s_n.length-1;(s_n[j]==0) && (j>0); j--);  //strip leading zeros
       
   449       for (zz=0,w=s_n[j]; w; (w>>=1),zz++);
       
   450       zz+=bpe*j;                             //zz=number of bits in s_n, ignoring leading zeros
       
   451       for (;;) {  //generate z-bit numbers until one falls in the range [0,s_n-1]
       
   452         randBigInt_(s_a,zz,0);
       
   453         if (greater(s_n,s_a))
       
   454           break;
       
   455       }                //now s_a is in the range [0,s_n-1]
       
   456       addInt_(s_n,3);  //now s_a is in the range [0,s_n-4]
       
   457       addInt_(s_a,2);  //now s_a is in the range [2,s_n-2]
       
   458       copy_(s_b,s_a);
       
   459       copy_(s_n1,s_n);
       
   460       addInt_(s_n1,-1);
       
   461       powMod_(s_b,s_n1,s_n);   //s_b=s_a^(s_n-1) modulo s_n
       
   462       addInt_(s_b,-1);
       
   463       if (isZero(s_b)) {
       
   464         copy_(s_b,s_a);
       
   465         powMod_(s_b,s_r2,s_n);
       
   466         addInt_(s_b,-1);
       
   467         copy_(s_aa,s_n);
       
   468         copy_(s_d,s_b);
       
   469         GCD_(s_d,s_n);  //if s_b and s_n are relatively prime, then s_n is a prime
       
   470         if (equalsInt(s_d,1)) {
       
   471           copy_(ans,s_aa);
       
   472           return;     //if we've made it this far, then s_n is absolutely guaranteed to be prime
       
   473         }
       
   474       }
       
   475     }
       
   476   }
       
   477 }
       
   478 
       
   479 //Return an n-bit random BigInt (n>=1).  If s=1, then the most significant of those n bits is set to 1.
       
   480 function randBigInt(n,s) {
       
   481   var a,b;
       
   482   a=Math.floor((n-1)/bpe)+2; //# array elements to hold the BigInt with a leading 0 element
       
   483   b=int2bigInt(0,0,a);
       
   484   randBigInt_(b,n,s);
       
   485   return b;
       
   486 }
       
   487 
       
   488 //Set b to an n-bit random BigInt.  If s=1, then the most significant of those n bits is set to 1.
       
   489 //Array b must be big enough to hold the result. Must have n>=1
       
   490 function randBigInt_(b,n,s) {
       
   491   var i,a;
       
   492   for (i=0;i<b.length;i++)
       
   493     b[i]=0;
       
   494   a=Math.floor((n-1)/bpe)+1; //# array elements to hold the BigInt
       
   495   for (i=0;i<a;i++) {
       
   496     b[i]=Math.floor(Math.random()*(1<<(bpe-1)));
       
   497   }
       
   498   b[a-1] &= (2<<((n-1)%bpe))-1;
       
   499   if (s==1)
       
   500     b[a-1] |= (1<<((n-1)%bpe));
       
   501 }
       
   502 
       
   503 //Return the greatest common divisor of bigInts x and y (each with same number of elements).
       
   504 function GCD(x,y) {
       
   505   var xc,yc;
       
   506   xc=dup(x);
       
   507   yc=dup(y);
       
   508   GCD_(xc,yc);
       
   509   return xc;
       
   510 }
       
   511 
       
   512 //set x to the greatest common divisor of bigInts x and y (each with same number of elements).
       
   513 //y is destroyed.
       
   514 function GCD_(x,y) {
       
   515   var i,xp,yp,A,B,C,D,q,sing;
       
   516   if (T.length!=x.length)
       
   517     T=dup(x);
       
   518 
       
   519   sing=1;
       
   520   while (sing) { //while y has nonzero elements other than y[0]
       
   521     sing=0;
       
   522     for (i=1;i<y.length;i++) //check if y has nonzero elements other than 0
       
   523       if (y[i]) {
       
   524         sing=1;
       
   525         break;
       
   526       }
       
   527     if (!sing) break; //quit when y all zero elements except possibly y[0]
       
   528 
       
   529     for (i=x.length;!x[i] && i>=0;i--);  //find most significant element of x
       
   530     xp=x[i];
       
   531     yp=y[i];
       
   532     A=1; B=0; C=0; D=1;
       
   533     while ((yp+C) && (yp+D)) {
       
   534       q =Math.floor((xp+A)/(yp+C));
       
   535       qp=Math.floor((xp+B)/(yp+D));
       
   536       if (q!=qp)
       
   537         break;
       
   538       t= A-q*C;   A=C;   C=t;    //  do (A,B,xp, C,D,yp) = (C,D,yp, A,B,xp) - q*(0,0,0, C,D,yp)      
       
   539       t= B-q*D;   B=D;   D=t;
       
   540       t=xp-q*yp; xp=yp; yp=t;
       
   541     }
       
   542     if (B) {
       
   543       copy_(T,x);
       
   544       linComb_(x,y,A,B); //x=A*x+B*y
       
   545       linComb_(y,T,D,C); //y=D*y+C*T
       
   546     } else {
       
   547       mod_(x,y);
       
   548       copy_(T,x);
       
   549       copy_(x,y);
       
   550       copy_(y,T);
       
   551     } 
       
   552   }
       
   553   if (y[0]==0)
       
   554     return;
       
   555   t=modInt(x,y[0]);
       
   556   copyInt_(x,y[0]);
       
   557   y[0]=t;
       
   558   while (y[0]) {
       
   559     x[0]%=y[0];
       
   560     t=x[0]; x[0]=y[0]; y[0]=t;
       
   561   }
       
   562 }
       
   563 
       
   564 //do x=x**(-1) mod n, for bigInts x and n.
       
   565 //If no inverse exists, it sets x to zero and returns 0, else it returns 1.
       
   566 //The x array must be at least as large as the n array.
       
   567 function inverseMod_(x,n) {
       
   568   var k=1+2*Math.max(x.length,n.length);
       
   569 
       
   570   if(!(x[0]&1)  && !(n[0]&1)) {  //if both inputs are even, then inverse doesn't exist
       
   571     copyInt_(x,0);
       
   572     return 0;
       
   573   }
       
   574 
       
   575   if (eg_u.length!=k) {
       
   576     eg_u=new Array(k);
       
   577     eg_v=new Array(k);
       
   578     eg_A=new Array(k);
       
   579     eg_B=new Array(k);
       
   580     eg_C=new Array(k);
       
   581     eg_D=new Array(k);
       
   582   }
       
   583 
       
   584   copy_(eg_u,x);
       
   585   copy_(eg_v,n);
       
   586   copyInt_(eg_A,1);
       
   587   copyInt_(eg_B,0);
       
   588   copyInt_(eg_C,0);
       
   589   copyInt_(eg_D,1);
       
   590   for (;;) {
       
   591     while(!(eg_u[0]&1)) {  //while eg_u is even
       
   592       halve_(eg_u);
       
   593       if (!(eg_A[0]&1) && !(eg_B[0]&1)) { //if eg_A==eg_B==0 mod 2
       
   594         halve_(eg_A);
       
   595         halve_(eg_B);      
       
   596       } else {
       
   597         add_(eg_A,n);  halve_(eg_A);
       
   598         sub_(eg_B,x);  halve_(eg_B);
       
   599       }
       
   600     }
       
   601 
       
   602     while (!(eg_v[0]&1)) {  //while eg_v is even
       
   603       halve_(eg_v);
       
   604       if (!(eg_C[0]&1) && !(eg_D[0]&1)) { //if eg_C==eg_D==0 mod 2
       
   605         halve_(eg_C);
       
   606         halve_(eg_D);      
       
   607       } else {
       
   608         add_(eg_C,n);  halve_(eg_C);
       
   609         sub_(eg_D,x);  halve_(eg_D);
       
   610       }
       
   611     }
       
   612 
       
   613     if (!greater(eg_v,eg_u)) { //eg_v <= eg_u
       
   614       sub_(eg_u,eg_v);
       
   615       sub_(eg_A,eg_C);
       
   616       sub_(eg_B,eg_D);
       
   617     } else {                   //eg_v > eg_u
       
   618       sub_(eg_v,eg_u);
       
   619       sub_(eg_C,eg_A);
       
   620       sub_(eg_D,eg_B);
       
   621     }
       
   622   
       
   623     if (equalsInt(eg_u,0)) {
       
   624       if (negative(eg_C)) //make sure answer is nonnegative
       
   625         add_(eg_C,n);
       
   626       copy_(x,eg_C);
       
   627 
       
   628       if (!equalsInt(eg_v,1)) { //if GCD_(x,n)!=1, then there is no inverse
       
   629         copyInt_(x,0);
       
   630         return 0;
       
   631       }
       
   632       return 1;
       
   633     }
       
   634   }
       
   635 }
       
   636 
       
   637 //return x**(-1) mod n, for integers x and n.  Return 0 if there is no inverse
       
   638 function inverseModInt(x,n) {
       
   639   var a=1,b=0,t;
       
   640   for (;;) {
       
   641     if (x==1) return a;
       
   642     if (x==0) return 0;
       
   643     b-=a*Math.floor(n/x);
       
   644     n%=x;
       
   645 
       
   646     if (n==1) return b; //to avoid negatives, change this b to n-b, and each -= to +=
       
   647     if (n==0) return 0;
       
   648     a-=b*Math.floor(x/n);
       
   649     x%=n;
       
   650   }
       
   651 }
       
   652 
       
   653 //this deprecated function is for backward compatibility only. 
       
   654 function inverseModInt_(x,n) {
       
   655    return inverseModInt(x,n);
       
   656 }
       
   657 
       
   658 
       
   659 //Given positive bigInts x and y, change the bigints v, a, and b to positive bigInts such that:
       
   660 //     v = GCD_(x,y) = a*x-b*y
       
   661 //The bigInts v, a, b, must have exactly as many elements as the larger of x and y.
       
   662 function eGCD_(x,y,v,a,b) {
       
   663   var g=0;
       
   664   var k=Math.max(x.length,y.length);
       
   665   if (eg_u.length!=k) {
       
   666     eg_u=new Array(k);
       
   667     eg_A=new Array(k);
       
   668     eg_B=new Array(k);
       
   669     eg_C=new Array(k);
       
   670     eg_D=new Array(k);
       
   671   }
       
   672   while(!(x[0]&1)  && !(y[0]&1)) {  //while x and y both even
       
   673     halve_(x);
       
   674     halve_(y);
       
   675     g++;
       
   676   }
       
   677   copy_(eg_u,x);
       
   678   copy_(v,y);
       
   679   copyInt_(eg_A,1);
       
   680   copyInt_(eg_B,0);
       
   681   copyInt_(eg_C,0);
       
   682   copyInt_(eg_D,1);
       
   683   for (;;) {
       
   684     while(!(eg_u[0]&1)) {  //while u is even
       
   685       halve_(eg_u);
       
   686       if (!(eg_A[0]&1) && !(eg_B[0]&1)) { //if A==B==0 mod 2
       
   687         halve_(eg_A);
       
   688         halve_(eg_B);      
       
   689       } else {
       
   690         add_(eg_A,y);  halve_(eg_A);
       
   691         sub_(eg_B,x);  halve_(eg_B);
       
   692       }
       
   693     }
       
   694 
       
   695     while (!(v[0]&1)) {  //while v is even
       
   696       halve_(v);
       
   697       if (!(eg_C[0]&1) && !(eg_D[0]&1)) { //if C==D==0 mod 2
       
   698         halve_(eg_C);
       
   699         halve_(eg_D);      
       
   700       } else {
       
   701         add_(eg_C,y);  halve_(eg_C);
       
   702         sub_(eg_D,x);  halve_(eg_D);
       
   703       }
       
   704     }
       
   705 
       
   706     if (!greater(v,eg_u)) { //v<=u
       
   707       sub_(eg_u,v);
       
   708       sub_(eg_A,eg_C);
       
   709       sub_(eg_B,eg_D);
       
   710     } else {                //v>u
       
   711       sub_(v,eg_u);
       
   712       sub_(eg_C,eg_A);
       
   713       sub_(eg_D,eg_B);
       
   714     }
       
   715     if (equalsInt(eg_u,0)) {
       
   716       if (negative(eg_C)) {   //make sure a (C)is nonnegative
       
   717         add_(eg_C,y);
       
   718         sub_(eg_D,x);
       
   719       }
       
   720       multInt_(eg_D,-1);  ///make sure b (D) is nonnegative
       
   721       copy_(a,eg_C);
       
   722       copy_(b,eg_D);
       
   723       leftShift_(v,g);
       
   724       return;
       
   725     }
       
   726   }
       
   727 }
       
   728 
       
   729 
       
   730 //is bigInt x negative?
       
   731 function negative(x) {
       
   732   return ((x[x.length-1]>>(bpe-1))&1);
       
   733 }
       
   734 
       
   735 
       
   736 //is (x << (shift*bpe)) > y?
       
   737 //x and y are nonnegative bigInts
       
   738 //shift is a nonnegative integer
       
   739 function greaterShift(x,y,shift) {
       
   740   var kx=x.length, ky=y.length;
       
   741   k=((kx+shift)<ky) ? (kx+shift) : ky;
       
   742   for (i=ky-1-shift; i<kx && i>=0; i++) 
       
   743     if (x[i]>0)
       
   744       return 1; //if there are nonzeros in x to the left of the first column of y, then x is bigger
       
   745   for (i=kx-1+shift; i<ky; i++)
       
   746     if (y[i]>0)
       
   747       return 0; //if there are nonzeros in y to the left of the first column of x, then x is not bigger
       
   748   for (i=k-1; i>=shift; i--)
       
   749     if      (x[i-shift]>y[i]) return 1;
       
   750     else if (x[i-shift]<y[i]) return 0;
       
   751   return 0;
       
   752 }
       
   753 
       
   754 //is x > y? (x and y both nonnegative)
       
   755 function greater(x,y) {
       
   756   var i;
       
   757   var k=(x.length<y.length) ? x.length : y.length;
       
   758 
       
   759   for (i=x.length;i<y.length;i++)
       
   760     if (y[i])
       
   761       return 0;  //y has more digits
       
   762 
       
   763   for (i=y.length;i<x.length;i++)
       
   764     if (x[i])
       
   765       return 1;  //x has more digits
       
   766 
       
   767   for (i=k-1;i>=0;i--)
       
   768     if (x[i]>y[i])
       
   769       return 1;
       
   770     else if (x[i]<y[i])
       
   771       return 0;
       
   772   return 0;
       
   773 }
       
   774 
       
   775 //divide x by y giving quotient q and remainder r.  (q=floor(x/y),  r=x mod y).  All 4 are bigints.
       
   776 //x must have at least one leading zero element.
       
   777 //y must be nonzero.
       
   778 //q and r must be arrays that are exactly the same length as x. (Or q can have more).
       
   779 //Must have x.length >= y.length >= 2.
       
   780 function divide_(x,y,q,r) {
       
   781   var kx, ky;
       
   782   var i,j,y1,y2,c,a,b;
       
   783   copy_(r,x);
       
   784   for (ky=y.length;y[ky-1]==0;ky--); //ky is number of elements in y, not including leading zeros
       
   785 
       
   786   //normalize: ensure the most significant element of y has its highest bit set  
       
   787   b=y[ky-1];
       
   788   for (a=0; b; a++)
       
   789     b>>=1;  
       
   790   a=bpe-a;  //a is how many bits to shift so that the high order bit of y is leftmost in its array element
       
   791   leftShift_(y,a);  //multiply both by 1<<a now, then divide both by that at the end
       
   792   leftShift_(r,a);
       
   793 
       
   794   //Rob Visser discovered a bug: the following line was originally just before the normalization.
       
   795   for (kx=r.length;r[kx-1]==0 && kx>ky;kx--); //kx is number of elements in normalized x, not including leading zeros
       
   796 
       
   797   copyInt_(q,0);                      // q=0
       
   798   while (!greaterShift(y,r,kx-ky)) {  // while (leftShift_(y,kx-ky) <= r) {
       
   799     subShift_(r,y,kx-ky);             //   r=r-leftShift_(y,kx-ky)
       
   800     q[kx-ky]++;                       //   q[kx-ky]++;
       
   801   }                                   // }
       
   802 
       
   803   for (i=kx-1; i>=ky; i--) {
       
   804     if (r[i]==y[ky-1])
       
   805       q[i-ky]=mask;
       
   806     else
       
   807       q[i-ky]=Math.floor((r[i]*radix+r[i-1])/y[ky-1]);	
       
   808 
       
   809     //The following for(;;) loop is equivalent to the commented while loop, 
       
   810     //except that the uncommented version avoids overflow.
       
   811     //The commented loop comes from HAC, which assumes r[-1]==y[-1]==0
       
   812     //  while (q[i-ky]*(y[ky-1]*radix+y[ky-2]) > r[i]*radix*radix+r[i-1]*radix+r[i-2])
       
   813     //    q[i-ky]--;    
       
   814     for (;;) {
       
   815       y2=(ky>1 ? y[ky-2] : 0)*q[i-ky];
       
   816       c=y2>>bpe;
       
   817       y2=y2 & mask;
       
   818       y1=c+q[i-ky]*y[ky-1];
       
   819       c=y1>>bpe;
       
   820       y1=y1 & mask;
       
   821 
       
   822       if (c==r[i] ? y1==r[i-1] ? y2>(i>1 ? r[i-2] : 0) : y1>r[i-1] : c>r[i]) 
       
   823         q[i-ky]--;
       
   824       else
       
   825         break;
       
   826     }
       
   827 
       
   828     linCombShift_(r,y,-q[i-ky],i-ky);    //r=r-q[i-ky]*leftShift_(y,i-ky)
       
   829     if (negative(r)) {
       
   830       addShift_(r,y,i-ky);         //r=r+leftShift_(y,i-ky)
       
   831       q[i-ky]--;
       
   832     }
       
   833   }
       
   834 
       
   835   rightShift_(y,a);  //undo the normalization step
       
   836   rightShift_(r,a);  //undo the normalization step
       
   837 }
       
   838 
       
   839 //do carries and borrows so each element of the bigInt x fits in bpe bits.
       
   840 function carry_(x) {
       
   841   var i,k,c,b;
       
   842   k=x.length;
       
   843   c=0;
       
   844   for (i=0;i<k;i++) {
       
   845     c+=x[i];
       
   846     b=0;
       
   847     if (c<0) {
       
   848       b=-(c>>bpe);
       
   849       c+=b*radix;
       
   850     }
       
   851     x[i]=c & mask;
       
   852     c=(c>>bpe)-b;
       
   853   }
       
   854 }
       
   855 
       
   856 //return x mod n for bigInt x and integer n.
       
   857 function modInt(x,n) {
       
   858   var i,c=0;
       
   859   for (i=x.length-1; i>=0; i--)
       
   860     c=(c*radix+x[i])%n;
       
   861   return c;
       
   862 }
       
   863 
       
   864 //convert the integer t into a bigInt with at least the given number of bits.
       
   865 //the returned array stores the bigInt in bpe-bit chunks, little endian (buff[0] is least significant word)
       
   866 //Pad the array with leading zeros so that it has at least minSize elements.
       
   867 //There will always be at least one leading 0 element.
       
   868 function int2bigInt(t,bits,minSize) {   
       
   869   var i,k;
       
   870   k=Math.ceil(bits/bpe)+1;
       
   871   k=minSize>k ? minSize : k;
       
   872   buff=new Array(k);
       
   873   copyInt_(buff,t);
       
   874   return buff;
       
   875 }
       
   876 
       
   877 //return the bigInt given a string representation in a given base.  
       
   878 //Pad the array with leading zeros so that it has at least minSize elements.
       
   879 //If base=-1, then it reads in a space-separated list of array elements in decimal.
       
   880 //The array will always have at least one leading zero, unless base=-1.
       
   881 function str2bigInt(s,base,minSize) {
       
   882   var d, i, j, x, y, kk;
       
   883   var k=s.length;
       
   884   if (base==-1) { //comma-separated list of array elements in decimal
       
   885     x=new Array(0);
       
   886     for (;;) {
       
   887       y=new Array(x.length+1);
       
   888       for (i=0;i<x.length;i++)
       
   889         y[i+1]=x[i];
       
   890       y[0]=parseInt(s,10);
       
   891       x=y;
       
   892       d=s.indexOf(',',0);
       
   893       if (d<1) 
       
   894         break;
       
   895       s=s.substring(d+1);
       
   896       if (s.length==0)
       
   897         break;
       
   898     }
       
   899     if (x.length<minSize) {
       
   900       y=new Array(minSize);
       
   901       copy_(y,x);
       
   902       return y;
       
   903     }
       
   904     return x;
       
   905   }
       
   906 
       
   907   x=int2bigInt(0,base*k,0);
       
   908   for (i=0;i<k;i++) {
       
   909     d=digitsStr.indexOf(s.substring(i,i+1),0);
       
   910     if (base<=36 && d>=36)  //convert lowercase to uppercase if base<=36
       
   911       d-=26;
       
   912     if (d<base && d>=0) {   //ignore illegal characters
       
   913       multInt_(x,base);
       
   914       addInt_(x,d);
       
   915     }
       
   916   }
       
   917 
       
   918   for (k=x.length;k>0 && !x[k-1];k--); //strip off leading zeros
       
   919   k=minSize>k+1 ? minSize : k+1;
       
   920   y=new Array(k);
       
   921   kk=k<x.length ? k : x.length;
       
   922   for (i=0;i<kk;i++)
       
   923     y[i]=x[i];
       
   924   for (;i<k;i++)
       
   925     y[i]=0;
       
   926   return y;
       
   927 }
       
   928 
       
   929 //is bigint x equal to integer y?
       
   930 //y must have less than bpe bits
       
   931 function equalsInt(x,y) {
       
   932   var i;
       
   933   if (x[0]!=y)
       
   934     return 0;
       
   935   for (i=1;i<x.length;i++)
       
   936     if (x[i])
       
   937       return 0;
       
   938   return 1;
       
   939 }
       
   940 
       
   941 //are bigints x and y equal?
       
   942 //this works even if x and y are different lengths and have arbitrarily many leading zeros
       
   943 function equals(x,y) {
       
   944   var i;
       
   945   var k=x.length<y.length ? x.length : y.length;
       
   946   for (i=0;i<k;i++)
       
   947     if (x[i]!=y[i])
       
   948       return 0;
       
   949   if (x.length>y.length) {
       
   950     for (;i<x.length;i++)
       
   951       if (x[i])
       
   952         return 0;
       
   953   } else {
       
   954     for (;i<y.length;i++)
       
   955       if (y[i])
       
   956         return 0;
       
   957   }
       
   958   return 1;
       
   959 }
       
   960 
       
   961 //is the bigInt x equal to zero?
       
   962 function isZero(x) {
       
   963   var i;
       
   964   for (i=0;i<x.length;i++)
       
   965     if (x[i])
       
   966       return 0;
       
   967   return 1;
       
   968 }
       
   969 
       
   970 //convert a bigInt into a string in a given base, from base 2 up to base 95.
       
   971 //Base -1 prints the contents of the array representing the number.
       
   972 function bigInt2str(x,base) {
       
   973   var i,t,s="";
       
   974 
       
   975   if (s6.length!=x.length) 
       
   976     s6=dup(x);
       
   977   else
       
   978     copy_(s6,x);
       
   979 
       
   980   if (base==-1) { //return the list of array contents
       
   981     for (i=x.length-1;i>0;i--)
       
   982       s+=x[i]+',';
       
   983     s+=x[0];
       
   984   }
       
   985   else { //return it in the given base
       
   986     while (!isZero(s6)) {
       
   987       t=divInt_(s6,base);  //t=s6 % base; s6=floor(s6/base);
       
   988       s=digitsStr.substring(t,t+1)+s;
       
   989     }
       
   990   }
       
   991   if (s.length==0)
       
   992     s="0";
       
   993   return s;
       
   994 }
       
   995 
       
   996 //returns a duplicate of bigInt x
       
   997 function dup(x) {
       
   998   var i;
       
   999   buff=new Array(x.length);
       
  1000   copy_(buff,x);
       
  1001   return buff;
       
  1002 }
       
  1003 
       
  1004 //do x=y on bigInts x and y.  x must be an array at least as big as y (not counting the leading zeros in y).
       
  1005 function copy_(x,y) {
       
  1006   var i;
       
  1007   var k=x.length<y.length ? x.length : y.length;
       
  1008   for (i=0;i<k;i++)
       
  1009     x[i]=y[i];
       
  1010   for (i=k;i<x.length;i++)
       
  1011     x[i]=0;
       
  1012 }
       
  1013 
       
  1014 //do x=y on bigInt x and integer y.  
       
  1015 function copyInt_(x,n) {
       
  1016   var i,c;
       
  1017   for (c=n,i=0;i<x.length;i++) {
       
  1018     x[i]=c & mask;
       
  1019     c>>=bpe;
       
  1020   }
       
  1021 }
       
  1022 
       
  1023 //do x=x+n where x is a bigInt and n is an integer.
       
  1024 //x must be large enough to hold the result.
       
  1025 function addInt_(x,n) {
       
  1026   var i,k,c,b;
       
  1027   x[0]+=n;
       
  1028   k=x.length;
       
  1029   c=0;
       
  1030   for (i=0;i<k;i++) {
       
  1031     c+=x[i];
       
  1032     b=0;
       
  1033     if (c<0) {
       
  1034       b=-(c>>bpe);
       
  1035       c+=b*radix;
       
  1036     }
       
  1037     x[i]=c & mask;
       
  1038     c=(c>>bpe)-b;
       
  1039     if (!c) return; //stop carrying as soon as the carry_ is zero
       
  1040   }
       
  1041 }
       
  1042 
       
  1043 //right shift bigInt x by n bits.  0 <= n < bpe.
       
  1044 function rightShift_(x,n) {
       
  1045   var i;
       
  1046   var k=Math.floor(n/bpe);
       
  1047   if (k) {
       
  1048     for (i=0;i<x.length-k;i++) //right shift x by k elements
       
  1049       x[i]=x[i+k];
       
  1050     for (;i<x.length;i++)
       
  1051       x[i]=0;
       
  1052     n%=bpe;
       
  1053   }
       
  1054   for (i=0;i<x.length-1;i++) {
       
  1055     x[i]=mask & ((x[i+1]<<(bpe-n)) | (x[i]>>n));
       
  1056   }
       
  1057   x[i]>>=n;
       
  1058 }
       
  1059 
       
  1060 //do x=floor(|x|/2)*sgn(x) for bigInt x in 2's complement
       
  1061 function halve_(x) {
       
  1062   var i;
       
  1063   for (i=0;i<x.length-1;i++) {
       
  1064     x[i]=mask & ((x[i+1]<<(bpe-1)) | (x[i]>>1));
       
  1065   }
       
  1066   x[i]=(x[i]>>1) | (x[i] & (radix>>1));  //most significant bit stays the same
       
  1067 }
       
  1068 
       
  1069 //left shift bigInt x by n bits.
       
  1070 function leftShift_(x,n) {
       
  1071   var i;
       
  1072   var k=Math.floor(n/bpe);
       
  1073   if (k) {
       
  1074     for (i=x.length; i>=k; i--) //left shift x by k elements
       
  1075       x[i]=x[i-k];
       
  1076     for (;i>=0;i--)
       
  1077       x[i]=0;  
       
  1078     n%=bpe;
       
  1079   }
       
  1080   if (!n)
       
  1081     return;
       
  1082   for (i=x.length-1;i>0;i--) {
       
  1083     x[i]=mask & ((x[i]<<n) | (x[i-1]>>(bpe-n)));
       
  1084   }
       
  1085   x[i]=mask & (x[i]<<n);
       
  1086 }
       
  1087 
       
  1088 //do x=x*n where x is a bigInt and n is an integer.
       
  1089 //x must be large enough to hold the result.
       
  1090 function multInt_(x,n) {
       
  1091   var i,k,c,b;
       
  1092   if (!n)
       
  1093     return;
       
  1094   k=x.length;
       
  1095   c=0;
       
  1096   for (i=0;i<k;i++) {
       
  1097     c+=x[i]*n;
       
  1098     b=0;
       
  1099     if (c<0) {
       
  1100       b=-(c>>bpe);
       
  1101       c+=b*radix;
       
  1102     }
       
  1103     x[i]=c & mask;
       
  1104     c=(c>>bpe)-b;
       
  1105   }
       
  1106 }
       
  1107 
       
  1108 //do x=floor(x/n) for bigInt x and integer n, and return the remainder
       
  1109 function divInt_(x,n) {
       
  1110   var i,r=0,s;
       
  1111   for (i=x.length-1;i>=0;i--) {
       
  1112     s=r*radix+x[i];
       
  1113     x[i]=Math.floor(s/n);
       
  1114     r=s%n;
       
  1115   }
       
  1116   return r;
       
  1117 }
       
  1118 
       
  1119 //do the linear combination x=a*x+b*y for bigInts x and y, and integers a and b.
       
  1120 //x must be large enough to hold the answer.
       
  1121 function linComb_(x,y,a,b) {
       
  1122   var i,c,k,kk;
       
  1123   k=x.length<y.length ? x.length : y.length;
       
  1124   kk=x.length;
       
  1125   for (c=0,i=0;i<k;i++) {
       
  1126     c+=a*x[i]+b*y[i];
       
  1127     x[i]=c & mask;
       
  1128     c>>=bpe;
       
  1129   }
       
  1130   for (i=k;i<kk;i++) {
       
  1131     c+=a*x[i];
       
  1132     x[i]=c & mask;
       
  1133     c>>=bpe;
       
  1134   }
       
  1135 }
       
  1136 
       
  1137 //do the linear combination x=a*x+b*(y<<(ys*bpe)) for bigInts x and y, and integers a, b and ys.
       
  1138 //x must be large enough to hold the answer.
       
  1139 function linCombShift_(x,y,b,ys) {
       
  1140   var i,c,k,kk;
       
  1141   k=x.length<ys+y.length ? x.length : ys+y.length;
       
  1142   kk=x.length;
       
  1143   for (c=0,i=ys;i<k;i++) {
       
  1144     c+=x[i]+b*y[i-ys];
       
  1145     x[i]=c & mask;
       
  1146     c>>=bpe;
       
  1147   }
       
  1148   for (i=k;c && i<kk;i++) {
       
  1149     c+=x[i];
       
  1150     x[i]=c & mask;
       
  1151     c>>=bpe;
       
  1152   }
       
  1153 }
       
  1154 
       
  1155 //do x=x+(y<<(ys*bpe)) for bigInts x and y, and integers a,b and ys.
       
  1156 //x must be large enough to hold the answer.
       
  1157 function addShift_(x,y,ys) {
       
  1158   var i,c,k,kk;
       
  1159   k=x.length<ys+y.length ? x.length : ys+y.length;
       
  1160   kk=x.length;
       
  1161   for (c=0,i=ys;i<k;i++) {
       
  1162     c+=x[i]+y[i-ys];
       
  1163     x[i]=c & mask;
       
  1164     c>>=bpe;
       
  1165   }
       
  1166   for (i=k;c && i<kk;i++) {
       
  1167     c+=x[i];
       
  1168     x[i]=c & mask;
       
  1169     c>>=bpe;
       
  1170   }
       
  1171 }
       
  1172 
       
  1173 //do x=x-(y<<(ys*bpe)) for bigInts x and y, and integers a,b and ys.
       
  1174 //x must be large enough to hold the answer.
       
  1175 function subShift_(x,y,ys) {
       
  1176   var i,c,k,kk;
       
  1177   k=x.length<ys+y.length ? x.length : ys+y.length;
       
  1178   kk=x.length;
       
  1179   for (c=0,i=ys;i<k;i++) {
       
  1180     c+=x[i]-y[i-ys];
       
  1181     x[i]=c & mask;
       
  1182     c>>=bpe;
       
  1183   }
       
  1184   for (i=k;c && i<kk;i++) {
       
  1185     c+=x[i];
       
  1186     x[i]=c & mask;
       
  1187     c>>=bpe;
       
  1188   }
       
  1189 }
       
  1190 
       
  1191 //do x=x-y for bigInts x and y.
       
  1192 //x must be large enough to hold the answer.
       
  1193 //negative answers will be 2s complement
       
  1194 function sub_(x,y) {
       
  1195   var i,c,k,kk;
       
  1196   k=x.length<y.length ? x.length : y.length;
       
  1197   for (c=0,i=0;i<k;i++) {
       
  1198     c+=x[i]-y[i];
       
  1199     x[i]=c & mask;
       
  1200     c>>=bpe;
       
  1201   }
       
  1202   for (i=k;c && i<x.length;i++) {
       
  1203     c+=x[i];
       
  1204     x[i]=c & mask;
       
  1205     c>>=bpe;
       
  1206   }
       
  1207 }
       
  1208 
       
  1209 //do x=x+y for bigInts x and y.
       
  1210 //x must be large enough to hold the answer.
       
  1211 function add_(x,y) {
       
  1212   var i,c,k,kk;
       
  1213   k=x.length<y.length ? x.length : y.length;
       
  1214   for (c=0,i=0;i<k;i++) {
       
  1215     c+=x[i]+y[i];
       
  1216     x[i]=c & mask;
       
  1217     c>>=bpe;
       
  1218   }
       
  1219   for (i=k;c && i<x.length;i++) {
       
  1220     c+=x[i];
       
  1221     x[i]=c & mask;
       
  1222     c>>=bpe;
       
  1223   }
       
  1224 }
       
  1225 
       
  1226 //do x=x*y for bigInts x and y.  This is faster when y<x.
       
  1227 function mult_(x,y) {
       
  1228   var i;
       
  1229   if (ss.length!=2*x.length)
       
  1230     ss=new Array(2*x.length);
       
  1231   copyInt_(ss,0);
       
  1232   for (i=0;i<y.length;i++)
       
  1233     if (y[i])
       
  1234       linCombShift_(ss,x,y[i],i);   //ss=1*ss+y[i]*(x<<(i*bpe))
       
  1235   copy_(x,ss);
       
  1236 }
       
  1237 
       
  1238 //do x=x mod n for bigInts x and n.
       
  1239 function mod_(x,n) {
       
  1240   if (s4.length!=x.length)
       
  1241     s4=dup(x);
       
  1242   else
       
  1243     copy_(s4,x);
       
  1244   if (s5.length!=x.length)
       
  1245     s5=dup(x);  
       
  1246   divide_(s4,n,s5,x);  //x = remainder of s4 / n
       
  1247 }
       
  1248 
       
  1249 //do x=x*y mod n for bigInts x,y,n.
       
  1250 //for greater speed, let y<x.
       
  1251 function multMod_(x,y,n) {
       
  1252   var i;
       
  1253   if (s0.length!=2*x.length)
       
  1254     s0=new Array(2*x.length);
       
  1255   copyInt_(s0,0);
       
  1256   for (i=0;i<y.length;i++)
       
  1257     if (y[i])
       
  1258       linCombShift_(s0,x,y[i],i);   //s0=1*s0+y[i]*(x<<(i*bpe))
       
  1259   mod_(s0,n);
       
  1260   copy_(x,s0);
       
  1261 }
       
  1262 
       
  1263 //do x=x*x mod n for bigInts x,n.
       
  1264 function squareMod_(x,n) {
       
  1265   var i,j,d,c,kx,kn,k;
       
  1266   for (kx=x.length; kx>0 && !x[kx-1]; kx--);  //ignore leading zeros in x
       
  1267   k=kx>n.length ? 2*kx : 2*n.length; //k=# elements in the product, which is twice the elements in the larger of x and n
       
  1268   if (s0.length!=k) 
       
  1269     s0=new Array(k);
       
  1270   copyInt_(s0,0);
       
  1271   for (i=0;i<kx;i++) {
       
  1272     c=s0[2*i]+x[i]*x[i];
       
  1273     s0[2*i]=c & mask;
       
  1274     c>>=bpe;
       
  1275     for (j=i+1;j<kx;j++) {
       
  1276       c=s0[i+j]+2*x[i]*x[j]+c;
       
  1277       s0[i+j]=(c & mask);
       
  1278       c>>=bpe;
       
  1279     }
       
  1280     s0[i+kx]=c;
       
  1281   }
       
  1282   mod_(s0,n);
       
  1283   copy_(x,s0);
       
  1284 }
       
  1285 
       
  1286 //return x with exactly k leading zero elements
       
  1287 function bigint_trim(x,k) {
       
  1288   var i,y;
       
  1289   for (i=x.length; i>0 && !x[i-1]; i--);
       
  1290   y=new Array(i+k);
       
  1291   copy_(y,x);
       
  1292   return y;
       
  1293 }
       
  1294 
       
  1295 //do x=x**y mod n, where x,y,n are bigInts and ** is exponentiation.  0**0=1.
       
  1296 //this is faster when n is odd.  x usually needs to have as many elements as n.
       
  1297 function powMod_(x,y,n) {
       
  1298   var k1,k2,kn,np;
       
  1299   if(s7.length!=n.length)
       
  1300     s7=dup(n);
       
  1301 
       
  1302   //for even modulus, use a simple square-and-multiply algorithm,
       
  1303   //rather than using the more complex Montgomery algorithm.
       
  1304   if ((n[0]&1)==0) {
       
  1305     copy_(s7,x);
       
  1306     copyInt_(x,1);
       
  1307     while(!equalsInt(y,0)) {
       
  1308       if (y[0]&1)
       
  1309         multMod_(x,s7,n);
       
  1310       divInt_(y,2);
       
  1311       squareMod_(s7,n); 
       
  1312     }
       
  1313     return;
       
  1314   }
       
  1315 
       
  1316   //calculate np from n for the Montgomery multiplications
       
  1317   copyInt_(s7,0);
       
  1318   for (kn=n.length;kn>0 && !n[kn-1];kn--);
       
  1319   np=radix-inverseModInt(modInt(n,radix),radix);
       
  1320   s7[kn]=1;
       
  1321   multMod_(x ,s7,n);   // x = x * 2**(kn*bp) mod n
       
  1322 
       
  1323   if (s3.length!=x.length)
       
  1324     s3=dup(x);
       
  1325   else
       
  1326     copy_(s3,x);
       
  1327 
       
  1328   for (k1=y.length-1;k1>0 & !y[k1]; k1--);  //k1=first nonzero element of y
       
  1329   if (y[k1]==0) {  //anything to the 0th power is 1
       
  1330     copyInt_(x,1);
       
  1331     return;
       
  1332   }
       
  1333   for (k2=1<<(bpe-1);k2 && !(y[k1] & k2); k2>>=1);  //k2=position of first 1 bit in y[k1]
       
  1334   for (;;) {
       
  1335     if (!(k2>>=1)) {  //look at next bit of y
       
  1336       k1--;
       
  1337       if (k1<0) {
       
  1338         mont_(x,one,n,np);
       
  1339         return;
       
  1340       }
       
  1341       k2=1<<(bpe-1);
       
  1342     }    
       
  1343     mont_(x,x,n,np);
       
  1344 
       
  1345     if (k2 & y[k1]) //if next bit is a 1
       
  1346       mont_(x,s3,n,np);
       
  1347   }
       
  1348 }    
       
  1349 
       
  1350 //do x=x*y*Ri mod n for bigInts x,y,n, 
       
  1351 //  where Ri = 2**(-kn*bpe) mod n, and kn is the 
       
  1352 //  number of elements in the n array, not 
       
  1353 //  counting leading zeros.  
       
  1354 //x must be large enough to hold the answer.
       
  1355 //It's OK if x and y are the same variable.
       
  1356 //must have:
       
  1357 //  x,y < n
       
  1358 //  n is odd
       
  1359 //  np = -(n^(-1)) mod radix
       
  1360 function mont_(x,y,n,np) {
       
  1361   var i,j,c,ui,t;
       
  1362   var kn=n.length;
       
  1363   var ky=y.length;
       
  1364 
       
  1365   if (sa.length!=kn)
       
  1366     sa=new Array(kn);
       
  1367 
       
  1368   for (;kn>0 && n[kn-1]==0;kn--); //ignore leading zeros of n
       
  1369   //this function sometimes gives wrong answers when the next line is uncommented
       
  1370   //for (;ky>0 && y[ky-1]==0;ky--); //ignore leading zeros of y
       
  1371 
       
  1372   copyInt_(sa,0);
       
  1373 
       
  1374   //the following loop consumes 95% of the runtime for randTruePrime_() and powMod_() for large keys
       
  1375   for (i=0; i<kn; i++) {
       
  1376     t=sa[0]+x[i]*y[0];
       
  1377     ui=((t & mask) * np) & mask;  //the inner "& mask" is needed on Macintosh MSIE, but not windows MSIE
       
  1378     c=(t+ui*n[0]) >> bpe;
       
  1379     t=x[i];
       
  1380 
       
  1381     //do sa=(sa+x[i]*y+ui*n)/b   where b=2**bpe
       
  1382     for (j=1;j<ky;j++) { 
       
  1383       c+=sa[j]+t*y[j]+ui*n[j];
       
  1384       sa[j-1]=c & mask;
       
  1385       c>>=bpe;
       
  1386     }    
       
  1387     for (;j<kn;j++) { 
       
  1388       c+=sa[j]+ui*n[j];
       
  1389       sa[j-1]=c & mask;
       
  1390       c>>=bpe;
       
  1391     }    
       
  1392     sa[j-1]=c & mask;
       
  1393   }
       
  1394 
       
  1395   if (!greater(n,sa))
       
  1396     sub_(sa,n);
       
  1397   copy_(x,sa);
       
  1398 }
       
  1399 
       
  1400