--- /n/sources/plan9/sys/include/mp.h	Thu Jun 27 00:10:11 2024
+++ /sys/include/mp.h	Thu Jun 27 00:16:49 2024
@@ -22,7 +22,9 @@
 
 enum
 {
-	MPstatic=	0x01,
+	MPstatic=	0x01,	/* static constant */
+	MPnorm=		0x02,	/* normalization status */
+	MPtimesafe=	0x04,	/* request time invariant computation */
 	Dbytes=		sizeof(mpdigit),	/* bytes per digit */
 	Dbits=		Dbytes*8		/* bits per digit */
 };
@@ -32,7 +34,7 @@
 mpint*	mpnew(int n);		/* create a new mpint with at least n bits */
 void	mpfree(mpint *b);
 void	mpbits(mpint *b, int n);	/* ensure that b has at least n bits */
-void	mpnorm(mpint *b);		/* dump leading zeros */
+mpint*	mpnorm(mpint *b);		/* dump leading zeros */
 mpint*	mpcopy(mpint *b);
 void	mpassign(mpint *old, mpint *new);
 
@@ -108,12 +110,14 @@
 /* prereq: p has room for n+1 digits */
 int	mpvecdigmulsub(mpdigit *b, int n, mpdigit m, mpdigit *p);
 
-/* p[0:alen*blen-1] = a[0:alen-1] * b[0:blen-1] */
+/* p[0:alen+blen-1] = a[0:alen-1] * b[0:blen-1] */
 /* prereq: alen >= blen, p has room for m*n digits */
 void	mpvecmul(mpdigit *a, int alen, mpdigit *b, int blen, mpdigit *p);
+void	mpvectsmul(mpdigit *a, int alen, mpdigit *b, int blen, mpdigit *p);
 
 /* sign of a - b or zero if the same */
 int	mpveccmp(mpdigit *a, int alen, mpdigit *b, int blen);
+int	mpvectscmp(mpdigit *a, int alen, mpdigit *b, int blen);
 
 /* divide the 2 digit dividend by the one digit divisor and stick in quotient */
 /* we assume that the result is one digit - overflow is all 1's */
--- /n/sources/plan9/sys/include/ape/mp.h	Thu Jun 27 00:10:33 2024
+++ /sys/include/ape/mp.h	Thu Jun 27 00:45:22 2024
@@ -31,7 +31,9 @@
 
 enum
 {
-	MPstatic=	0x01,
+	MPstatic=	0x01,	/* static constant */
+	MPnorm=		0x02,	/* normalization status */
+	MPtimesafe=	0x04,	/* request time invariant computation */
 	Dbytes=		sizeof(mpdigit),	/* bytes per digit */
 	Dbits=		Dbytes*8		/* bits per digit */
 };
@@ -41,7 +43,7 @@
 mpint*	mpnew(int n);		/* create a new mpint with at least n bits */
 void	mpfree(mpint *b);
 void	mpbits(mpint *b, int n);	/* ensure that b has at least n bits */
-void	mpnorm(mpint *b);		/* dump leading zeros */
+mpint*	mpnorm(mpint *b);		/* dump leading zeros */
 mpint*	mpcopy(mpint *b);
 void	mpassign(mpint *old, mpint *new);
 
@@ -117,12 +119,14 @@
 /* prereq: p has room for n+1 digits */
 int	mpvecdigmulsub(mpdigit *b, int n, mpdigit m, mpdigit *p);
 
-/* p[0:alen*blen-1] = a[0:alen-1] * b[0:blen-1] */
+/* p[0:alen+blen-1] = a[0:alen-1] * b[0:blen-1] */
 /* prereq: alen >= blen, p has room for m*n digits */
 void	mpvecmul(mpdigit *a, int alen, mpdigit *b, int blen, mpdigit *p);
+void	mpvectsmul(mpdigit *a, int alen, mpdigit *b, int blen, mpdigit *p);
 
 /* sign of a - b or zero if the same */
 int	mpveccmp(mpdigit *a, int alen, mpdigit *b, int blen);
+int	mpvectscmp(mpdigit *a, int alen, mpdigit *b, int blen);
 
 /* divide the 2 digit dividend by the one digit divisor and stick in quotient */
 /* we assume that the result is one digit - overflow is all 1's */
--- /n/sources/plan9/sys/man/2/mp	Thu Jun 27 00:11:01 2024
+++ /sys/man/2/mp	Thu Jun 27 00:16:50 2024
@@ -22,7 +22,7 @@
 void	mpbits(mpint *b, int n)
 .PP
 .B
-void	mpnorm(mpint *b)
+mpint*	mpnorm(mpint *b)
 .PP
 .B
 mpint*	mpcopy(mpint *b)
@@ -588,8 +588,8 @@
 -1 if negative.
 .TP
 .I mpvecmul
