148370Sbostic /*-
248370Sbostic * %sccs.include.proprietary.c%
348370Sbostic */
448370Sbostic
519819Sdist #ifndef lint
6*61321Sbostic static char sccsid[] = "@(#)mdiv.c 8.1 (Berkeley) 06/04/93";
748370Sbostic #endif /* not lint */
819819Sdist
99946Ssam #include <mp.h>
mdiv(a,b,q,r)109946Ssam mdiv(a,b,q,r) MINT *a,*b,*q,*r;
119946Ssam { MINT x,y;
129946Ssam int sign;
139946Ssam sign=1;
149946Ssam x.val=a->val;
159946Ssam y.val=b->val;
169946Ssam x.len=a->len;
179946Ssam if(x.len<0) {sign= -1; x.len= -x.len;}
189946Ssam y.len=b->len;
199946Ssam if(y.len<0) {sign= -sign; y.len= -y.len;}
209946Ssam xfree(q);
219946Ssam xfree(r);
229946Ssam m_div(&x,&y,q,r);
239946Ssam if(sign==-1)
249946Ssam { q->len= -q->len;
259946Ssam r->len = - r->len;
269946Ssam }
279946Ssam return;
289946Ssam }
m_dsb(q,n,a,b)299946Ssam m_dsb(q,n,a,b) short *a,*b;
309946Ssam { long int x,qx;
3114332Ssam int borrow,j;
3214332Ssam short u;
339946Ssam qx=q;
349946Ssam borrow=0;
359946Ssam for(j=0;j<n;j++)
369946Ssam { x=borrow-a[j]*qx+b[j];
379946Ssam b[j]=x&077777;
389946Ssam borrow=x>>15;
399946Ssam }
409946Ssam x=borrow+b[j];
419946Ssam b[j]=x&077777;
429946Ssam if(x>>15 ==0) { return(0);}
439946Ssam borrow=0;
449946Ssam for(j=0;j<n;j++)
459946Ssam { u=a[j]+b[j]+borrow;
469946Ssam if(u<0) borrow=1;
479946Ssam else borrow=0;
489946Ssam b[j]=u&077777;
499946Ssam }
509946Ssam { return(1);}
519946Ssam }
m_trq(v1,v2,u1,u2,u3)529946Ssam m_trq(v1,v2,u1,u2,u3)
539946Ssam { long int d;
549946Ssam long int x1;
559946Ssam if(u1==v1) d=077777;
569946Ssam else d=(u1*0100000L+u2)/v1;
579946Ssam while(1)
589946Ssam { x1=u1*0100000L+u2-v1*d;
599946Ssam x1=x1*0100000L+u3-v2*d;
609946Ssam if(x1<0) d=d-1;
619946Ssam else {return(d);}
629946Ssam }
639946Ssam }
m_div(a,b,q,r)649946Ssam m_div(a,b,q,r) MINT *a,*b,*q,*r;
659946Ssam { MINT u,v,x,w;
669946Ssam short d,*qval;
679946Ssam int qq,n,v1,v2,j;
689946Ssam u.len=v.len=x.len=w.len=0;
699946Ssam if(b->len==0) { fatal("mdiv divide by zero"); return;}
709946Ssam if(b->len==1)
719946Ssam { r->val=xalloc(1,"m_div1");
729946Ssam sdiv(a,b->val[0],q,r->val);
739946Ssam if(r->val[0]==0)
749946Ssam { shfree(r->val);
759946Ssam r->len=0;
769946Ssam }
779946Ssam else r->len=1;
789946Ssam return;
799946Ssam }
809946Ssam if(a->len < b->len)
819946Ssam { q->len=0;
829946Ssam r->len=a->len;
839946Ssam r->val=xalloc(r->len,"m_div2");
849946Ssam for(qq=0;qq<r->len;qq++) r->val[qq]=a->val[qq];
859946Ssam return;
869946Ssam }
879946Ssam x.len=1;
889946Ssam x.val = &d;
899946Ssam n=b->len;
909946Ssam d=0100000L/(b->val[n-1]+1L);
919946Ssam mult(a,&x,&u); /*subtle: relies on fact that mult allocates extra space */
929946Ssam mult(b,&x,&v);
939946Ssam v1=v.val[n-1];
949946Ssam v2=v.val[n-2];
959946Ssam qval=xalloc(a->len-n+1,"m_div3");
969946Ssam for(j=a->len-n;j>=0;j--)
979946Ssam { qq=m_trq(v1,v2,u.val[j+n],u.val[j+n-1],u.val[j+n-2]);
989946Ssam if(m_dsb(qq,n,v.val,&(u.val[j]))) qq -= 1;
999946Ssam qval[j]=qq;
1009946Ssam }
1019946Ssam x.len=n;
1029946Ssam x.val=u.val;
1039946Ssam mcan(&x);
1049946Ssam sdiv(&x,d,&w,(short *)&qq);
1059946Ssam r->len=w.len;
1069946Ssam r->val=w.val;
1079946Ssam q->val=qval;
1089946Ssam qq=a->len-n+1;
1099946Ssam if(qq>0 && qval[qq-1]==0) qq -= 1;
1109946Ssam q->len=qq;
1119946Ssam if(qq==0) shfree(qval);
1129946Ssam if(x.len!=0) xfree(&u);
1139946Ssam xfree(&v);
1149946Ssam return;
1159946Ssam }
116