| 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 | } |