Bell 32V development

Work on file usr/src/libmp/compall
Work on file usr/src/libmp/mult.c
Work on file usr/src/libmp/mdiv.c
Work on file usr/src/libmp/mout.c
Work on file usr/src/libmp/msqrt.c
Work on file usr/src/libmp/gcd.c
Work on file usr/src/libmp/madd.c
Work on file usr/src/libmp/mklib
Work on file usr/src/libmp/pow.c

Co-Authored-By: John Reiser <jfr@research.uucp>
Synthesized-from: 32v

diff --git a/usr/src/libmp/compall b/usr/src/libmp/compall
new file mode 100755 (executable)
index 0000000..ab8eaf9
--- /dev/null
+++ b/usr/src/libmp/compall
@@ -0,0 +1 @@
+cc -c -O pow.c gcd.c msqrt.c mdiv.c mout.c mult.c madd.c util.c

diff --git a/usr/src/libmp/gcd.c b/usr/src/libmp/gcd.c
new file mode 100644 (file)
index 0000000..cedad0a
--- /dev/null
+++ b/usr/src/libmp/gcd.c
@@ -0,0 +1,46 @@
+#include <mp.h>
+gcd(a,b,c) MINT *a,*b,*c;
+{      MINT x,y,z,w;
+       x.len=y.len=z.len=w.len=0;
+       move(a,&x);
+       move(b,&y);
+       while(y.len!=0)
+       {       mdiv(&x,&y,&w,&z);
+               move(&y,&x);
+               move(&z,&y);
+       }
+       move(&x,c);
+       xfree(&x);
+       xfree(&y);
+       xfree(&z);
+       xfree(&w);
+       return;
+}
+invert(a, b, c) MINT *a, *b, *c;
+{      MINT x, y, z, w, Anew, Aold;
+       int i = 0;
+       x.len = y.len = z.len = w.len = Aold.len = 0;
+       Anew.len = 1;
+       Anew.val = xalloc(1);
+       *Anew.val = 1;
+       move(b, &x);
+       move(a, &y);
+       while(y.len != 0)
+       {       mdiv(&x, &y, &w, &z);
+               move(&Anew, &x);
+               mult(&w, &Anew, &Anew);
+               madd(&Anew, &Aold, &Anew);
+               move(&x, &Aold);
+               move(&y, &x);
+               move(&z, &y);
+               i++;
+       }
+       move(&Aold, c);
+       if( (i&01) == 0) msub(b, c, c);
+       xfree(&x);
+       xfree(&y);
+       xfree(&z);
+       xfree(&w);
+       xfree(&Aold);
+       xfree(&Anew);
+}