-.BR "p[0:alen*blen] = a[0:alen-1] * b[0:blen-1]" .
-We assume that p has room for alen*blen+1 digits.
+.BR "p[0:alen+blen] = a[0:alen-1] * b[0:blen-1]" .
+We assume that p has room for alen+blen+1 digits.
 .TP
 .I mpveccmp
 This returns -1, 0, or +1 as a - b is negative, 0, or positive.
@@ -600,6 +600,17 @@
 and
 .I mpzero
 are the constants 2, 1 and 0.  These cannot be freed.
+.SS "Time invariant computation"
+.PP
+In the field of cryptography, it is sometimes necessary to implement
+algorithms such that the runtime of the algorithm is not dependent on
+the input data. This library provides partial support for time
+invariant computation with the
+.I MPtimesafe
+flag that can be set on input or destination operands to request timing
+safe operation. The result of a timing safe operation will also have the
+.I MPtimesafe
+flag set and is not normalized.
 .SS "Chinese remainder theorem
 .PP
 When computing in a non-prime modulus, 
diff -Nru /n/sources/plan9/sys/src/ape/lib/mp/port/libc.h /sys/src/ape/lib/mp/port/libc.h
--- /n/sources/plan9/sys/src/ape/lib/mp/port/libc.h	Mon Jan 28 17:24:15 2013
+++ /sys/src/ape/lib/mp/port/libc.h	Thu Jun 27 01:11:12 2024
@@ -4,6 +4,7 @@
 #include <string.h>
 #include <utf.h>
 #include <fmt.h>
+#include <assert.h>
 
 typedef unsigned int u32int;
 typedef unsigned long long u64int;
@@ -22,5 +23,7 @@
 extern int  enc64(char *, int, uchar *, int);
 
 extern	vlong	nsec(void);
+
+extern int	_tas(int*);
 
 extern void sysfatal(char*, ...);
diff -Nru /n/sources/plan9/sys/src/ape/lib/mp/port/mkfile /sys/src/ape/lib/mp/port/mkfile
--- /n/sources/plan9/sys/src/ape/lib/mp/port/mkfile	Thu Jun 27 00:04:08 2024
+++ /sys/src/ape/lib/mp/port/mkfile	Thu Jun 27 00:35:51 2024
@@ -24,6 +24,7 @@
 	mpvecsub\
 	mpvecdigmuladd\
 	mpveccmp\
+	mpvectscmp\
 	mpdigdiv\
 	mpdiv\
 	mpexp\
diff -Nru /n/sources/plan9/sys/src/libmp/port/betomp.c /sys/src/libmp/port/betomp.c
--- /n/sources/plan9/sys/src/libmp/port/betomp.c	Mon Sep 19 07:19:11 2005
+++ /sys/src/libmp/port/betomp.c	Thu Jun 27 00:16:50 2024
@@ -13,19 +13,12 @@
 		b = mpnew(0);
 		setmalloctag(b, getcallerpc(&p));
 	}
-
-	// dump leading zeros
-	while(*p == 0 && n > 1){
-		p++;
-		n--;
-	}
-
-	// get the space
 	mpbits(b, n*8);
-	b->top = DIGITS(n*8);
-	m = b->top-1;
 
-	// first digit might not be Dbytes long
+	m = DIGITS(n*8);
+	b->top = m--;
+	b->sign = 1;
+
 	s = ((n-1)*8)%Dbits;
 	x = 0;
 	for(; n > 0; n--){
@@ -37,6 +30,5 @@
 			x = 0;
 		}
 	}
-
-	return b;
+	return mpnorm(b);
 }
diff -Nru /n/sources/plan9/sys/src/libmp/port/letomp.c /sys/src/libmp/port/letomp.c
--- /n/sources/plan9/sys/src/libmp/port/letomp.c	Wed Jun 26 23:28:43 2024
+++ /sys/src/libmp/port/letomp.c	Thu Jun 27 00:16:50 2024
@@ -26,5 +26,6 @@
 	if(i > 0)
 		b->p[m++] = x;
 	b->top = m;
-	return b;
+	b->sign = 1;
+	return mpnorm(b);
 }
diff -Nru /n/sources/plan9/sys/src/libmp/port/mkfile /sys/src/libmp/port/mkfile
--- /n/sources/plan9/sys/src/libmp/port/mkfile	Wed Jun 26 23:59:55 2024
+++ /sys/src/libmp/port/mkfile	Thu Jun 27 00:16:50 2024
@@ -22,6 +22,7 @@
 	mpvecsub\
 	mpvecdigmuladd\
 	mpveccmp\
+	mpvectscmp\
 	mpdigdiv\
 	mpdiv\
 	mpexp\
diff -Nru /n/sources/plan9/sys/src/libmp/port/mpadd.c /sys/src/libmp/port/mpadd.c
--- /n/sources/plan9/sys/src/libmp/port/mpadd.c	Thu Mar 15 21:43:44 2001
+++ /sys/src/libmp/port/mpadd.c	Thu Jun 27 00:16:52 2024
@@ -9,6 +9,8 @@
 	int m, n;
 	mpint *t;
 
