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