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