+	sum->flags |= (b1->flags | b2->flags) & MPtimesafe;
+
 	// get the sizes right
 	if(b2->top > b1->top){
 		t = b1;
@@ -41,6 +43,7 @@
 	int sign;
 
 	if(b1->sign != b2->sign){
+		assert(((b1->flags | b2->flags | sum->flags) & MPtimesafe) == 0);
 		if(b1->sign < 0)
 			mpmagsub(b2, b1, sum);
 		else
diff -Nru /n/sources/plan9/sys/src/libmp/port/mpaux.c /sys/src/libmp/port/mpaux.c
--- /n/sources/plan9/sys/src/libmp/port/mpaux.c	Wed Jun 26 23:28:52 2024
+++ /sys/src/libmp/port/mpaux.c	Thu Jun 27 00:30:49 2024
@@ -9,7 +9,7 @@
 	1,
 	1,
 	_mptwodata,
-	MPstatic
+	MPstatic|MPnorm
 };
 mpint *mptwo = &_mptwo;
 
@@ -20,7 +20,7 @@
 	1,
 	1,
 	_mponedata,
-	MPstatic
+	MPstatic|MPnorm
 };
 mpint *mpone = &_mpone;
 
@@ -31,7 +31,7 @@
 	1,
 	0,
 	_mpzerodata,
-	MPstatic
+	MPstatic|MPnorm
 };
 mpint *mpzero = &_mpzero;
 
@@ -57,18 +57,17 @@
 	if(n < 0)
 		sysfatal("mpsetminbits: n < 0");
 
-	b = mallocz(sizeof(mpint), 1);
-	if(b == nil)
-		sysfatal("mpnew: %r");
-	setmalloctag(b, getcallerpc(&n));
 	n = DIGITS(n);
 	if(n < mpmindigits)
 		n = mpmindigits;
-	b->p = (mpdigit*)mallocz(n*Dbytes, 1);
-	if(b->p == nil)
+	b = mallocz(sizeof(mpint) + n*Dbytes, 1);
+	if(b == nil)
 		sysfatal("mpnew: %r");
+	setmalloctag(b, getcallerpc(&n));
+	b->p = (mpdigit*)&b[1];
 	b->size = n;
 	b->sign = 1;
+	b->flags = MPnorm;
 
 	return b;
 }
@@ -83,16 +82,23 @@
 	if(b->size >= n){
 		if(b->top >= n)
 			return;
-		memset(&b->p[b->top], 0, Dbytes*(n - b->top));
-		b->top = n;
-		return;
+	} else {
+		if(b->p == (mpdigit*)&b[1]){
+			b->p = (mpdigit*)mallocz(n*Dbytes, 0);
+			if(b->p == nil)
+				sysfatal("mpbits: %r");
+			memmove(b->p, &b[1], Dbytes*b->top);
+			memset(&b[1], 0, Dbytes*b->size);
+		} else {
+			b->p = (mpdigit*)realloc(b->p, n*Dbytes);
+			if(b->p == nil)
+				sysfatal("mpbits: %r");
+		}
+		b->size = n;
 	}
-	b->p = (mpdigit*)realloc(b->p, n*Dbytes);
-	if(b->p == nil)
-		sysfatal("mpbits: %r");
 	memset(&b->p[b->top], 0, Dbytes*(n - b->top));
-	b->size = n;
 	b->top = n;
+	b->flags &= ~MPnorm;
 }
 
 void
@@ -103,21 +109,29 @@
 	if(b->flags & MPstatic)
 		sysfatal("freeing mp constant");
 	memset(b->p, 0, b->size*Dbytes);	// information hiding
-	free(b->p);
+	if(b->p != (mpdigit*)&b[1])
+		free(b->p);
 	free(b);
 }
 
-void
+mpint*
 mpnorm(mpint *b)
 {
 	int i;
 
+	if(b->flags & MPtimesafe){
+		assert(b->sign == 1);
+		b->flags &= ~MPnorm;
+		return b;
+	}
 	for(i = b->top-1; i >= 0; i--)
 		if(b->p[i] != 0)
 			break;
 	b->top = i+1;
 	if(b->top == 0)
 		b->sign = 1;
+	b->flags |= MPnorm;
+	return b;
 }
 
 mpint*
@@ -129,6 +143,7 @@
 	setmalloctag(new, getcallerpc(&old));
 	new->top = old->top;
 	new->sign = old->sign;
+	new->flags = old->flags & ~MPstatic;
 	memmove(new->p, old->p, Dbytes*old->top);
 	return new;
 }
@@ -136,9 +151,14 @@
 void
 mpassign(mpint *old, mpint *new)
 {
+	if(new == nil || old == new)
+		return;
+	new->top = 0;
 	mpbits(new, Dbits*old->top);
 	new->sign = old->sign;
 	new->top = old->top;
+	new->flags &= ~MPnorm;
+	new->flags |= old->flags & ~MPstatic;
 	memmove(new->p, old->p, Dbytes*old->top);
 }
 
