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