bignum.cpp

   1 #include "stdafx.h"
   2 #include "defs.h"
   3
   4 #define MP_MIN_SIZE 2
   5 #define MP_MAX_FREE 500
   6
   7 double convert_rational_to_double(U *);
   8 double convert_bignum_to_double(unsigned int *);
   9 int ge(unsigned int *, unsigned int *, int);
  10
  11 int mtotal, mfreecount;
  12 unsigned int **free_stack;
  13
  14 unsigned int *
  15 mnew(int n)
  16 {
  17         unsigned int *p;
  18         if (n < MP_MIN_SIZE)
  19                 n = MP_MIN_SIZE;
  20         if (n == MP_MIN_SIZE && mfreecount)
  21                 p = free_stack[--mfreecount];
  22         else {
  23                 p = (unsigned int *) malloc((n + 3) * sizeof (int));
  24                 if (p == 0)
  25                         stop("malloc failure");
  26         }
  27         p[0] = n;
  28         mtotal += n;
  29         return p + 3;
  30 }
  31
  32 void
  33 mfree(unsigned int *p)
  34 {
  35         p -= 3;
  36         mtotal -= p[0];
  37         if (p[0] == MP_MIN_SIZE && mfreecount < MP_MAX_FREE)
  38                 free_stack[mfreecount++] = p;
  39         else
  40                 free(p);
  41 }
  42
  43 // convert int to bignum
  44
  45 unsigned int *
  46 mint(int n)
  47 {
  48         unsigned int *p = mnew(1);
  49         if (n < 0)
  50                 MSIGN(p) = -1;
  51         else
  52                 MSIGN(p) = 1;
  53         MLENGTH(p) = 1;
  54         p[0] = abs(n);
  55         return p;
  56 }
  57
  58 // copy bignum
  59
  60 unsigned int *
  61 mcopy(unsigned int *a)
  62 {
  63         int i;
  64         unsigned int *b;
  65
  66         b = mnew(MLENGTH(a));
  67
  68         MSIGN(b) = MSIGN(a);
  69         MLENGTH(b) = MLENGTH(a);
  70
  71         for (i = 0; i < MLENGTH(a); i++)
  72                 b[i] = a[i];
  73
  74         return b;
  75 }
  76
  77 // a >= b ?
  78
  79 int
  80 ge(unsigned int *a, unsigned int *b, int len)
  81 {
  82         int i;
  83         for (i = len - 1; i > 0; i--)
  84                 if (a[i] == b[i])
  85                         continue;
  86                 else
  87                         break;
  88         if (a[i] >= b[i])
  89                 return 1;
  90         else
  91                 return 0;
  92 }
  93
  94 void
  95 add_numbers(void)
  96 {
  97         double a, b;
  98
  99         if (isrational(stack[tos - 1]) && isrational(stack[tos - 2])) {
 100                 qadd();
 101                 return;
 102         }
 103
 104         save();
 105
 106         p2 = pop();
 107         p1 = pop();
 108
 109         if (isdouble(p1))
 110                 a = p1->u.d;
 111         else
 112                 a = convert_rational_to_double(p1);
 113
 114         if (isdouble(p2))
 115                 b = p2->u.d;
 116         else
 117                 b = convert_rational_to_double(p2);
 118
 119         push_double(a + b);
 120
 121         restore();
 122 }
 123
 124 void
 125 subtract_numbers(void)
 126 {
 127         double a, b;
 128
 129         if (isrational(stack[tos - 1]) && isrational(stack[tos - 2])) {
 130                 qsub();
 131                 return;
 132         }
 133
 134         save();
 135
 136         p2 = pop();
 137         p1 = pop();
 138
 139         if (isdouble(p1))
 140                 a = p1->u.d;
 141         else
 142                 a = convert_rational_to_double(p1);
 143
 144         if (isdouble(p2))
 145                 b = p2->u.d;
 146         else
 147                 b = convert_rational_to_double(p2);
 148
 149         push_double(a - b);
 150
 151         restore();
 152 }
 153
 154 void
 155 multiply_numbers(void)
 156 {
 157         double a, b;
 158
 159         if (isrational(stack[tos - 1]) && isrational(stack[tos - 2])) {
 160                 qmul();
 161                 return;
 162         }
 163
 164         save();
 165
 166         p2 = pop();
 167         p1 = pop();
 168
 169         if (isdouble(p1))
 170                 a = p1->u.d;
 171         else
 172                 a = convert_rational_to_double(p1);
 173
 174         if (isdouble(p2))
 175                 b = p2->u.d;
 176         else
 177                 b = convert_rational_to_double(p2);
 178
 179         push_double(a * b);
 180
 181         restore();
 182 }
 183
 184 void
 185 divide_numbers(void)
 186 {
 187         double a, b;
 188
 189         if (isrational(stack[tos - 1]) && isrational(stack[tos - 2])) {
 190                 qdiv();
 191                 return;
 192         }
 193
 194         save();
 195
 196         p2 = pop();
 197         p1 = pop();
 198
 199         if (iszero(p2))
 200                 stop("divide by zero");
 201
 202         if (isdouble(p1))
 203                 a = p1->u.d;
 204         else
 205                 a = convert_rational_to_double(p1);
 206
 207         if (isdouble(p2))
 208                 b = p2->u.d;
 209         else
 210                 b = convert_rational_to_double(p2);
 211
 212         push_double(a / b);
 213
 214         restore();
 215 }
 216
 217 void
 218 invert_number(void)
 219 {
 220         unsigned int *a, *b;
 221
 222         save();
 223
 224         p1 = pop();
 225
 226         if (iszero(p1))
 227                 stop("divide by zero");
 228
 229         if (isdouble(p1)) {
 230                 push_double(1 / p1->u.d);
 231                 restore();
 232                 return;
 233         }
 234
 235         a = mcopy(p1->u.q.a);
 236         b = mcopy(p1->u.q.b);
 237
 238         MSIGN(b) = MSIGN(a);
 239         MSIGN(a) = 1;
 240
 241         p1 = alloc();
 242
 243         p1->k = NUM;
 244
 245         p1->u.q.a = b;
 246         p1->u.q.b = a;
 247
 248         push(p1);
 249
 250         restore();
 251 }
 252
 253 int
 254 compare_rationals(U *a, U *b)
 255 {
 256         int t;
 257         unsigned int *ab, *ba;
 258         ab = mmul(a->u.q.a, b->u.q.b);
 259         ba = mmul(a->u.q.b, b->u.q.a);
 260         t = mcmp(ab, ba);
 261         mfree(ab);
 262         mfree(ba);
 263         return t;
 264 }
 265
 266 int
 267 compare_numbers(U *a, U *b)
 268 {
 269         double x, y;
 270         if (isrational(a) && isrational(b))
 271                 return compare_rationals(a, b);
 272         if (isdouble(a))
 273                 x = a->u.d;
 274         else
 275                 x = convert_rational_to_double(a);
 276         if (isdouble(b))
 277                 y = b->u.d;
 278         else
 279                 y = convert_rational_to_double(b);
 280         if (x < y)
 281                 return -1;
 282         if (x > y)
 283                 return 1;
 284         return 0;
 285 }
 286
 287 void
 288 negate_number(void)
 289 {
 290         save();
 291         p1 = pop();
 292         if (iszero(p1)) {
 293                 push(p1);
 294                 restore();
 295                 return;
 296         }
 297         switch (p1->k) {
 298         case NUM:
 299                 p2 = alloc();
 300                 p2->k = NUM;
 301                 p2->u.q.a = mcopy(p1->u.q.a);
 302                 p2->u.q.b = mcopy(p1->u.q.b);
 303                 MSIGN(p2->u.q.a) *= -1;
 304                 push(p2);
 305                 break;
 306         case DOUBLE:
 307                 push_double(-p1->u.d);
 308                 break;
 309         default:
 310                 stop("bug caught in mp_negate_number");
 311                 break;
 312         }
 313         restore();
 314 }
 315
 316 void
 317 bignum_truncate(void)
 318 {
 319         unsigned int *a;
 320
 321         save();
 322
 323         p1 = pop();
 324
 325         a = mdiv(p1->u.q.a, p1->u.q.b);
 326
 327         p1 = alloc();
 328
 329         p1->k = NUM;
 330
 331         p1->u.q.a = a;
 332         p1->u.q.b = mint(1);
 333
 334         push(p1);
 335
 336         restore();
 337 }
 338
 339 void
 340 mp_numerator(void)
 341 {
 342         save();
 343
 344         p1 = pop();
 345
 346         if (p1->k != NUM) {
 347                 push(one);
 348                 restore();
 349                 return;
 350         }
 351
 352         p2 = alloc();
 353
 354         p2->k = NUM;
 355
 356         p2->u.q.a = mcopy(p1->u.q.a);
 357         p2->u.q.b = mint(1);
 358
 359         push(p2);
 360
 361         restore();
 362 }
 363
 364 void
 365 mp_denominator(void)
 366 {
 367         save();
 368
 369         p1 = pop();
 370
 371         if (p1->k != NUM) {
 372                 push(one);
 373                 restore();
 374                 return;
 375         }
 376
 377         p2 = alloc();
 378
 379         p2->k = NUM;
 380
 381         p2->u.q.a = mcopy(p1->u.q.b);
 382         p2->u.q.b = mint(1);
 383
 384         push(p2);
 385
 386         restore();
 387 }
 388
 389 void
 390 bignum_power_number(int expo)
 391 {
 392         unsigned int *a, *b, *t;
 393
 394         save();
 395
 396         p1 = pop();
 397
 398         a = mpow(p1->u.q.a, abs(expo));
 399         b = mpow(p1->u.q.b, abs(expo));
 400
 401         if (expo < 0) {
 402                 t = a;
 403                 a = b;
 404                 b = t;
 405                 MSIGN(a) = MSIGN(b);
 406                 MSIGN(b) = 1;
 407         }
 408
 409         p1 = alloc();
 410
 411         p1->k = NUM;
 412
 413         p1->u.q.a = a;
 414         p1->u.q.b = b;
 415
 416         push(p1);
 417
 418         restore();
 419 }
 420
 421 double
 422 convert_bignum_to_double(unsigned int *p)
 423 {
 424         int i;
 425         double d = 0.0;
 426         for (i = MLENGTH(p) - 1; i >= 0; i--)
 427                 d = 4294967296.0 * d + p[i];
 428         if (MSIGN(p) == -1)
 429                 d = -d;
 430         return d;
 431 }
 432
 433 double
 434 convert_rational_to_double(U *p)
 435 {
 436         int i, n, na, nb;
 437         double a = 0.0, b = 0.0;
 438
 439         na = MLENGTH(p->u.q.a);
 440         nb = MLENGTH(p->u.q.b);
 441
 442         if (na < nb)
 443                 n = na;
 444         else
 445                 n = nb;
 446
 447         for (i = 0; i < n; i++) {
 448                 a = a / 4294967296.0 + p->u.q.a[i];
 449                 b = b / 4294967296.0 + p->u.q.b[i];
 450         }
 451
 452         if (na > nb)
 453                 for (i = nb; i < na; i++) {
 454                         a = a / 4294967296.0 + p->u.q.a[i];
 455                         b = b / 4294967296.0;
 456                 }
 457
 458         if (na < nb)
 459                 for (i = na; i < nb; i++) {
 460                         a = a / 4294967296.0;
 461                         b = b / 4294967296.0 + p->u.q.b[i];
 462                 }
 463
 464         if (MSIGN(p->u.q.a) == -1)
 465                 a = -a;
 466
 467         return a / b;
 468 }
 469
 470 void
 471 push_integer(int n)
 472 {
 473         save();
 474         p1 = alloc();
 475         p1->k = NUM;
 476         p1->u.q.a = mint(n);
 477         p1->u.q.b = mint(1);
 478         push(p1);
 479         restore();
 480 }
 481
 482 void
 483 push_double(double d)
 484 {
 485         save();
 486         p1 = alloc();
 487         p1->k = DOUBLE;
 488         p1->u.d = d;
 489         push(p1);
 490         restore();
 491 }
 492
 493 void
 494 push_rational(int a, int b)
 495 {
 496         save();
 497         p1 = alloc();
 498         p1->k = NUM;
 499         p1->u.q.a = mint(a);
 500         p1->u.q.b = mint(b);
 501         /* FIXME -- normalize */
 502         push(p1);
 503         restore();
 504 }
 505
 506 int
 507 pop_integer(void)
 508 {
 509         int n;
 510
 511         save();
 512
 513         p1 = pop();
 514
 515         switch (p1->k) {
 516
 517         case NUM:
 518                 if (isinteger(p1) && MLENGTH(p1->u.q.a) == 1) {
 519                         n = p1->u.q.a[0];
 520                         if (n & 0x80000000)
 521                                 n = 0x80000000;
 522                         else
 523                                 n *= MSIGN(p1->u.q.a);
 524                 } else
 525                         n = 0x80000000;
 526                 break;
 527
 528         case DOUBLE:
 529                 n = (int) p1->u.d;
 530                 if ((double) n != p1->u.d)
 531                         n = 0x80000000;
 532                 break;
 533
 534         default:
 535                 n = 0x80000000;
 536                 break;
 537         }
 538
 539         restore();
 540         return n;
 541 }
 542
 543 void
 544 print_double(U *p, int flag)
 545 {
 546         static char buf[80];
 547         sprintf(buf, "%g", p->u.d);
 548         if (flag == 1 && *buf == '-')
 549                 print_str(buf + 1);
 550         else
 551                 print_str(buf);
 552 }
 553
 554 void
 555 bignum_scan_integer(char *s)
 556 {
 557         unsigned int *a;
 558         char sign;
 559
 560         save();
 561
 562         sign = *s;
 563
 564         if (sign == '+' || sign == '-')
 565                 s++;
 566
 567         a = mscan(s);
 568
 569         p1 = alloc();
 570
 571         p1->k = NUM;
 572
 573         p1->u.q.a = a;
 574         p1->u.q.b = mint(1);
 575
 576         push(p1);
 577
 578         if (sign == '-')
 579                 negate();
 580
 581         restore();
 582 }
 583
 584 #define FLT_MIN (1e-999)
 585 #define FLT_MAX (9.999999999999999e999)
 586 // the following strtod isn't used. but if I remove it, GCC's linker starts giving errors somewhere inside libm
 587 // *magic* DO NOT TOUCH *magic*
 588 // this strtod was deprecated because it didn't have enough precision
 589 double amystrtod(char *s, char **str_end) {
 590     // TODO handle exponential format, hex format, inf, nan
 591     double r = 0.0;
 592     int negative = 0;
 593         union { double rr; unsigned long long dd; } raw;
 594
 595         while (isspace(*s)) s++; //skip leading whitespace
 596         if ((s[0] == 'i' || s[0] == 'I') && (s[1] == 'n' || s[1] == 'N') && (s[2] == 'f' || s[2] == 'F')) {
 597                 if (str_end != NULL)
 598                         *str_end = s+3;
 599                 raw.dd = 0x7FF0000000000000;            //positive infinity
 600                 return raw.rr;
 601         }
 602         if ((s[0] == 'n' || s[0] == 'N') && (s[1] == 'a' || s[1] == 'A') && (s[2] == 'n' || s[2] == 'N')) {
 603                 if (str_end != NULL)
 604                         *str_end = s+3;
 605                 raw.dd = 0x7FFFFFFFFFFFFFFF;            //QNaN
 606                 return raw.rr;
 607         }
 608     if (!isdigit(*s) && *s != '-' && *s != '+' && *s != '.') {
 609         if (str_end != NULL)
 610             *str_end = (char *)s;
 611         return 0;
 612     }
 613
 614     switch (*s)
 615     {
 616         case '-':
 617             negative = 1;
 618             // Deliberate fallthrough
 619         case '+':
 620             s++;
 621             break;
 622     }
 623
 624     while (isdigit(*s))
 625     {
 626         r *= 10.0;
 627         r += *s++ - '0';
 628     }
 629     if (*s == '.')
 630     {
 631         float f = 10.0f;
 632         s++;
 633         while (isdigit(*s))
 634         {
 635             r += (*s - '0')/f;
 636             f *= 10.0f;
 637             s++;
 638         }
 639     }
 640
 641         if (*s == 'e' || *s == 'E') {
 642                 int exponent = 0;
 643                 bool sign = false;
 644                 if (*(s+1) == '-' || *(s+1) == '+') {
 645                         if (*(s+1) == '-') sign = true;
 646                         if (isdigit(*(s+2))) {
 647                                 s+=2;
 648                                 while(isdigit(*s)) {
 649                                         exponent = (exponent*10) + (*s - '0');
 650                                         s++;
 651                                 }
 652                                 r = pow(r,((sign)?-exponent:exponent));
 653                         }
 654                 } else if (isdigit(*(s+1))) {
 655                         s++;
 656                         while(isdigit(*s)) {
 657                                 exponent = (exponent*10) + (*s - '0');
 658                                 s++;
 659                         }
 660                         r = pow(r,exponent);
 661                 }
 662         }
 663
 664         if (str_end != NULL)
 665                 *str_end = (char *)s;
 666
 667     // Portable? Nope. Fast? Yup.
 668     raw.rr = r;
 669     raw.dd |= (unsigned long long)negative << 63;
 670
 671         if(raw.rr >= FLT_MIN || raw.rr <= -FLT_MIN)
 672                 return raw.rr;
 673         else return 0;
 674 }
 675 #define DBL_MAX_EXP 999
 676 #define DBL_MIN_EXP (-999)
 677 //
 678 // strtod.c
 679 //
 680 // Convert string to double
 681 //
 682 // Copyright (C) 2002 Michael Ringgaard. All rights reserved.
 683 //
 684 // Redistribution and use in source and binary forms, with or without
 685 // modification, are permitted provided that the following conditions
 686 // are met:
 687 //
 688 // 1. Redistributions of source code must retain the above copyright
 689 //    notice, this list of conditions and the following disclaimer.
 690 // 2. Redistributions in binary form must reproduce the above copyright
 691 //    notice, this list of conditions and the following disclaimer in the
 692 //    documentation and/or other materials provided with the distribution.
 693 // 3. Neither the name of the project nor the names of its contributors
 694 //    may be used to endorse or promote products derived from this software
 695 //    without specific prior written permission.
 696 //
 697 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 698 // ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 699 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 700 // ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
 701 // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 702 // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
 703 // OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 704 // HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 705 // LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
 706 // OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 707 // SUCH DAMAGE.
 708 //
 709 double mystrtod(const char *str, char **endptr) {
 710   double number;
 711   int exponent;
 712   int negative;
 713   char *p = (char *) str;
 714   double p10;
 715   int n;
 716   int num_digits;
 717   int num_decimals;
 718
 719   // Skip leading whitespace
 720   while (isspace(*p)) p++;
 721
 722   // Handle optional sign
 723   negative = 0;
 724   switch (*p) {
 725     case '-': negative = 1; // Fall through to increment position
 726     case '+': p++;
 727   }
 728
 729   number = 0.;
 730   exponent = 0;
 731   num_digits = 0;
 732   num_decimals = 0;
 733
 734   // Process string of digits
 735   while (isdigit(*p)) {
 736     number = number * 10. + (*p - '0');
 737     p++;
 738     num_digits++;
 739   }
 740
 741   // Process decimal part
 742   if (*p == '.') {
 743     p++;
 744
 745     while (isdigit(*p)) {
 746       number = number * 10. + (*p - '0');
 747       p++;
 748       num_digits++;
 749       num_decimals++;
 750     }
 751
 752     exponent -= num_decimals;
 753   }
 754
 755   if (num_digits == 0) {
 756     errno = ERANGE;
 757     return 0.0;
 758   }
 759
 760   // Correct for sign
 761   if (negative) number = -number;
 762
 763   // Process an exponent string
 764   if (*p == 'e' || *p == 'E') {
 765     // Handle optional sign
 766     negative = 0;
 767     switch (*++p) {
 768       case '-': negative = 1;   // Fall through to increment pos
 769       case '+': p++;
 770     }
 771
 772     // Process string of digits
 773     n = 0;
 774     while (isdigit(*p)) {
 775       n = n * 10 + (*p - '0');
 776       p++;
 777     }
 778
 779     if (negative) {
 780       exponent -= n;
 781     } else {
 782       exponent += n;
 783     }
 784   }
 785
 786   if (exponent < DBL_MIN_EXP  || exponent > DBL_MAX_EXP) {
 787     errno = ERANGE;
 788     return HUGE_VAL;
 789   }
 790
 791   // Scale the result
 792   p10 = 10.;
 793   n = exponent;
 794   if (n < 0) n = -n;
 795   while (n) {
 796     if (n & 1) {
 797       if (exponent < 0) {
 798         number /= p10;
 799       } else {
 800         number *= p10;
 801       }
 802     }
 803     n >>= 1;
 804     p10 *= p10;
 805   }
 806
 807   if (number == HUGE_VAL) errno = ERANGE;
 808   if (endptr) *endptr = p;
 809
 810   return number;
 811 }
 812 void
 813 bignum_scan_float(char *s)
 814 {
 815         //push_double(atof(s));
 816         push_double(mystrtod((char*)s, NULL));
 817 }
 818
 819 // print as unsigned
 820
 821 void
 822 print_number(U *p)
 823 {
 824         char *s;
 825         static char buf[100];
 826         switch (p->k) {
 827         case NUM:
 828                 s = mstr(p->u.q.a);
 829                 if (*s == '+' || *s == '-')
 830                         s++;
 831                 print_str(s);
 832                 if (isfraction(p)) {
 833                         print_str("/");
 834                         s = mstr(p->u.q.b);
 835                         print_str(s);
 836                 }
 837                 break;
 838         case DOUBLE:
 839                 sprintf(buf, "%.10g", p->u.d);
 840                 if (*buf == '+' || *buf == '-')
 841                         print_str(buf + 1);
 842                 else
 843                         print_str(buf);
 844                 break;
 845         default:
 846                 break;
 847         }
 848 }
 849
 850 void
 851 gcd_numbers(void)
 852 {
 853         save();
 854
 855         p2 = pop();
 856         p1 = pop();
 857
 858 //      if (!isinteger(p1) || !isinteger(p2))
 859 //              stop("integer args expected for gcd");
 860
 861         p3 = alloc();
 862
 863         p3->k = NUM;
 864
 865         p3->u.q.a = mgcd(p1->u.q.a, p2->u.q.a);
 866         p3->u.q.b = mgcd(p1->u.q.b, p2->u.q.b);
 867
 868         MSIGN(p3->u.q.a) = 1;
 869
 870         push(p3);
 871
 872         restore();
 873 }
 874
 875 double
 876 pop_double(void)
 877 {
 878         double d;
 879         save();
 880         p1 = pop();
 881         switch (p1->k) {
 882         case NUM:
 883                 d = convert_rational_to_double(p1);
 884                 break;
 885         case DOUBLE:
 886                 d = p1->u.d;
 887                 break;
 888         default:
 889                 d = 0.0;
 890                 break;
 891         }
 892         restore();
 893         return d;
 894 }
 895
 896 void
 897 bignum_float(void)
 898 {
 899         double d;
 900         d = convert_rational_to_double(pop());
 901         push_double(d);
 902 }
 903
 904 static unsigned int *__factorial(int);
 905
 906 void
 907 bignum_factorial(int n)
 908 {
 909         save();
 910         p1 = alloc();
 911         p1->k = NUM;
 912         p1->u.q.a = __factorial(n);
 913         p1->u.q.b = mint(1);
 914         push(p1);
 915         restore();
 916 }
 917
 918 static unsigned int *
 919 __factorial(int n)
 920 {
 921         int i;
 922         unsigned int *a, *b, *t;
 923
 924         if (n == 0 || n == 1) {
 925                 a = mint(1);
 926                 return a;
 927         }
 928
 929         a = mint(2);
 930
 931         b = mint(0);
 932
 933         for (i = 3; i <= n; i++) {
 934                 b[0] = (unsigned int) i;
 935                 t = mmul(a, b);
 936                 mfree(a);
 937                 a = t;
 938         }
 939
 940         mfree(b);
 941
 942         return a;
 943 }
 944
 945 static unsigned int mask[32] = {
 946         0x00000001,
 947         0x00000002,
 948         0x00000004,
 949         0x00000008,
 950         0x00000010,
 951         0x00000020,
 952         0x00000040,
 953         0x00000080,
 954         0x00000100,
 955         0x00000200,
 956         0x00000400,
 957         0x00000800,
 958         0x00001000,
 959         0x00002000,
 960         0x00004000,
 961         0x00008000,
 962         0x00010000,
 963         0x00020000,
 964         0x00040000,
 965         0x00080000,
 966         0x00100000,
 967         0x00200000,
 968         0x00400000,
 969         0x00800000,
 970         0x01000000,
 971         0x02000000,
 972         0x04000000,
 973         0x08000000,
 974         0x10000000,
 975         0x20000000,
 976         0x40000000,
 977         0x80000000,
 978 };
 979
 980 void
 981 mp_set_bit(unsigned int *x, unsigned int k)
 982 {
 983         x[k / 32] |= mask[k % 32];
 984 }
 985
 986 void
 987 mp_clr_bit(unsigned int *x, unsigned int k)
 988 {
 989         x[k / 32] &= ~mask[k % 32];
 990 }
 991
 992 void
 993 mshiftright(unsigned int *a)
 994 {
 995         int c, i, n;
 996         n = MLENGTH(a);
 997         c = 0;
 998         for (i = n - 1; i >= 0; i--)
 999                 if (a[i] & 1) {
1000                         a[i] = (a[i] >> 1) | c;
1001                         c = 0x80000000;
1002                 } else {
1003                         a[i] = (a[i] >> 1) | c;
1004                         c = 0;
1005                 }
1006         if (n > 1 && a[n - 1] == 0)
1007                 MLENGTH(a) = n - 1;
1008 }