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