From: Tom London <tbl@research.uucp>
Date: Wed, 24 Jan 1979 13:25:21 +0000 (-0500)
Subject: Bell 32V development
X-Git-Tag: Bell-32V~274
X-Git-Url: https://git.subgeniuskitty.com/unix-history/.git/commitdiff_plain/9cedae4652fbcf3cb0b08b0b73cb31f0f12c6f52

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
index 0000000000..ab8eaf9a41
--- /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
index 0000000000..cedad0ab58
--- /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
index 0000000000..155d8c06ac
--- /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
index 0000000000..1f8a0fc950
--- /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
index 0000000000..20c16a6dc6
--- /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
index 0000000000..afb30a82c5
--- /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
index 0000000000..a9545938be
--- /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
index 0000000000..d4d8de924b
--- /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
index 0000000000..4a745e4cb8
--- /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;
+}