diff --git a/usr/src/libmp/madd.c b/usr/src/libmp/madd.c
new file mode 100644 (file)
index 0000000..155d8c0
--- /dev/null
+++ b/usr/src/libmp/madd.c
@@ -0,0 +1,146 @@
+#include <mp.h>
+m_add(a,b,c) struct mint *a,*b,*c;
+{      int carry,i;
+       int x;
+       short *cval;
+       cval=xalloc(a->len+1,"m_add");
+       carry=0;
+       for(i=0;i<b->len;i++)
+       {       x=carry+a->val[i]+b->val[i];
+               if(x&0100000)
+               {       carry=1;
+                       cval[i]=x&077777;
+               }
+               else
+               {       carry=0;
+                       cval[i]=x;
+               }
+       }
+       for(;i<a->len;i++)
+       {       x=carry+a->val[i];
+               if(x&0100000) cval[i]=x&077777;
+               else
+               {       carry=0;
+                       cval[i]=x;
+               }
+
+       }
+       if(carry==1)
+       {       cval[i]=1;
+               c->len=i+1;
+       }
+       else c->len=a->len;
+       c->val=cval;
+       if(c->len==0) shfree(cval);
+       return;
+}
+madd(a,b,c) struct mint *a,*b,*c;
+{      struct mint x,y,z;
+       int sign;
+       x.len=a->len;
+       x.val=a->val;
+       y.len=b->len;
+       y.val=b->val;
+       z.len=0;
+       sign=1;
+       if(x.len>=0)
+               if(y.len>=0)
+                       if(x.len>=y.len) m_add(&x,&y,&z);
+                       else m_add(&y,&x,&z);
+               else
+               {       y.len= -y.len;
+                       msub(&x,&y,&z);
+               }
+       else    if(y.len<=0)
+               {       x.len = -x.len;
+                       y.len= -y.len;
+                       sign= -1;
+                       madd(&x,&y,&z);
+               }
+               else
+               {       x.len= -x.len;
+                       msub(&y,&x,&z);
+               }
+        xfree(c);
+       c->val=z.val;
+       c->len=sign*z.len;
+       return;
+}
+m_sub(a,b,c) struct mint *a,*b,*c;
+{      int x,i;
+       int borrow;
+       short one;
+       struct mint mone;
+       one=1; mone.len= 1; mone.val= &one;
+       c->val=xalloc(a->len,"m_sub");
+       borrow=0;
+       for(i=0;i<b->len;i++)
+       {       x=borrow+a->val[i]-b->val[i];
+               if(x&0100000)
+               {       borrow= -1;
+                       c->val[i]=x&077777;
+               }
+               else
+               {       borrow=0;
+                       c->val[i]=x;
+               }
+       }
+       for(;i<a->len;i++)
+       {       x=borrow+a->val[i];
+               if(x&0100000) c->val[i]=x&077777;
+               else
+               {       borrow=0;
+                       c->val[i]=x;
+               }
+       }
+       if(borrow<0)
+       {       for(i=0;i<a->len;i++) c->val[i] ^= 077777;
+               c->len=a->len;
+               madd(c,&mone,c);
+       }
+       for(i=a->len-1;i>=0;--i) if(c->val[i]>0)
+                               {       if(borrow==0) c->len=i+1;
+                                       else c->len= -i-1;
+                                       return;
+                               }
+       shfree(c->val);
+       return;
+}
+msub(a,b,c) struct mint *a,*b,*c;
+{      struct mint x,y,z;
+       int sign;
+       x.len=a->len;
+       y.len=b->len;
+       x.val=a->val;
+       y.val=b->val;
+       z.len=0;
+       sign=1;
+       if(x.len>=0)
+               if(y.len>=0)
+                       if(x.len>=y.len) m_sub(&x,&y,&z);
+                       else
+                       {       sign= -1;
+                               msub(&y,&x,&z);
+                       }
+               else
+               {       y.len= -y.len;
+                       madd(&x,&y,&z);
+               }
+       else    if(y.len<=0)
+               {       sign= -1;
+                       x.len= -x.len;
+                       y.len= -y.len;
+                       msub(&y,&x,&z);
+               }
+               else
+               {       x.len= -x.len;
+                       madd(&x,&y,&z);
+                       sign= -1;
+               }
+       if(a==c && x.len!=0) xfree(a);
+       else if(b==c && y.len!=0) xfree(b);
+       else xfree(c);
+       c->val=z.val;
+       c->len=sign*z.len;
+       return;
+}

