ghc/interpreter/sainteger.c

   1
   2 /* --------------------------------------------------------------------------
   3  * Yet another implementation of Integer
   4  *
   5  * Copyright (c) Glasgow University, 1999.
   6  * All rights reserved. See NOTICE for details and conditions of use etc...
   7  * ------------------------------------------------------------------------*/
   8
   9 #include <stdlib.h>
  10 #include <stdio.h>
  11 #include <assert.h>
  12 #include <ctype.h>
  13
  14 #include "sainteger.h"
  15
  16
  17 /* --------------------------------------------------------------------------
  18  * Local fns
  19  * ------------------------------------------------------------------------*/
  20
  21 typedef unsigned char uchar;
  22 typedef unsigned short ush;
  23
  24
  25 static int maxused_add ( B*, B* );
  26 static int maxused_sub ( B*, B* );
  27 static int maxused_mul ( B*, B* );
  28 static int maxused_qrm ( B*, B* );
  29 static int maxused_neg ( B* );
  30
  31 static int  ucmp ( B*, B* );
  32 static void uadd ( B*, B*, B* );
  33 static void usub ( B*, B*, B* );
  34 static void umul ( B*, B*, B* );
  35 static void uqrm ( B*, B*, B*, B* );
  36
  37 /*#define DEBUG_SAINTEGER*/
  38 /*#define DEBUG_SAINTEGER_UQRM*/
  39
  40
  41 #ifdef DEBUG_SAINTEGER
  42 #define myassert(zzzz) assert(zzzz)
  43 #else
  44 #define myassert(zzzz) /* */
  45 #endif
  46
  47
  48 /* --------------------------------------------------------------------------
  49  * Basics
  50  * ------------------------------------------------------------------------*/
  51
  52 void pp ( B* x )
  53 {
  54    int i;
  55    printf ( "sign=%2d  used=%d  size=%d   ", x->sign, x->used, x->size );
  56    for (i = x->used-1; i >= 0; i--)
  57       printf ( "%2x ", (int)(x->stuff[i]) );
  58    printf ( "\n" );
  59 }
  60
  61
  62 static int sane ( B* x )
  63 {
  64    int i;
  65
  66    if (x->sign == 0 && x->used != 0) return 0;
  67    if (x->sign != -1 && x->sign != 0 && x->sign != 1) return 0;
  68
  69    if (x->used < 0) return 0;
  70    if (x->size < 0) return 0;
  71    if (x->used > x->size) return 0;
  72    if (x->used == 0) return 1;
  73    if (x->stuff[x->used-1] == 0) return 0;
  74    for (i = 0; i < x->used; i++)
  75       if (x->stuff[i] >= B_BASE) return 0;
  76    return 1;
  77 }
  78
  79
  80 int is_sane ( B* x )
  81 {
  82    return sane(x);
  83 }
  84
  85
  86 static void u_renormalise ( B* b )
  87 {
  88    while (b->used > 0 && b->stuff[b->used-1] == 0) b->used--;
  89    if (b->used == 0) b->sign = 0; else b->sign = 1;
  90 }
  91
  92
  93 void do_renormalise ( B* b )
  94 {
  95    while (b->used > 0 && b->stuff[b->used-1] == 0) b->used--;
  96    if (b->used == 0) b->sign = 0;
  97 }
  98
  99 /* --------------------------------------------------------------------------
 100  * Size of things
 101  * ------------------------------------------------------------------------*/
 102
 103 static int maxused_add ( B* x, B* y )
 104 {
 105    myassert(sane(x));
 106    myassert(sane(y));
 107    return 1 + (x->used > y->used ? x->used : y->used);
 108 }
 109
 110 static int maxused_sub ( B* x, B* y )
 111 {
 112    myassert(sane(x));
 113    myassert(sane(y));
 114    return 1 + (x->used > y->used ? x->used : y->used);
 115 }
 116
 117 static int maxused_mul ( B* x, B* y )
 118 {
 119    myassert(sane(x));
 120    myassert(sane(y));
 121    return x->used + y->used;
 122 }
 123
 124 static int maxused_qrm ( B* x, B* y )
 125 {
 126    myassert(sane(x));
 127    myassert(sane(y));
 128    return (x->used > y->used ? x->used : y->used);
 129 }
 130
 131 static int maxused_neg ( B* x )
 132 {
 133    myassert(sane(x));
 134    return x->used;
 135 }
 136
 137
 138 /* quick, safe approx */
 139 static int maxused_fromInt ( int sizeof_int )
 140 {
 141    if (B_BASE == 256)  return     sizeof_int;
 142    if (B_BASE >= 16)   return 2 * sizeof_int;
 143    if (B_BASE >= 4)    return 4 * sizeof_int;
 144    /* (B_BASE >= 2) */ return 8 * sizeof_int;
 145 }
 146
 147 /* ditto */
 148 static int maxused_fromStr ( char* str )
 149 {
 150    int nd = 0;
 151    if (*str == '-') str++;
 152    while (isdigit((int)(*str))) { str++; nd++; };
 153
 154    if (B_BASE >= 100) return ((nd+1) / 2);
 155    if (B_BASE >= 10)  return nd;
 156    /* (B_BASE >= 2)*/ return 4 * nd;
 157 }
 158
 159
 160 int size_add ( B* x, B* y )
 161 {
 162    return sizeof(B) + maxused_add(x,y);
 163 }
 164
 165 int size_sub ( B* x, B* y )
 166 {
 167    return sizeof(B) + maxused_sub(x,y);
 168 }
 169
 170 int size_mul ( B* x, B* y )
 171 {
 172    return sizeof(B) + maxused_mul(x,y);
 173 }
 174
 175 int size_qrm ( B* x, B* y )
 176 {
 177    return sizeof(B) + maxused_qrm(x,y);
 178 }
 179
 180 int size_neg ( B* x )
 181 {
 182    return sizeof(B) + maxused_neg(x);
 183 }
 184
 185 int size_fromInt ( void )
 186 {
 187    int sizeof_int = sizeof(int);
 188    return sizeof(B) + maxused_fromInt ( sizeof_int );
 189 }
 190
 191 int size_fromWord ( void )
 192 {
 193    int sizeof_word = sizeof(unsigned int);
 194    return sizeof(B) + maxused_fromInt ( sizeof_word );
 195 }
 196
 197 int size_fromStr ( char* str )
 198 {
 199    return sizeof(B) + maxused_fromStr ( str );
 200 }
 201
 202 int size_fltmantissa ( void )
 203 {
 204    return sizeof(B) + sizeof(float);
 205 }
 206
 207 int size_dblmantissa ( void )
 208 {
 209    return sizeof(B) + sizeof(double);
 210 }
 211
 212
 213 /* --------------------------------------------------------------------------
 214  * Conversions
 215  * ------------------------------------------------------------------------*/
 216
 217 void do_fromInt  ( int n, int sizeRes, B* res )
 218 {
 219
 220    res->size = sizeRes - sizeof(B);
 221    res->sign = res->used = 0;
 222    if (n == 0) { myassert(sane(res)); return; };
 223    if (n < 0) res->sign = -1; else res->sign = 1;
 224    if (n < 0) n = -n;
 225
 226    while (n != 0) {
 227       res->stuff[res->used] = (uchar)(n % B_BASE);
 228       n /= B_BASE;
 229       res->used++;
 230    }
 231    myassert(sane(res));
 232 }
 233
 234 void do_fromWord  ( unsigned int n, int sizeRes, B* res )
 235 {
 236
 237    res->size = sizeRes - sizeof(B);
 238    res->sign = res->used = 0;
 239    if (n == 0) { myassert(sane(res)); return; };
 240    res->sign = 1;
 241
 242    while (n != 0) {
 243       res->stuff[res->used] = (uchar)(n % B_BASE);
 244       n /= B_BASE;
 245       res->used++;
 246    }
 247    myassert(sane(res));
 248 }
 249
 250 /* NOTE: This only works currectly if B_BASE >= 10 */
 251 void do_fromStr ( char* str, int sizeRes, B* res )
 252 {
 253    int sign, d, t, j, carry;
 254
 255    res->size = sizeRes - sizeof(B);
 256    res->sign = res->used = 0;
 257    sign = 1;
 258    if (*str == '-') { str++; sign = -1; };
 259
 260    while (isdigit((int)(*str))) {
 261
 262       /* multiply res by 10 */
 263       carry = 0;
 264       for (j = 0; j < res->used; j++) {
 265          t = 10 * res->stuff[j] + carry;
 266          res->stuff[j] = t % B_BASE;
 267          carry = t / B_BASE;
 268       }
 269       myassert(carry < B_BASE);
 270       if (carry > 0)
 271          res->stuff[res->used++] = carry;
 272
 273       /* add a digit on */
 274       d = *str - '0';
 275       str++;
 276
 277       carry = d;
 278       for (j = 0; j < res->used; j++) {
 279          carry += res->stuff[j];
 280          res->stuff[j] = carry % B_BASE;
 281          carry /= B_BASE;
 282          if (carry == 0) break;
 283       }
 284       if (carry > 0)
 285          res->stuff[res->used++] = carry;
 286    }
 287
 288    res->sign = sign;
 289    myassert(sane(res));
 290 }
 291
 292 int do_toInt ( B* x )
 293 {
 294    int i, d, res;
 295    if (x->sign == 0) return 0;
 296    res = 0;
 297    for (i = x->used-1; i >= 0; i--) {
 298       d = x->stuff[i];
 299       res = res * B_BASE + d;
 300    }
 301    if (x->sign < 0) res = -res;
 302    return res;
 303 }
 304
 305 unsigned int do_toWord ( B* x )
 306 {
 307    int i, d;
 308    unsigned int res;
 309    if (x->sign == 0) return 0;
 310    res = 0;
 311    for (i = x->used-1; i >= 0; i--) {
 312       d = x->stuff[i];
 313       res = res * B_BASE + d;
 314    }
 315    return res;
 316 }
 317
 318 float do_toFloat ( B* x )
 319 {
 320    int i, d;
 321    float res;
 322    if (x->sign == 0) return 0.0;
 323    res = 0.0;
 324    for (i = x->used-1; i >= 0; i--) {
 325       d = x->stuff[i];
 326       res = res * B_BASE_FLT + d;
 327    }
 328    if (x->sign < 0) res = -res;
 329    return res;
 330 }
 331
 332 double do_toDouble ( B* x )
 333 {
 334    int i, d;
 335    double res;
 336    if (x->sign == 0) return 0.0;
 337    res = 0.0;
 338    for (i = x->used-1; i >= 0; i--) {
 339       d = x->stuff[i];
 340       res = res * B_BASE_FLT + d;
 341    }
 342    if (x->sign < 0) res = -res;
 343    return res;
 344 }
 345
 346
 347 /* --------------------------------------------------------------------------
 348  * Signed ops
 349  * ------------------------------------------------------------------------*/
 350
 351 /* A helper for signed + and -.  sdiff(x,y) ignores the signs of x and y
 352    sets p to the signed value abs(x)-abs(y).
 353 */
 354 static void sdiff ( B* x, B* y, B* res )
 355 {
 356    int t;
 357    myassert(sane(x));
 358    myassert(sane(y));
 359    myassert(res->size == maxused_sub(x,y));
 360    t = ucmp(x,y);
 361    if (t == 0) { res->sign = res->used = 0; return; }
 362    if (t == -1) {
 363       /* x < y */
 364       usub(y,x,res);
 365       res->sign = -1;
 366    } else {
 367       /* x > y */
 368       usub(x,y,res);
 369       res->sign = 1;
 370    }
 371    myassert(sane(res));
 372 }
 373
 374 int do_getsign ( B* x )
 375 {
 376    myassert(sane(x));
 377    return x->sign;
 378 }
 379
 380 void do_neg ( B* x, int sizeRes, B* res )
 381 {
 382    int i;
 383    myassert(sane(x));
 384    res->size = sizeRes - sizeof(B);
 385    res->used = x->used;
 386    for (i = 0; i < x->used; i++)
 387       res->stuff[i] = x->stuff[i];
 388    res->sign = - (x->sign);
 389 }
 390
 391 void do_add ( B* x, B* y, int sizeRes, B* res )
 392 {
 393    myassert(sane(x));
 394    myassert(sane(y));
 395    res->size = sizeRes - sizeof(B);
 396    res->used = res->sign = 0;
 397
 398    if ( (x->sign >= 0 && y->sign >= 0) ||
 399         (x->sign < 0  && y->sign < 0)) {
 400       /* same sign; add magnitude and clone sign */
 401       uadd(x,y,res);
 402       if (x->sign < 0 && res->sign != 0) res->sign = -1;
 403    }
 404    else
 405    /* signs differ; employ sdiff */
 406    if (x->sign >= 0 && y->sign < 0) {
 407       sdiff(x,y,res);
 408    } else {
 409       myassert(x->sign < 0 && y->sign >= 0);
 410       sdiff(y,x,res);
 411    }
 412    myassert(sane(res));
 413 }
 414
 415 void do_sub ( B* x, B* y, int sizeRes, B* res )
 416 {
 417    myassert(sane(x));
 418    myassert(sane(y));
 419    res->size = sizeRes - sizeof(B);
 420    res->used = res->sign = 0;
 421
 422    if ( (x->sign >= 0 && y->sign < 0) ||
 423         (x->sign < 0  && y->sign >= 0)) {
 424       /* opposite signs; add magnitudes and clone sign of x */
 425       uadd(x,y,res);
 426       myassert(res->sign != 0);
 427       if (x->sign < 0) res->sign = -1;
 428    }
 429    else
 430    /* signs are the same; employ sdiff */
 431    if (x->sign >= 0 && y->sign >= 0) {
 432       sdiff(x,y,res);
 433    } else {
 434       myassert(x->sign < 0 && y->sign < 0);
 435       sdiff(y,x,res);
 436    }
 437    myassert(sane(res));
 438 }
 439
 440
 441 void do_mul ( B* x, B* y, int sizeRes, B* res )
 442 {
 443    myassert(sane(x));
 444    myassert(sane(y));
 445    res->size = sizeRes - sizeof(B);
 446    res->used = res->sign = 0;
 447
 448    if (x->sign == 0 || y->sign == 0) {
 449       res->sign = res->used = 0;
 450       myassert(sane(res));
 451       return;
 452    }
 453    umul(x,y,res);
 454    if (x->sign != y->sign) res->sign = -1;
 455    myassert(sane(res));
 456 }
 457
 458
 459 void do_qrm ( B* x, B* y, int sizeRes, B* q, B* r )
 460 {
 461    myassert(sane(x));
 462    myassert(sane(y));
 463
 464    q->size = r->size = sizeRes - sizeof(B);
 465    q->used = r->used = q->sign = r->sign = 0;
 466
 467    if (y->sign == 0) {
 468       fprintf(stderr, "do_qrm: division by zero -- exiting now!\n");
 469       exit(1);
 470       return;
 471    }
 472
 473    if (x->sign == 0) {
 474       q->used = r->used = q->sign = r->sign = 0;
 475       myassert(sane(q)); myassert(sane(r));
 476       return;
 477    }
 478
 479    uqrm ( x, y, q, r );
 480    if (x->sign != y->sign && q->sign != 0) q->sign = -1;
 481    if (x->sign == -1 && r->sign != 0) r->sign = -1;
 482
 483    myassert(sane(q)); myassert(sane(r));
 484 }
 485
 486 int do_cmp ( B* x, B* y )
 487 {
 488    if (!sane(x))
 489       pp(x);
 490    myassert(sane(x));
 491    myassert(sane(y));
 492    if (x->sign < y->sign) return -1;
 493    if (x->sign > y->sign) return 1;
 494    myassert(x->sign == y->sign);
 495    if (x->sign == 0) return 0;
 496    if (x->sign == 1) return ucmp(x,y); else return ucmp(y,x);
 497 }
 498
 499
 500 /* --------------------------------------------------------------------------
 501  * Unsigned ops
 502  * ------------------------------------------------------------------------*/
 503
 504 static int ucmp ( B* x, B* y )
 505 {
 506    int i;
 507    myassert(sane(x));
 508    myassert(sane(y));
 509    if (x->used < y->used) return -1;
 510    if (x->used > y->used) return 1;
 511    for (i = x->used-1; i >= 0; i--) {
 512       if (x->stuff[i] < y->stuff[i]) return -1;
 513       if (x->stuff[i] > y->stuff[i]) return 1;
 514    }
 515    return 0;
 516 }
 517
 518
 519
 520 static void uadd ( B* x, B* y, B* res )
 521 {
 522    int c, i, t, n;
 523    B* longer;
 524
 525    myassert(sane(x));
 526    myassert(sane(y));
 527    myassert (res->size == maxused_add(x,y));
 528    res->used = res->size;
 529    res->stuff[res->used-1] = 0;
 530
 531    if (x->used > y->used) {
 532       n = y->used;
 533       longer = x;
 534    } else {
 535       n = x->used;
 536       longer = y;
 537    }
 538
 539    c = 0;
 540    for (i = 0; i < n; i++) {
 541       t = x->stuff[i] + y->stuff[i] + c;
 542       if (t >= B_BASE) {
 543          res->stuff[i] = t-B_BASE;
 544          c = 1;
 545       } else {
 546          res->stuff[i] = t;
 547          c = 0;
 548       }
 549    }
 550
 551    for (i = n; i < longer->used; i++) {
 552       t = longer->stuff[i] + c;
 553       if (t >= B_BASE) {
 554          res->stuff[i] = t-B_BASE;
 555       } else {
 556          res->stuff[i] = t;
 557          c = 0;
 558       }
 559    }
 560    if (c > 0) {
 561       myassert(res->used == longer->used+1);
 562       res->stuff[longer->used] = c;
 563    }
 564
 565    u_renormalise(res);
 566    myassert(sane(res));
 567 }
 568
 569
 570 static void usub ( B* x, B* y, B* res )
 571 {
 572    int b, i, t;
 573    myassert(sane(x));
 574    myassert(sane(y));
 575    myassert (x->used >= y->used);
 576    myassert (res->size == maxused_sub(x,y));
 577
 578    b = 0;
 579    for (i = 0; i < y->used; i++) {
 580       t = x->stuff[i] - y->stuff[i] - b;
 581       if (t < 0) {
 582          res->stuff[i] = t + B_BASE;
 583          b = 1;
 584       } else {
 585          res->stuff[i] = t;
 586          b = 0;
 587       }
 588    }
 589
 590    for (i = y->used; i < x->used; i++) {
 591       t = x->stuff[i] - b;
 592       if (t < 0) {
 593          res->stuff[i] = t + B_BASE;
 594       } else {
 595          res->stuff[i] = t;
 596          b = 0;
 597       }
 598    }
 599    myassert (b == 0);
 600
 601    res->used = x->used;
 602    u_renormalise(res);
 603    myassert(sane(res));
 604 }
 605
 606
 607 void umul ( B* x, B* y, B* res )
 608 {
 609    int i, j, carry;
 610
 611    myassert(sane(x));
 612    myassert(sane(y));
 613    myassert(res->size == maxused_mul(x,y));
 614
 615    for (j = 0; j < y->used; j++) res->stuff[j] = 0;
 616
 617    for (i = 0; i < x->used; i++) {
 618       carry = 0;
 619       for (j = 0; j < y->used; j++) {
 620          carry += res->stuff[i+j] + x->stuff[i]*y->stuff[j];
 621          res->stuff[i+j] = carry % B_BASE;
 622          carry /= B_BASE;
 623          myassert (carry < B_BASE);
 624       }
 625       res->stuff[i+y->used] = carry;
 626    }
 627
 628    res->used = x->used+y->used;
 629    u_renormalise(res);
 630    myassert(sane(res));
 631 }
 632
 633
 634 static void uqrm ( B* dend, B* isor, B* dres, B* mres )
 635 {
 636    int i, j, t, vh, toolarge, delta, carry, scaleup;
 637    uchar *dend_stuff, *isor_stuff, *tmp;
 638
 639    myassert(sane(isor));
 640    myassert(sane(dend));
 641    myassert(isor->used > 0);  // against division by zero
 642
 643    myassert(dres->size == maxused_qrm(isor,dend));
 644    myassert(mres->size == maxused_qrm(isor,dend));
 645
 646    if (dend->used < isor->used) {
 647       // Result of division must be zero, since dividend has
 648       // fewer digits than the divisor.  Remainder is the
 649       // original dividend.
 650       dres->used = 0;
 651       mres->used = dend->used;
 652       for (j = 0; j < mres->used; j++) mres->stuff[j] = dend->stuff[j];
 653       u_renormalise(dres); u_renormalise(mres);
 654       myassert(sane(dres));
 655       myassert(sane(mres));
 656       return;
 657    }
 658
 659    if (isor->used == 1) {
 660
 661       // Simple case; divisor is a single digit
 662       carry = 0;
 663       for (j = dend->used-1; j >= 0; j--) {
 664          carry += dend->stuff[j];
 665          dres->stuff[j] = carry/isor->stuff[0];
 666          carry = B_BASE*(carry%isor->stuff[0]);
 667       }
 668       carry /= B_BASE;
 669       dres->used = dend->used;
 670       u_renormalise(dres);
 671
 672       // Remainder is the final carry value
 673       mres->used = 0;
 674       if (carry > 0) {
 675          mres->used = 1;
 676          mres->stuff[0] = carry;
 677       }
 678       u_renormalise(dres); u_renormalise(mres);
 679       myassert(sane(dres));
 680       myassert(sane(mres));
 681       return;
 682
 683    } else {
 684
 685       // Complex case: both dividend and divisor have two or more digits.
 686       myassert(isor->used >= 2);
 687       myassert(dend->used >= 2);
 688
 689       // Allocate space for a copy of both dividend and divisor, since we
 690       // need to mess with them.  Also allocate tmp as a place to hold
 691       // values of the form   quotient_digit * divisor.
 692       dend_stuff = malloc ( sizeof(uchar)*(dend->used+1) );
 693       isor_stuff = malloc ( sizeof(uchar)*isor->used     );
 694       tmp        = malloc ( sizeof(uchar)*(isor->used+1) );
 695       myassert (dend_stuff && isor_stuff && tmp);
 696
 697       // Calculate a scaling-up factor, and multiply both divisor and
 698       // dividend by it.  Doing this reduces the number of corrections
 699       // needed to the quotient-digit-estimates made in the loop below,
 700       // and thus speeds up division, but is not actually needed to
 701       // get the correct results.  The scaleup factor should not increase
 702       // the number of digits needed to represent either the divisor
 703       // (since the factor is derived from it) or the dividend (since
 704       // we already gave it a new leading zero).
 705       scaleup = B_BASE / (1 + isor->stuff[isor->used-1]);
 706       myassert (1 <= scaleup && scaleup <= B_BASE/2);
 707
 708       if (scaleup == 1) {
 709          // Don't bother to multiply; just copy.
 710          for (j = 0; j < dend->used; j++) dend_stuff[j] = dend->stuff[j];
 711          for (j = 0; j < isor->used; j++) isor_stuff[j] = isor->stuff[j];
 712
 713          // Extend dividend with leading zero.
 714          dend_stuff[dend->used] = tmp[isor->used] = 0;
 715
 716       } else {
 717          carry = 0;
 718          for (j = 0; j < isor->used; j++) {
 719             t = scaleup * isor->stuff[j] + carry;
 720             isor_stuff[j] = t % B_BASE;
 721             carry = t / B_BASE;
 722          }
 723          myassert (carry == 0);
 724
 725          carry = 0;
 726          for (j = 0; j < dend->used; j++) {
 727             t = scaleup * dend->stuff[j] + carry;
 728             dend_stuff[j] = t % B_BASE;
 729             carry = t / B_BASE;
 730          }
 731          dend_stuff[dend->used] = carry;
 732          tmp[isor->used] = 0;
 733       }
 734
 735       // For each quotient digit ...
 736       for (i = dend->used; i >= isor->used; i--) {
 737          myassert (i-2 >= 0);
 738          myassert (i <= dend->used);
 739          myassert (isor->used >= 2);
 740
 741 #if DEBUG_SAINTEGER_UQRM
 742          printf("\n---------\nqdigit %d\n", i );
 743          printf("dend_stuff is ");
 744          for (j = dend->used; j>= 0; j--) printf("%d ",dend_stuff[j]);
 745          printf("\n");
 746 #endif
 747          // Make a guess vh of the quotient digit
 748          vh = (B_BASE*B_BASE*dend_stuff[i] + B_BASE*dend_stuff[i-1] + dend_stuff[i-2])
 749               /
 750               (B_BASE*isor_stuff[isor->used-1] + isor_stuff[isor->used-2]);
 751          if (vh > B_BASE-1) vh = B_BASE-1;
 752 #if DEBUG_SAINTEGER_UQRM
 753          printf("guess formed from %d %d %d   %d %d\n",
 754                  dend_stuff[i], dend_stuff[i-1] , dend_stuff[i-2],
 755                  isor_stuff[isor->used-1], isor_stuff[isor->used-2]);
 756          printf("guess is %d\n", vh );
 757 #endif
 758          // Check if vh is too large (by 1).  Calculate vh * isor into tmp
 759          // and see if it exceeds the same length prefix of dend.  If so,
 760          // vh needs to be decremented.
 761          carry = 0;
 762          for (j = 0; j < isor->used; j++) {
 763             t = vh * isor_stuff[j] + carry;
 764             tmp[j] = t % B_BASE;
 765             carry = t / B_BASE;
 766          }
 767          tmp[isor->used] = carry;
 768          delta = i - isor->used;
 769 #if DEBUG_SAINTEGER_UQRM
 770          printf("final carry is %d\n", carry);
 771          printf("vh * isor is " );
 772          for (j = isor->used; j >=0; j--) printf("%d ",tmp[j]);printf("\n");
 773          printf("delta = %d\n", delta );
 774 #endif
 775          toolarge = 0;
 776          for (j = isor->used; j >= 0; j--) {
 777 #if DEBUG_SAINTEGER_UQRM
 778             printf ( "(%d,%d)  ", (int)(tmp[j]), (int)(dend_stuff[j+delta]) );
 779 #endif
 780             if (tmp[j] > dend_stuff[j+delta]) {toolarge=1; break;};
 781             if (tmp[j] < dend_stuff[j+delta]) break;
 782          }
 783
 784          // If we did guess too large, decrement vh and subtract a copy of
 785          // isor from tmp.  This had better not go negative!
 786          if (toolarge) {
 787 #if DEBUG_SAINTEGER_UQRM
 788             printf ( "guess too large\n" );
 789 #endif
 790             vh--;
 791             carry = 0;
 792             for (j = 0; j < isor->used; j++) {
 793                if (carry + isor_stuff[j] > tmp[j]) {
 794                   tmp[j] = (B_BASE + tmp[j]) - isor_stuff[j] - carry;
 795                   carry = 1;
 796                } else {
 797                   tmp[j] = tmp[j] - isor_stuff[j] - carry;
 798                   carry = 0;
 799                }
 800             }
 801             //if (carry > 0) {pp(isor);pp(dend);};
 802             //myassert(carry == 0);
 803             if (carry > 0) {
 804                myassert(tmp[isor->used] > 0);
 805                tmp[isor->used]--;
 806             }
 807 #if DEBUG_SAINTEGER_UQRM
 808             printf("after adjustment of tmp ");
 809             for (j = isor->used; j >=0; j--) printf("%d ",tmp[j]);
 810             printf("\n");
 811 #endif
 812          }
 813
 814          // Now vh really is the i'th quotient digit.
 815          // Subtract (tmp << delta) from
 816          // the dividend.
 817          carry = 0;
 818          for (j = 0; j <= isor->used; j++) {
 819             if (carry + tmp[j] > dend_stuff[j+delta]) {
 820                dend_stuff[j+delta] = (B_BASE+dend_stuff[j+delta]) - tmp[j] - carry;
 821                carry = 1;
 822             } else {
 823                dend_stuff[j+delta] = dend_stuff[j+delta] - tmp[j] - carry;
 824                carry = 0;
 825             }
 826          }
 827          myassert(carry==0);
 828
 829 #if DEBUG_SAINTEGER_UQRM
 830          printf("after final sub ");
 831          for(j=dend->used; j>=0; j--) printf("%d ", dend_stuff[j]);
 832          printf("\n");
 833 #endif
 834
 835          // park vh in the result array
 836 #if DEBUG_SAINTEGER_UDIV
 837          printf("[%d] <- %d\n", i-isor->used, vh );
 838 #endif
 839          dres->stuff[i-isor->used] = vh;
 840       }
 841    }
 842
 843    // Now we've got all the quotient digits.  Zap leading zeroes.
 844    dres->used = dend->used - isor->used + 1;
 845    u_renormalise(dres);
 846    myassert(sane(dres));
 847
 848    // The remainder is in dend_stuff.  Copy, divide by the original scaling
 849    // factor, and zap leading zeroes.
 850    mres->used = dend->used;
 851    for (j = 0; j < dend->used; j++) mres->stuff[j] = dend_stuff[j];
 852    u_renormalise(mres);
 853    myassert(sane(mres));
 854
 855    if (scaleup > 1) {
 856       carry = 0;
 857       for (j = mres->used-1; j >= 0; j--) {
 858          carry += mres->stuff[j];
 859          mres->stuff[j] = carry/scaleup;
 860          carry = B_BASE*(carry%scaleup);
 861       }
 862       myassert (carry == 0);
 863       u_renormalise(mres);
 864       myassert(sane(mres));
 865    }
 866
 867    free(tmp);
 868    free(isor_stuff);
 869    free(dend_stuff);
 870 }
 871
 872
 873 /* --------------------------------------------------------------------------
 874  * Test framework
 875  * ------------------------------------------------------------------------*/
 876
 877 #if 0
 878 int main ( int argc, char** argv )
 879 {
 880    int i, j, t, k, m;
 881    B *bi, *bj, *bk, *bm;
 882
 883    for (i = -10007; i <= 10007; i++) {
 884       printf ( "i = %d\n", i );
 885
 886       t = size_fromInt(); bi = malloc(t); myassert(bi);
 887       do_fromInt(i, t, bi);
 888
 889       t = do_toInt(bi); myassert(i == t);
 890
 891       for (j = -10007; j <= 10007; j++) {
 892
 893          t = size_fromInt(); bj = malloc(t); myassert(bj);
 894          do_fromInt(j, t, bj);
 895
 896          t = do_toInt(bj); myassert(j == t);
 897
 898          if (1) {
 899             t = size_add(bi,bj); bk = malloc(t); myassert(bk);
 900             do_add(bi,bj,t,bk);
 901             k = do_toInt(bk);
 902             if (i+j != k) {
 903                pp(bi); pp(bj); pp(bk);
 904                myassert(i+j == k);
 905             }
 906             free(bk);
 907          }
 908
 909          if (1) {
 910             t = size_sub(bi,bj); bk = malloc(t); myassert(bk);
 911             do_sub(bi,bj,t,bk);
 912             k = do_toInt(bk);
 913             if (i-j != k) {
 914                pp(bi); pp(bj); pp(bk);
 915                myassert(i-j == k);
 916             }
 917             free(bk);
 918          }
 919
 920          if (1) {
 921             t = size_mul(bi,bj); bk = malloc(t); myassert(bk);
 922             do_mul(bi,bj,t,bk);
 923             k = do_toInt(bk);
 924             if (i*j != k) {
 925                pp(bi); pp(bj); pp(bk);
 926                myassert(i*j == k);
 927             }
 928             free(bk);
 929          }
 930
 931          if (j != 0) {
 932             t = size_qrm(bi,bj);
 933             bk = malloc(t); myassert(bk);
 934             bm = malloc(t); myassert(bm);
 935             do_qrm(bi,bj,t,bk,bm);
 936             k = do_toInt(bk);
 937             m = do_toInt(bm);
 938             myassert(k == i/j);
 939             myassert(m == i%j);
 940             free(bk); free(bm);
 941          }
 942
 943          free(bj);
 944       }
 945       free(bi);
 946
 947    }
 948    printf("done\n");
 949    return 0;
 950 }
 951 #endif
 952
 953 #if 0
 954 int main ( int argc, char** argv )
 955 {
 956    B *a, *b, *c, *d, *e;
 957    a = fromInt(1); b=fromInt(9); pp(a); pp(b);
 958    c = mkB( maxused_uqrm(a,b) );
 959    d = mkB( maxused_uqrm(a,b) );
 960    e = mkB( maxused_uadd(a,b) );
 961    uadd(a,b,e); pp(e);
 962    //uqrm(a,b,c,d); pp(c); pp(d);
 963
 964    return 0;
 965 }
 966 #endif
 967
 968 /*-------------------------------------------------------------------------*/