--- /n/sources/plan9/sys/include/mp.h	Thu Jun 27 02:14:12 2024
+++ /sys/include/mp.h	Thu Jun 27 02:17:28 2024
@@ -75,6 +75,11 @@
 void	mpexp(mpint *b, mpint *e, mpint *m, mpint *res);	/* res = b**e mod m */
 void	mpmod(mpint *b, mpint *m, mpint *remainder);	/* remainder = b mod m */
 
+/* modular arithmetic, time invariant when 0≤b1≤m-1 and 0≤b2≤m-1 */
+void	mpmodadd(mpint *b1, mpint *b2, mpint *m, mpint *sum);	/* sum = b1+b2 % m */
+void	mpmodsub(mpint *b1, mpint *b2, mpint *m, mpint *diff);	/* diff = b1-b2 % m */
+void	mpmodmul(mpint *b1, mpint *b2, mpint *m, mpint *prod);	/* prod = b1*b2 % m */
+
 /* quotient = dividend/divisor, remainder = dividend % divisor */
 void	mpdiv(mpint *dividend, mpint *divisor,  mpint *quotient, mpint *remainder);
 
--- /n/sources/plan9/sys/include/ape/mp.h	Thu Jun 27 02:14:12 2024
+++ /sys/include/ape/mp.h	Thu Jun 27 02:21:36 2024
@@ -84,6 +84,11 @@
 void	mpexp(mpint *b, mpint *e, mpint *m, mpint *res);	/* res = b**e mod m */
 void	mpmod(mpint *b, mpint *m, mpint *remainder);	/* remainder = b mod m */
 
+/* modular arithmetic, time invariant when 0≤b1≤m-1 and 0≤b2≤m-1 */
+void	mpmodadd(mpint *b1, mpint *b2, mpint *m, mpint *sum);	/* sum = b1+b2 % m */
+void	mpmodsub(mpint *b1, mpint *b2, mpint *m, mpint *diff);	/* diff = b1-b2 % m */
+void	mpmodmul(mpint *b1, mpint *b2, mpint *m, mpint *prod);	/* prod = b1*b2 % m */
+
 /* quotient = dividend/divisor, remainder = dividend % divisor */
 void	mpdiv(mpint *dividend, mpint *divisor,  mpint *quotient, mpint *remainder);
 
--- /n/sources/plan9/sys/man/2/mp	Thu Jun 27 02:14:12 2024
+++ /sys/man/2/mp	Thu Jun 27 02:17:28 2024
@@ -1,6 +1,6 @@
 .TH MP 2
 .SH NAME
-mpsetminbits, mpnew, mpfree, mpbits, mpnorm, mpcopy, mpassign, mprand, mpnrand, strtomp, mpfmt,mptoa, betomp, mptobe, mptober, letomp, mptole, mptolel, mptoui, uitomp, mptoi, itomp, uvtomp, mptouv, vtomp, mptov, mpdigdiv, mpadd, mpsub, mpleft, mpright, mpmul, mpexp, mpmod, mpdiv, mpcmp¸ mpsel, mpextendedgcd, mpinvert, mpsignif, mplowbits0, mpvecdigmuladd, mpvecdigmulsub, mpvecadd, mpvecsub, mpveccmp, mpvecmul, mpmagcmp, mpmagadd, mpmagsub, crtpre, crtin, crtout, crtprefree, crtresfree \- extended precision arithmetic
+mpsetminbits, mpnew, mpfree, mpbits, mpnorm, mpcopy, mpassign, mprand, mpnrand, strtomp, mpfmt,mptoa, betomp, mptobe, mptober, letomp, mptole, mptolel, mptoui, uitomp, mptoi, itomp, uvtomp, mptouv, vtomp, mptov, mpdigdiv, mpadd, mpsub, mpleft, mpright, mpmul, mpexp, mpmod, mpmodadd, mpmodsub, mpmodmul, mpdiv, mpcmp¸ mpsel, mpextendedgcd, mpinvert, mpsignif, mplowbits0, mpvecdigmuladd, mpvecdigmulsub, mpvecadd, mpvecsub, mpveccmp, mpvecmul, mpmagcmp, mpmagadd, mpmagsub, crtpre, crtin, crtout, crtprefree, crtresfree \- extended precision arithmetic
 .SH SYNOPSIS
 .B #include <u.h>
 .br
@@ -120,6 +120,15 @@
 .B
 	mpint *remainder)
 .PP
+void	mpmodadd(mpint *b1, mpint *b2, mpint *m, mpint *sum)
+.PP
+.B
+void	mpmodsub(mpint *b1, mpint *b2, mpint *m, mpint *diff)
+.PP
+.B
+void	mpmodmul(mpint *b1, mpint *b2, mpint *m, mpint *prod)
+.PP
+.B
 .B
 int	mpcmp(mpint *b1, mpint *b2)
 .PP
