![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to b1nuQ.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 b1nuQ.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: b1nuQ.c,v 1.4 85/08/22 16:51:40 timo Exp $
5: */
6:
7: #include "b.h"
8: #include "b1obj.h"
9: #include "b0con.h"
10: #include "b1num.h"
11:
12:
13: /* Product of integer and single "digit" */
14:
15: Visible integer int1mul(v, n1) integer v; digit n1; {
16: integer a;
17: digit save, bigcarry, carry = 0;
18: double z, zz, n = n1;
19: register int i;
20: struct integer vv;
21:
22: FreezeSmallInt(v, vv);
23:
24: a = (integer) grab_num(Length(v)+2);
25:
26: for (i = 0; i < Length(v); ++i) {
27: z = Digit(v,i) * n;
28: bigcarry = zz = floor(z/BASE);
29: carry += z - zz*BASE;
30: Digit(a,i) = save = Modulo(carry, BASE);
31: carry = (carry-save)/BASE + bigcarry;
32: }
33:
34: Digit(a,i) = save = Modulo(carry, BASE);
35: Digit(a,i+1) = (carry-save)/BASE;
36:
37: return int_canon(a);
38: }
39:
40:
41: /* Quotient of positive integer and single "digit" > 0 */
42:
43: Hidden integer int1div(v, n1, prem) integer v; digit n1, *prem; {
44: integer q;
45: double r_over_n, r = 0, n = n1;
46: register int i;
47: struct integer vv;
48:
49: FreezeSmallInt(v, vv);
50:
51: q = (integer) grab_num(Length(v));
52: for (i = Length(v)-1; i >= 0; --i) {
53: r = r*BASE + Digit(v,i);
54: Digit(q,i) = r_over_n = floor(r/n);
55: r -= r_over_n * n;
56: }
57: if (prem)
58: *prem = r;
59: return int_canon(q);
60: }
61:
62:
63: /* Long division routine, gives access to division algorithm. */
64:
65: Visible digit int_ldiv(v1, w1, pquot, prem) integer v1, w1, *pquot, *prem; {
66: integer a;
67: int sign = 1, rel_v = 0, rel_w = 0;
68: digit div, rem;
69: struct integer vv1, ww1;
70:
71: if (w1 == int_0) syserr(MESS(1100, "zero division (int_ldiv)"));
72:
73: /* Make v, w positive */
74: if (Msd(v1) < 0) {
75: sign = -1;
76: ++rel_v;
77: v1 = int_neg(v1);
78: }
79:
80: if (Msd(w1) < 0) {
81: sign *= -1;
82: ++rel_w;
83: w1 = int_neg(w1);
84: }
85:
86: FreezeSmallInt(v1, vv1);
87: FreezeSmallInt(w1, ww1);
88:
89: div = sign;
90:
91: /* Check v << w or single-digit w */
92: if (Length(v1) < Length(w1)
93: || Length(v1) == Length(w1)
94: && Digit(v1, Length(v1)-1) < Digit(w1, Length(w1)-1)) {
95: a = int_0;
96: if (prem) {
97: if (v1 == &vv1) *prem= (integer) MkSmallInt(Digit(v1,0));
98: else *prem = (integer) Copy(v1);
99: }
100: }
101: else if (Length(w1) == 1) {
102: /* Single-precision division */
103: a = int1div(v1, Digit(w1,0), &rem);
104: if (prem) *prem = mk_int((double)rem);
105: }
106: else {
107: /* Multi-precision division */
108: /* Cf. Knuth II Sec. 4.3.1. Algorithm D */
109: /* Note that we count in the reverse direction (not easier!) */
110:
111: double z, zz;
112: digit carry, save, bigcarry;
113: double q, d = BASE/(Digit(w1, Length(w1)-1)+1);
114: register int i, j, k;
115: integer v, w;
116: digit vj;
117:
118: /* Normalize: make Msd(w) >= BASE/2 by multiplying
119: both v and w by d */
120:
121: v = int1mul(v1, (digit)d);
122: /* v is used as accumulator, must make a copy */
123: /* v cannot be int_1 */
124: /* (then it would be one of the cases above) */
125:
126: if (d == 1) w = (integer) Copy(w1);
127: else w = int1mul(w1, (digit)d);
128:
129: a = (integer) grab_num(Length(v1)-Length(w)+1);
130:
131: /* Division loop */
132:
133: for (j = Length(v1), k = Length(a)-1; k >= 0; --j, --k) {
134: vj = j >= Length(v) ? 0 : Digit(v,j);
135:
136: /* Find trial digit */
137:
138: if (vj == Digit(w, Length(w)-1)) q = BASE-1;
139: else q = floor( ((double)vj*BASE
140: + Digit(v,j-1)) / Digit(w, Length(w)-1) );
141:
142: /* Correct trial digit */
143:
144: while (Digit(w,Length(w)-2) * q >
145: ((double)vj*BASE + Digit(v,j-1)
146: - q*Digit(w, Length(w)-1)) *BASE + Digit(v,j-2))
147: --q;
148:
149: /* Subtract q*w from v */
150:
151: carry = 0;
152: for (i = 0; i < Length(w) && i+k < Length(v); ++i) {
153: z = Digit(w,i) * q;
154: bigcarry = zz = floor(z/BASE);
155: carry += Digit(v,i+k) - z + zz*BASE;
156: Digit(v,i+k) =
157: save = Modulo(carry, BASE);
158: carry = (carry-save)/BASE - bigcarry;
159: }
160:
161: if (i+k < Length(v))
162: carry += Digit(v, i+k), Digit(v, i+k) = 0;
163:
164: /* Add back necessary? */
165:
166: /* It is very difficult to find test cases
167: where add back is necessary if BASE is large.
168: Thanks to Arjen Lenstra, we have v=n*n-1, w=n,
169: where n = 8109636009903000000 (the last six
170: digits are not important). */
171:
172: if (carry == 0) /* No */
173: Digit(a,k) = q;
174: else { /* Yes, add back */
175: if (carry != -1) syserr(MESS(1101, "int_ldiv internal failure"));
176: Digit(a,k) = q-1;
177: carry = 0;
178: for (i = 0; i < Length(w) && i+k < Length(v); ++i) {
179: carry += Digit(v, i+k) + Digit(w,i);
180: Digit(v,i+k) =
181: save = Modulo(carry, BASE);
182: carry = (carry-save)/BASE;
183: }
184: }
185: } /* End for(j) */
186:
187: if (prem) *prem = int_canon(v); /* Store remainder */
188: else release((value) v);
189: div = sign*d; /* Store normalization factor */
190: release((value) w);
191: a = int_canon(a);
192: }
193:
194: if (rel_v) release((value) v1);
195: if (rel_w) release((value) w1);
196:
197: if (sign < 0) {
198: integer temp = a;
199: a = int_neg(a);
200: release((value) temp);
201: }
202:
203: if (pquot) *pquot = a;
204: else release((value) a);
205: return div;
206: }
207:
208:
209: Visible integer int_quot(v, w) integer v, w; {
210: integer quo;
211: VOID int_ldiv(v, w, &quo, (integer*)0);
212: return quo;
213: }
214:
215: Visible integer int_mod(v, w) integer v, w; {
216: integer rem;
217: digit div;
218: bool flag;
219: div = int_ldiv(v, w, (integer*)0, &rem); /* Rem. is always positive */
220: if (rem == int_0)
221: return rem; /* v mod w = 0 */
222: flag = (div < 0);
223: if (flag || Msd(w) < 0) div = -div;
224: if (div != 1) { /* Divide by div to get proper remainder back */
225: v = int1div(rem, div, (digit*)0);
226: release((value) rem);
227: rem = v;
228: }
229: if (flag) { /* Make same sign as w */
230: if (Msd(w) < 0) v = int_sum(w, rem);
231: else v = int_diff(w, rem);
232: release((value) rem);
233: rem = v;
234: }
235: return rem;
236: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.