diff --git a/usr/src/libmp/mdiv.c b/usr/src/libmp/mdiv.c
new file mode 100644 (file)
index 0000000..1f8a0fc
--- /dev/null
+++ b/usr/src/libmp/mdiv.c
@@ -0,0 +1,106 @@
+#include <mp.h>
+mdiv(a,b,q,r) MINT *a,*b,*q,*r;
+{      MINT x,y;
+       int sign;
+       sign=1;
+       x.val=a->val;
+       y.val=b->val;
+       x.len=a->len;
+       if(x.len<0) {sign= -1; x.len= -x.len;}
+       y.len=b->len;
+       if(y.len<0) {sign= -sign; y.len= -y.len;}
+       xfree(q);
+       xfree(r);
+       m_div(&x,&y,q,r);
+       if(sign==-1)
+       {       q->len= -q->len;
+               r->len = - r->len;
+       }
+       return;
+}
+m_dsb(q,n,a,b) short *a,*b;
+{      long int x,qx;
+       int borrow,j,u;
+       qx=q;
+       borrow=0;
+       for(j=0;j<n;j++)
+       {       x=borrow-a[j]*qx+b[j];
+               b[j]=x&077777;
+               borrow=x>>15;
+       }
+       x=borrow+b[j];
+       b[j]=x&077777;
+       if(x>>15 ==0) { return(0);}
+       borrow=0;
+       for(j=0;j<n;j++)
+       {       u=a[j]+b[j]+borrow;
+               if(u<0) borrow=1;
+               else borrow=0;
+               b[j]=u&077777;
+       }
+       { return(1);}
+}
+m_trq(v1,v2,u1,u2,u3)
+{      long int d;
+       long int x1;
+       if(u1==v1) d=077777;
+       else d=(u1*0100000L+u2)/v1;
+       while(1)
+       {       x1=u1*0100000L+u2-v1*d;
+               x1=x1*0100000L+u3-v2*d;
+               if(x1<0) d=d-1;
+               else {return(d);}
+       }
+}
+m_div(a,b,q,r) MINT *a,*b,*q,*r;
+{      MINT u,v,x,w;
+       short d,*qval;
+       int qq,n,v1,v2,j;
+       u.len=v.len=x.len=w.len=0;
+       if(b->len==0) { fatal("mdiv divide by zero"); return;}
+       if(b->len==1)
+       {       r->val=xalloc(1,"m_div1");
+               sdiv(a,b->val[0],q,r->val);
+               if(r->val[0]==0)
+               {       shfree(r->val);
+                       r->len=0;
+               }
+               else r->len=1;
+               return;
+       }
+       if(a->len < b->len)
+       {       q->len=0;
+               r->len=a->len;
+               r->val=xalloc(r->len,"m_div2");
+               for(qq=0;qq<r->len;qq++) r->val[qq]=a->val[qq];
+               return;
+       }
+       x.len=1;
+       x.val = &d;
+       n=b->len;
+       d=0100000L/(b->val[n-1]+1L);
+       mult(a,&x,&u); /*subtle: relies on fact that mult allocates extra space */
+       mult(b,&x,&v);
+       v1=v.val[n-1];
+       v2=v.val[n-2];
+       qval=xalloc(a->len-n+1,"m_div3");
+       for(j=a->len-n;j>=0;j--)
+       {       qq=m_trq(v1,v2,u.val[j+n],u.val[j+n-1],u.val[j+n-2]);
+               if(m_dsb(qq,n,v.val,&(u.val[j]))) qq -= 1;
+               qval[j]=qq;
+       }
+       x.len=n;
+       x.val=u.val;
+       mcan(&x);
+       sdiv(&x,d,&w,(short *)&qq);
+       r->len=w.len;
+       r->val=w.val;
+       q->val=qval;
+       qq=a->len-n+1;
+       if(qq>0 && qval[qq-1]==0) qq -= 1;
+       q->len=qq;
+       if(qq==0) shfree(qval);
+       if(x.len!=0) xfree(&u);
+       xfree(&v);
+       return;
+}

diff --git a/usr/src/libmp/mklib b/usr/src/libmp/mklib
new file mode 100755 (executable)
index 0000000..20c16a6
--- /dev/null
+++ b/usr/src/libmp/mklib
@@ -0,0 +1 @@
+ar cr libmp.a pow.o gcd.o msqrt.o mdiv.o mout.o mult.o madd.o util.o

