OpenCores
URL https://opencores.org/ocsvn/openrisc/openrisc/trunk

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [gnu-dev/] [or1k-gcc/] [libjava/] [classpath/] [gnu/] [java/] [math/] [MPN.java] - Blame information for rev 769

Details | Compare with Previous | View Log

Line No. Rev Author Line
1 769 jeremybenn
/* gnu.java.math.MPN
2
   Copyright (C) 1999, 2000, 2001, 2004  Free Software Foundation, Inc.
3
 
4
This file is part of GNU Classpath.
5
 
6
GNU Classpath is free software; you can redistribute it and/or modify
7
it under the terms of the GNU General Public License as published by
8
the Free Software Foundation; either version 2, or (at your option)
9
any later version.
10
 
11
GNU Classpath is distributed in the hope that it will be useful, but
12
WITHOUT ANY WARRANTY; without even the implied warranty of
13
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
14
General Public License for more details.
15
 
16
You should have received a copy of the GNU General Public License
17
along with GNU Classpath; see the file COPYING.  If not, write to the
18
Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
19
02110-1301 USA.
20
 
21
Linking this library statically or dynamically with other modules is
22
making a combined work based on this library.  Thus, the terms and
23
conditions of the GNU General Public License cover the whole
24
combination.
25
 
26
As a special exception, the copyright holders of this library give you
27
permission to link this library with independent modules to produce an
28
executable, regardless of the license terms of these independent
29
modules, and to copy and distribute the resulting executable under
30
terms of your choice, provided that you also meet, for each linked
31
independent module, the terms and conditions of the license of that
32
module.  An independent module is a module which is not derived from
33
or based on this library.  If you modify this library, you may extend
34
this exception to your version of the library, but you are not
35
obligated to do so.  If you do not wish to do so, delete this
36
exception statement from your version. */
37
 
38
// Included from Kawa 1.6.62 with permission of the author,
39
// Per Bothner <per@bothner.com>.
40
 
41
package gnu.java.math;
42
 
43
/** This contains various low-level routines for unsigned bigints.
44
 * The interfaces match the mpn interfaces in gmp,
45
 * so it should be easy to replace them with fast native functions
46
 * that are trivial wrappers around the mpn_ functions in gmp
47
 * (at least on platforms that use 32-bit "limbs").
48
 */
49
 
