fix from the net for sign extension problem
[unix-history] / usr / src / lib / libmp / mdiv.c
CommitLineData
a4b0a236 1/* @(#)mdiv.c 4.2 %G% */
6a73aa09
SL
2
3#include <mp.h>
4mdiv(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}
23m_dsb(q,n,a,b) short *a,*b;
24{ long int x,qx;
a4b0a236
SL
25 int borrow,j;
26 short u;
6a73aa09
SL
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}
46m_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}
58m_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}