![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to b1nuG.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [CSRG BSD Unix] / 43BSD / contrib / B / src / bint|
Return to b1nuG.c CVS log | Up to [CSRG BSD Unix] / 43BSD / contrib / B / src / bint |
1.1 root 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: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.