@@ -168,6 +188,7 @@
 	int k, bit, digit;
 	mpdigit d;
 
+	assert(n->flags & MPnorm);
 	if(n->top==0)
 		return 0;
 	k = 0;
@@ -188,4 +209,3 @@
 	}
 	return k;
 }
-
diff -Nru /n/sources/plan9/sys/src/libmp/port/mpcmp.c /sys/src/libmp/port/mpcmp.c
--- /n/sources/plan9/sys/src/libmp/port/mpcmp.c	Sun Jun 15 14:03:30 2003
+++ /sys/src/libmp/port/mpcmp.c	Thu Jun 27 00:16:52 2024
@@ -8,10 +8,14 @@
 {
 	int i;
 
-	i = b1->top - b2->top;
-	if(i)
-		return i;
-
+	i = b1->flags | b2->flags;
+	if(i & MPtimesafe)
+		return mpvectscmp(b1->p, b1->top, b2->p, b2->top);
+	if(i & MPnorm){
+		i = b1->top - b2->top;
+		if(i)
+			return i;
+	}
 	return mpveccmp(b1->p, b1->top, b2->p, b2->top);
 }
 
@@ -19,10 +23,8 @@
 int
 mpcmp(mpint *b1, mpint *b2)
 {
-	if(b1->sign != b2->sign)
-		return b1->sign - b2->sign;
-	if(b1->sign < 0)
-		return mpmagcmp(b2, b1);
-	else
-		return mpmagcmp(b1, b2);
+	int sign;
+
+	sign = (b1->sign - b2->sign) >> 1;	// -1, 0, 1
+	return sign | (sign&1)-1 & mpmagcmp(b1, b2)*b1->sign;
 }
diff -Nru /n/sources/plan9/sys/src/libmp/port/mpdiv.c /sys/src/libmp/port/mpdiv.c
--- /n/sources/plan9/sys/src/libmp/port/mpdiv.c	Tue Jul 25 11:59:24 2000
+++ /sys/src/libmp/port/mpdiv.c	Thu Jun 27 00:16:52 2024
@@ -13,10 +13,29 @@
 	mpdigit qd, *up, *vp, *qp;
 	mpint *u, *v, *t;
 
+	assert(quotient != remainder);
+	assert(divisor->flags & MPnorm);
+
 	// divide bv zero
 	if(divisor->top == 0)
 		abort();
 
+	// division by one or small powers of two
+	if(divisor->top == 1 && (divisor->p[0] & divisor->p[0]-1) == 0){
+		vlong r = (vlong)dividend->sign * (dividend->p[0] & divisor->p[0]-1);
+		if(quotient != nil){
+			for(s = 0; ((divisor->p[0] >> s) & 1) == 0; s++)
+				;
+			mpright(dividend, s, quotient);
+		}
+		if(remainder != nil){
+			remainder->flags |= dividend->flags & MPtimesafe;
+			vtomp(r, remainder);
+		}
+		return;
+	}
+	assert((dividend->flags & MPtimesafe) == 0);
+
 	// quick check
 	if(mpmagcmp(dividend, divisor) < 0){
 		if(remainder != nil)
@@ -95,12 +114,14 @@
 		*up-- = 0;
 	}
 	if(qp != nil){
+		assert((quotient->flags & MPtimesafe) == 0);
 		mpnorm(quotient);
 		if(dividend->sign != divisor->sign)
 			quotient->sign = -1;
 	}
 
 	if(remainder != nil){
+		assert((remainder->flags & MPtimesafe) == 0);
 		mpright(u, s, remainder);	// u is the remainder shifted
 		remainder->sign = dividend->sign;
 	}
