# 2.11BSD#431 - /usr/src/new/PORT/B/src/bint/b1nuQ.c

##### Functions

    int1div           X
    int1mul           X
    int_ldiv          X
    int_mod           X
    int_quot          X

       1: /* Copyright (c) Stichting Mathematisch Centrum, Amsterdam, 1985. */
       2:
       3: /*
       4:   $Header: b1nuQ.c,v 1.4 85/08/22 16:51:40 timo Exp $
       5: */
       6:
       7: #include "b.h"
       8: #include "b1obj.h"
       9: #include "b0con.h"
      10: #include "b1num.h"
      11:
      12:
      13: /* Product of integer and single "digit" */
      14:
      15: Visible integer int1mul(v, n1) integer v; digit n1; {
      16:     integer a;
      17:     digit save, bigcarry, carry = 0;
      18:     double z, zz, n = n1;
      19:     register int i;
      20:     struct integer vv;
      21:
      22:     FreezeSmallInt(v, vv);
      23:
      24:     a = (integer) grab_num(Length(v)+2);
      25:
      26:     for (i = 0; i < Length(v); ++i) {
      27:         z = Digit(v,i) * n;
      28:         bigcarry = zz = floor(z/BASE);
      29:         carry += z - zz*BASE;
      30:         Digit(a,i) = save = Modulo(carry, BASE);
      31:         carry = (carry-save)/BASE + bigcarry;
      32:     }
      33:
      34:     Digit(a,i) = save = Modulo(carry, BASE);
      35:     Digit(a,i+1) = (carry-save)/BASE;
      36:
      37:     return int_canon(a);
      38: }
      39:
      40:
      41: /* Quotient of positive integer and single "digit" > 0 */
      42:
      43: Hidden integer int1div(v, n1, prem) integer v; digit n1, *prem; {
      44:     integer q;
      45:     double r_over_n, r = 0, n = n1;
      46:     register int i;
      47:     struct integer vv;
      48:
      49:     FreezeSmallInt(v, vv);
      50:
      51:     q = (integer) grab_num(Length(v));
      52:     for (i = Length(v)-1; i >= 0; --i) {
      53:         r = r*BASE + Digit(v,i);
      54:         Digit(q,i) = r_over_n = floor(r/n);
      55:         r -= r_over_n * n;
      56:     }
      57:     if (prem)
      58:         *prem = r;
      59:     return int_canon(q);
      60: }
      61:
      62:
      63: /* Long division routine, gives access to division algorithm. */
      64:
      65: Visible digit int_ldiv(v1, w1, pquot, prem) integer v1, w1, *pquot, *prem; {
      66:     integer a;
      67:     int sign = 1, rel_v = 0, rel_w = 0;
      68:     digit div, rem;
      69:     struct integer vv1, ww1;
      70:
      71:     if (w1 == int_0) syserr(MESS(1100, "zero division (int_ldiv)"));
      72:
      73:     /* Make v, w positive */
      74:     if (Msd(v1) < 0) {
      75:         sign = -1;
      76:         ++rel_v;
      77:         v1 = int_neg(v1);
      78:     }
      79:
      80:     if (Msd(w1) < 0) {
      81:         sign *= -1;
      82:         ++rel_w;
      83:         w1 = int_neg(w1);
      84:     }
      85:
      86:     FreezeSmallInt(v1, vv1);
      87:     FreezeSmallInt(w1, ww1);
      88:
      89:     div = sign;
      90:
      91:     /* Check v << w or single-digit w */
      92:     if (Length(v1) < Length(w1)
      93:         || Length(v1) == Length(w1)
      94:             && Digit(v1, Length(v1)-1) < Digit(w1, Length(w1)-1)) {
      95:         a = int_0;
      96:         if (prem) {
      97:             if (v1 == &vv1) *prem= (integer) MkSmallInt(Digit(v1,0));
      98:             else *prem = (integer) Copy(v1);
      99:         }
     100:     }
     101:     else if (Length(w1) == 1) {
     102:         /* Single-precision division */
     103:         a = int1div(v1, Digit(w1,0), &rem);
     104:         if (prem) *prem = mk_int((double)rem);
     105:     }
     106:     else {
     107:         /* Multi-precision division */
     108:         /* Cf. Knuth II Sec. 4.3.1. Algorithm D */
     109:         /* Note that we count in the reverse direction (not easier!) */
     110:
     111:         double z, zz;
     112:         digit carry, save, bigcarry;
     113:         double q, d = BASE/(Digit(w1, Length(w1)-1)+1);
     114:         register int i, j, k;
     115:         integer v, w;
     116:         digit vj;
     117:
     118:         /* Normalize: make Msd(w) >= BASE/2 by multiplying
     119: 		   both v and w by d */
     120:
     121:         v = int1mul(v1, (digit)d);
     122:             /* v is used as accumulator, must make a copy */
     123:             /* v cannot be int_1 */
     124:             /* (then it would be one of the cases above) */
     125:
     126:         if (d == 1) w = (integer) Copy(w1);
     127:         else w = int1mul(w1, (digit)d);
     128:
     129:         a = (integer) grab_num(Length(v1)-Length(w)+1);
     130:
     131:         /* Division loop */
     132:
     133:         for (j = Length(v1), k = Length(a)-1; k >= 0; --j, --k) {
     134:             vj = j >= Length(v) ? 0 : Digit(v,j);
     135:
     136:             /* Find trial digit */
     137:
     138:             if (vj == Digit(w, Length(w)-1)) q = BASE-1;
     139:             else q = floor( ((double)vj*BASE
     140:                     + Digit(v,j-1)) / Digit(w, Length(w)-1) );
     141:
     142:             /* Correct trial digit */
     143:
     144:             while (Digit(w,Length(w)-2) * q >
     145:                 ((double)vj*BASE + Digit(v,j-1)
     146:                     - q*Digit(w, Length(w)-1)) *BASE + Digit(v,j-2))
     147:                 --q;
     148:
     149:             /* Subtract q*w from v */
     150:
     151:             carry = 0;
     152:             for (i = 0; i < Length(w) && i+k < Length(v); ++i) {
     153:                 z = Digit(w,i) * q;
     154:                 bigcarry = zz = floor(z/BASE);
     155:                 carry += Digit(v,i+k) - z + zz*BASE;
     156:                 Digit(v,i+k) =
     157:                     save = Modulo(carry, BASE);
     158:                 carry = (carry-save)/BASE - bigcarry;
     159:             }
     160:
     161:             if (i+k < Length(v))
     162:                 carry += Digit(v, i+k), Digit(v, i+k) = 0;
     163:
     164:             /* Add back necessary? */
     165:
     166:                 /* It is very difficult to find test cases
     167: 				   where add back is necessary if BASE is large.
     168: 				   Thanks to Arjen Lenstra, we have v=n*n-1, w=n,
     169: 				   where n = 8109636009903000000 (the last six
     170: 				   digits are not important). */
     171:
     172:             if (carry == 0)     /* No */
     173:                 Digit(a,k) = q;
     174:             else {      /* Yes, add back */
     175:                 if (carry != -1) syserr(MESS(1101, "int_ldiv internal failure"));
     176:                 Digit(a,k) = q-1;
     177:                 carry = 0;
     178:                 for (i = 0; i < Length(w) && i+k < Length(v); ++i) {
     179:                     carry += Digit(v, i+k) + Digit(w,i);
     180:                     Digit(v,i+k) =
     181:                         save = Modulo(carry, BASE);
     182:                     carry = (carry-save)/BASE;
     183:                 }
     184:             }
     185:         }   /* End for(j) */
     186:
     187:         if (prem) *prem = int_canon(v); /* Store remainder */
     188:         else release((value) v);
     189:         div = sign*d;   /* Store normalization factor */
     190:         release((value) w);
     191:         a = int_canon(a);
     192:     }
     193:
     194:     if (rel_v) release((value) v1);
     195:     if (rel_w) release((value) w1);
     196:
     197:     if (sign < 0) {
     198:         integer temp = a;
     199:         a = int_neg(a);
     200:         release((value) temp);
     201:     }
     202:
     203:     if (pquot) *pquot = a;
     204:     else release((value) a);
     205:     return div;
     206: }
     207:
     208:
     209: Visible integer int_quot(v, w) integer v, w; {
     210:     integer quo;
     211:     VOID int_ldiv(v, w, &quo, (integer*)0);
     212:     return quo;
     213: }
     214:
     215: Visible integer int_mod(v, w) integer v, w; {
     216:     integer rem;
     217:     digit div;
     218:     bool flag;
     219:     div = int_ldiv(v, w, (integer*)0, &rem); /* Rem. is always positive */
     220:     if (rem == int_0)
     221:         return rem; /* v mod w = 0 */
     222:     flag = (div < 0);
     223:     if (flag || Msd(w) < 0) div = -div;
     224:     if (div != 1) { /* Divide by div to get proper remainder back */
     225:         v = int1div(rem, div, (digit*)0);
     226:         release((value) rem);
     227:         rem = v;
     228:     }
     229:     if (flag) { /* Make same sign as w */
     230:         if (Msd(w) < 0) v = int_sum(w, rem);
     231:         else v = int_diff(w, rem);
     232:         release((value) rem);
     233:         rem = v;
     234:     }
     235:     return rem;
     236: }

#### Defined functions

int1div defined in line 43; used 2 times

- in line 103, 225

int1mul defined in line 15; used 3 times

- in line 121, 127
- in /usr/src/new/PORT/B/src/bint/b1num.h line 50

int_ldiv defined in line 65; used 4 times

- in line 211, 219
- in /usr/src/new/PORT/B/src/bint/b1fun.c line 288
- in /usr/src/new/PORT/B/src/bint/b1num.h line 54

int_mod defined in line 215; used 2 times

- in /usr/src/new/PORT/B/src/bint/b1fun.c line 179
- in /usr/src/new/PORT/B/src/bint/b1num.h line 53

int_quot defined in line 209; used 3 times

- in /usr/src/new/PORT/B/src/bint/b1nuR.c line 68-69(2)
- in /usr/src/new/PORT/B/src/bint/b1num.h line 47

Last modified: 1985-08-27 Generated: 2016-12-26 | Generated by src2html V0.67 page hit count: 731 |