KickJava   Java API By Example, From Geeks To Geeks.

Java > Open Source Codes > gnu > math > MPN


1 // Copyright (c) 1997 Per M.A. Bothner.
2
// This is free software; for terms and warranty disclaimer see ./COPYING.
3

4 package gnu.math;
5
6 /** This contains various low-level routines for unsigned bigints.
7  * The interfaces match the mpn interfaces in gmp,
8  * so it should be easy to replace them with fast native functions
9  * that are trivial wrappers around the mpn_ functions in gmp
10  * (at least on platforms that use 32-bit "limbs").
11  */

12
13 class MPN
14 {
15   /** Add x[0:size-1] and y, and write the size least
16    * significant words of the result to dest.
17    * Return carry, either 0 or 1.
18    * All values are unsigned.
19    * This is basically the same as gmp's mpn_add_1. */

20   public static int add_1 (int[] dest, int[] x, int size, int y)
21   {
22     long carry = (long) y & 0xffffffffL;
23     for (int i = 0; i < size; i++)
24       {
25     carry += ((long) x[i] & 0xffffffffL);
26     dest[i] = (int) carry;
27     carry >>= 32;
28       }
29     return (int) carry;
30   }
31
32   /** Add x[0:len-1] and y[0:len-1] and write the len least
33    * significant words of the result to dest[0:len-1].
34    * All words are treated as unsigned.
35    * @return the carry, either 0 or 1
36    * This function is basically the same as gmp's mpn_add_n.
37    */

38   public static int add_n (int dest[], int[] x, int[] y, int len)
39   {
40     long carry = 0;
41     for (int i = 0; i < len; i++)
42       {
43     carry += ((long) x[i] & 0xffffffffL)
44       + ((long) y[i] & 0xffffffffL);
45     dest[i] = (int) carry;
46     carry >>>= 32;
47       }
48     return (int) carry;
49   }
50
51   /** Subtract Y[0:size-1] from X[0:size-1], and write
52    * the size least significant words of the result to dest[0:size-1].
53    * Return borrow, either 0 or 1.
54    * This is basically the same as gmp's mpn_sub_n function.
55    */

56
57   public static int sub_n (int[] dest, int[] X, int[] Y, int size)
58   {
59     int cy = 0;
60     for (int i = 0; i < size; i++)
61       {
62     int y = Y[i];
63     int x = X[i];
64     y += cy; /* add previous carry to subtrahend */
65     // Invert the high-order bit, because: (unsigned) X > (unsigned) Y
66
// iff: (int) (X^0x80000000) > (int) (Y^0x80000000).
67
cy = (y^0x80000000) < (cy^0x80000000) ? 1 : 0;
68     y = x - y;
69     cy += (y^0x80000000) > (x ^ 0x80000000) ? 1 : 0;
70     dest[i] = y;
71       }
72     return cy;
73   }
74
75   /** Multiply x[0:len-1] by y, and write the len least
76    * significant words of the product to dest[0:len-1].
77    * Return the most significant word of the product.
78    * All values are treated as if they were unsigned
79    * (i.e. masked with 0xffffffffL).
80    * OK if dest==x (not sure if this is guaranteed for mpn_mul_1).
81    * This function is basically the same as gmp's mpn_mul_1.
82    */

83
84   public static int mul_1 (int[] dest, int[] x, int len, int y)
85   {
86     long yword = (long) y & 0xffffffffL;
87     long carry = 0;
88     for (int j = 0; j < len; j++)
89       {
90         carry += ((long) x[j] & 0xffffffffL) * yword;
91         dest[j] = (int) carry;
92         carry >>>= 32;
93       }
94     return (int) carry;
95   }
96
97   /**
98    * Multiply x[0:xlen-1] and y[0:ylen-1], and
99    * write the result to dest[0:xlen+ylen-1].
100    * The destination has to have space for xlen+ylen words,
101    * even if the result might be one limb smaller.
102    * This function requires that xlen >= ylen.
103    * The destination must be distinct from either input operands.
104    * All operands are unsigned.
105    * This function is basically the same gmp's mpn_mul. */

106
107   public static void mul (int[] dest,
108               int[] x, int xlen,
109               int[] y, int ylen)
110   {
111     dest[xlen] = MPN.mul_1 (dest, x, xlen, y[0]);
112
113     for (int i = 1; i < ylen; i++)
114       {
115     long yword = (long) y[i] & 0xffffffffL;
116     long carry = 0;
117     for (int j = 0; j < xlen; j++)
118       {
119         carry += ((long) x[j] & 0xffffffffL) * yword
120           + ((long) dest[i+j] & 0xffffffffL);
121         dest[i+j] = (int) carry;
122         carry >>>= 32;
123       }
124     dest[i+xlen] = (int) carry;
125       }
126   }
127
128   /* Divide (unsigned long) N by (unsigned int) D.
129    * Returns (remainder << 32)+(unsigned int)(quotient).
130    * Assumes (unsigned int)(N>>32) < (unsigned int)D.
131    * Code transcribed from gmp-2.0's mpn_udiv_w_sdiv function.
132    */

133   public static long udiv_qrnnd (long N, int D)
134   {
135     long q, r;
136     long a1 = N >>> 32;
137     long a0 = N & 0xffffffffL;
138     if (D >= 0)
139       {
140     if (a1 < ((D - a1 - (a0 >>> 31)) & 0xffffffffL))
141       {
142         /* dividend, divisor, and quotient are nonnegative */
143         q = N / D;
144         r = N % D;
145       }
146     else
147       {
148         /* Compute c1*2^32 + c0 = a1*2^32 + a0 - 2^31*d */
149         long c = N - ((long) D << 31);
150         /* Divide (c1*2^32 + c0) by d */
151         q = c / D;
152         r = c % D;
153         /* Add 2^31 to quotient */
154         q += 1 << 31;
155       }
156       }
157     else
158       {
159     long b1 = D >>> 1; /* d/2, between 2^30 and 2^31 - 1 */
160     //long c1 = (a1 >> 1); /* A/2 */
161
//int c0 = (a1 << 31) + (a0 >> 1);
162
long c = N >>> 1;
163     if (a1 < b1 || (a1 >> 1) < b1)
164       {
165         if (a1 < b1)
166           {
167         q = c / b1;
168         r = c % b1;
169           }
170         else /* c1 < b1, so 2^31 <= (A/2)/b1 < 2^32 */
171           {
172         c = ~(c - (b1 << 32));
173         q = c / b1; /* (A/2) / (d/2) */
174         r = c % b1;
175         q = (~q) & 0xffffffffL; /* (A/2)/b1 */
176         r = (b1 - 1) - r; /* r < b1 => new r >= 0 */
177           }
178         r = 2 * r + (a0 & 1);
179         if ((D & 1) != 0)
180           {
181         if (r >= q) {
182                 r = r - q;
183         } else if (q - r <= ((long) D & 0xffffffffL)) {
184                        r = r - q + D;
185                 q -= 1;
186         } else {
187                        r = r - q + D + D;
188                 q -= 2;
189         }
190           }
191       }
192     else /* Implies c1 = b1 */
193       { /* Hence a1 = d - 1 = 2*b1 - 1 */
194         if (a0 >= ((long)(-D) & 0xffffffffL))
195           {
196         q = -1;
197             r = a0 + D;
198           }
199         else
200           {
201         q = -2;
202             r = a0 + D + D;
203           }
204       }
205       }
206
207     return (r << 32) | (q & 0xFFFFFFFFl);
208   }
209
210     /** Divide divident[0:len-1] by (unsigned int)divisor.
211      * Write result into quotient[0:len-1].
212      * Return the one-word (unsigned) remainder.
213      * OK for quotient==dividend.
214      */

215
216   public static int divmod_1 (int[] quotient, int[] dividend,
217                   int len, int divisor)
218   {
219     int i = len - 1;
220     long r = dividend[i];
221     if ((r & 0xffffffffL) >= ((long)divisor & 0xffffffffL))
222       r = 0;
223     else
224       {
225     quotient[i--] = 0;
226     r <<= 32;
227       }
228
229     for (; i >= 0; i--)
230       {
231     int n0 = dividend[i];
232     r = (r & ~0xffffffffL) | (n0 & 0xffffffffL);
233     r = udiv_qrnnd (r, divisor);
234     quotient[i] = (int) r;
235       }
236     return (int)(r >> 32);
237   }
238
239   /* Subtract x[0:len-1]*y from dest[offset:offset+len-1].
240    * All values are treated as if unsigned.
241    * @return the most significant word of
242    * the product, minus borrow-out from the subtraction.
243    */

244   public static int submul_1 (int[] dest, int offset, int[] x, int len, int y)
245   {
246     long yl = (long) y & 0xffffffffL;
247     int carry = 0;
248     int j = 0;
249     do
250       {
251     long prod = ((long) x[j] & 0xffffffffL) * yl;
252     int prod_low = (int) prod;
253     int prod_high = (int) (prod >> 32);
254     prod_low += carry;
255     // Invert the high-order bit, because: (unsigned) X > (unsigned) Y
256
// iff: (int) (X^0x80000000) > (int) (Y^0x80000000).
257
carry = ((prod_low ^ 0x80000000) < (carry ^ 0x80000000) ? 1 : 0)
258       + prod_high;
259     int x_j = dest[offset+j];
260     prod_low = x_j - prod_low;
261     if ((prod_low ^ 0x80000000) > (x_j ^ 0x80000000))
262       carry++;
263     dest[offset+j] = prod_low;
264       }
265     while (++j < len);
266     return carry;
267   }
268
269   /** Divide zds[0:nx] by y[0:ny-1].
270    * The remainder ends up in zds[0:ny-1].
271    * The quotient ends up in zds[ny:nx].
272    * Assumes: nx>ny.
273    * (int)y[ny-1] < 0 (i.e. most significant bit set)
274    */

275
276   public static void divide (int[] zds, int nx, int[] y, int ny)
277   {
278     // This is basically Knuth's formulation of the classical algorithm,
279
// but translated from in scm_divbigbig in Jaffar's SCM implementation.
280

281     // Correspondance with Knuth's notation:
282
// Knuth's u[0:m+n] == zds[nx:0].
283
// Knuth's v[1:n] == y[ny-1:0]
284
// Knuth's n == ny.
285
// Knuth's m == nx-ny.
286
// Our nx == Knuth's m+n.
287

288     // Could be re-implemented using gmp's mpn_divrem:
289
// zds[nx] = mpn_divrem (&zds[ny], 0, zds, nx, y, ny).
290

291     int j = nx;
292     do
293       { // loop over digits of quotient
294
// Knuth's j == our nx-j.
295
// Knuth's u[j:j+n] == our zds[j:j-ny].
296
int qhat; // treated as unsigned
297
if (zds[j]==y[ny-1])
298       qhat = -1; // 0xffffffff
299
else
300       {
301         long w = (((long)(zds[j])) << 32) + ((long)zds[j-1] & 0xffffffffL);
302         qhat = (int) udiv_qrnnd (w, y[ny-1]);
303       }
304     if (qhat != 0)
305       {
306         int borrow = submul_1 (zds, j - ny, y, ny, qhat);
307         int save = zds[j];
308         long num = ((long)save&0xffffffffL) - ((long)borrow&0xffffffffL);
309             while (num != 0)
310           {
311         qhat--;
312         long carry = 0;
313         for (int i = 0; i < ny; i++)
314           {
315             carry += ((long) zds[j-ny+i] & 0xffffffffL)
316               + ((long) y[i] & 0xffffffffL);
317             zds[j-ny+i] = (int) carry;
318             carry >>>= 32;
319           }
320         zds[j] += carry;
321         num = carry - 1;
322           }
323       }
324     zds[j] = qhat;
325       } while (--j >= ny);
326   }
327
328   /** Number of digits in the conversion base that always fits in a word.
329    * For example, for base 10 this is 9, since 10**9 is the
330    * largest number that fits into a words (assuming 32-bit words).
331    * This is the same as gmp's __mp_bases[radix].chars_per_limb.
332    * @param radix the base
333    * @return number of digits */

334   public static int chars_per_word (int radix)
335   {
336     if (radix < 10)
337       {
338     if (radix < 8)
339       {
340         if (radix <= 2)
341           return 32;
342         else if (radix == 3)
343           return 20;
344         else if (radix == 4)
345           return 16;
346         else
347           return 18 - radix;
348       }
349     else
350       return 10;
351       }
352     else if (radix < 12)
353       return 9;
354     else if (radix <= 16)
355       return 8;
356     else if (radix <= 23)
357       return 7;
358     else if (radix <= 40)
359       return 6;
360     // The following are conservative, but we don't care.
361
else if (radix <= 256)
362       return 4;
363     else
364       return 1;
365   }
366
367   /** Count the number of leading zero bits in an int. */
368   public static int count_leading_zeros (int i)
369   {
370     if (i == 0)
371       return 32;
372     int count = 0;
373     for (int k = 16; k > 0; k = k >> 1) {
374       int j = i >>> k;
375       if (j == 0)
376     count += k;
377       else
378     i = j;
379     }
380     return count;
381   }
382
383   public static int set_str (int dest[], byte[] str, int str_len, int base)
384   {
385     int size = 0;
386     if ((base & (base - 1)) == 0)
387       {
388     // The base is a power of 2. Read the input string from
389
// least to most significant character/digit. */
390

391     int next_bitpos = 0;
392     int bits_per_indigit = 0;
393     for (int i = base; (i >>= 1) != 0; ) bits_per_indigit++;
394     int res_digit = 0;
395
396     for (int i = str_len; --i >= 0; )
397       {
398         int inp_digit = str[i];
399         res_digit |= inp_digit << next_bitpos;
400         next_bitpos += bits_per_indigit;
401         if (next_bitpos >= 32)
402           {
403         dest[size++] = res_digit;
404         next_bitpos -= 32;
405         res_digit = inp_digit >> (bits_per_indigit - next_bitpos);
406           }
407       }
408
409     if (res_digit != 0)
410       dest[size++] = res_digit;
411       }
412     else
413       {
414     // General case. The base is not a power of 2.
415
int indigits_per_limb = MPN.chars_per_word (base);
416     int str_pos = 0;
417
418     while (str_pos < str_len)
419       {
420         int chunk = str_len - str_pos;
421         if (chunk > indigits_per_limb)
422           chunk = indigits_per_limb;
423         int res_digit = str[str_pos++];
424         int big_base = base;
425
426         while (--chunk > 0)
427           {
428         res_digit = res_digit * base + str[str_pos++];
429         big_base *= base;
430           }
431
432         int cy_limb;
433         if (size == 0)
434           cy_limb = res_digit;
435         else
436           {
437         cy_limb = MPN.mul_1 (dest, dest, size, big_base);
438         cy_limb += MPN.add_1 (dest, dest, size, res_digit);
439           }
440         if (cy_limb != 0)
441           dest[size++] = cy_limb;
442       }
443        }
444     return size;
445   }
446
447   /** Compare {@code x[0:size-1]} with {@code y[0:size-1]}, treating them as unsigned integers.
448    * @return -1, 0, or 1 depending on if {@code x<y}, {@code x==y}, or {@code x>y}.
449    * This is basically the same as gmp's {@code mpn_cmp} function.
450    */

451   public static int cmp (int[] x, int[] y, int size)
452   {
453     while (--size >= 0)
454       {
455     int x_word = x[size];
456     int y_word = y[size];
457     if (x_word != y_word)
458       {
459         // Invert the high-order bit, because:
460
// (unsigned) X > (unsigned) Y iff
461
// (int) (X^0x80000000) > (int) (Y^0x80000000).
462
return (x_word ^ 0x80000000) > (y_word ^0x80000000) ? 1 : -1;
463       }
464       }
465     return 0;
466   }
467
468   /** Compare {@code x[0:xlen-1]} with {@code y[0:ylen-1]}, treating them as unsigned integers.
469    * @return -1, 0, or 1 depending on
470    * whether {@code x<y}, {@code x==y}, or {@code x>y}.
471    */

472   public static int cmp (int[] x, int xlen, int[] y, int ylen)
473   {
474     return xlen > ylen ? 1 : xlen < ylen ? -1 : cmp (x, y, xlen);
475   }
476
477   /* Shift x[x_start:x_start+len-1] count bits to the "right"
478    * (i.e. divide by 2**count).
479    * Store the len least significant words of the result at dest.
480    * The bits shifted out to the right are returned.
481    * OK if dest==x.
482    * Assumes: 0 < count < 32
483    */

484
485   public static int rshift (int[] dest, int[] x, int x_start,
486                 int len, int count)
487   {
488     int count_2 = 32 - count;
489     int low_word = x[x_start];
490     int retval = low_word << count_2;
491     int i = 1;
492     for (; i < len; i++)
493       {
494     int high_word = x[x_start+i];
495     dest[i-1] = (low_word >>> count) | (high_word << count_2);
496     low_word = high_word;
497       }
498     dest[i-1] = low_word >>> count;
499     return retval;
500   }
501
502   /* Shift x[x_start:x_start+len-1] count bits to the "right"
503    * (i.e. divide by 2**count).
504    * Store the len least significant words of the result at dest.
505    * OK if dest==x.
506    * Assumes: 0 <= count < 32
507    * Same as rshift, but handles count==0 (and has no return value).
508    */

509   public static void rshift0 (int[] dest, int[] x, int x_start,
510                   int len, int count)
511   {
512     if (count > 0)
513       rshift(dest, x, x_start, len, count);
514     else
515       for (int i = 0; i < len; i++)
516     dest[i] = x[i + x_start];
517   }
518
519   /** Return the long-truncated value of right shifting.
520   * @param x a two's-complement "bignum"
521   * @param len the number of significant words in x
522   * @param count the shift count
523   * @return (long)(x[0..len-1] >> count).
524   */

525   public static long rshift_long (int[] x, int len, int count)
526   {
527     int wordno = count >> 5;
528     count &= 31;
529     int sign = x[len-1] < 0 ? -1 : 0;
530     int w0 = wordno >= len ? sign : x[wordno];
531     wordno++;
532     int w1 = wordno >= len ? sign : x[wordno];
533     if (count != 0)
534       {
535     wordno++;
536     int w2 = wordno >= len ? sign : x[wordno];
537     w0 = (w0 >>> count) | (w1 << (32-count));
538     w1 = (w1 >>> count) | (w2 << (32-count));
539       }
540     return ((long)w1 << 32) | ((long)w0 & 0xffffffffL);
541   }
542
543   /* Shift x[0:len-1] left by count bits, and store the len least
544    * significant words of the result in dest[d_offset:d_offset+len-1].
545    * Return the bits shifted out from the most significant digit.
546    * Assumes 0 < count < 32.
547    * OK if dest==x.
548    */

549
550   public static int lshift (int[] dest, int d_offset,
551                 int[] x, int len, int count)
552   {
553     int count_2 = 32 - count;
554     int i = len - 1;
555     int high_word = x[i];
556     int retval = high_word >>> count_2;
557     d_offset++;
558     while (--i >= 0)
559       {
560     int low_word = x[i];
561     dest[d_offset+i] = (high_word << count) | (low_word >>> count_2);
562     high_word = low_word;
563       }
564     dest[d_offset+i] = high_word << count;
565     return retval;
566   }
567
568   /** Return least i such that word&(1<<i). Assumes word!=0. */
569
570   static int findLowestBit (int word)
571   {
572     int i = 0;
573     while ((word & 0xF) == 0)
574       {
575     word >>= 4;
576     i += 4;
577       }
578     if ((word & 3) == 0)
579       {
580     word >>= 2;
581     i += 2;
582       }
583     if ((word & 1) == 0)
584       i += 1;
585     return i;
586   }
587
588   /** Return least i such that words & (1<<i). Assumes there is such an i. */
589
590   static int findLowestBit (int[] words)
591   {
592     for (int i = 0; ; i++)
593       {
594     if (words[i] != 0)
595       return 32 * i + findLowestBit (words[i]);
596       }
597   }
598
599   /** Calculate Greatest Common Divisior of x[0:len-1] and y[0:len-1].
600     * Assumes both arguments are non-zero.
601     * Leaves result in x, and returns len of result.
602     * Also destroys y (actually sets it to a copy of the result). */

603
604   public static int gcd (int[] x, int[] y, int len)
605   {
606     int i, word;
607     // Find sh such that both x and y are divisible by 2**sh.
608
for (i = 0; ; i++)
609       {
610     word = x[i] | y[i];
611     if (word != 0)
612       {
613         // Must terminate, since x and y are non-zero.
614
break;
615       }
616       }
617     int initShiftWords = i;
618     int initShiftBits = findLowestBit (word);
619     // Logically: sh = initShiftWords * 32 + initShiftBits
620

621     // Temporarily devide both x and y by 2**sh.
622
len -= initShiftWords;
623     MPN.rshift0 (x, x, initShiftWords, len, initShiftBits);
624     MPN.rshift0 (y, y, initShiftWords, len, initShiftBits);
625
626     int[] odd_arg; /* One of x or y which is odd. */
627     int[] other_arg; /* The other one can be even or odd. */
628     if ((x[0] & 1) != 0)
629       {
630     odd_arg = x;
631     other_arg = y;
632       }
633     else
634       {
635     odd_arg = y;
636     other_arg = x;
637       }
638
639     for (;;)
640       {
641     // Shift other_arg until it is odd; this doesn't
642
// affect the gcd, since we divide by 2**k, which does not
643
// divide odd_arg.
644
for (i = 0; other_arg[i] == 0; ) i++;
645     if (i > 0)
646       {
647         int j;
648         for (j = 0; j < len-i; j++)
649         other_arg[j] = other_arg[j+i];
650         for ( ; j < len; j++)
651           other_arg[j] = 0;
652       }
653     i = findLowestBit(other_arg[0]);
654     if (i > 0)
655       MPN.rshift (other_arg, other_arg, 0, len, i);
656
657     // Now both odd_arg and other_arg are odd.
658

659     // Subtract the smaller from the larger.
660
// This does not change the result, since gcd(a-b,b)==gcd(a,b).
661
i = MPN.cmp(odd_arg, other_arg, len);
662     if (i == 0)
663         break;
664     if (i > 0)
665       { // odd_arg > other_arg
666
MPN.sub_n (odd_arg, odd_arg, other_arg, len);
667         // Now odd_arg is even, so swap with other_arg;
668
int[] tmp = odd_arg; odd_arg = other_arg; other_arg = tmp;
669       }
670     else
671       { // other_arg > odd_arg
672
MPN.sub_n (other_arg, other_arg, odd_arg, len);
673     }
674     while (odd_arg[len-1] == 0 && other_arg[len-1] == 0)
675       len--;
676     }
677     if (initShiftWords + initShiftBits > 0)
678       {
679     if (initShiftBits > 0)
680       {
681         int sh_out = MPN.lshift (x, initShiftWords, x, len, initShiftBits);
682         if (sh_out != 0)
683           x[(len++)+initShiftWords] = sh_out;
684       }
685     else
686       {
687         for (i = len; --i >= 0;)
688           x[i+initShiftWords] = x[i];
689       }
690     for (i = initShiftWords; --i >= 0; )
691       x[i] = 0;
692     len += initShiftWords;
693       }
694     return len;
695   }
696
697   public static int intLength (int i)
698   {
699     return 32 - count_leading_zeros (i < 0 ? ~i : i);
700   }
701
702   /** Calcaulte the Common Lisp "integer-length" function.
703    * Assumes input is canonicalized: len==IntNum.wordsNeeded(words,len) */

704   public static int intLength (int[] words, int len)
705   {
706     len--;
707     return intLength (words[len]) + 32 * len;
708   }
709
710   /* DEBUGGING:
711   public static void dprint (IntNum x)
712   {
713     if (x.words == null)
714       System.err.print(Long.toString((long) x.ival & 0xffffffffL, 16));
715     else
716       dprint (System.err, x.words, x.ival);
717   }
718   public static void dprint (int[] x) { dprint (System.err, x, x.length); }
719   public static void dprint (int[] x, int len) { dprint (System.err, x, len); }
720   public static void dprint (java.io.PrintStream ps, int[] x, int len)
721   {
722     ps.print('(');
723     for (int i = 0; i < len; i++)
724       {
725     if (i > 0)
726       ps.print (' ');
727     ps.print ("#x" + Long.toString ((long) x[i] & 0xffffffffL, 16));
728       }
729     ps.print(')');
730   }
731   */

732 }
733
Popular Tags