tmpdiv.c - plan9port - [fork] Plan 9 from user space
  HTML git clone git://src.adamsgaard.dk/plan9port
   DIR Log
   DIR Files
   DIR Refs
   DIR README
   DIR LICENSE
       ---
       tmpdiv.c (2480B)
       ---
            1 #include "os.h"
            2 #include <mp.h>
            3 #include "dat.h"
            4 
            5 /* division ala knuth, seminumerical algorithms, pp 237-238 */
            6 /* the numbers are stored backwards to what knuth expects so j */
            7 /* counts down rather than up. */
            8 
            9 void
           10 mpdiv(mpint *dividend, mpint *divisor, mpint *quotient, mpint *remainder)
           11 {
           12         int j, s, vn, sign;
           13         mpdigit qd, *up, *vp, *qp;
           14         mpint *u, *v, *t;
           15 
           16         /* divide bv zero */
           17         if(divisor->top == 0)
           18                 abort();
           19 
           20         /* quick check */
           21         if(mpmagcmp(dividend, divisor) < 0){
           22                 if(remainder != nil)
           23                         mpassign(dividend, remainder);
           24                 if(quotient != nil)
           25                         mpassign(mpzero, quotient);
           26                 return;
           27         }
           28 
           29         /* D1: shift until divisor, v, has hi bit set (needed to make trial */
           30         /*     divisor accurate) */
           31         qd = divisor->p[divisor->top-1];
           32         for(s = 0; (qd & mpdighi) == 0; s++)
           33                 qd <<= 1;
           34         u = mpnew((dividend->top+2)*Dbits + s);
           35         if(s == 0 && divisor != quotient && divisor != remainder) {
           36                 mpassign(dividend, u);
           37                 v = divisor;
           38         } else {
           39                 mpleft(dividend, s, u);
           40                 v = mpnew(divisor->top*Dbits);
           41                 mpleft(divisor, s, v);
           42         }
           43         up = u->p+u->top-1;
           44         vp = v->p+v->top-1;
           45         vn = v->top;
           46 
           47         /* D1a: make sure high digit of dividend is less than high digit of divisor */
           48         if(*up >= *vp){
           49                 *++up = 0;
           50                 u->top++;
           51         }
           52 
           53         /* storage for multiplies */
           54         t = mpnew(4*Dbits);
           55 
           56         qp = nil;
           57         if(quotient != nil){
           58                 mpbits(quotient, (u->top - v->top)*Dbits);
           59                 quotient->top = u->top - v->top;
           60                 qp = quotient->p+quotient->top-1;
           61         }
           62 
           63         /* D2, D7: loop on length of dividend */
           64         for(j = u->top; j > vn; j--){
           65 
           66                 /* D3: calculate trial divisor */
           67                 mpdigdiv(up-1, *vp, &qd);
           68 
           69                 /* D3a: rule out trial divisors 2 greater than real divisor */
           70                 if(vn > 1) for(;;){
           71                         memset(t->p, 0, 3*Dbytes);        /* mpvecdigmuladd adds to what's there */
           72                         mpvecdigmuladd(vp-1, 2, qd, t->p);
           73                         if(mpveccmp(t->p, 3, up-2, 3) > 0)
           74                                 qd--;
           75                         else
           76                                 break;
           77                 }
           78 
           79                 /* D4: u -= v*qd << j*Dbits */
           80                 sign = mpvecdigmulsub(v->p, vn, qd, up-vn);
           81                 if(sign < 0){
           82 
           83                         /* D6: trial divisor was too high, add back borrowed */
           84                         /*     value and decrease divisor */
           85                         mpvecadd(up-vn, vn+1, v->p, vn, up-vn);
           86                         qd--;
           87                 }
           88 
           89                 /* D5: save quotient digit */
           90                 if(qp != nil)
           91                         *qp-- = qd;
           92 
           93                 /* push top of u down one */
           94                 u->top--;
           95                 *up-- = 0;
           96         }
           97         if(qp != nil){
           98                 mpnorm(quotient);
           99                 if(dividend->sign != divisor->sign)
          100                         quotient->sign = -1;
          101         }
          102 
          103         if(remainder != nil){
          104                 mpright(u, s, remainder);        /* u is the remainder shifted */
          105                 remainder->sign = dividend->sign;
          106         }
          107 
          108         mpfree(t);
          109         mpfree(u);
          110         if(v != divisor)
          111                 mpfree(v);
          112 }