2.11BSD#431 - /usr/src/new/PORT/B/src/bint/b1nuG.c

   1: /* Copyright (c) Stichting Mathematisch Centrum, Amsterdam, 1985. */
   2: 
   3: /*
   4:   $Header: b1nuG.c,v 1.4 85/08/22 16:50:59 timo Exp $
   5: */
   6: 
   7: #include "b.h"
   8: #include "b0con.h"
   9: #include "b1obj.h"
  10: #include "b1num.h"
  11: 
  12: 
  13: /*
  14:  * Routines for greatest common divisor calculation
  15:  * Cf. Knuth II Sec. 4.5.2. Algorithm B
  16:  * "Binary gcd algorithm"
  17:  * The labels correspond with those in the book
  18:  *
  19:  * Assumptions about built-in arithmetic:
  20:  * x>>1 == x/2  (if x >= 0)
  21:  * 1<<k == 2**k (if it fits in a word)
  22:  */
  23: 
  24: /* Single-precision gcd for integers > 0 */
  25: 
  26: Hidden digit dig_gcd(u, v) register digit u, v; {
  27:     register digit t;
  28:     register int k = 0;
  29: 
  30:     if (u <= 0 || v <= 0) syserr(MESS(900, "dig_gcd of number(s) <= 0"));
  31: 
  32:  /*B1*/ while (Even(u) && Even(v)) ++k, u >>= 1, v >>= 1;
  33: 
  34:  /*B2*/ if (Even(u)) t = u;
  35:     else {  t = -v;
  36:         goto B4;
  37:     }
  38: 
  39:     do {
  40:  /*B3*/     do {
  41:             t /= 2;
  42:  B4:            ;
  43:         } while (Even(t));
  44: 
  45:  /*B5*/     if (t > 0) u = t;
  46:         else v = -t;
  47: 
  48:  /*B6*/     t = u-v;
  49:     } while (t);
  50: 
  51:     return u * (1<<k);
  52: }
  53: 
  54: 
  55: Visible integer int_half(v) integer v; {
  56:     register int i;
  57:     register long carry;
  58: 
  59:     if (IsSmallInt(v))  return (integer) MkSmallInt(SmallIntVal(v) / 2);
  60: 
  61:     if (Msd(v) < 0) {
  62:         i = Length(v)-2;
  63:         if (i < 0) {
  64:             release((value) v);
  65:             return int_0;
  66:         }
  67:         carry = BASE;
  68:     }
  69:     else {
  70:         carry = 0;
  71:         i = Length(v)-1;
  72:     }
  73: 
  74:     if (Refcnt(v) > 1) uniql((value *) &v);
  75: 
  76:     for (; i >= 0; --i) {
  77:         carry += Digit(v,i);
  78:         Digit(v,i) = carry/2;
  79:         carry = carry&1 ? BASE : 0;
  80:     }
  81: 
  82:     return int_canon(v);
  83: }
  84: 
  85: 
  86: Hidden integer twice(v) integer v; {
  87:     digit carry = 0;
  88:     int i;
  89: 
  90:     if (IsSmallInt(v)) return mk_int(2.0 * SmallIntVal(v));
  91: 
  92:     if (Refcnt(v) > 1) uniql((value *) &v);
  93: 
  94:     for (i = 0; i < Length(v); ++i) {
  95:         carry += Digit(v,i) * 2;
  96:         if (carry >= BASE)
  97:             Digit(v,i) = carry-BASE, carry = 1;
  98:         else
  99:             Digit(v,i) = carry, carry = 0;
 100:     }
 101: 
 102:     if (carry) {    /* Enlarge the number */
 103:         v = (integer) regrab_num((value) v, Length(v)+1);
 104:         Digit(v, Length(v)-1) = carry;
 105:     }
 106: 
 107:     return v;
 108: }
 109: 
 110: 
 111: Hidden bool even(u) integer u; {
 112:     if (IsSmallInt(u))
 113:         return Even(SmallIntVal(u));
 114:     return Even(Lsd(u));
 115: }
 116: 
 117: 
 118: /* Multi-precision gcd of integers > 0 */
 119: 
 120: Visible integer int_gcd(u1, v1) integer u1, v1; {
 121:     struct integer uu1, vv1;
 122: 
 123:     if (Msd(u1) <= 0 || Msd(v1) <= 0)
 124:         syserr(MESS(901, "gcd of number(s) <= 0"));
 125: 
 126:     if (u1==int_1 || v1==int_1) return int_1;
 127:         /* Speed-up for e.g. 1E-100000 */
 128: 
 129:     if (IsSmallInt(u1) && IsSmallInt(v1))
 130:         return (integer)
 131:             MkSmallInt(dig_gcd(SmallIntVal(u1), SmallIntVal(v1)));
 132: 
 133:     FreezeSmallInt(u1, uu1);
 134:     FreezeSmallInt(v1, vv1);
 135: 
 136:     /* Multi-precision binary gcd algorithm */
 137: 
 138:     {   long k = 0;
 139:         integer t, u, v;
 140: 
 141:         u = (integer) Copy((value) u1);
 142:         v = (integer) Copy((value) v1);
 143: 
 144:  /*B1*/     while (even(u) && even(v)) {
 145:             u = int_half(u);
 146:             v = int_half(v);
 147:             if (++k < 0) {
 148:                 /*It's a number we can't cope with,
 149: 				  with too many common factors 2.
 150: 				  Though the user can't help it,
 151: 				  the least we can do is to allow
 152: 				  continuation of the session.
 153: 				*/
 154:                 error(MESS(902, "exceptionally large rational number"));
 155:                 k = 0;
 156:             }
 157:         }
 158: 
 159:  /*B2*/     if (even(u)) t = (integer) Copy(u);
 160:         else {
 161:             t = int_neg(v);
 162:             goto B4;
 163:         }
 164: 
 165:         do {
 166:  /*B3*/         do {
 167:                 t = int_half(t);
 168:  B4:                ;
 169:             } while (even(t));
 170: 
 171:  /*B5*/         if (Msd(t) >= 0) {
 172:                 release((value) u);
 173:                 u = t;
 174:             }
 175:             else {
 176:                 release((value) v);
 177:                 v = int_neg(t);
 178:                 release((value) t);
 179:             }
 180: 
 181:  /*B6*/         t = int_diff(u, v);
 182:             /* t cannot be int_1 since both u and v are odd! */
 183:         } while (t != int_0);
 184: 
 185:         release((value) t);
 186:         release((value) v);
 187: 
 188:         while (--k >= 0) u = twice(u);
 189: 
 190:         return u;
 191:     }
 192: }

Defined functions