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

##### Functions

    dig_gcd           X
    even              X
    int_gcd           X
    int_half          X
    twice             X

       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

dig_gcd defined in line 26; used 1 times

- in line 131

even defined in line 111; used 4 times

- in line 144(2), 159, 169

int_gcd defined in line 120; used 2 times

- in /usr/src/new/PORT/B/src/bint/b1nuR.c line 64
- in /usr/src/new/PORT/B/src/bint/b1num.h line 48

int_half defined in line 55; used 5 times

- in line 145-146(2), 167
- in /usr/src/new/PORT/B/src/bint/b1nuR.c line 229
- in /usr/src/new/PORT/B/src/bint/b1num.h line 52

twice defined in line 86; used 1 times

- in line 188

Last modified: 1985-08-27 Generated: 2016-12-26 | Generated by src2html V0.67 page hit count: 693 |