50
public class MPN
51
{
52
  /** Add x[0:size-1] and y, and write the size least
53
   * significant words of the result to dest.
54
   * Return carry, either 0 or 1.
55
   * All values are unsigned.
56
   * This is basically the same as gmp's mpn_add_1. */
57
  public static int add_1 (int[] dest, int[] x, int size, int y)
58
  {
59
    long carry = (long) y & 0xffffffffL;
60
    for (int i = 0;  i < size;  i++)
61
      {
62
        carry += ((long) x[i] & 0xffffffffL);
63
        dest[i] = (int) carry;
64
        carry >>= 32;
65
      }
66
    return (int) carry;
67
  }
68
 
69
  /** Add x[0:len-1] and y[0:len-1] and write the len least
70
   * significant words of the result to dest[0:len-1].
71
   * All words are treated as unsigned.
72
   * @return the carry, either 0 or 1
73
   * This function is basically the same as gmp's mpn_add_n.
74
   */
75
  public static int add_n (int dest[], int[] x, int[] y, int len)
76
  {
77
    long carry = 0;
78
    for (int i = 0; i < len;  i++)
79
      {
80
        carry += ((long) x[i] & 0xffffffffL)
81
          + ((long) y[i] & 0xffffffffL);
82
        dest[i] = (int) carry;
83
        carry >>>= 32;
84
      }
85
    return (int) carry;
86
  }
87
 
88
  /** Subtract Y[0:size-1] from X[0:size-1], and write
89
   * the size least significant words of the result to dest[0:size-1].
90
   * Return borrow, either 0 or 1.
91
   * This is basically the same as gmp's mpn_sub_n function.
92
   */
93
 
94
  public static int sub_n (int[] dest, int[] X, int[] Y, int size)
95
  {
96
    int cy = 0;
97
    for (int i = 0;  i < size;  i++)
98
      {
99
        int y = Y[i];
100
        int x = X[i];
101
        y += cy;        /* add previous carry to subtrahend */
102
        // Invert the high-order bit, because: (unsigned) X > (unsigned) Y
103
        // iff: (int) (X^0x80000000) > (int) (Y^0x80000000).
104
        cy = (y^0x80000000) < (cy^0x80000000) ? 1 : 0;
105
        y = x - y;
106
        cy += (y^0x80000000) > (x ^ 0x80000000) ? 1 : 0;
107
        dest[i] = y;
108
      }
109
    return cy;
110
  }
111
 
112
  /** Multiply x[0:len-1] by y, and write the len least
113
   * significant words of the product to dest[0:len-1].
114
   * Return the most significant word of the product.
115
   * All values are treated as if they were unsigned
116
   * (i.e. masked with 0xffffffffL).
117
   * OK if dest==x (not sure if this is guaranteed for mpn_mul_1).
118
   * This function is basically the same as gmp's mpn_mul_1.
119
   */
120
 
121
  public static int mul_1 (int[] dest, int[] x, int len, int y)
122
  {
123
    long yword = (long) y & 0xffffffffL;
124
    long carry = 0;
125
    for (int j = 0;  j < len; j++)
126
      {
127
        carry += ((long) x[j] & 0xffffffffL) * yword;
128
        dest[j] = (int) carry;
129
        carry >>>= 32;
130
      }
131
    return (int) carry;
132
  }
133
 
134
  /**
135
   * Multiply x[0:xlen-1] and y[0:ylen-1], and
136
   * write the result to dest[0:xlen+ylen-1].
137
   * The destination has to have space for xlen+ylen words,
138
   * even if the result might be one limb smaller.
139
   * This function requires that xlen >= ylen.
140
   * The destination must be distinct from either input operands.
141
   * All operands are unsigned.
142
   * This function is basically the same gmp's mpn_mul. */
143
 
144
  public static void mul (int[] dest,
145
                          int[] x, int xlen,
146
                          int[] y, int ylen)
147
  {
148
    dest[xlen] = MPN.mul_1 (dest, x, xlen, y[0]);
149
 
150
    for (int i = 1;  i < ylen; i++)
151
      {
152
        long yword = (long) y[i] & 0xffffffffL;
153
        long carry = 0;
154
        for (int j = 0;  j < xlen; j++)
155
          {
156
            carry += ((long) x[j] & 0xffffffffL) * yword
157
              + ((long) dest[i+j] & 0xffffffffL);
158
            dest[i+j] = (int) carry;
159
            carry >>>= 32;
160
          }
161
        dest[i+xlen] = (int) carry;
162
      }
163
  }
164
 
165
  /* Divide (unsigned long) N by (unsigned int) D.
166
   * Returns (remainder << 32)+(unsigned int)(quotient).
167
   * Assumes (unsigned int)(N>>32) < (unsigned int)D.
168
   * Code transcribed from gmp-2.0's mpn_udiv_w_sdiv function.
169
   */
170
  public static long udiv_qrnnd (long N, int D)
171
  {
172
    long q, r;
173
    long a1 = N >>> 32;
174
    long a0 = N & 0xffffffffL;
175
    if (D >= 0)
176
      {
177
        if (a1 < ((D - a1 - (a0 >>> 31)) & 0xffffffffL))
178
          {
179
            /* dividend, divisor, and quotient are nonnegative */
180
            q = N / D;
181
            r = N % D;
182
          }
183
        else
184
          {
185
            /* Compute c1*2^32 + c0 = a1*2^32 + a0 - 2^31*d */
186
            long c = N - ((long) D << 31);
187
            /* Divide (c1*2^32 + c0) by d */
188
            q = c / D;
189
            r = c % D;
190
            /* Add 2^31 to quotient */
191
            q += 1 << 31;
192
          }
193
      }
194
    else
195
      {
196
        long b1 = D >>> 1;      /* d/2, between 2^30 and 2^31 - 1 */
197
        //long c1 = (a1 >> 1); /* A/2 */
198
        //int c0 = (a1 << 31) + (a0 >> 1);
199
        long c = N >>> 1;
200
        if (a1 < b1 || (a1 >> 1) < b1)
201
          {
202
            if (a1 < b1)
203
              {
204
                q = c / b1;
205
                r = c % b1;
206
              }
207
            else /* c1 < b1, so 2^31 <= (A/2)/b1 < 2^32 */
208
              {
209
                c = ~(c - (b1 << 32));
210
                q = c / b1;  /* (A/2) / (d/2) */
211
                r = c % b1;
212
                q = (~q) & 0xffffffffL;    /* (A/2)/b1 */
213
                r = (b1 - 1) - r; /* r < b1 => new r >= 0 */
214
              }
215
            r = 2 * r + (a0 & 1);
216
            if ((D & 1) != 0)
217
              {
218
                if (r >= q) {
219
                        r = r - q;
220
                } else if (q - r <= ((long) D & 0xffffffffL)) {
221
                       r = r - q + D;
222
                        q -= 1;
223
                } else {
224
                       r = r - q + D + D;
225
                        q -= 2;
226
                }
227
              }
228
          }
229
        else                            /* Implies c1 = b1 */
230
          {                             /* Hence a1 = d - 1 = 2*b1 - 1 */
231
            if (a0 >= ((long)(-D) & 0xffffffffL))
232
              {
233
                q = -1;
234
                r = a0 + D;
235
              }
236
            else
237
              {
238
                q = -2;
239
                r = a0 + D + D;
240
              }
241
          }
242
      }
243
 
244
    return (r << 32) | (q & 0xFFFFFFFFl);
245
  }
246
 
247
    /** Divide divident[0:len-1] by (unsigned int)divisor.
248
     * Write result into quotient[0:len-1.
249
     * Return the one-word (unsigned) remainder.
250
     * OK for quotient==dividend.
251
     */
252
 
253
  public static int divmod_1 (int[] quotient, int[] dividend,
254
                              int len, int divisor)
255
  {
256
    int i = len - 1;
257
    long r = dividend[i];
258
    if ((r & 0xffffffffL) >= ((long)divisor & 0xffffffffL))
259
      r = 0;
260
    else
261
      {
262
        quotient[i--] = 0;
263
        r <<= 32;
264
      }
265
 
266
    for (;  i >= 0;  i--)
267
      {
268
        int n0 = dividend[i];
269
        r = (r & ~0xffffffffL) | (n0 & 0xffffffffL);
270
        r = udiv_qrnnd (r, divisor);
271
        quotient[i] = (int) r;
272
      }
273
    return (int)(r >> 32);
274
  }
275
 
276
  /* Subtract x[0:len-1]*y from dest[offset:offset+len-1].
277
   * All values are treated as if unsigned.
278
   * @return the most significant word of
279
   * the product, minus borrow-out from the subtraction.
280
   */
281
  public static int submul_1 (int[] dest, int offset, int[] x, int len, int y)
282
  {
283
    long yl = (long) y & 0xffffffffL;
284
    int carry = 0;
285
    int j = 0;
286
    do
287
      {
288
        long prod = ((long) x[j] & 0xffffffffL) * yl;
289
        int prod_low = (int) prod;
290
        int prod_high = (int) (prod >> 32);
291
        prod_low += carry;
292
        // Invert the high-order bit, because: (unsigned) X > (unsigned) Y
293
        // iff: (int) (X^0x80000000) > (int) (Y^0x80000000).
294
        carry = ((prod_low ^ 0x80000000) < (carry ^ 0x80000000) ? 1 : 0)
295
          + prod_high;
296
        int x_j = dest[offset+j];
297
        prod_low = x_j - prod_low;
298
        if ((prod_low ^ 0x80000000) > (x_j ^ 0x80000000))
299
          carry++;
300
        dest[offset+j] = prod_low;
301
      }
302
    while (++j < len);
303
    return carry;
304
  }
305
 
306
  /** Divide zds[0:nx] by y[0:ny-1].
307
   * The remainder ends up in zds[0:ny-1].
308
   * The quotient ends up in zds[ny:nx].
309
   * Assumes:  nx>ny.
310
   * (int)y[ny-1] < 0  (i.e. most significant bit set)
311
   */
312
 
313
  public static void divide (int[] zds, int nx, int[] y, int ny)
314
  {
315
    // This is basically Knuth's formulation of the classical algorithm,
316
    // but translated from in scm_divbigbig in Jaffar's SCM implementation.
317
 
318
    // Correspondance with Knuth's notation:
319
    // Knuth's u[0:m+n] == zds[nx:0].
320
    // Knuth's v[1:n] == y[ny-1:0]
321
    // Knuth's n == ny.
322
    // Knuth's m == nx-ny.
323
    // Our nx == Knuth's m+n.
324
 
325
    // Could be re-implemented using gmp's mpn_divrem:
326
    // zds[nx] = mpn_divrem (&zds[ny], 0, zds, nx, y, ny).
327
 
328
    int j = nx;
329
    do
330
      {                          // loop over digits of quotient
331
        // Knuth's j == our nx-j.
332
        // Knuth's u[j:j+n] == our zds[j:j-ny].
333
        int qhat;  // treated as unsigned
334
        if (zds[j]==y[ny-1])
335
          qhat = -1;  // 0xffffffff
336
        else
337
          {
338
            long w = (((long)(zds[j])) << 32) + ((long)zds[j-1] & 0xffffffffL);
339
            qhat = (int) udiv_qrnnd (w, y[ny-1]);
340
          }
341
        if (qhat != 0)
342
          {
343
            int borrow = submul_1 (zds, j - ny, y, ny, qhat);
344
            int save = zds[j];
345
            long num = ((long)save&0xffffffffL) - ((long)borrow&0xffffffffL);
346
            while (num != 0)
347
              {
348
                qhat--;
349
                long carry = 0;
350
                for (int i = 0;  i < ny; i++)
351
                  {
352
                    carry += ((long) zds[j-ny+i] & 0xffffffffL)
353
                      + ((long) y[i] & 0xffffffffL);
354
                    zds[j-ny+i] = (int) carry;
355
                    carry >>>= 32;
356
                  }
357
                zds[j] += carry;
358
                num = carry - 1;
359
              }
360
          }
361
        zds[j] = qhat;
362
      } while (--j >= ny);
363
  }
364
 
365
  /** Number of digits in the conversion base that always fits in a word.
366
   * For example, for base 10 this is 9, since 10**9 is the
367
   * largest number that fits into a words (assuming 32-bit words).
368
   * This is the same as gmp's __mp_bases[radix].chars_per_limb.
369
   * @param radix the base
370
   * @return number of digits */
371
  public static int chars_per_word (int radix)
372
  {
373
    if (radix < 10)
374
      {
375
        if (radix < 8)
376
          {
377
            if (radix <= 2)
378
              return 32;
379
            else if (radix == 3)
380
              return 20;
381
            else if (radix == 4)
382
              return 16;
383
            else
384
              return 18 - radix;
385
          }
386
        else
387
          return 10;
388
      }
389
    else if (radix < 12)
390
      return 9;
391
    else if (radix <= 16)
392
      return 8;
393
    else if (radix <= 23)
394
      return 7;
395
    else if (radix <= 40)
396
      return 6;
397
    // The following are conservative, but we don't care.
398
    else if (radix <= 256)
399
      return 4;
400
    else
401
      return 1;
402
  }
403
 
404
  /** Count the number of leading zero bits in an int. */
405
  public static int count_leading_zeros (int i)
406
  {
407
    if (i == 0)
408
      return 32;
409
    int count = 0;
410
    for (int k = 16;  k > 0;  k = k >> 1) {
411
      int j = i >>> k;
412
      if (j == 0)
413
        count += k;
414
      else
415
        i = j;
416
    }
417
    return count;
418
  }
419
 
420
  public static int set_str (int dest[], byte[] str, int str_len, int base)
421
  {
422
    int size = 0;
423
    if ((base & (base - 1)) == 0)
424
      {
425
        // The base is a power of 2.  Read the input string from
426
        // least to most significant character/digit.  */
427
 
428
        int next_bitpos = 0;
429
        int bits_per_indigit = 0;
430
        for (int i = base; (i >>= 1) != 0; ) bits_per_indigit++;
431
        int res_digit = 0;
432
 
433
        for (int i = str_len;  --i >= 0; )
434
          {
435
            int inp_digit = str[i];
436
            res_digit |= inp_digit << next_bitpos;
437
            next_bitpos += bits_per_indigit;
438
            if (next_bitpos >= 32)
439
              {
440
                dest[size++] = res_digit;
441
                next_bitpos -= 32;
442
                res_digit = inp_digit >> (bits_per_indigit - next_bitpos);
443
              }
444
          }
445
 
446
        if (res_digit != 0)
447
          dest[size++] = res_digit;
448
      }
449
    else
450
      {
451
        // General case.  The base is not a power of 2.
452
        int indigits_per_limb = MPN.chars_per_word (base);
453
        int str_pos = 0;
454
 
455
        while (str_pos < str_len)
456
          {
457
            int chunk = str_len - str_pos;
458
            if (chunk > indigits_per_limb)
459
              chunk = indigits_per_limb;
460
            int res_digit = str[str_pos++];
461
            int big_base = base;
462
 
463
            while (--chunk > 0)
464
              {
465
                res_digit = res_digit * base + str[str_pos++];
466
                big_base *= base;
467
              }
468
 
469
            int cy_limb;
470
            if (size == 0)
471
              cy_limb = res_digit;
472
            else
473
              {
474
                cy_limb = MPN.mul_1 (dest, dest, size, big_base);
475
                cy_limb += MPN.add_1 (dest, dest, size, res_digit);
476
              }
477
            if (cy_limb != 0)
478
              dest[size++] = cy_limb;
479
          }
480
       }
481
    return size;
482
  }
483
 
484
  /** Compare x[0:size-1] with y[0:size-1], treating them as unsigned integers.
485
   * @result -1, 0, or 1 depending on if x&lt;y, x==y, or x&gt;y.
486
   * This is basically the same as gmp's mpn_cmp function.
487
   */
488
  public static int cmp (int[] x, int[] y, int size)
489
  {
490
    while (--size >= 0)
491
      {
492
        int x_word = x[size];
493
        int y_word = y[size];
494
        if (x_word != y_word)
495
          {
496
            // Invert the high-order bit, because:
497
            // (unsigned) X > (unsigned) Y iff
498
            // (int) (X^0x80000000) > (int) (Y^0x80000000).
499
            return (x_word ^ 0x80000000) > (y_word ^0x80000000) ? 1 : -1;
500
          }
501
      }
502
    return 0;
503
  }
504
 
505
  /**
506
   * Compare x[0:xlen-1] with y[0:ylen-1], treating them as unsigned integers.
507
   *
508
   * @return -1, 0, or 1 depending on if x&lt;y, x==y, or x&gt;y.
509
   */
510
  public static int cmp (int[] x, int xlen, int[] y, int ylen)
511
  {
512
    return xlen > ylen ? 1 : xlen < ylen ? -1 : cmp (x, y, xlen);
513
  }
514
 
515
  /**
516
   * Shift x[x_start:x_start+len-1] count bits to the "right"
517
   * (i.e. divide by 2**count).
518
   * Store the len least significant words of the result at dest.
519
   * The bits shifted out to the right are returned.
520
   * OK if dest==x.
521
   * Assumes: 0 &lt; count &lt; 32
522
   */
523
  public static int rshift (int[] dest, int[] x, int x_start,
524
                            int len, int count)
525
  {
526
    int count_2 = 32 - count;
527
    int low_word = x[x_start];
528
    int retval = low_word << count_2;
529
    int i = 1;
530
    for (; i < len;  i++)
531
      {
532
        int high_word = x[x_start+i];
533
        dest[i-1] = (low_word >>> count) | (high_word << count_2);
534
        low_word = high_word;
535
      }
536
    dest[i-1] = low_word >>> count;
537
    return retval;
538
  }
539
 
540
  /**
541
   * Shift x[x_start:x_start+len-1] count bits to the "right"
542
   * (i.e. divide by 2**count).
543
   * Store the len least significant words of the result at dest.
544
   * OK if dest==x.
545
   * Assumes: 0 &lt;= count &lt; 32
546
   * Same as rshift, but handles count==0 (and has no return value).
547
   */
548
  public static void rshift0 (int[] dest, int[] x, int x_start,
549
                              int len, int count)
550
  {
551
    if (count > 0)
552
      rshift(dest, x, x_start, len, count);
553
    else
554
      for (int i = 0;  i < len;  i++)
555
        dest[i] = x[i + x_start];
556
  }
557
 
558
  /** Return the long-truncated value of right shifting.
559
  * @param x a two's-complement "bignum"
560
  * @param len the number of significant words in x
561
  * @param count the shift count
562
  * @return (long)(x[0..len-1] &gt;&gt; count).
563
  */
564
  public static long rshift_long (int[] x, int len, int count)
565
  {
566
    int wordno = count >> 5;
567
    count &= 31;
568
    int sign = x[len-1] < 0 ? -1 : 0;
569
    int w0 = wordno >= len ? sign : x[wordno];
570
    wordno++;
571
    int w1 = wordno >= len ? sign : x[wordno];
572
    if (count != 0)
573
      {
574
        wordno++;
575
        int w2 = wordno >= len ? sign : x[wordno];
576
        w0 = (w0 >>> count) | (w1 << (32-count));
577
        w1 = (w1 >>> count) | (w2 << (32-count));
578
      }
579
    return ((long)w1 << 32) | ((long)w0 & 0xffffffffL);
580
  }
581
 
582
  /* Shift x[0:len-1] left by count bits, and store the len least
583
   * significant words of the result in dest[d_offset:d_offset+len-1].
584
   * Return the bits shifted out from the most significant digit.
585
   * Assumes 0 &lt; count &lt; 32.
586
   * OK if dest==x.
587
   */
588
 
589
  public static int lshift (int[] dest, int d_offset,
590
                            int[] x, int len, int count)
591
  {
592
    int count_2 = 32 - count;
593
    int i = len - 1;
594
    int high_word = x[i];
595
    int retval = high_word >>> count_2;
596
    d_offset++;
597
    while (--i >= 0)
598
      {
599
        int low_word = x[i];
600
        dest[d_offset+i] = (high_word << count) | (low_word >>> count_2);
601
        high_word = low_word;
602
      }
603
    dest[d_offset+i] = high_word << count;
604
    return retval;
605
  }
606
 
607
  /** Return least i such that word &amp; (1&lt;&lt;i). Assumes word!=0. */
608
 
609
  public static int findLowestBit (int word)
610
  {
611
    int i = 0;
612
    while ((word & 0xF) == 0)
613
      {
614
        word >>= 4;
615
        i += 4;
616
      }
617
    if ((word & 3) == 0)
618
      {
619
        word >>= 2;
620
        i += 2;
621
      }
622
    if ((word & 1) == 0)
623
      i += 1;
624
    return i;
625
  }
626
 
627
  /** Return least i such that words &amp; (1&lt;&lt;i). Assumes there is such an i. */
628
 
629
  public static int findLowestBit (int[] words)
630
  {
631
    for (int i = 0;  ; i++)
632
      {
633
        if (words[i] != 0)
634
          return 32 * i + findLowestBit (words[i]);
635
      }
636
  }
637
 
638
  /** Calculate Greatest Common Divisior of x[0:len-1] and y[0:len-1].
639
    * Assumes both arguments are non-zero.
640
    * Leaves result in x, and returns len of result.
641
    * Also destroys y (actually sets it to a copy of the result). */
642
 
643
  public static int gcd (int[] x, int[] y, int len)
644
  {
645
    int i, word;
646
    // Find sh such that both x and y are divisible by 2**sh.
647
    for (i = 0; ; i++)
648
      {
649
        word = x[i] | y[i];
650
        if (word != 0)
651
          {
652
            // Must terminate, since x and y are non-zero.
653
            break;
654
          }
655
      }
656
    int initShiftWords = i;
657
    int initShiftBits = findLowestBit (word);
658
    // Logically: sh = initShiftWords * 32 + initShiftBits
659
 
660
    // Temporarily devide both x and y by 2**sh.
661
    len -= initShiftWords;
662
    MPN.rshift0 (x, x, initShiftWords, len, initShiftBits);
663
    MPN.rshift0 (y, y, initShiftWords, len, initShiftBits);
664
 
665
    int[] odd_arg; /* One of x or y which is odd. */
666
    int[] other_arg; /* The other one can be even or odd. */
667
    if ((x[0] & 1) != 0)
668
      {
669
        odd_arg = x;
670
        other_arg = y;
671
      }
672
    else
673
      {
674
        odd_arg = y;
675
        other_arg = x;
676
      }
677
 
678
    for (;;)
679
      {
680
        // Shift other_arg until it is odd; this doesn't
681
        // affect the gcd, since we divide by 2**k, which does not
682
        // divide odd_arg.
683
        for (i = 0; other_arg[i] == 0; ) i++;
684
        if (i > 0)
685
          {
686
            int j;
687
            for (j = 0; j < len-i; j++)
688
                other_arg[j] = other_arg[j+i];
689
            for ( ; j < len; j++)
690
              other_arg[j] = 0;
691
          }
692
        i = findLowestBit(other_arg[0]);
693
        if (i > 0)
694
          MPN.rshift (other_arg, other_arg, 0, len, i);
695
 
696
        // Now both odd_arg and other_arg are odd.
697
 
698
        // Subtract the smaller from the larger.
699
        // This does not change the result, since gcd(a-b,b)==gcd(a,b).
700
        i = MPN.cmp(odd_arg, other_arg, len);
701
        if (i == 0)
702
            break;
703
        if (i > 0)
704
          { // odd_arg > other_arg
705
            MPN.sub_n (odd_arg, odd_arg, other_arg, len);
706
            // Now odd_arg is even, so swap with other_arg;
707
            int[] tmp = odd_arg; odd_arg = other_arg; other_arg = tmp;
708
          }
709
        else
710
          { // other_arg > odd_arg
711
            MPN.sub_n (other_arg, other_arg, odd_arg, len);
712
        }
713
        while (odd_arg[len-1] == 0 && other_arg[len-1] == 0)
714
          len--;
715
    }
716
    if (initShiftWords + initShiftBits > 0)
717
      {
718
        if (initShiftBits > 0)
719
          {
720
            int sh_out = MPN.lshift (x, initShiftWords, x, len, initShiftBits);
721
            if (sh_out != 0)
722
              x[(len++)+initShiftWords] = sh_out;
723
          }
724
        else
725
          {
726
            for (i = len; --i >= 0;)
727
              x[i+initShiftWords] = x[i];
728
          }
729
        for (i = initShiftWords;  --i >= 0; )
730
          x[i] = 0;
731
        len += initShiftWords;
732
      }
733
    return len;
734
  }
735
 
736
  public static int intLength (int i)
737
  {
738
    return 32 - count_leading_zeros (i < 0 ? ~i : i);
739
  }
740
 
741
  /** Calcaulte the Common Lisp "integer-length" function.
742
   * Assumes input is canonicalized:  len==BigInteger.wordsNeeded(words,len) */
743
  public static int intLength (int[] words, int len)
744
  {
745
    len--;
746
    return intLength (words[len]) + 32 * len;
747
  }
748
 
749
  /* DEBUGGING:
750
  public static void dprint (BigInteger x)
751
  {
752
    if (x.words == null)
753
      System.err.print(Long.toString((long) x.ival & 0xffffffffL, 16));
754
    else
755
      dprint (System.err, x.words, x.ival);
756
  }
757
  public static void dprint (int[] x) { dprint (System.err, x, x.length); }
758
  public static void dprint (int[] x, int len) { dprint (System.err, x, len); }
759
  public static void dprint (java.io.PrintStream ps, int[] x, int len)
760
  {
761
    ps.print('(');
762
    for (int i = 0;  i < len; i++)
763
      {
764
        if (i > 0)
765
          ps.print (' ');
766
        ps.print ("#x" + Long.toString ((long) x[i] & 0xffffffffL, 16));
767
      }
768
    ps.print(')');
769
  }
770
  */
771
}

powered by: WebSVN 2.1.0

© copyright 1999-2024 OpenCores.org, equivalent to Oliscience, all rights reserved. OpenCores®, registered trademark.