diff --git a/usr/src/libmp/mout.c b/usr/src/libmp/mout.c
new file mode 100644 (file)
index 0000000..afb30a8
--- /dev/null
+++ b/usr/src/libmp/mout.c
@@ -0,0 +1,141 @@
+#include <stdio.h>
+#include <mp.h>
+m_in(a,b,f) MINT *a; FILE *f;
+{      MINT x,y,ten;
+       int sign,c;
+       short qten,qy;
+       xfree(a);
+       sign=1;
+       ten.len=1;
+       ten.val= &qten;
+       qten=b;
+       x.len=0;
+       y.len=1;
+       y.val= &qy;
+       while((c=getc(f))!=EOF)
+       switch(c)
+       {
+       case '\\':      getc(f);
+               continue;
+       case '\t':
+       case '\n': a->len *= sign;
+               xfree(&x);
+               return(0);
+       case ' ':
+               continue;
+       case '-': sign = -sign;
+               continue;
+       default: if(c>='0' && c<= '9')
+               {       qy=c-'0';
+                       mult(&x,&ten,a);
+                       madd(a,&y,a);
+                       move(a,&x);
+                       continue;
+               }
+               else
+               {       VOID ungetc(c,stdin);
+                       a->len *= sign;
+                       return(0);
+               }
+       }
+       return(EOF);
+}
+m_out(a,b,f) MINT *a; FILE *f;
+{      int sign,xlen,i;
+       short r;
+       MINT x;
+       char *obuf;
+       register char *bp;
+       sign=1;
+       xlen=a->len;
+       if(xlen<0)
+       {       xlen= -xlen;
+               sign= -1;
+       }
+       if(xlen==0)
+       {       fprintf(f,"0\n");
+               return;
+       }
+       x.len=xlen;
+       x.val=xalloc(xlen,"m_out");
+       for(i=0;i<xlen;i++) x.val[i]=a->val[i];
+       obuf=(char *)malloc(7*xlen);
+       bp=obuf+7*xlen-1;
+       *bp--=0;
+       while(x.len>0)
+       {       for(i=0;i<10&&x.len>0;i++)
+               {       sdiv(&x,b,&x,&r);
+                       *bp--=r+'0';
+               }
+               if(x.len>0) *bp--=' ';
+       }
+       if(sign==-1) *bp--='-';
+       fprintf(f,"%s\n",bp+1);
+       free(obuf);
+       FREE(x)
+       return;
+}
+sdiv(a,n,q,r) MINT *a,*q; short *r;
+{      MINT x,y;
+       int sign;
+       sign=1;
+       x.len=a->len;
+       x.val=a->val;
+       if(n<0)
+       {       sign= -sign;
+               n= -n;
+       }
+       if(x.len<0)
+       {       sign = -sign;
+               x.len= -x.len;
+       }
+       s_div(&x,n,&y,r);
+       xfree(q);
+       q->val=y.val;
+       q->len=sign*y.len;
+       *r = *r*sign;
+       return;
+}
+s_div(a,n,q,r) MINT *a,*q; short *r;
+{      int qlen,i;
+       long int x;
+       short *qval;
+       x=0;
+       qlen=a->len;
+       qval=xalloc(qlen,"s_div");
+       for(i=qlen-1;i>=0;i--)
+       {
+               x=x*0100000L+a->val[i];
+               qval[i]=x/n;
+               x=x%n;
+       }
+       *r=x;
+       if(qval[qlen-1]==0) qlen--;
+       q->len=qlen;
+       q->val=qval;
+       if(qlen==0) shfree(qval);
+       return;
+}
+min(a) MINT *a;
+{
+       return(m_in(a,10,stdin));
+}
+omin(a) MINT *a;
+{
+       return(m_in(a,8,stdin));
+}
+mout(a) MINT *a;
+{
+       m_out(a,10,stdout);
+}
+omout(a) MINT *a;
+{
+       m_out(a,8,stdout);
+}
+fmout(a,f) MINT *a; FILE *f;
+{      m_out(a,10,f);
+}
+fmin(a,f) MINT *a; FILE *f;
+{
+       return(m_in(a,10,f));
+}

diff --git a/usr/src/libmp/msqrt.c b/usr/src/libmp/msqrt.c
new file mode 100644 (file)
index 0000000..a954593
--- /dev/null
+++ b/usr/src/libmp/msqrt.c
@@ -0,0 +1,37 @@
+#include <mp.h>
+msqrt(a,b,r) MINT *a,*b,*r;
+{      MINT x,junk,y;
+       int j;
+       x.len=junk.len=y.len=0;
+       if(a->len<0) fatal("msqrt: neg arg");
+       if(a->len==0)
+       {       b->len=0;
+               r->len=0;
+               return(0);
+       }
+       if(a->len%2==1) x.len=(1+a->len)/2;
+       else x.len=1+a->len/2;
+       x.val=xalloc(x.len,"msqrt");
+       for(j=0;j<x.len;x.val[j++]=0);
+       if(a->len%2==1) x.val[x.len-1]=0400;
+       else x.val[x.len-1]=1;
+       xfree(b);
+       xfree(r);
+loop:
+       mdiv(a,&x,&y,&junk);
+       xfree(&junk);
+       madd(&x,&y,&y);
+       sdiv(&y,2,&y,(short *)&j);
+       if(mcmp(&x,&y)>0)
+       {       xfree(&x);
+               move(&y,&x);
+               xfree(&y);
+               goto loop;
+       }
+       xfree(&y);
+       move(&x,b);
+       mult(&x,&x,&x);
+       msub(a,&x,r);
+       xfree(&x);
+       return(r->len);
+}