diff -Nru /n/sources/plan9/sys/src/libmp/port/mpeuclid.c /sys/src/libmp/port/mpeuclid.c
--- /n/sources/plan9/sys/src/libmp/port/mpeuclid.c	Wed Sep 17 22:56:35 2003
+++ /sys/src/libmp/port/mpeuclid.c	Thu Jun 27 00:16:52 2024
@@ -13,6 +13,9 @@
 {
 	mpint *tmp, *x0, *x1, *x2, *y0, *y1, *y2, *q, *r;
 
+	assert((a->flags&b->flags) & MPnorm);
+	assert(((a->flags|b->flags|d->flags|x->flags|y->flags) & MPtimesafe) == 0);
+
 	if(a->sign<0 || b->sign<0)
 		sysfatal("mpeuclid: negative arg");
 
diff -Nru /n/sources/plan9/sys/src/libmp/port/mpexp.c /sys/src/libmp/port/mpexp.c
--- /n/sources/plan9/sys/src/libmp/port/mpexp.c	Sun Jun 15 14:04:21 2003
+++ /sys/src/libmp/port/mpexp.c	Thu Jun 27 00:16:53 2024
@@ -22,6 +22,10 @@
 	mpdigit d, bit;
 	int i, j;
 
+	assert(m->flags & MPnorm);
+	assert((e->flags & MPtimesafe) == 0);
+	res->flags |= b->flags & MPtimesafe;
+
 	i = mpcmp(e,mpzero);
 	if(i==0){
 		mpassign(mpone, res);
diff -Nru /n/sources/plan9/sys/src/libmp/port/mpextendedgcd.c /sys/src/libmp/port/mpextendedgcd.c
--- /n/sources/plan9/sys/src/libmp/port/mpextendedgcd.c	Fri Dec 12 09:21:48 2003
+++ /sys/src/libmp/port/mpextendedgcd.c	Thu Jun 27 00:16:54 2024
@@ -15,6 +15,9 @@
 	mpint *u, *A, *B, *C, *D;
 	int g;
 
+	assert((a->flags&b->flags) & MPnorm);
+	assert(((a->flags|b->flags|v->flags|x->flags|y->flags) & MPtimesafe) == 0);
+
 	if(a->sign < 0 || b->sign < 0){
 		mpassign(mpzero, v);
 		mpassign(mpzero, y);
diff -Nru /n/sources/plan9/sys/src/libmp/port/mpfmt.c /sys/src/libmp/port/mpfmt.c
--- /n/sources/plan9/sys/src/libmp/port/mpfmt.c	Mon Feb 17 12:16:26 2003
+++ /sys/src/libmp/port/mpfmt.c	Thu Jun 27 00:16:54 2024
@@ -102,6 +102,7 @@
 		return -1;
 
 	d = mpcopy(b);
+	mpnorm(d);
 	r = mpnew(0);
 	billion = uitomp(1000000000, nil);
 	out = buf+len;
@@ -128,14 +129,19 @@
 mpfmt(Fmt *fmt)
 {
 	mpint *b;
-	char *p;
+	char *p, f;
 
 	b = va_arg(fmt->args, mpint*);
 	if(b == nil)
 		return fmtstrcpy(fmt, "*");
-	
+
+	f = b->flags;
+	b->flags &= ~MPtimesafe;
+
 	p = mptoa(b, fmt->prec, nil, 0);
 	fmt->flags &= ~FmtPrec;
+
+	b->flags = f;
 
 	if(p == nil)
 		return fmtstrcpy(fmt, "*");
diff -Nru /n/sources/plan9/sys/src/libmp/port/mpleft.c /sys/src/libmp/port/mpleft.c
--- /n/sources/plan9/sys/src/libmp/port/mpleft.c	Sun Jun 15 14:04:22 2003
+++ /sys/src/libmp/port/mpleft.c	Thu Jun 27 00:16:54 2024
@@ -15,8 +15,8 @@
 		return;
 	}
 
-	// a negative left shift is a right shift
-	if(shift < 0){
+	// a zero or negative left shift is a right shift
+	if(shift <= 0){
 		mpright(b, -shift, res);
 		return;
 	}
@@ -46,7 +46,6 @@
 	for(i = 0; i < d; i++)
 		res->p[i] = 0;
 
-	// normalize
-	while(res->top > 0 && res->p[res->top-1] == 0)
-		res->top--;
+	res->flags |= b->flags & MPtimesafe;
+	mpnorm(res);
 }
diff -Nru /n/sources/plan9/sys/src/libmp/port/mpmod.c /sys/src/libmp/port/mpmod.c
--- /n/sources/plan9/sys/src/libmp/port/mpmod.c	Wed Feb  9 08:58:45 2000
+++ /sys/src/libmp/port/mpmod.c	Thu Jun 27 00:16:54 2024
@@ -2,14 +2,100 @@
 #include <mp.h>
 #include "dat.h"
 
-// remainder = b mod m
-//
-// knuth, vol 2, pp 398-400
-
 void
-mpmod(mpint *b, mpint *m, mpint *remainder)
+mpmod(mpint *x, mpint *n, mpint *r)
 {
-	mpdiv(b, m, nil, remainder);
-	if(remainder->sign < 0)
-		mpadd(m, remainder, remainder);
+	static int busy;
+	static mpint *p, *m, *c, *v;
+	mpdigit q[32], t[64], d;
+	int sign, k, s, qn, tn;
+
+	sign = x->sign;
+
+	assert(n->flags & MPnorm);
+	if(n->top < 2 || n->top > nelem(q) || (x->top-n->top) > nelem(q))
+		goto hard;
+
+	/*
+	 * check if n = 2**k - c where c has few power of two factors
+	 * above the lowest digit.
+	 */
+	for(k = n->top-1; k > 0; k--){
+		d = n->p[k] >> 1;
+		if((d+1 & d) != 0)
+			goto hard;
+	}
+
+	d = n->p[n->top-1];
+	for(s = 0; (d & (mpdigit)1<<Dbits-1) == 0; s++)
+		d <<= 1;
+
+	/* lo(x) = x[0:k-1], hi(x) = x[k:xn-1] */
+	k = n->top;
+
+	while(_tas(&busy))
+		;
+
+	if(p == nil || mpmagcmp(n, p) != 0){
+		if(m == nil){
+			m = mpnew(0);
+			c = mpnew(0);
+			p = mpnew(0);
+		}
+		mpassign(n, p);
+
+		mpleft(n, s, m);
+		mpleft(mpone, k*Dbits, c);
+		mpsub(c, m, c);
+	}
+
+	mpleft(x, s, r);
+	if(r->top <= k){
+		mpbits(r, (k+1)*Dbits);
+		r->top = k+1;
+	}
+
+	/* q = hi(r) */
+	qn = r->top - k;
+	memmove(q, r->p+k, qn*Dbytes);
+
+	/* r = lo(r) */
+	r->top = k;
+
+	do {
+		/* t = q*c */
+		tn = qn + c->top;
+		memset(t, 0, tn*Dbytes);
+		mpvecmul(q, qn, c->p, c->top, t);
+
+		/* q = hi(t) */
+		qn = tn - k;
+		if(qn <= 0) qn = 0;
+		else memmove(q, t+k, qn*Dbytes);
+
+		/* r += lo(t) */
+		if(tn > k)
+			tn = k;
+		mpvecadd(r->p, k, t, tn, r->p);
+
+		/* if(r >= m) r -= m */
+		mpvecsub(r->p, k+1, m->p, k, t), d = t[k];
+		for(tn = 0; tn < k; tn++)
+			r->p[tn] = (r->p[tn] & d) | (t[tn] & ~d);
+	} while(qn > 0);
+
+	busy = 0;
+
+	if(s != 0)
+		mpright(r, s, r);
+	else
+		mpnorm(r);
+	goto done;
+
+hard:
+	mpdiv(x, n, nil, r);
+
+done:
+	if(sign < 0)
+		mpmagsub(n, r, r);
 }
diff -Nru /n/sources/plan9/sys/src/libmp/port/mpmul.c /sys/src/libmp/port/mpmul.c
--- /n/sources/plan9/sys/src/libmp/port/mpmul.c	Thu Feb 28 16:08:53 2002
+++ /sys/src/libmp/port/mpmul.c	Thu Jun 27 00:16:55 2024
@@ -113,10 +113,6 @@
 		a = b;
 		b = t;
 	}
-	if(blen == 0){
-		memset(p, 0, Dbytes*(alen+blen));
-		return;
-	}
 
 	if(alen >= KARATSUBAMIN && blen > 1){
 		// O(n^1.585)
@@ -132,24 +128,48 @@
 }
 
 void
+mpvectsmul(mpdigit *a, int alen, mpdigit *b, int blen, mpdigit *p)
+{
+	int i;
+	mpdigit *t;
+
+	if(alen < blen){
+		i = alen;
+		alen = blen;
+		blen = i;
+		t = a;
+		a = b;
+		b = t;
+	}
+	if(blen == 0)
+		return;
+	for(i = 0; i < blen; i++)
+		mpvecdigmuladd(a, alen, b[i], &p[i]);
+}
+
+void
 mpmul(mpint *b1, mpint *b2, mpint *prod)
 {
 	mpint *oprod;
 
-	oprod = nil;
+	oprod = prod;
 	if(prod == b1 || prod == b2){
-		oprod = prod;
 		prod = mpnew(0);
+		prod->flags = oprod->flags;
 	}
+	prod->flags |= (b1->flags | b2->flags) & MPtimesafe;
 
 	prod->top = 0;
 	mpbits(prod, (b1->top+b2->top+1)*Dbits);
-	mpvecmul(b1->p, b1->top, b2->p, b2->top, prod->p);
+	if(prod->flags & MPtimesafe)
+		mpvectsmul(b1->p, b1->top, b2->p, b2->top, prod->p);
+	else
+		mpvecmul(b1->p, b1->top, b2->p, b2->top, prod->p);
 	prod->top = b1->top+b2->top+1;
 	prod->sign = b1->sign*b2->sign;
 	mpnorm(prod);
 
-	if(oprod != nil){
+	if(oprod != prod){
 		mpassign(prod, oprod);
 		mpfree(prod);
 	}
diff -Nru /n/sources/plan9/sys/src/libmp/port/mprand.c /sys/src/libmp/port/mprand.c
--- /n/sources/plan9/sys/src/libmp/port/mprand.c	Wed Jun 26 23:28:55 2024
+++ /sys/src/libmp/port/mprand.c	Thu Jun 27 00:16:55 2024
@@ -6,8 +6,8 @@
 mpint*
 mprand(int bits, void (*gen)(uchar*, int), mpint *b)
 {
-	int n, m;
 	mpdigit mask;
+	int n, m;
 	uchar *p;
 
 	n = DIGITS(bits);
@@ -26,18 +26,12 @@
 
 	// make sure we don't give too many bits
 	m = bits%Dbits;
-	n--;
-	if(m > 0){
-		mask = 1;
-		mask <<= m;
-		mask--;
-		b->p[n] &= mask;
-	}
+	if(m == 0)
+		return b;
 
-	for(; n >= 0; n--)
-		if(b->p[n] != 0)
-			break;
-	b->top = n+1;
-	b->sign = 1;
-	return b;
+	mask = 1;
+	mask <<= m;
+	mask--;
+	b->p[n-1] &= mask;
+	return mpnorm(b);
 }
diff -Nru /n/sources/plan9/sys/src/libmp/port/mpright.c /sys/src/libmp/port/mpright.c
--- /n/sources/plan9/sys/src/libmp/port/mpright.c	Sun Jun 15 14:04:22 2003
+++ /sys/src/libmp/port/mpright.c	Thu Jun 27 00:16:56 2024
@@ -23,12 +23,16 @@
 
 	if(res != b)
 		mpbits(res, b->top*Dbits - shift);
+	else if(shift == 0)
+		return;
+
 	d = shift/Dbits;
 	r = shift - d*Dbits;
 	l = Dbits - r;
 
 	//  shift all the bits out == zero
 	if(d>=b->top){
+		res->sign = 1;
 		res->top = 0;
 		return;
 	}
@@ -46,9 +50,8 @@
 		}
 		res->p[i++] = last>>r;
 	}
-	while(i > 0 && res->p[i-1] == 0)
-		i--;
+
 	res->top = i;
-	if(i==0)
-		res->sign = 1;
+	res->flags |= b->flags & MPtimesafe;
+	mpnorm(res);
 }
diff -Nru /n/sources/plan9/sys/src/libmp/port/mpsub.c /sys/src/libmp/port/mpsub.c
--- /n/sources/plan9/sys/src/libmp/port/mpsub.c	Thu Mar 15 21:43:45 2001
+++ /sys/src/libmp/port/mpsub.c	Thu Jun 27 00:16:56 2024
@@ -11,12 +11,15 @@
 
 	// get the sizes right
 	if(mpmagcmp(b1, b2) < 0){
+		assert(((b1->flags | b2->flags | diff->flags) & MPtimesafe) == 0);
 		sign = -1;
 		t = b1;
 		b1 = b2;
 		b2 = t;
-	} else
+	} else {
+		diff->flags |= (b1->flags | b2->flags) & MPtimesafe;
 		sign = 1;
+	}
 	n = b1->top;
 	m = b2->top;
 	if(m == 0){
@@ -39,6 +42,7 @@
 	int sign;
 
 	if(b1->sign != b2->sign){
+		assert(((b1->flags | b2->flags | diff->flags) & MPtimesafe) == 0);
 		sign = b1->sign;
 		mpmagadd(b1, b2, diff);
 		diff->sign = sign;
diff -Nru /n/sources/plan9/sys/src/libmp/port/mptoi.c /sys/src/libmp/port/mptoi.c
--- /n/sources/plan9/sys/src/libmp/port/mptoi.c	Wed Jun 26 23:29:01 2024
+++ /sys/src/libmp/port/mptoi.c	Thu Jun 27 00:16:56 2024
@@ -14,15 +14,11 @@
 		b = mpnew(0);
 		setmalloctag(b, getcallerpc(&i));
 	}
-	mpassign(mpzero, b);
-	if(i != 0)
-		b->top = 1;
-	if(i < 0){
-		b->sign = -1;
-		*b->p = -i;
-	} else
-		*b->p = i;
-	return b;
+	b->sign = (i >> (sizeof(i)*8 - 1)) | 1;
+	i *= b->sign;
+	*b->p = i;
+	b->top = 1;
+	return mpnorm(b);
 }
 
 int
