|
|
1.1 ! root 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: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.