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