001: // Copyright (c) 1997 Per M.A. Bothner.
002: // This is free software; for terms and warranty disclaimer see ./COPYING.
003:
004: package gnu.math;
005:
006: /** This contains various low-level routines for unsigned bigints.
007: * The interfaces match the mpn interfaces in gmp,
008: * so it should be easy to replace them with fast native functions
009: * that are trivial wrappers around the mpn_ functions in gmp
010: * (at least on platforms that use 32-bit "limbs").
011: */
012:
013: class MPN {
014: /** Add x[0:size-1] and y, and write the size least
015: * significant words of the result to dest.
016: * Return carry, either 0 or 1.
017: * All values are unsigned.
018: * This is basically the same as gmp's mpn_add_1. */
019: public static int add_1(int[] dest, int[] x, int size, int y) {
020: long carry = (long) y & 0xffffffffL;
021: for (int i = 0; i < size; i++) {
022: carry += ((long) x[i] & 0xffffffffL);
023: dest[i] = (int) carry;
024: carry >>= 32;
025: }
026: return (int) carry;
027: }
028:
029: /** Add x[0:len-1] and y[0:len-1] and write the len least
030: * significant words of the result to dest[0:len-1].
031: * All words are treated as unsigned.
032: * @return the carry, either 0 or 1
033: * This function is basically the same as gmp's mpn_add_n.
034: */
035: public static int add_n(int dest[], int[] x, int[] y, int len) {
036: long carry = 0;
037: for (int i = 0; i < len; i++) {
038: carry += ((long) x[i] & 0xffffffffL)
039: + ((long) y[i] & 0xffffffffL);
040: dest[i] = (int) carry;
041: carry >>>= 32;
042: }
043: return (int) carry;
044: }
045:
046: /** Subtract Y[0:size-1] from X[0:size-1], and write
047: * the size least significant words of the result to dest[0:size-1].
048: * Return borrow, either 0 or 1.
049: * This is basically the same as gmp's mpn_sub_n function.
050: */
051:
052: public static int sub_n(int[] dest, int[] X, int[] Y, int size) {
053: int cy = 0;
054: for (int i = 0; i < size; i++) {
055: int y = Y[i];
056: int x = X[i];
057: y += cy; /* add previous carry to subtrahend */
058: // Invert the high-order bit, because: (unsigned) X > (unsigned) Y
059: // iff: (int) (X^0x80000000) > (int) (Y^0x80000000).
060: cy = (y ^ 0x80000000) < (cy ^ 0x80000000) ? 1 : 0;
061: y = x - y;
062: cy += (y ^ 0x80000000) > (x ^ 0x80000000) ? 1 : 0;
063: dest[i] = y;
064: }
065: return cy;
066: }
067:
068: /** Multiply x[0:len-1] by y, and write the len least
069: * significant words of the product to dest[0:len-1].
070: * Return the most significant word of the product.
071: * All values are treated as if they were unsigned
072: * (i.e. masked with 0xffffffffL).
073: * OK if dest==x (not sure if this is guaranteed for mpn_mul_1).
074: * This function is basically the same as gmp's mpn_mul_1.
075: */
076:
077: public static int mul_1(int[] dest, int[] x, int len, int y) {
078: long yword = (long) y & 0xffffffffL;
079: long carry = 0;
080: for (int j = 0; j < len; j++) {
081: carry += ((long) x[j] & 0xffffffffL) * yword;
082: dest[j] = (int) carry;
083: carry >>>= 32;
084: }
085: return (int) carry;
086: }
087:
088: /**
089: * Multiply x[0:xlen-1] and y[0:ylen-1], and
090: * write the result to dest[0:xlen+ylen-1].
091: * The destination has to have space for xlen+ylen words,
092: * even if the result might be one limb smaller.
093: * This function requires that xlen >= ylen.
094: * The destination must be distinct from either input operands.
095: * All operands are unsigned.
096: * This function is basically the same gmp's mpn_mul. */
097:
098: public static void mul(int[] dest, int[] x, int xlen, int[] y,
099: int ylen) {
100: dest[xlen] = MPN.mul_1(dest, x, xlen, y[0]);
101:
102: for (int i = 1; i < ylen; i++) {
103: long yword = (long) y[i] & 0xffffffffL;
104: long carry = 0;
105: for (int j = 0; j < xlen; j++) {
106: carry += ((long) x[j] & 0xffffffffL) * yword
107: + ((long) dest[i + j] & 0xffffffffL);
108: dest[i + j] = (int) carry;
109: carry >>>= 32;
110: }
111: dest[i + xlen] = (int) carry;
112: }
113: }
114:
115: /* Divide (unsigned long) N by (unsigned int) D.
116: * Returns (remainder << 32)+(unsigned int)(quotient).
117: * Assumes (unsigned int)(N>>32) < (unsigned int)D.
118: * Code transcribed from gmp-2.0's mpn_udiv_w_sdiv function.
119: */
120: public static long udiv_qrnnd(long N, int D) {
121: long q, r;
122: long a1 = N >>> 32;
123: long a0 = N & 0xffffffffL;
124: if (D >= 0) {
125: if (a1 < ((D - a1 - (a0 >>> 31)) & 0xffffffffL)) {
126: /* dividend, divisor, and quotient are nonnegative */
127: q = N / D;
128: r = N % D;
129: } else {
130: /* Compute c1*2^32 + c0 = a1*2^32 + a0 - 2^31*d */
131: long c = N - ((long) D << 31);
132: /* Divide (c1*2^32 + c0) by d */
133: q = c / D;
134: r = c % D;
135: /* Add 2^31 to quotient */
136: q += 1 << 31;
137: }
138: } else {
139: long b1 = D >>> 1; /* d/2, between 2^30 and 2^31 - 1 */
140: //long c1 = (a1 >> 1); /* A/2 */
141: //int c0 = (a1 << 31) + (a0 >> 1);
142: long c = N >>> 1;
143: if (a1 < b1 || (a1 >> 1) < b1) {
144: if (a1 < b1) {
145: q = c / b1;
146: r = c % b1;
147: } else /* c1 < b1, so 2^31 <= (A/2)/b1 < 2^32 */
148: {
149: c = ~(c - (b1 << 32));
150: q = c / b1; /* (A/2) / (d/2) */
151: r = c % b1;
152: q = (~q) & 0xffffffffL; /* (A/2)/b1 */
153: r = (b1 - 1) - r; /* r < b1 => new r >= 0 */
154: }
155: r = 2 * r + (a0 & 1);
156: if ((D & 1) != 0) {
157: if (r >= q) {
158: r = r - q;
159: } else if (q - r <= ((long) D & 0xffffffffL)) {
160: r = r - q + D;
161: q -= 1;
162: } else {
163: r = r - q + D + D;
164: q -= 2;
165: }
166: }
167: } else /* Implies c1 = b1 */
168: { /* Hence a1 = d - 1 = 2*b1 - 1 */
169: if (a0 >= ((long) (-D) & 0xffffffffL)) {
170: q = -1;
171: r = a0 + D;
172: } else {
173: q = -2;
174: r = a0 + D + D;
175: }
176: }
177: }
178:
179: return (r << 32) | (q & 0xFFFFFFFFl);
180: }
181:
182: /** Divide divident[0:len-1] by (unsigned int)divisor.
183: * Write result into quotient[0:len-1].
184: * Return the one-word (unsigned) remainder.
185: * OK for quotient==dividend.
186: */
187:
188: public static int divmod_1(int[] quotient, int[] dividend, int len,
189: int divisor) {
190: int i = len - 1;
191: long r = dividend[i];
192: if ((r & 0xffffffffL) >= ((long) divisor & 0xffffffffL))
193: r = 0;
194: else {
195: quotient[i--] = 0;
196: r <<= 32;
197: }
198:
199: for (; i >= 0; i--) {
200: int n0 = dividend[i];
201: r = (r & ~0xffffffffL) | (n0 & 0xffffffffL);
202: r = udiv_qrnnd(r, divisor);
203: quotient[i] = (int) r;
204: }
205: return (int) (r >> 32);
206: }
207:
208: /* Subtract x[0:len-1]*y from dest[offset:offset+len-1].
209: * All values are treated as if unsigned.
210: * @return the most significant word of
211: * the product, minus borrow-out from the subtraction.
212: */
213: public static int submul_1(int[] dest, int offset, int[] x,
214: int len, int y) {
215: long yl = (long) y & 0xffffffffL;
216: int carry = 0;
217: int j = 0;
218: do {
219: long prod = ((long) x[j] & 0xffffffffL) * yl;
220: int prod_low = (int) prod;
221: int prod_high = (int) (prod >> 32);
222: prod_low += carry;
223: // Invert the high-order bit, because: (unsigned) X > (unsigned) Y
224: // iff: (int) (X^0x80000000) > (int) (Y^0x80000000).
225: carry = ((prod_low ^ 0x80000000) < (carry ^ 0x80000000) ? 1
226: : 0)
227: + prod_high;
228: int x_j = dest[offset + j];
229: prod_low = x_j - prod_low;
230: if ((prod_low ^ 0x80000000) > (x_j ^ 0x80000000))
231: carry++;
232: dest[offset + j] = prod_low;
233: } while (++j < len);
234: return carry;
235: }
236:
237: /** Divide zds[0:nx] by y[0:ny-1].
238: * The remainder ends up in zds[0:ny-1].
239: * The quotient ends up in zds[ny:nx].
240: * Assumes: nx>ny.
241: * (int)y[ny-1] < 0 (i.e. most significant bit set)
242: */
243:
244: public static void divide(int[] zds, int nx, int[] y, int ny) {
245: // This is basically Knuth's formulation of the classical algorithm,
246: // but translated from in scm_divbigbig in Jaffar's SCM implementation.
247:
248: // Correspondance with Knuth's notation:
249: // Knuth's u[0:m+n] == zds[nx:0].
250: // Knuth's v[1:n] == y[ny-1:0]
251: // Knuth's n == ny.
252: // Knuth's m == nx-ny.
253: // Our nx == Knuth's m+n.
254:
255: // Could be re-implemented using gmp's mpn_divrem:
256: // zds[nx] = mpn_divrem (&zds[ny], 0, zds, nx, y, ny).
257:
258: int j = nx;
259: do { // loop over digits of quotient
260: // Knuth's j == our nx-j.
261: // Knuth's u[j:j+n] == our zds[j:j-ny].
262: int qhat; // treated as unsigned
263: if (zds[j] == y[ny - 1])
264: qhat = -1; // 0xffffffff
265: else {
266: long w = (((long) (zds[j])) << 32)
267: + ((long) zds[j - 1] & 0xffffffffL);
268: qhat = (int) udiv_qrnnd(w, y[ny - 1]);
269: }
270: if (qhat != 0) {
271: int borrow = submul_1(zds, j - ny, y, ny, qhat);
272: int save = zds[j];
273: long num = ((long) save & 0xffffffffL)
274: - ((long) borrow & 0xffffffffL);
275: while (num != 0) {
276: qhat--;
277: long carry = 0;
278: for (int i = 0; i < ny; i++) {
279: carry += ((long) zds[j - ny + i] & 0xffffffffL)
280: + ((long) y[i] & 0xffffffffL);
281: zds[j - ny + i] = (int) carry;
282: carry >>>= 32;
283: }
284: zds[j] += carry;
285: num = carry - 1;
286: }
287: }
288: zds[j] = qhat;
289: } while (--j >= ny);
290: }
291:
292: /** Number of digits in the conversion base that always fits in a word.
293: * For example, for base 10 this is 9, since 10**9 is the
294: * largest number that fits into a words (assuming 32-bit words).
295: * This is the same as gmp's __mp_bases[radix].chars_per_limb.
296: * @param radix the base
297: * @return number of digits */
298: public static int chars_per_word(int radix) {
299: if (radix < 10) {
300: if (radix < 8) {
301: if (radix <= 2)
302: return 32;
303: else if (radix == 3)
304: return 20;
305: else if (radix == 4)
306: return 16;
307: else
308: return 18 - radix;
309: } else
310: return 10;
311: } else if (radix < 12)
312: return 9;
313: else if (radix <= 16)
314: return 8;
315: else if (radix <= 23)
316: return 7;
317: else if (radix <= 40)
318: return 6;
319: // The following are conservative, but we don't care.
320: else if (radix <= 256)
321: return 4;
322: else
323: return 1;
324: }
325:
326: /** Count the number of leading zero bits in an int. */
327: public static int count_leading_zeros(int i) {
328: if (i == 0)
329: return 32;
330: int count = 0;
331: for (int k = 16; k > 0; k = k >> 1) {
332: int j = i >>> k;
333: if (j == 0)
334: count += k;
335: else
336: i = j;
337: }
338: return count;
339: }
340:
341: public static int set_str(int dest[], byte[] str, int str_len,
342: int base) {
343: int size = 0;
344: if ((base & (base - 1)) == 0) {
345: // The base is a power of 2. Read the input string from
346: // least to most significant character/digit. */
347:
348: int next_bitpos = 0;
349: int bits_per_indigit = 0;
350: for (int i = base; (i >>= 1) != 0;)
351: bits_per_indigit++;
352: int res_digit = 0;
353:
354: for (int i = str_len; --i >= 0;) {
355: int inp_digit = str[i];
356: res_digit |= inp_digit << next_bitpos;
357: next_bitpos += bits_per_indigit;
358: if (next_bitpos >= 32) {
359: dest[size++] = res_digit;
360: next_bitpos -= 32;
361: res_digit = inp_digit >> (bits_per_indigit - next_bitpos);
362: }
363: }
364:
365: if (res_digit != 0)
366: dest[size++] = res_digit;
367: } else {
368: // General case. The base is not a power of 2.
369: int indigits_per_limb = MPN.chars_per_word(base);
370: int str_pos = 0;
371:
372: while (str_pos < str_len) {
373: int chunk = str_len - str_pos;
374: if (chunk > indigits_per_limb)
375: chunk = indigits_per_limb;
376: int res_digit = str[str_pos++];
377: int big_base = base;
378:
379: while (--chunk > 0) {
380: res_digit = res_digit * base + str[str_pos++];
381: big_base *= base;
382: }
383:
384: int cy_limb;
385: if (size == 0)
386: cy_limb = res_digit;
387: else {
388: cy_limb = MPN.mul_1(dest, dest, size, big_base);
389: cy_limb += MPN.add_1(dest, dest, size, res_digit);
390: }
391: if (cy_limb != 0)
392: dest[size++] = cy_limb;
393: }
394: }
395: return size;
396: }
397:
398: /** Compare {@code x[0:size-1]} with {@code y[0:size-1]}, treating them as unsigned integers.
399: * @return -1, 0, or 1 depending on if {@code x<y}, {@code x==y}, or {@code x>y}.
400: * This is basically the same as gmp's {@code mpn_cmp} function.
401: */
402: public static int cmp(int[] x, int[] y, int size) {
403: while (--size >= 0) {
404: int x_word = x[size];
405: int y_word = y[size];
406: if (x_word != y_word) {
407: // Invert the high-order bit, because:
408: // (unsigned) X > (unsigned) Y iff
409: // (int) (X^0x80000000) > (int) (Y^0x80000000).
410: return (x_word ^ 0x80000000) > (y_word ^ 0x80000000) ? 1
411: : -1;
412: }
413: }
414: return 0;
415: }
416:
417: /** Compare {@code x[0:xlen-1]} with {@code y[0:ylen-1]}, treating them as unsigned integers.
418: * @return -1, 0, or 1 depending on
419: * whether {@code x<y}, {@code x==y}, or {@code x>y}.
420: */
421: public static int cmp(int[] x, int xlen, int[] y, int ylen) {
422: return xlen > ylen ? 1 : xlen < ylen ? -1 : cmp(x, y, xlen);
423: }
424:
425: /* Shift x[x_start:x_start+len-1] count bits to the "right"
426: * (i.e. divide by 2**count).
427: * Store the len least significant words of the result at dest.
428: * The bits shifted out to the right are returned.
429: * OK if dest==x.
430: * Assumes: 0 < count < 32
431: */
432:
433: public static int rshift(int[] dest, int[] x, int x_start, int len,
434: int count) {
435: int count_2 = 32 - count;
436: int low_word = x[x_start];
437: int retval = low_word << count_2;
438: int i = 1;
439: for (; i < len; i++) {
440: int high_word = x[x_start + i];
441: dest[i - 1] = (low_word >>> count) | (high_word << count_2);
442: low_word = high_word;
443: }
444: dest[i - 1] = low_word >>> count;
445: return retval;
446: }
447:
448: /* Shift x[x_start:x_start+len-1] count bits to the "right"
449: * (i.e. divide by 2**count).
450: * Store the len least significant words of the result at dest.
451: * OK if dest==x.
452: * Assumes: 0 <= count < 32
453: * Same as rshift, but handles count==0 (and has no return value).
454: */
455: public static void rshift0(int[] dest, int[] x, int x_start,
456: int len, int count) {
457: if (count > 0)
458: rshift(dest, x, x_start, len, count);
459: else
460: for (int i = 0; i < len; i++)
461: dest[i] = x[i + x_start];
462: }
463:
464: /** Return the long-truncated value of right shifting.
465: * @param x a two's-complement "bignum"
466: * @param len the number of significant words in x
467: * @param count the shift count
468: * @return (long)(x[0..len-1] >> count).
469: */
470: public static long rshift_long(int[] x, int len, int count) {
471: int wordno = count >> 5;
472: count &= 31;
473: int sign = x[len - 1] < 0 ? -1 : 0;
474: int w0 = wordno >= len ? sign : x[wordno];
475: wordno++;
476: int w1 = wordno >= len ? sign : x[wordno];
477: if (count != 0) {
478: wordno++;
479: int w2 = wordno >= len ? sign : x[wordno];
480: w0 = (w0 >>> count) | (w1 << (32 - count));
481: w1 = (w1 >>> count) | (w2 << (32 - count));
482: }
483: return ((long) w1 << 32) | ((long) w0 & 0xffffffffL);
484: }
485:
486: /* Shift x[0:len-1] left by count bits, and store the len least
487: * significant words of the result in dest[d_offset:d_offset+len-1].
488: * Return the bits shifted out from the most significant digit.
489: * Assumes 0 < count < 32.
490: * OK if dest==x.
491: */
492:
493: public static int lshift(int[] dest, int d_offset, int[] x,
494: int len, int count) {
495: int count_2 = 32 - count;
496: int i = len - 1;
497: int high_word = x[i];
498: int retval = high_word >>> count_2;
499: d_offset++;
500: while (--i >= 0) {
501: int low_word = x[i];
502: dest[d_offset + i] = (high_word << count)
503: | (low_word >>> count_2);
504: high_word = low_word;
505: }
506: dest[d_offset + i] = high_word << count;
507: return retval;
508: }
509:
510: /** Return least i such that word&(1<<i). Assumes word!=0. */
511:
512: static int findLowestBit(int word) {
513: int i = 0;
514: while ((word & 0xF) == 0) {
515: word >>= 4;
516: i += 4;
517: }
518: if ((word & 3) == 0) {
519: word >>= 2;
520: i += 2;
521: }
522: if ((word & 1) == 0)
523: i += 1;
524: return i;
525: }
526:
527: /** Return least i such that words & (1<<i). Assumes there is such an i. */
528:
529: static int findLowestBit(int[] words) {
530: for (int i = 0;; i++) {
531: if (words[i] != 0)
532: return 32 * i + findLowestBit(words[i]);
533: }
534: }
535:
536: /** Calculate Greatest Common Divisior of x[0:len-1] and y[0:len-1].
537: * Assumes both arguments are non-zero.
538: * Leaves result in x, and returns len of result.
539: * Also destroys y (actually sets it to a copy of the result). */
540:
541: public static int gcd(int[] x, int[] y, int len) {
542: int i, word;
543: // Find sh such that both x and y are divisible by 2**sh.
544: for (i = 0;; i++) {
545: word = x[i] | y[i];
546: if (word != 0) {
547: // Must terminate, since x and y are non-zero.
548: break;
549: }
550: }
551: int initShiftWords = i;
552: int initShiftBits = findLowestBit(word);
553: // Logically: sh = initShiftWords * 32 + initShiftBits
554:
555: // Temporarily devide both x and y by 2**sh.
556: len -= initShiftWords;
557: MPN.rshift0(x, x, initShiftWords, len, initShiftBits);
558: MPN.rshift0(y, y, initShiftWords, len, initShiftBits);
559:
560: int[] odd_arg; /* One of x or y which is odd. */
561: int[] other_arg; /* The other one can be even or odd. */
562: if ((x[0] & 1) != 0) {
563: odd_arg = x;
564: other_arg = y;
565: } else {
566: odd_arg = y;
567: other_arg = x;
568: }
569:
570: for (;;) {
571: // Shift other_arg until it is odd; this doesn't
572: // affect the gcd, since we divide by 2**k, which does not
573: // divide odd_arg.
574: for (i = 0; other_arg[i] == 0;)
575: i++;
576: if (i > 0) {
577: int j;
578: for (j = 0; j < len - i; j++)
579: other_arg[j] = other_arg[j + i];
580: for (; j < len; j++)
581: other_arg[j] = 0;
582: }
583: i = findLowestBit(other_arg[0]);
584: if (i > 0)
585: MPN.rshift(other_arg, other_arg, 0, len, i);
586:
587: // Now both odd_arg and other_arg are odd.
588:
589: // Subtract the smaller from the larger.
590: // This does not change the result, since gcd(a-b,b)==gcd(a,b).
591: i = MPN.cmp(odd_arg, other_arg, len);
592: if (i == 0)
593: break;
594: if (i > 0) { // odd_arg > other_arg
595: MPN.sub_n(odd_arg, odd_arg, other_arg, len);
596: // Now odd_arg is even, so swap with other_arg;
597: int[] tmp = odd_arg;
598: odd_arg = other_arg;
599: other_arg = tmp;
600: } else { // other_arg > odd_arg
601: MPN.sub_n(other_arg, other_arg, odd_arg, len);
602: }
603: while (odd_arg[len - 1] == 0 && other_arg[len - 1] == 0)
604: len--;
605: }
606: if (initShiftWords + initShiftBits > 0) {
607: if (initShiftBits > 0) {
608: int sh_out = MPN.lshift(x, initShiftWords, x, len,
609: initShiftBits);
610: if (sh_out != 0)
611: x[(len++) + initShiftWords] = sh_out;
612: } else {
613: for (i = len; --i >= 0;)
614: x[i + initShiftWords] = x[i];
615: }
616: for (i = initShiftWords; --i >= 0;)
617: x[i] = 0;
618: len += initShiftWords;
619: }
620: return len;
621: }
622:
623: public static int intLength(int i) {
624: return 32 - count_leading_zeros(i < 0 ? ~i : i);
625: }
626:
627: /** Calcaulte the Common Lisp "integer-length" function.
628: * Assumes input is canonicalized: len==IntNum.wordsNeeded(words,len) */
629: public static int intLength(int[] words, int len) {
630: len--;
631: return intLength(words[len]) + 32 * len;
632: }
633:
634: /* DEBUGGING:
635: public static void dprint (IntNum x)
636: {
637: if (x.words == null)
638: System.err.print(Long.toString((long) x.ival & 0xffffffffL, 16));
639: else
640: dprint (System.err, x.words, x.ival);
641: }
642: public static void dprint (int[] x) { dprint (System.err, x, x.length); }
643: public static void dprint (int[] x, int len) { dprint (System.err, x, len); }
644: public static void dprint (java.io.PrintStream ps, int[] x, int len)
645: {
646: ps.print('(');
647: for (int i = 0; i < len; i++)
648: {
649: if (i > 0)
650: ps.print (' ');
651: ps.print ("#x" + Long.toString ((long) x[i] & 0xffffffffL, 16));
652: }
653: ps.print(')');
654: }
655: */
656: }
|