diff -Nru /n/sources/plan9/sys/src/libmp/port/mptoui.c /sys/src/libmp/port/mptoui.c
--- /n/sources/plan9/sys/src/libmp/port/mptoui.c	Wed Jun 26 23:29:08 2024
+++ /sys/src/libmp/port/mptoui.c	Thu Jun 27 00:16:57 2024
@@ -14,11 +14,10 @@
 		b = mpnew(0);
 		setmalloctag(b, getcallerpc(&i));
 	}
-	mpassign(mpzero, b);
-	if(i != 0)
-		b->top = 1;
 	*b->p = i;
-	return b;
+	b->top = 1;
+	b->sign = 1;
+	return mpnorm(b);
 }
 
 uint
diff -Nru /n/sources/plan9/sys/src/libmp/port/mptouv.c /sys/src/libmp/port/mptouv.c
--- /n/sources/plan9/sys/src/libmp/port/mptouv.c	Wed Jun 26 23:29:12 2024
+++ /sys/src/libmp/port/mptouv.c	Thu Jun 27 00:16:57 2024
@@ -13,20 +13,17 @@
 {
 	int s;
 
-	if(b == nil){
+	if(b == nil)
 		b = mpnew(VLDIGITS*sizeof(mpdigit));
-		setmalloctag(b, getcallerpc(&v));
-	}else
+	else
 		mpbits(b, VLDIGITS*sizeof(mpdigit));
-	mpassign(mpzero, b);
-	if(v == 0)
-		return b;
-	for(s = 0; s < VLDIGITS && v != 0; s++){
+	b->sign = 1;
+	for(s = 0; s < VLDIGITS; s++){
 		b->p[s] = v;
 		v >>= sizeof(mpdigit)*8;
 	}
 	b->top = s;
-	return b;
+	return mpnorm(b);
 }
 
 uvlong
@@ -38,7 +35,6 @@
 	if(b->top == 0)
 		return 0LL;
 
-	mpnorm(b);
 	if(b->top > VLDIGITS)
 		return MAXVLONG;
 
diff -Nru /n/sources/plan9/sys/src/libmp/port/mptov.c /sys/src/libmp/port/mptov.c
--- /n/sources/plan9/sys/src/libmp/port/mptov.c	Wed Jun 26 23:29:13 2024
+++ /sys/src/libmp/port/mptov.c	Thu Jun 27 00:16:58 2024
@@ -19,20 +19,14 @@
 		setmalloctag(b, getcallerpc(&v));
 	}else
 		mpbits(b, VLDIGITS*sizeof(mpdigit));
