```   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
even defined in line 111; used 4 times
int_half defined in line 55; used 5 times
twice defined in line 86; used 1 times
 Last modified: 1985-08-27 Generated: 2016-12-26 Generated by src2html V0.67 page hit count: 2128