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