@@ -525,6 +534,19 @@
 .I b2
 is assigned to
 .IR res .
+.PD
+.PP
+Modular arithmetic:
+.TF mpmodmul_
+.TP
+.I mpmodadd
+.BR "sum = b1+b2 mod m" .
+.TP
+.I mpmodsub
+.BR "diff = b1-b2 mod m" .
+.TP
+.I mpmodmul
+.BR "prod = b1*b2 mod m" .
 .PD
 .PP
 .I Mpextendedgcd
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 02:04:34 2024
+++ /sys/src/ape/lib/mp/port/mkfile	Thu Jun 27 02:18:17 2024
@@ -30,6 +30,7 @@
 	mpdiv\
 	mpexp\
 	mpmod\
+	mpmodop\
 	mpextendedgcd\
 	mpinvert\
 	mprand\
diff -Nru /n/sources/plan9/sys/src/libmp/port/mkfile /sys/src/libmp/port/mkfile
--- /n/sources/plan9/sys/src/libmp/port/mkfile	Thu Jun 27 02:02:06 2024
+++ /sys/src/libmp/port/mkfile	Thu Jun 27 02:17:30 2024
@@ -28,6 +28,7 @@
 	mpdiv\
 	mpexp\
 	mpmod\
+	mpmodop\
 	mpextendedgcd\
 	mpinvert\
 	mprand\
diff -Nru /n/sources/plan9/sys/src/libmp/port/mpmodop.c /sys/src/libmp/port/mpmodop.c
--- /n/sources/plan9/sys/src/libmp/port/mpmodop.c	Wed Dec 31 19:00:00 1969
+++ /sys/src/libmp/port/mpmodop.c	Thu Jun 27 02:17:31 2024
@@ -0,0 +1,96 @@
+#include <u.h>
+#include <libc.h>
+#include <mp.h>
+
+/* operands need to have m->top+1 digits of space and satisfy 0 ≤ a ≤ m-1 */
+static mpint*
+modarg(mpint *a, mpint *m)
+{
+	if(a->size <= m->top || a->sign < 0 || mpmagcmp(a, m) >= 0){
+		a = mpcopy(a);
+		mpmod(a, m, a);
+		mpbits(a, Dbits*(m->top+1));
+		a->top = m->top;
+	} else if(a->top < m->top){
+		memset(&a->p[a->top], 0, (m->top - a->top)*Dbytes);
+	}
+	return a;
+}
+
+void
+mpmodadd(mpint *b1, mpint *b2, mpint *m, mpint *sum)
+{
+	mpint *a, *b;
+	mpdigit d;
+	int i, j;
+
+	a = modarg(b1, m);
+	b = modarg(b2, m);
+
+	sum->flags |= (a->flags | b->flags) & MPtimesafe;
+	mpbits(sum, Dbits*2*(m->top+1));
+
+	mpvecadd(a->p, m->top, b->p, m->top, sum->p);
+	mpvecsub(sum->p, m->top+1, m->p, m->top, sum->p+m->top+1);
+
+	d = sum->p[2*m->top+1];
+	for(i = 0, j = m->top+1; i < m->top; i++, j++)
+		sum->p[i] = (sum->p[i] & d) | (sum->p[j] & ~d);
+
+	sum->top = m->top;
+	sum->sign = 1;
+	mpnorm(sum);
+
+	if(a != b1)
+		mpfree(a);
+	if(b != b2)
+		mpfree(b);
+}
+
+void
+mpmodsub(mpint *b1, mpint *b2, mpint *m, mpint *diff)
+{
+	mpint *a, *b;
+	mpdigit d;
+	int i, j;
+
+	a = modarg(b1, m);
+	b = modarg(b2, m);
+
+	diff->flags |= (a->flags | b->flags) & MPtimesafe;
+	mpbits(diff, Dbits*2*(m->top+1));
+
+	a->p[m->top] = 0;
+	mpvecsub(a->p, m->top+1, b->p, m->top, diff->p);
+	mpvecadd(diff->p, m->top, m->p, m->top, diff->p+m->top+1);
+
+	d = ~diff->p[m->top];
+	for(i = 0, j = m->top+1; i < m->top; i++, j++)
+		diff->p[i] = (diff->p[i] & d) | (diff->p[j] & ~d);
+
+	diff->top = m->top;
+	diff->sign = 1;
+	mpnorm(diff);
+
+	if(a != b1)
+		mpfree(a);
+	if(b != b2)
+		mpfree(b);
+}
+
+void
+mpmodmul(mpint *b1, mpint *b2, mpint *m, mpint *prod)
+{
+	mpint *a, *b;
+
+	a = modarg(b1, m);
+	b = modarg(b2, m);
+
+	mpmul(a, b, prod);
+	mpmod(prod, m, prod);
+
+	if(a != b1)
+		mpfree(a);
+	if(b != b2)
+		mpfree(b);
+}