2.9BSD - /usr/src/lib/ape/mdiv.c

   1: #include "lint.h"
   2: #ifndef lint
   3: static char sccs_id[] = "@(#)mdiv.c	2.1	7/6/82";
   4: #endif lint
   5: 
   6: #include <ape.h>
   7: 
   8: sdiv(a,n,q,r)  /* q is quotient, r remainder of a/n for n int */
   9: MINT *a,*q; short *r;
  10: {   MINT x,y;
  11:     int sign;
  12: 
  13:     if (n==0) aperror("sdiv: zero divisor");
  14:     sign=1;
  15:     x.len=a->len;
  16:     x.val=a->val;
  17:     if (n<0)
  18:     {   sign= -sign;
  19:         n= -n;
  20:     }
  21:     if (x.len<0)
  22:     {   sign = -sign;
  23:         x.len= -x.len;
  24:     }
  25:     s_div(&x,n,&y,r);
  26:     xfree(q);
  27:     q->val=y.val;
  28:     q->len=sign*y.len;
  29:     *r = *r*sign;
  30:     return;
  31: }
  32: 
  33: s_div(a,n,q,r)  /* dirtywork function for sdiv */
  34: MINT *a,*q; short *r;
  35: {   int qlen,i;
  36:     long int x;
  37:     short *qval;
  38: 
  39:     if (n==0) aperror("sdiv: zero divisor");
  40:     x=0L;
  41:     qlen=a->len;
  42:     qval=xalloc(qlen,"s_div");  /* I guess qlen and qval are only used
  43: 				  	 * to cut down on the amount of typing
  44: 				  	 * required. */
  45:     for(i=qlen-1;i>=0;i--)
  46:     {
  47:         x=x*LONGCARRY+a->val[i];
  48:         qval[i]=x/n;
  49:         x=x%n;
  50:     }
  51:     *r=x;
  52:     if (qval[qlen-1]==0) qlen--;
  53:     q->len=qlen;
  54:     q->val=qval;
  55:     if (qlen==0) shfree(qval);
  56:     return;
  57: }
  58: 
  59: mdiv(a,b,q,r)  /* q = quotient, r = remainder, of a/b */
  60: MINT *a,*b,*q,*r;
  61: {   MINT x,y,w,z;
  62:     int sign;
  63: 
  64:     if (b->len == 0) aperror("mdiv: zero divisor");
  65:     sign=1;
  66:     x.val=a->val;
  67:     y.val=b->val;
  68:     x.len=a->len;
  69:     if(x.len<0) {sign= -1; x.len= -x.len;}
  70:     y.len=b->len;
  71:     if(y.len<0) {sign= -sign; y.len= -y.len;}
  72:     z.len = w.len = 0;
  73:     m_div(&x,&y,&z,&w);
  74:     xfree(q);
  75:     xfree(r);
  76:     q->val = z.val; q->len = z.len;
  77:     r->val = w.val; r->len = w.len;
  78:     if(sign==-1)
  79:     {   q->len = -q->len;
  80:         r->len = -r->len;
  81:     }
  82:     return;
  83: }
  84: 
  85: m_dsb(q,n,a,b)
  86: short *a,*b;
  87: {   long int x,qx;
  88:     int borrow,j,u;
  89: 
  90:     qx=q;
  91:     /* qx is used merely to force arithmetic
  92: 	 * to be done long (double integral precision) */
  93:     borrow=0;
  94:     for(j=0;j<n;j++)
  95:     {   x=borrow-a[j]*qx+b[j];
  96:         b[j]=x&TOPSHORT;
  97:         borrow=x>>15;
  98:     }
  99:     x=borrow+b[j];
 100:     b[j]=x&TOPSHORT;
 101:     if (x>>15 ==0)  return(0);
 102:     borrow=0;
 103:     for(j=0;j<n;j++)
 104:     {   u=a[j]+b[j]+borrow;
 105:         if(u<0) borrow=1;
 106:         else borrow=0;
 107:         b[j]=u&TOPSHORT;
 108:     }
 109:     return(1);
 110: }
 111: 
 112: m_trq(v1,v2,u1,u2,u3)
 113: {   long int d;
 114:     long int x1;
 115: 
 116:     if(u1==v1) d=TOPSHORT;
 117:     else d=(u1*LONGCARRY+u2)/v1;
 118:     while(1)
 119:     {   x1=u1*LONGCARRY+u2-v1*d;
 120:         x1=x1*LONGCARRY+u3-v2*d;
 121:         if(x1<0) d=d-1;
 122:         else return(d);
 123:     }
 124: }
 125: 
 126: m_div(a,b,q,r)  /* basic "long division" routine */
 127: MINT *a,*b,*q,*r;
 128: {   MINT u,v,x,w;
 129:     short d,*qval;
 130:     int qq,n,v1,v2,j;
 131: 
 132:     u.len=v.len=x.len=w.len=0;
 133:     if (b->len==0) { aperror("mdiv divide by zero"); return;}
 134:     if (b->len==1)
 135:     {   r->val=xalloc(1,"m_div1");
 136:         sdiv(a,b->val[0],q,r->val);
 137:         if (r->val[0]==0)
 138:             {
 139:             xfree(r);
 140:             return;
 141:             }
 142:         else r->len=1;
 143:         return;
 144:     }
 145:     if (mcmp(a,b) <0)
 146:     {   xfree(q);
 147:         move(a,r);
 148:         return;
 149:     }
 150:     x.len=1;
 151:     x.val = &d;
 152:     n=b->len;
 153:     d=LONGCARRY/(b->val[n-1]+1L);
 154:     mult(a,&x,&u); /*subtle: relies on fact that mult allocates extra space */
 155:     mult(b,&x,&v);
 156:     v1=v.val[n-1];
 157:     v2=v.val[n-2];
 158:     qval=xalloc(a->len-n+1,"m_div3");
 159:     for(j=a->len-n;j>=0;j--)
 160:     {   qq=m_trq(v1,v2,u.val[j+n],u.val[j+n-1],u.val[j+n-2]);
 161:         if(m_dsb(qq,n,v.val,&(u.val[j]))) qq -= 1;
 162:         qval[j]=qq;
 163:     }
 164:     x.len=n;
 165:     x.val=u.val;
 166:     mcan(&x);
 167:     sdiv(&x,d,&w,(short *)&qq);
 168:     r->len=w.len;
 169:     r->val=w.val;
 170:     q->val=qval;
 171:     qq=a->len-n+1;
 172:     if(qq>0 && qval[qq-1]==0) qq -= 1;
 173:     q->len=qq;
 174:     if(qq==0) shfree(qval);
 175:     if(x.len!=0) xfree(&u);
 176:     xfree(&v);
 177:     return;
 178: }

Defined functions

Defined variables