diff --git a/usr/src/libmp/mult.c b/usr/src/libmp/mult.c
new file mode 100644 (file)
index 0000000..d4d8de9
--- /dev/null
+++ b/usr/src/libmp/mult.c
@@ -0,0 +1,87 @@
+#include <mp.h>
+mult(a,b,c) struct mint *a,*b,*c;
+{      struct mint x,y,z;
+       int sign;
+       sign = 1;
+       x.val=a->val;
+       y.val=b->val;
+       z.len=0;
+       if(a->len<0)
+       {       x.len= -a->len;
+               sign= -sign;
+       }
+       else    x.len=a->len;
+       if(b->len<0)
+       {       y.len= -b->len;
+               sign= -sign;
+       }
+       else    y.len=b->len;
+       if(x.len<y.len) m_mult(&y,&x,&z);
+       else m_mult(&x,&y,&z);
+       xfree(c);
+       if(sign<0) c->len= -z.len;
+       else c->len=z.len;
+       if(c->len==0) shfree(z.val);
+       else c->val=z.val;
+       return;
+}
+#define S2 x=a->val[j];
+#define S3 x=x*b->val[i-j];
+#define S4 tradd(&carry,&sum,x);
+#define S5 c->val[i]=sum.yy.low&077777;
+#define S6 sum.xx=sum.xx>>15;
+#define S7 sum.yy.high=carry;
+m_mult(a,b,c) struct mint *a,*b,*c;
+{      long x;
+       union {long xx; struct half yy;} sum;
+       int carry;
+       int i,j;
+       c->val=xalloc(a->len+b->len,"m_mult");
+       sum.xx=0;
+       for(i=0;i<b->len;i++)
+       {       carry=0;
+               for(j=0;j<i+1;j++)
+               {       S2
+                       S3
+                       S4
+               }
+               S5
+               S6
+               S7
+       }
+       for(;i<a->len;i++)
+       {       carry=0;
+               for(j=i-b->len+1;j<i+1;j++)
+               {       S2
+                       S3
+                       S4
+               }
+               S5
+               S6
+               S7
+       }
+       for(;i<a->len+b->len;i++)
+       {       carry=0;
+               for(j=i-b->len+1;j<a->len;j++)
+               {       S2
+                       S3
+                       S4
+               }
+               S5
+               S6
+               S7
+       }
+       if(c->val[i-1]!=0)
+               c->len=a->len+b->len;
+       else    c->len=a->len+b->len-1;
+       return;
+}
+tradd(a,b,c) long c; int *a; union g {long xx; struct half yy;} *b;
+{
+       b->xx= b->xx+c;
+       if(b->yy.high&0100000)
+       {       b->yy.high= b->yy.high&077777;
+               *a += 1;
+       }
+       return;
+}

diff --git a/usr/src/libmp/pow.c b/usr/src/libmp/pow.c
new file mode 100644 (file)
index 0000000..4a745e4
--- /dev/null
+++ b/usr/src/libmp/pow.c
@@ -0,0 +1,40 @@
+#include <mp.h>
+pow(a,b,c,d) MINT *a,*b,*c,*d;
+{      int i,j,n;
+       MINT x,y;
+       x.len=y.len=0;
+       xfree(d);
+       d->len=1;
+       d->val=xalloc(1,"pow");
+       *d->val=1;
+       for(j=0;j<b->len;j++)
+       {       n=b->val[b->len-j-1];
+               for(i=0;i<15;i++)
+               {       mult(d,d,&x);
+                       mdiv(&x,c,&y,d);
+                       if((n=n<<1)&0100000)
+                       {       mult(a,d,&x);
+                               mdiv(&x,c,&y,d);
+                       }
+               }
+       }
+       xfree(&x);
+       xfree(&y);
+       return;
+}
+rpow(a,n,b) MINT *a,*b;
+{      MINT x,y;
+       int i;
+       x.len=1;
+       x.val=xalloc(1,"rpow");
+       *x.val=n;
+       y.len=n*a->len+4;
+       y.val=xalloc(y.len,"rpow2");
+       for(i=0;i<y.len;i++) y.val[i]=0;
+       y.val[y.len-1]=010000;
+       xfree(b);
+       pow(a,&x,&y,b);
+       xfree(&x);
+       xfree(&y);
+       return;
+}