-	mpassign(mpzero, b);
-	if(v == 0)
-		return b;
-	if(v < 0){
-		b->sign = -1;
-		uv = -v;
-	} else
-		uv = v;
-	for(s = 0; s < VLDIGITS && uv != 0; s++){
+	b->sign = (v >> (sizeof(v)*8 - 1)) | 1;
+	uv = v * b->sign;
+	for(s = 0; s < VLDIGITS; s++){
 		b->p[s] = uv;
 		uv >>= sizeof(mpdigit)*8;
 	}
 	b->top = s;
-	return b;
+	return mpnorm(b);
 }
 
 vlong
@@ -44,7 +38,6 @@
 	if(b->top == 0)
 		return 0LL;
 
-	mpnorm(b);
 	if(b->top > VLDIGITS){
 		if(b->sign > 0)
 			return (vlong)MAXVLONG;
diff -Nru /n/sources/plan9/sys/src/libmp/port/mpvectscmp.c /sys/src/libmp/port/mpvectscmp.c
--- /n/sources/plan9/sys/src/libmp/port/mpvectscmp.c	Wed Dec 31 19:00:00 1969
+++ /sys/src/libmp/port/mpvectscmp.c	Thu Jun 27 00:17:00 2024
@@ -0,0 +1,34 @@
+#include "os.h"
+#include <mp.h>
+#include "dat.h"
+
+int
+mpvectscmp(mpdigit *a, int alen, mpdigit *b, int blen)
+{
+	mpdigit x, y, z, v;
+	int m, p;
+
+	if(alen > blen){
+		v = 0;
+		while(alen > blen)
+			v |= a[--alen];
+		m = p = (-v^v|v)>>Dbits-1;
+	} else if(blen > alen){
+		v = 0;
+		while(blen > alen)
+			v |= b[--blen];
+		m = (-v^v|v)>>Dbits-1;
+		p = m^1;
+	} else
+		m = p = 0;
+	while(alen-- > 0){
+		x = a[alen];
+		y = b[alen];
+		z = x - y;
+		x = ~x;
+		v = ((-z^z|z)>>Dbits-1) & ~m;
+		p = ((~(x&y|x&z|y&z)>>Dbits-1) & v) | (p & ~v);
+		m |= v;
+	}
+	return (p-m) | m;
+}
diff -Nru /n/sources/plan9/sys/src/libmp/port/strtomp.c /sys/src/libmp/port/strtomp.c
--- /n/sources/plan9/sys/src/libmp/port/strtomp.c	Wed Jun 26 23:29:15 2024
+++ /sys/src/libmp/port/strtomp.c	Thu Jun 27 00:17:00 2024
@@ -50,7 +50,6 @@
 	int i;
 	mpdigit x;
 
-	b->top = 0;
 	for(p = a; *p; p++)
 		if(tab.t16[*(uchar*)p] == INVAL)
 			break;
@@ -144,12 +143,9 @@
 	mpbits(b, n*5);
 	p = malloc(n);
 	if(p == nil)
-		return buf;
+		return a;
 	m = dec32(p, n, buf, n);
-	if(m == -1)
-		a = buf;
-	else
-		betomp(p, m, b);
+	betomp(p, m, b);
 	free(p);
 	return a;
 }
@@ -201,10 +197,9 @@
 	if(e == a)
 		return nil;
 
-	mpnorm(b);
-	b->sign = sign;
 	if(pp != nil)
 		*pp = e;
 
-	return b;
+	b->sign = sign;
+	return mpnorm(b);
 }