[BACK]Return to gcdext.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib / gmp / mpn / generic

Annotation of OpenXM_contrib/gmp/mpn/generic/gcdext.c, Revision 1.1.1.1

1.1       maekawa     1: /* mpn_gcdext -- Extended Greatest Common Divisor.
                      2:
                      3: Copyright (C) 1996 Free Software Foundation, Inc.
                      4:
                      5: This file is part of the GNU MP Library.
                      6:
                      7: The GNU MP Library is free software; you can redistribute it and/or modify
                      8: it under the terms of the GNU Library General Public License as published by
                      9: the Free Software Foundation; either version 2 of the License, or (at your
                     10: option) any later version.
                     11:
                     12: The GNU MP Library is distributed in the hope that it will be useful, but
                     13: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
                     14: or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
                     15: License for more details.
                     16:
                     17: You should have received a copy of the GNU Library General Public License
                     18: along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
                     19: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
                     20: MA 02111-1307, USA. */
                     21:
                     22: #include "gmp.h"
                     23: #include "gmp-impl.h"
                     24: #include "longlong.h"
                     25:
                     26: #ifndef EXTEND
                     27: #define EXTEND 1
                     28: #endif
                     29:
                     30: #if STAT
                     31: int arr[BITS_PER_MP_LIMB];
                     32: #endif
                     33:
                     34: #define SGN(A) (((A) < 0) ? -1 : ((A) > 0))
                     35:
                     36: /* Idea 1: After we have performed a full division, don't shift operands back,
                     37:           but instead account for the extra factors-of-2 thus introduced.
                     38:    Idea 2: Simple generalization to use divide-and-conquer would give us an
                     39:           algorithm that runs faster than O(n^2).
                     40:    Idea 3: The input numbers need less space as the computation progresses,
                     41:           while the s0 and s1 variables need more space.  To save space, we
                     42:           could make them share space, and have the latter variables grow
                     43:           into the former.  */
                     44:
                     45: /* Precondition: U >= V.  */
                     46:
                     47: mp_size_t
                     48: #if EXTEND
                     49: #if __STDC__
                     50: mpn_gcdext (mp_ptr gp, mp_ptr s0p,
                     51:            mp_ptr up, mp_size_t size, mp_ptr vp, mp_size_t vsize)
                     52: #else
                     53: mpn_gcdext (gp, s0p, up, size, vp, vsize)
                     54:      mp_ptr gp;
                     55:      mp_ptr s0p;
                     56:      mp_ptr up;
                     57:      mp_size_t size;
                     58:      mp_ptr vp;
                     59:      mp_size_t vsize;
                     60: #endif
                     61: #else
                     62: #if __STDC__
                     63: mpn_gcd (mp_ptr gp,
                     64:         mp_ptr up, mp_size_t size, mp_ptr vp, mp_size_t vsize)
                     65: #else
                     66: mpn_gcd (gp, up, size, vp, vsize)
                     67:      mp_ptr gp;
                     68:      mp_ptr up;
                     69:      mp_size_t size;
                     70:      mp_ptr vp;
                     71:      mp_size_t vsize;
                     72: #endif
                     73: #endif
                     74: {
                     75:   mp_limb_t uh, vh;
                     76:   mp_limb_signed_t A, B, C, D;
                     77:   int cnt;
                     78:   mp_ptr tp, wp;
                     79: #if RECORD
                     80:   mp_limb_signed_t min = 0, max = 0;
                     81: #endif
                     82: #if EXTEND
                     83:   mp_ptr s1p;
                     84:   mp_ptr orig_s0p = s0p;
                     85:   mp_size_t ssize, orig_size = size;
                     86:   TMP_DECL (mark);
                     87:
                     88:   TMP_MARK (mark);
                     89:
                     90:   tp = (mp_ptr) TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB);
                     91:   wp = (mp_ptr) TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB);
                     92:   s1p = (mp_ptr) TMP_ALLOC (size * BYTES_PER_MP_LIMB);
                     93:
                     94:   MPN_ZERO (s0p, size);
                     95:   MPN_ZERO (s1p, size);
                     96:
                     97:   s0p[0] = 1;
                     98:   s1p[0] = 0;
                     99:   ssize = 1;
                    100: #endif
                    101:
                    102:   if (size > vsize)
                    103:     {
                    104:       /* Normalize V (and shift up U the same amount).  */
                    105:       count_leading_zeros (cnt, vp[vsize - 1]);
                    106:       if (cnt != 0)
                    107:        {
                    108:          mp_limb_t cy;
                    109:          mpn_lshift (vp, vp, vsize, cnt);
                    110:          cy = mpn_lshift (up, up, size, cnt);
                    111:          up[size] = cy;
                    112:          size += cy != 0;
                    113:        }
                    114:
                    115:       mpn_divmod (up + vsize, up, size, vp, vsize);
                    116: #if EXTEND
                    117:       /* This is really what it boils down to in this case... */
                    118:       s0p[0] = 0;
                    119:       s1p[0] = 1;
                    120: #endif
                    121:       size = vsize;
                    122:       if (cnt != 0)
                    123:        {
                    124:          mpn_rshift (up, up, size, cnt);
                    125:          mpn_rshift (vp, vp, size, cnt);
                    126:        }
                    127:       {
                    128:        mp_ptr xp;
                    129:        xp = up; up = vp; vp = xp;
                    130:       }
                    131:     }
                    132:
                    133:   for (;;)
                    134:     {
                    135:       /* Figure out exact size of V.  */
                    136:       vsize = size;
                    137:       MPN_NORMALIZE (vp, vsize);
                    138:       if (vsize <= 1)
                    139:        break;
                    140:
                    141:       /* Make UH be the most significant limb of U, and make VH be
                    142:         corresponding bits from V.  */
                    143:       uh = up[size - 1];
                    144:       vh = vp[size - 1];
                    145:       count_leading_zeros (cnt, uh);
                    146:       if (cnt != 0)
                    147:        {
                    148:          uh = (uh << cnt) | (up[size - 2] >> (BITS_PER_MP_LIMB - cnt));
                    149:          vh = (vh << cnt) | (vp[size - 2] >> (BITS_PER_MP_LIMB - cnt));
                    150:        }
                    151:
                    152: #if 0
                    153:       /* For now, only handle BITS_PER_MP_LIMB-1 bits.  This makes
                    154:         room for sign bit.  */
                    155:       uh >>= 1;
                    156:       vh >>= 1;
                    157: #endif
                    158:       A = 1;
                    159:       B = 0;
                    160:       C = 0;
                    161:       D = 1;
                    162:
                    163:       for (;;)
                    164:        {
                    165:          mp_limb_signed_t q, T;
                    166:          if (vh + C == 0 || vh + D == 0)
                    167:            break;
                    168:
                    169:          q = (uh + A) / (vh + C);
                    170:          if (q != (uh + B) / (vh + D))
                    171:            break;
                    172:
                    173:          T = A - q * C;
                    174:          A = C;
                    175:          C = T;
                    176:          T = B - q * D;
                    177:          B = D;
                    178:          D = T;
                    179:          T = uh - q * vh;
                    180:          uh = vh;
                    181:          vh = T;
                    182:        }
                    183:
                    184: #if RECORD
                    185:       min = MIN (A, min);  min = MIN (B, min);
                    186:       min = MIN (C, min);  min = MIN (D, min);
                    187:       max = MAX (A, max);  max = MAX (B, max);
                    188:       max = MAX (C, max);  max = MAX (D, max);
                    189: #endif
                    190:
                    191:       if (B == 0)
                    192:        {
                    193:          mp_limb_t qh;
                    194:          mp_size_t i;
                    195:
                    196:          /* This is quite rare.  I.e., optimize something else!  */
                    197:
                    198:          /* Normalize V (and shift up U the same amount).  */
                    199:          count_leading_zeros (cnt, vp[vsize - 1]);
                    200:          if (cnt != 0)
                    201:            {
                    202:              mp_limb_t cy;
                    203:              mpn_lshift (vp, vp, vsize, cnt);
                    204:              cy = mpn_lshift (up, up, size, cnt);
                    205:              up[size] = cy;
                    206:              size += cy != 0;
                    207:            }
                    208:
                    209:          qh = mpn_divmod (up + vsize, up, size, vp, vsize);
                    210: #if EXTEND
                    211:          MPN_COPY (tp, s0p, ssize);
                    212:          for (i = 0; i < size - vsize; i++)
                    213:            {
                    214:              mp_limb_t cy;
                    215:              cy = mpn_addmul_1 (tp + i, s1p, ssize, up[vsize + i]);
                    216:              if (cy != 0)
                    217:                tp[ssize++] = cy;
                    218:            }
                    219:          if (qh != 0)
                    220:            {
                    221:              mp_limb_t cy;
                    222:              abort ();
                    223:              /* XXX since qh == 1, mpn_addmul_1 is overkill */
                    224:              cy = mpn_addmul_1 (tp + size - vsize, s1p, ssize, qh);
                    225:              if (cy != 0)
                    226:                tp[ssize++] = cy;
                    227:            }
                    228: #if 0
                    229:          MPN_COPY (s0p, s1p, ssize); /* should be old ssize, kind of */
                    230:          MPN_COPY (s1p, tp, ssize);
                    231: #else
                    232:          {
                    233:            mp_ptr xp;
                    234:            xp = s0p; s0p = s1p; s1p = xp;
                    235:            xp = s1p; s1p = tp; tp = xp;
                    236:          }
                    237: #endif
                    238: #endif
                    239:          size = vsize;
                    240:          if (cnt != 0)
                    241:            {
                    242:              mpn_rshift (up, up, size, cnt);
                    243:              mpn_rshift (vp, vp, size, cnt);
                    244:            }
                    245:
                    246:          {
                    247:            mp_ptr xp;
                    248:            xp = up; up = vp; vp = xp;
                    249:          }
                    250:          MPN_NORMALIZE (up, size);
                    251:        }
                    252:       else
                    253:        {
                    254:          /* T = U*A + V*B
                    255:             W = U*C + V*D
                    256:             U = T
                    257:             V = W         */
                    258:
                    259:          if (SGN(A) == SGN(B)) /* should be different sign */
                    260:            abort ();
                    261:          if (SGN(C) == SGN(D)) /* should be different sign */
                    262:            abort ();
                    263: #if STAT
                    264:          { mp_limb_t x;
                    265:            x = ABS (A) | ABS (B) | ABS (C) | ABS (D);
                    266:            count_leading_zeros (cnt, x);
                    267:            arr[BITS_PER_MP_LIMB - cnt]++; }
                    268: #endif
                    269:          if (A == 0)
                    270:            {
                    271:              if (B != 1) abort ();
                    272:              MPN_COPY (tp, vp, size);
                    273:            }
                    274:          else
                    275:            {
                    276:              if (A < 0)
                    277:                {
                    278:                  mpn_mul_1 (tp, vp, size, B);
                    279:                  mpn_submul_1 (tp, up, size, -A);
                    280:                }
                    281:              else
                    282:                {
                    283:                  mpn_mul_1 (tp, up, size, A);
                    284:                  mpn_submul_1 (tp, vp, size, -B);
                    285:                }
                    286:            }
                    287:          if (C < 0)
                    288:            {
                    289:              mpn_mul_1 (wp, vp, size, D);
                    290:              mpn_submul_1 (wp, up, size, -C);
                    291:            }
                    292:          else
                    293:            {
                    294:              mpn_mul_1 (wp, up, size, C);
                    295:              mpn_submul_1 (wp, vp, size, -D);
                    296:            }
                    297:
                    298:          {
                    299:            mp_ptr xp;
                    300:            xp = tp; tp = up; up = xp;
                    301:            xp = wp; wp = vp; vp = xp;
                    302:          }
                    303:
                    304: #if EXTEND
                    305:          { mp_limb_t cy;
                    306:          MPN_ZERO (tp, orig_size);
                    307:          if (A == 0)
                    308:            {
                    309:              if (B != 1) abort ();
                    310:              MPN_COPY (tp, s1p, ssize);
                    311:            }
                    312:          else
                    313:            {
                    314:              if (A < 0)
                    315:                {
                    316:                  cy = mpn_mul_1 (tp, s1p, ssize, B);
                    317:                  cy += mpn_addmul_1 (tp, s0p, ssize, -A);
                    318:                }
                    319:              else
                    320:                {
                    321:                  cy = mpn_mul_1 (tp, s0p, ssize, A);
                    322:                  cy += mpn_addmul_1 (tp, s1p, ssize, -B);
                    323:                }
                    324:              if (cy != 0)
                    325:                tp[ssize++] = cy;
                    326:            }
                    327:          MPN_ZERO (wp, orig_size);
                    328:          if (C < 0)
                    329:            {
                    330:              cy = mpn_mul_1 (wp, s1p, ssize, D);
                    331:              cy += mpn_addmul_1 (wp, s0p, ssize, -C);
                    332:            }
                    333:          else
                    334:            {
                    335:              cy = mpn_mul_1 (wp, s0p, ssize, C);
                    336:              cy += mpn_addmul_1 (wp, s1p, ssize, -D);
                    337:            }
                    338:          if (cy != 0)
                    339:            wp[ssize++] = cy;
                    340:          }
                    341:          {
                    342:            mp_ptr xp;
                    343:            xp = tp; tp = s0p; s0p = xp;
                    344:            xp = wp; wp = s1p; s1p = xp;
                    345:          }
                    346: #endif
                    347: #if 0  /* Is it a win to remove multiple zeros here? */
                    348:          MPN_NORMALIZE (up, size);
                    349: #else
                    350:          if (up[size - 1] == 0)
                    351:            size--;
                    352: #endif
                    353:        }
                    354:     }
                    355:
                    356: #if RECORD
                    357:   printf ("min: %ld\n", min);
                    358:   printf ("max: %ld\n", max);
                    359: #endif
                    360:
                    361:   if (vsize == 0)
                    362:     {
                    363:       if (gp != up)
                    364:        MPN_COPY (gp, up, size);
                    365: #if EXTEND
                    366:       if (orig_s0p != s0p)
                    367:        MPN_COPY (orig_s0p, s0p, ssize);
                    368: #endif
                    369:       TMP_FREE (mark);
                    370:       return size;
                    371:     }
                    372:   else
                    373:     {
                    374:       mp_limb_t vl, ul, t;
                    375: #if EXTEND
                    376:       mp_limb_t cy;
                    377:       mp_size_t i;
                    378: #endif
                    379:       vl = vp[0];
                    380: #if EXTEND
                    381:       t = mpn_divmod_1 (wp, up, size, vl);
                    382:       MPN_COPY (tp, s0p, ssize);
                    383:       for (i = 0; i < size; i++)
                    384:        {
                    385:          cy = mpn_addmul_1 (tp + i, s1p, ssize, wp[i]);
                    386:          if (cy != 0)
                    387:            tp[ssize++] = cy;
                    388:        }
                    389: #if 0
                    390:       MPN_COPY (s0p, s1p, ssize);
                    391:       MPN_COPY (s1p, tp, ssize);
                    392: #else
                    393:       {
                    394:        mp_ptr xp;
                    395:        xp = s0p; s0p = s1p; s1p = xp;
                    396:        xp = s1p; s1p = tp; tp = xp;
                    397:       }
                    398: #endif
                    399: #else
                    400:       t = mpn_mod_1 (up, size, vl);
                    401: #endif
                    402:       ul = vl;
                    403:       vl = t;
                    404:       while (vl != 0)
                    405:        {
                    406:          mp_limb_t t;
                    407: #if EXTEND
                    408:          mp_limb_t q, cy;
                    409:          q = ul / vl;
                    410:          t = ul - q*vl;
                    411:
                    412:          MPN_COPY (tp, s0p, ssize);
                    413:          cy = mpn_addmul_1 (tp, s1p, ssize, q);
                    414:          if (cy != 0)
                    415:            tp[ssize++] = cy;
                    416: #if 0
                    417:          MPN_COPY (s0p, s1p, ssize);
                    418:          MPN_COPY (s1p, tp, ssize);
                    419: #else
                    420:          {
                    421:            mp_ptr xp;
                    422:            xp = s0p; s0p = s1p; s1p = xp;
                    423:            xp = s1p; s1p = tp; tp = xp;
                    424:          }
                    425: #endif
                    426:
                    427: #else
                    428:          t = ul % vl;
                    429: #endif
                    430:          ul = vl;
                    431:          vl = t;
                    432:        }
                    433:       gp[0] = ul;
                    434: #if EXTEND
                    435:       if (orig_s0p != s0p)
                    436:        MPN_COPY (orig_s0p, s0p, ssize);
                    437: #endif
                    438:       TMP_FREE (mark);
                    439:       return 1;
                    440:     }
                    441: }

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>