![[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.