home *** CD-ROM | disk | FTP | other *** search
/ minnie.tuhs.org / unixen.tar / unixen / PDP-11 / Trees / V7 / usr / src / libmp / mdiv.c < prev    next >
Encoding:
C/C++ Source or Header  |  1979-01-10  |  1.9 KB  |  107 lines
  1. #include <mp.h>
  2. mdiv(a,b,q,r) MINT *a,*b,*q,*r;
  3. {    MINT x,y;
  4.     int sign;
  5.     sign=1;
  6.     x.val=a->val;
  7.     y.val=b->val;
  8.     x.len=a->len;
  9.     if(x.len<0) {sign= -1; x.len= -x.len;}
  10.     y.len=b->len;
  11.     if(y.len<0) {sign= -sign; y.len= -y.len;}
  12.     xfree(q);
  13.     xfree(r);
  14.     m_div(&x,&y,q,r);
  15.     if(sign==-1)
  16.     {    q->len= -q->len;
  17.         r->len = - r->len;
  18.     }
  19.     return;
  20. }
  21. m_dsb(q,n,a,b) short *a,*b;
  22. {    long int x,qx;
  23.     int borrow,j,u;
  24.     qx=q;
  25.     borrow=0;
  26.     for(j=0;j<n;j++)
  27.     {    x=borrow-a[j]*qx+b[j];
  28.         b[j]=x&077777;
  29.         borrow=x>>15;
  30.     }
  31.     x=borrow+b[j];
  32.     b[j]=x&077777;
  33.     if(x>>15 ==0) { return(0);}
  34.     borrow=0;
  35.     for(j=0;j<n;j++)
  36.     {    u=a[j]+b[j]+borrow;
  37.         if(u<0) borrow=1;
  38.         else borrow=0;
  39.         b[j]=u&077777;
  40.     }
  41.     { return(1);}
  42. }
  43. m_trq(v1,v2,u1,u2,u3)
  44. {    long int d;
  45.     long int x1;
  46.     if(u1==v1) d=077777;
  47.     else d=(u1*0100000L+u2)/v1;
  48.     while(1)
  49.     {    x1=u1*0100000L+u2-v1*d;
  50.         x1=x1*0100000L+u3-v2*d;
  51.         if(x1<0) d=d-1;
  52.         else {return(d);}
  53.     }
  54. }
  55. m_div(a,b,q,r) MINT *a,*b,*q,*r;
  56. {    MINT u,v,x,w;
  57.     short d,*qval;
  58.     int qq,n,v1,v2,j;
  59.     u.len=v.len=x.len=w.len=0;
  60.     if(b->len==0) { fatal("mdiv divide by zero"); return;}
  61.     if(b->len==1)
  62.     {    r->val=xalloc(1,"m_div1");
  63.         sdiv(a,b->val[0],q,r->val);
  64.         if(r->val[0]==0)
  65.         {    shfree(r->val);
  66.             r->len=0;
  67.         }
  68.         else r->len=1;
  69.         return;
  70.     }
  71.     if(a->len < b->len)
  72.     {    q->len=0;
  73.         r->len=a->len;
  74.         r->val=xalloc(r->len,"m_div2");
  75.         for(qq=0;qq<r->len;qq++) r->val[qq]=a->val[qq];
  76.         return;
  77.     }
  78.     x.len=1;
  79.     x.val = &d;
  80.     n=b->len;
  81.     d=0100000L/(b->val[n-1]+1L);
  82.     mult(a,&x,&u); /*subtle: relies on fact that mult allocates extra space */
  83.     mult(b,&x,&v);
  84.     v1=v.val[n-1];
  85.     v2=v.val[n-2];
  86.     qval=xalloc(a->len-n+1,"m_div3");
  87.     for(j=a->len-n;j>=0;j--)
  88.     {    qq=m_trq(v1,v2,u.val[j+n],u.val[j+n-1],u.val[j+n-2]);
  89.         if(m_dsb(qq,n,v.val,&(u.val[j]))) qq -= 1;
  90.         qval[j]=qq;
  91.     }
  92.     x.len=n;
  93.     x.val=u.val;
  94.     mcan(&x);
  95.     sdiv(&x,d,&w,(short *)&qq);
  96.     r->len=w.len;
  97.     r->val=w.val;
  98.     q->val=qval;
  99.     qq=a->len-n+1;
  100.     if(qq>0 && qval[qq-1]==0) qq -= 1;
  101.     q->len=qq;
  102.     if(qq==0) shfree(qval);
  103.     if(x.len!=0) xfree(&u);
  104.     xfree(&v);
  105.     return;
  106. }
  107.