![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to mdiv.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [MW Coherent from dump] / coherent / d / usr / lib / libmp|
Return to mdiv.c CVS log | Up to [MW Coherent from dump] / coherent / d / usr / lib / libmp |
1.1 root 1: #include "mprec.h"
2: #include <assert.h>
3: #include <sys/mdata.h>
4:
5:
6: /*
7: * Mdiv sets the mints pointed to by "q" and "r" to the quotient
8: * and remainder (respectively) of dividing "a" by "b". If "b"
9: * is zero, then it calls mperr with the appropriate message.
10: * The division is performed such that the following two properties
11: * hold :
12: * 1. "r" + "q" * "b" = "a".
13: * 2. The sign of "r" = the sign of "q".
14: * 3. The abs. value of "r" < abs. value of "b".
15: * Note that no assumption is made as to the distinctness of "a",
16: * "b", "q" and "r".
17: */
18:
19: void
20: mdiv(a, b, q, r)
21: register mint *a, *b, *q, *r;
22: {
23: int apos, bpos;
24: mint al, bl;
25: int ispos();
26: void shortd(), longd();
27:
28: /* take care of small denominators */
29: if (b->len <= 2 && (
30: b->len == 1 /* 0 to BASE - 2 */
31: || b->val[1] == NEFL /* -BASE to -2 */
32: || b->val[1] == 0 /* BASE - 1 */
33: || (b->val[1] == 1 && b->val[0] == 0) /* BASE */
34: )
35: ) {
36: shortd(a, b, q, r);
37: return;
38: }
39:
40: /* make denominator negative */
41: if (bpos = ispos(b)) {
42: minit(&bl);
43: mneg(b, &bl);
44: b = &bl;
45: }
46: assert(b->len > 2);
47:
48: /* copy numerator, makeing it positive */
49: minit(&al);
50: if (apos = ispos(a))
51: mcopy(a, &al);
52: else
53: mneg(a, &al);
54: a = &al;
55:
56: /*
57: * BASE ^ (a->len - 1) <= a < BASE ^ a->len
58: * BASE ^ (b->len - 2) < -b <= BASE ^ (b->len - 1)
59: *
60: * do divide.
61: */
62:
63: if (a->len <= b->len - 2)
64: mcopy(mzero, q);
65: else
66: longd(a, b, q);
67:
68: /* adjust signs */
69: if (apos ? !bpos : bpos) {
70: mneg(q, q);
71: mneg(a, r);
72: } else
73: mcopy(a, r);
74:
75: /* throw away garbage */
76: mpfree(al.val);
77: if (bpos)
78: mpfree(bl.val);
79: }
80:
81:
82: /*
83: * Shortd divides the mints pointed to by "a" and "b".
84: * It sets the mints pointed to by "q" and "r" to the quotient
85: * and remainder respectively. It assumes that "b" is between
86: * -BASE and BASE.
87: */
88:
89: static void
90: shortd(a, b, q, r)
91: register mint *a, *b, *q, *r;
92: {
93: int rem;
94: int denom;
95:
96: assert(b->len <= 2);
97: if (ispos(b)) {
98: denom = b->val[0];
99: if (b->len != 1 && denom == 0)
100: denom = BASE;
101: } else
102: denom = b->val[0] - BASE;
103: assert(-BASE <= denom && denom <= BASE);
104: sdiv(a, denom, q, &rem);
105: mitom(rem, r);
106: }
107:
108:
109: /*
110: * Longd does the actual long divide. It sets "q" and "a" to
111: * the quotient and remainder (respectively) of "a" / (- "b").
112: * It assumes the following:
113: * 1. "b" is negative
114: * 2. "a" is positive
115: * 3. "a" is a destroyable copy
116: * 4. ("a"->len + 2) > ("b"->len) >= 3
117: */
118:
119: static void
120: longd(a, b, q)
121: register mint *a, *b;
122: mint *q;
123: {
124: register char *rp;
125: unsigned shift;
126: int tden;
127: unsigned tbias;
128: mint res;
129:
130:
131: /*
132: * For "a" we have the range
133: * (BASE - 1) * BASE ^ (a->len - 1) > a
134: * a >= (BASE - 1) * BASE ^ (a->len - 2)
135: * and for "b" we want
136: * tden * BASE ^ tbias >= -b > (tden - 1) * BASE ^ tbias
137: */
138:
139: tden = est(b);
140: tbias = b->len - 3;
141:
142: assert(a->len >= tbias + 2);
143: res.len = a->len - tbias;
144: res.val = (char *)mpalc(res.len);
145: rp = &res.val[res.len - 1];
146: *rp-- = 0;
147: shift = a->len - (tbias + 2);
148:
149: do {
150: *rp = guess(tden, tbias + shift, a);
151: msma(b, *rp, shift, a);
152: if (snc(b, shift, a)) {
153: ++*rp;
154: msma(b, 1, shift, a);
155: }
156: --rp;
157: } while (shift-- != 0);
158: norm(a);
159: norm(&res);
160: mpfree(q->val);
161: *q = res;
162: }
163:
164:
165:
166:
167: /*
168: * Est estimates the mint pointed to by "a". Its result
169: * satisfies the following:
170: * 1. result * BASE ^ (a->len - 3) >= -a
171: * && -a > (result - 1) * BASE ^ (a->len - 3)
172: * 2. BASE ^ 2 >= mant > BASE
173: * Note that "a" is assumed to be negative and have length atleast
174: * three.
175: */
176: static int
177: est(a)
178: register mint *a;
179: {
180: register int result;
181:
182: assert(!ispos(a) && a->len >= 3);
183: result = (NEFL - a->val[a->len - 2]) << L2BASE;
184: result += NEFL - a->val[a->len - 3] + 1;
185: assert(BASE * BASE >= result && result > BASE);
186: return (result);
187: }
188:
189:
190: /*
191: * Guess returns an integer which is the integral part of
192: * "a" / ("den" * BASE ^ "shift")
193: * Note that it is assumed that the following hold:
194: * 1. "a" is positive.
195: * 2. BASE ^ 2 >= "den" > BASE
196: * 3. the quotient is a single "digit" (ie. NEFL > quotient >= 0)
197: */
198:
199: static int
200: guess(den, shift, a)
201: int den;
202: unsigned shift;
203: mint *a;
204: {
205: register char *ap = &a->val[a->len - 1];
206: long result;
207:
208: if (*ap == 0)
209: --ap;
210: assert(ap + 1 - a->val <= shift + 3);
211: if (ap + 1 - a->val < shift + 2)
212: return (0);
213: result = (ap[0] << L2BASE) + ap[-1];
214: if (ap + 1 - a->val == shift + 3)
215: result = (result << L2BASE) + ap[-2];
216: assert(result < den * (long)BASE);
217: return (result / den);
218: }
219:
220:
221: /*
222: * Snc does a shift-negate-compare. Specifically, it assumes that
223: * "b" is negative, "a" is positive and it returns wether or not
224: * a >= -b * BASE ^ shift
225: * (ie. a + b * BASE ^ shift is positive (ie. >= 0).)
226: */
227:
228: static
229: snc(b, shift, a)
230: mint *a, *b;
231: unsigned shift;
232: {
233: register char *ap;
234: register char *bp;
235: register char *limit;
236: int aleng, bleng;
237: int temp;
238:
239: /*
240: * correct length of a
241: */
242: limit = a->val;
243: for (ap = &limit[a->len - 1]; *ap == 0 && ap > limit; --ap)
244: ;
245: a->len = ap - limit + 1;
246: if (*ap == NEFL)
247: ++a->len;
248: aleng = ap - limit + 1;
249: bleng = b->len + shift - 1;
250:
251: /*
252: * BASE ^ aleng > a >= BASE ^ (aleng - 1)
253: * BASE ^ (b->len - 1) >= -b > BASE ^ (b->len - 2)
254: * so multiplying by BASE ^ shift gives
255: * BASE ^ bleng >= -b * BASE ^ shift > BASE ^ (bleng - 1)
256: * This means that
257: * IF bleng > aleng
258: * THEN bleng - 1 >= aleng
259: * SO -b * BASE ^ shift > BASE ^ (bleng - 1)
260: * >= BASE ^ aleng
261: * > a
262: * and thus we return FALSE.
263: * Also we have
264: * IF aleng > bleng
265: * THEN aleng - 1 >= bleng
266: * SO a >= BASE ^ (aleng - 1)
267: * >= BASE ^ bleng
268: * >= -b * BASE ^ shift
269: * and thus we return TRUE.
270: * Finally
271: * If aleng = bleng then the two numbers line up and
272: * it is just a matter of returning wether or not their
273: * sum would cause a carry.
274: */
275: if (aleng != bleng)
276: return (aleng > bleng);
277: limit = b->val;
278: for (bp = &limit[b->len - 2]; bp >= limit; --bp, --ap) {
279: temp = *ap + *bp;
280: if (temp >= BASE)
281: return (TRUE);
282: if (temp < NEFL)
283: return (FALSE);
284: }
285: return (FALSE);
286: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.