mmul.cpp

   1 // Bignum multiplication and division
   2
   3 #include "stdafx.h"
   4 #include "defs.h"
   5
   6 extern int ge(unsigned int *, unsigned int *, int);
   7
   8 static void mulf(unsigned int *, unsigned int *, int, unsigned int);
   9 static void addf(unsigned int *, unsigned int *, int);
  10 static void subf(unsigned int *, unsigned int *, int);
  11
  12 unsigned int *
  13 mmul(unsigned int *a, unsigned int *b)
  14 {
  15         int alen, blen, i, n;
  16         unsigned int *t, *x;
  17
  18         if (MZERO(a) || MZERO(b))
  19                 return mint(0);
  20
  21         if (MLENGTH(a) == 1 && a[0] == 1) {
  22                 t = mcopy(b);
  23                 MSIGN(t) *= MSIGN(a);
  24                 return t;
  25         }
  26
  27         if (MLENGTH(b) == 1 && b[0] == 1) {
  28                 t = mcopy(a);
  29                 MSIGN(t) *= MSIGN(b);
  30                 return t;
  31         }
  32
  33         alen = MLENGTH(a);
  34         blen = MLENGTH(b);
  35
  36         n = alen + blen;
  37
  38         x = mnew(n);
  39
  40         t = mnew(alen + 1);
  41
  42         for (i = 0; i < n; i++)
  43                 x[i] = 0;
  44
  45         /* sum of partial products */
  46
  47         for (i = 0; i < blen; i++) {
  48                 mulf(t, a, alen, b[i]);
  49                 addf(x + i, t, alen + 1);
  50         }
  51
  52         mfree(t);
  53
  54         /* length of product */
  55
  56         for (i = n - 1; i > 0; i--)
  57                 if (x[i])
  58                         break;
  59
  60         MLENGTH(x) = i + 1;
  61
  62         MSIGN(x) = MSIGN(a) * MSIGN(b);
  63
  64         return x;
  65 }
  66
  67 unsigned int *
  68 mdiv(unsigned int *a, unsigned int *b)
  69 {
  70         int alen, blen, i, n;
  71         unsigned int c, *t, *x, *y;
  72         unsigned long long jj, kk;
  73
  74         if (MZERO(b))
  75                 stop("divide by zero");
  76
  77         if (MZERO(a))
  78                 return mint(0);
  79
  80         alen = MLENGTH(a);
  81         blen = MLENGTH(b);
  82
  83         n = alen - blen;
  84
  85         if (n < 0)
  86                 return mint(0);
  87
  88         x = mnew(alen + 1);
  89
  90         for (i = 0; i < alen; i++)
  91                 x[i] = a[i];
  92
  93         x[i] = 0;
  94
  95         y = mnew(n + 1);
  96
  97         t = mnew(blen + 1);
  98
  99         /* Add 1 here to round up in case the remaining words are non-zero. */
 100
 101         kk = (unsigned long long) b[blen - 1] + 1;
 102
 103         for (i = 0; i <= n; i++) {
 104
 105                 y[n - i] = 0;
 106
 107                 for (;;) {
 108
 109                         /* estimate the partial quotient */
 110
 111                         if (little_endian()) {
 112                                 ((unsigned int *) &jj)[0] = x[alen - i - 1];
 113                                 ((unsigned int *) &jj)[1] = x[alen - i - 0];
 114                         } else {
 115                                 ((unsigned int *) &jj)[1] = x[alen - i - 1];
 116                                 ((unsigned int *) &jj)[0] = x[alen - i - 0];
 117                         }
 118
 119                         c = (unsigned int) (jj / kk);
 120
 121                         if (c == 0) {
 122                                 if (ge(x + n - i, b, blen)) { /* see note 1 */
 123                                         y[n - i]++;
 124                                         subf(x + n - i, b, blen);
 125                                 }
 126                                 break;
 127                         }
 128
 129                         y[n - i] += c;
 130                         mulf(t, b, blen, c);
 131                         subf(x + n - i, t, blen + 1);
 132                 }
 133         }
 134
 135         mfree(t);
 136
 137         mfree(x);
 138
 139         /* length of quotient */
 140
 141         for (i = n; i > 0; i--)
 142                 if (y[i])
 143                         break;
 144
 145         if (i == 0 && y[0] == 0) {
 146                 mfree(y);
 147                 y = mint(0);
 148         } else {
 149                 MLENGTH(y) = i + 1;
 150                 MSIGN(y) = MSIGN(a) * MSIGN(b);
 151         }
 152
 153         return y;
 154 }
 155
 156 // a = a + b
 157
 158 static void
 159 addf(unsigned int *a, unsigned int *b, int len)
 160 {
 161         int i;
 162         long long t = 0; /* can be signed or unsigned */
 163         for (i = 0; i < len; i++) {
 164                 t += (long long) a[i] + b[i];
 165                 a[i] = (unsigned int) t;
 166                 t >>= 32;
 167         }
 168 }
 169
 170 // a = a - b
 171
 172 static void
 173 subf(unsigned int *a, unsigned int *b, int len)
 174 {
 175         int i;
 176         long long t = 0; /* must be signed */
 177         for (i = 0; i < len; i++) {
 178                 t += (long long) a[i] - b[i];
 179                 a[i] = (unsigned int) t;
 180                 t >>= 32;
 181         }
 182 }
 183
 184 // a = b * c
 185
 186 // 0xffffffff + 0xffffffff * 0xffffffff == 0xffffffff00000000
 187
 188 static void
 189 mulf(unsigned int *a, unsigned int *b, int len, unsigned int c)
 190 {
 191         int i;
 192         unsigned long long t = 0; /* must be unsigned */
 193         for (i = 0; i < len; i++) {
 194                 t += (unsigned long long) b[i] * c;
 195                 a[i] = (unsigned int) t;
 196                 t >>= 32;
 197         }
 198         a[i] = (unsigned int) t;
 199 }
 200
 201 unsigned int *
 202 mmod(unsigned int *a, unsigned int *b)
 203 {
 204         int alen, blen, i, n;
 205         unsigned int c, *t, *x, *y;
 206         unsigned long long jj, kk;
 207
 208         if (MZERO(b))
 209                 stop("divide by zero");
 210
 211         if (MZERO(a))
 212                 return mint(0);
 213
 214         alen = MLENGTH(a);
 215         blen = MLENGTH(b);
 216
 217         n = alen - blen;
 218
 219         if (n < 0)
 220                 return mcopy(a);
 221
 222         x = mnew(alen + 1);
 223
 224         for (i = 0; i < alen; i++)
 225                 x[i] = a[i];
 226
 227         x[i] = 0;
 228
 229         y = mnew(n + 1);
 230
 231         t = mnew(blen + 1);
 232
 233         kk = (unsigned long long) b[blen - 1] + 1;
 234
 235         for (i = 0; i <= n; i++) {
 236
 237                 y[n - i] = 0;
 238
 239                 for (;;) {
 240
 241                         /* estimate the partial quotient */
 242
 243                         if (little_endian()) {
 244                                 ((unsigned int *) &jj)[0] = x[alen - i - 1];
 245                                 ((unsigned int *) &jj)[1] = x[alen - i - 0];
 246                         } else {
 247                                 ((unsigned int *) &jj)[1] = x[alen - i - 1];
 248                                 ((unsigned int *) &jj)[0] = x[alen - i - 0];
 249                         }
 250
 251                         c = (int) (jj / kk);
 252
 253                         if (c == 0) {
 254                                 if (ge(x + n - i, b, blen)) { /* see note 1 */
 255                                         y[n - i]++;
 256                                         subf(x + n - i, b, blen);
 257                                 }
 258                                 break;
 259                         }
 260
 261                         y[n - i] += c;
 262                         mulf(t, b, blen, c);
 263                         subf(x + n - i, t, blen + 1);
 264                 }
 265         }
 266
 267         mfree(t);
 268
 269         mfree(y);
 270
 271         /* length of remainder */
 272
 273         for (i = blen - 1; i > 0; i--)
 274                 if (x[i])
 275                         break;
 276
 277         if (i == 0 && x[0] == 0) {
 278                 mfree(x);
 279                 x = mint(0);
 280         } else {
 281                 MLENGTH(x) = i + 1;
 282                 MSIGN(x) = MSIGN(a);
 283         }
 284
 285         return x;
 286 }
 287
 288 // return both quotient and remainder of a/b
 289
 290 void
 291 mdivrem(unsigned int **q, unsigned int **r, unsigned int *a, unsigned int *b)
 292 {
 293         int alen, blen, i, n;
 294         unsigned int c, *t, *x, *y;
 295         unsigned long long jj, kk;
 296
 297         if (MZERO(b))
 298                 stop("divide by zero");
 299
 300         if (MZERO(a)) {
 301                 *q = mint(0);
 302                 *r = mint(0);
 303                 return;
 304         }
 305
 306         alen = MLENGTH(a);
 307         blen = MLENGTH(b);
 308
 309         n = alen - blen;
 310
 311         if (n < 0) {
 312                 *q = mint(0);
 313                 *r = mcopy(a);
 314                 return;
 315         }
 316
 317         x = mnew(alen + 1);
 318
 319         for (i = 0; i < alen; i++)
 320                 x[i] = a[i];
 321
 322         x[i] = 0;
 323
 324         y = mnew(n + 1);
 325
 326         t = mnew(blen + 1);
 327
 328         kk = (unsigned long long) b[blen - 1] + 1;
 329
 330         for (i = 0; i <= n; i++) {
 331
 332                 y[n - i] = 0;
 333
 334                 for (;;) {
 335
 336                         /* estimate the partial quotient */
 337
 338                         if (little_endian()) {
 339                                 ((unsigned int *) &jj)[0] = x[alen - i - 1];
 340                                 ((unsigned int *) &jj)[1] = x[alen - i - 0];
 341                         } else {
 342                                 ((unsigned int *) &jj)[1] = x[alen - i - 1];
 343                                 ((unsigned int *) &jj)[0] = x[alen - i - 0];
 344                         }
 345
 346                         c = (int) (jj / kk);
 347
 348                         if (c == 0) {
 349                                 if (ge(x + n - i, b, blen)) { /* see note 1 */
 350                                         y[n - i]++;
 351                                         subf(x + n - i, b, blen);
 352                                 }
 353                                 break;
 354                         }
 355
 356                         y[n - i] += c;
 357                         mulf(t, b, blen, c);
 358                         subf(x + n - i, t, blen + 1);
 359                 }
 360         }
 361
 362         mfree(t);
 363
 364         /* length of quotient */
 365
 366         for (i = n; i > 0; i--)
 367                 if (y[i])
 368                         break;
 369
 370         if (i == 0 && y[0] == 0) {
 371                 mfree(y);
 372                 y = mint(0);
 373         } else {
 374                 MLENGTH(y) = i + 1;
 375                 MSIGN(y) = MSIGN(a) * MSIGN(b);
 376         }
 377
 378         /* length of remainder */
 379
 380         for (i = blen - 1; i > 0; i--)
 381                 if (x[i])
 382                         break;
 383
 384         if (i == 0 && x[0] == 0) {
 385                 mfree(x);
 386                 x = mint(0);
 387         } else {
 388                 MLENGTH(x) = i + 1;
 389                 MSIGN(x) = MSIGN(a);
 390         }
 391
 392         *q = y;
 393         *r = x;
 394 }