1 /* Copyright (c) 1984, 1986, 1987, 1988, 1989 AT&T */ 2 /* All Rights Reserved */ 3 4 5 /* 6 * Copyright (c) 1980 Regents of the University of California. 7 * All rights reserved. The Berkeley software License Agreement 8 * specifies the terms and conditions for redistribution. 9 */ 10 /* Portions Copyright(c) 1988, Sun Microsystems Inc. */ 11 /* All Rights Reserved */ 12 13 /* 14 * Copyright (c) 1997, by Sun Microsystems, Inc. 15 * All rights reserved. 16 */ 17 18 #ident "%Z%%M% %I% %E% SMI" /* SVr4.0 1.1 */ 19 20 /* LINTLIBRARY */ 21 22 #include <mp.h> 23 #include <stdio.h> 24 #include <stdlib.h> 25 #include <sys/types.h> 26 #include "libmp.h" 27 28 static void m_div(MINT *, MINT *, MINT *, MINT *); 29 30 void 31 mp_mdiv(MINT *a, MINT *b, MINT *q, MINT *r) 32 { 33 MINT x, y; 34 int sign; 35 36 sign = 1; 37 x.len = y.len = 0; 38 _mp_move(a, &x); 39 _mp_move(b, &y); 40 if (x.len < 0) { 41 sign = -1; 42 x.len = -x.len; 43 } 44 if (y.len < 0) { 45 sign = -sign; 46 y.len = -y.len; 47 } 48 _mp_xfree(q); 49 _mp_xfree(r); 50 m_div(&x, &y, q, r); 51 if (sign == -1) { 52 q->len = -q->len; 53 r->len = -r->len; 54 } 55 _mp_xfree(&x); 56 _mp_xfree(&y); 57 } 58 59 static int 60 m_dsb(int qx, int n, short *a, short *b) 61 { 62 int borrow; 63 int s3b2shit; 64 int j; 65 short fifteen = 15; 66 short *aptr, *bptr; 67 #ifdef DEBUGDSB 68 (void) printf("m_dsb %d %d %d %d\n", qx, n, *a, *b); 69 #endif 70 71 borrow = 0; 72 aptr = a; 73 bptr = b; 74 for (j = n; j > 0; j--) { 75 #ifdef DEBUGDSB 76 (void) printf("1 borrow=%x %d %d %d\n", borrow, (*aptr * qx), 77 *bptr, *aptr); 78 #endif 79 borrow -= (*aptr++) * qx - *bptr; 80 #ifdef DEBUGDSB 81 (void) printf("2 borrow=%x %d %d %d\n", borrow, (*aptr * qx), 82 *bptr, *aptr); 83 #endif 84 *bptr++ = (short)(borrow & 077777); 85 #ifdef DEBUGDSB 86 (void) printf("3 borrow=%x %d %d %d\n", borrow, (*aptr * qx), 87 *bptr, *aptr); 88 #endif 89 if (borrow >= 0) borrow >>= fifteen; /* 3b2 */ 90 else borrow = 0xfffe0000 | (borrow >> fifteen); 91 #ifdef DEBUGDSB 92 (void) printf("4 borrow=%x %d %d %d\n", borrow, (*aptr * qx), 93 *bptr, *aptr); 94 #endif 95 } 96 borrow += *bptr; 97 *bptr = (short)(borrow & 077777); 98 if (borrow >= 0) s3b2shit = borrow >> fifteen; /* 3b2 */ 99 else s3b2shit = 0xfffe0000 | (borrow >> fifteen); 100 if (s3b2shit == 0) { 101 #ifdef DEBUGDSB 102 (void) printf("mdsb 0\n"); 103 #endif 104 return (0); 105 } 106 borrow = 0; 107 aptr = a; 108 bptr = b; 109 for (j = n; j > 0; j--) { 110 borrow += *aptr++ + *bptr; 111 *bptr++ = (short)(borrow & 077777); 112 if (borrow >= 0) borrow >>= fifteen; /* 3b2 */ 113 else borrow = 0xfffe0000 | (borrow >>fifteen); 114 } 115 #ifdef DEBUGDSB 116 (void) printf("mdsb 1\n"); 117 #endif 118 return (1); 119 } 120 121 static int 122 m_trq(short v1, short v2, short u1, short u2, short u3) 123 { 124 short d; 125 int x1; 126 int c1; 127 128 c1 = u1 * 0100000 + u2; 129 if (u1 == v1) { 130 d = 077777; 131 } else { 132 d = (short)(c1 / v1); 133 } 134 do { 135 x1 = c1 - v1 * d; 136 x1 = x1 * 0100000 + u3 - v2 * d; 137 --d; 138 } while (x1 < 0); 139 #ifdef DEBUGMTRQ 140 (void) printf("mtrq %d %d %d %d %d %d\n", v1, v2, u1, u2, u3, (d+1)); 141 #endif 142 return ((int)d + 1); 143 } 144 145 static void 146 m_div(MINT *a, MINT *b, MINT *q, MINT *r) 147 { 148 MINT u, v, x, w; 149 short d; 150 short *qval; 151 short *uval; 152 int j; 153 int qq; 154 int n; 155 short v1; 156 short v2; 157 158 u.len = v.len = x.len = w.len = 0; 159 if (b->len == 0) { 160 _mp_fatal("mdiv divide by zero"); 161 return; 162 } 163 if (b->len == 1) { 164 r->val = _mp_xalloc(1, "m_div1"); 165 mp_sdiv(a, b->val[0], q, r->val); 166 if (r->val[0] == 0) { 167 free(r->val); 168 r->len = 0; 169 } else { 170 r->len = 1; 171 } 172 return; 173 } 174 if (a -> len < b -> len) { 175 q->len = 0; 176 r->len = a->len; 177 r->val = _mp_xalloc(r->len, "m_div2"); 178 for (qq = 0; qq < r->len; qq++) { 179 r->val[qq] = a->val[qq]; 180 } 181 return; 182 } 183 x.len = 1; 184 x.val = &d; 185 n = b->len; 186 d = 0100000 / (b->val[n - 1] + 1); 187 mp_mult(a, &x, &u); /* subtle: relies on mult allocing extra space */ 188 mp_mult(b, &x, &v); 189 #ifdef DEBUG_MDIV 190 (void) printf(" u=%s\n", mtox(&u)); 191 (void) printf(" v=%s\n", mtox(&v)); 192 #endif 193 v1 = v.val[n - 1]; 194 v2 = v.val[n - 2]; 195 qval = _mp_xalloc(a -> len - n + 1, "m_div3"); 196 uval = u.val; 197 for (j = a->len - n; j >= 0; j--) { 198 qq = m_trq(v1, v2, uval[j + n], uval[j + n - 1], 199 uval[j + n - 2]); 200 if (m_dsb(qq, n, v.val, uval + j)) 201 qq -= 1; 202 qval[j] = (short)qq; 203 } 204 x.len = n; 205 x.val = u.val; 206 _mp_mcan(&x); 207 #ifdef DEBUG_MDIV 208 (void) printf(" x=%s\n", mtox(&x)); 209 (void) printf(" d(in)=%d\n", (d)); 210 #endif 211 mp_sdiv(&x, d, &w, &d); 212 #ifdef DEBUG_MDIV 213 (void) printf(" w=%s\n", mtox(&w)); 214 (void) printf(" d(out)=%d\n", (d)); 215 #endif 216 r->len = w.len; 217 r->val = w.val; 218 q->val = qval; 219 qq = a->len - n + 1; 220 if (qq > 0 && qval[qq - 1] == 0) 221 qq -= 1; 222 q->len = qq; 223 if (qq == 0) 224 free(qval); 225 if (x.len != 0) 226 _mp_xfree(&u); 227 _mp_xfree(&v); 228 } 229