Add ResidueClass(k, n), an element of ℤ/nℤ, and list ℤ/nℤ through it

CHANGELOG.md

@@ -71,6 +71,33 @@
   does for `sin(x) = 2`. So `x·eˣ = 3` still solves to `W(3)` only, and
   `x·eˣ = -1` has no root.
 
+### New Features
+
+- **`ResidueClass(k, n)` is an element of ℤ/nℤ**, and `QuotientRing(Integers, n)`
+  lists them (#399, contributed by [enumeratio](https://github.com/enumeratio)).
+  The canonical form reduces `k` into `0…n−1`, so `ResidueClass(7, 5)` is
+  `ResidueClass(2, 5)` and `ResidueClass(-1, 7)` is `ResidueClass(6, 7)`, and
+  `ResidueClass(7, 5) == ResidueClass(2, 5)` is `True`. Sums, differences,
+  products and integer powers of classes are classes:
+  `ResidueClass(5, 7) + ResidueClass(4, 7)` is `ResidueClass(2, 7)` and
+  `ResidueClass(3, 7)^6` is `ResidueClass(1, 7)`. An integer next to a class is
+  read in its ring: `ResidueClass(5, 7) + 3` is `ResidueClass(1, 7)`, and
+  `ResidueClass(1/3, 7)` is `ResidueClass(5, 7)`. An inverse, a division and a
+  negative power need `gcd(k, n) = 1`: `1/ResidueClass(3, 7)` is
+  `ResidueClass(5, 7)`, and `1/ResidueClass(2, 4)` stays unevaluated. Classes of
+  two moduli meet in ℤ/gcd(m, n): `ResidueClass(2, 4) + ResidueClass(1, 6)` is
+  `ResidueClass(1, 2)`. The class has type `value`, and the element type of
+  `QuotientRing(Integers, n)` is `value`, not `unknown`. ℤ/5ℤ lists
+  `ResidueClass(0, 5)` … `ResidueClass(4, 5)` (`ListFrom`, `Union`, `Tally`…
+  walk them), and `Element(ResidueClass(7, 5), ℤ/5ℤ)` is `True`; the
+  membership of a bare integer, `Element(7, ℤ/5ℤ)`, is still not decided. In
+  LaTeX the class is `\overline{k}_{n}` for integer literals `k` and `n ≥ 1`:
+  `\overline{7}_{5}` was `Subscript(Conjugate(7), 5)` and is now
+  `ResidueClass(2, 5)`, while `\overline{7}`, `\overline{z}_1` and the
+  repeating decimal `0.\overline{3}` are unchanged. `Mod` stays the remainder,
+  classes are not ordered, a modulus past 2^53 is exact, and compiled code does
+  not lower a class (it falls back, as for any operator without a lowering).
+
 ### Improvements
 
 - **More accurate `erf`, `erfc` and `erfi` in doubles.** The machine kernels
@@ -125,14 +152,14 @@
   same element of ℤ/5ℤ. For a positive integer literal `n`, the collection
   now has the count `n`, is finite and is not empty, and `Count` evaluates:
   `Count(\mathbb{Z}/5\mathbb{Z})` is `5`. A symbolic, zero or negative
-  modulus, or a base other than `Integers`, stays inert. The type is
-  `set<unknown>`, as for an `Adjoin` adjunct that the engine cannot type. The
-  engine has no value for a residue class, so the classes are not listed and
-  membership is not decided.
+  modulus, or a base other than `Integers`, stays inert. Its elements are
+  `ResidueClass` values (see New Features), so the type is `set<value>` for
+  the base `Integers`, and `set<unknown>`, as for an `Adjoin` adjunct that the
+  engine cannot type, for another base.
 
 - **A finite collection whose elements cannot be computed no longer reads as
   empty.** `Linspace(a, 1, 3)` with a symbolic `a` has the count 3 but no
-  elements that can be computed, and so does ℤ/5ℤ. Many operators checked
+  elements that can be computed. Many operators checked
   only that a collection was finite and then walked its elements, so they
   read it as empty and gave a wrong answer with no diagnostic:
   `Union(Linspace(a, 1, 3), {1})` was `Set(1)`, `Unique(…)` was `[]`,

src/compute-engine/boxed-expression/residue-class.ts

@@ -0,0 +1,216 @@
+/**
+ * Residue classes: the elements of ℤ/nℤ, `ResidueClass(k, n)`.
+ *
+ * A class is held as a pair of integers `k` in [0, n) and `n ≥ 1`, in
+ * `bigint`, so a modulus past 2^53 stays exact. Every function here works on
+ * that pair; the arithmetic handlers of `Add`, `Multiply`, `Divide`, `Power`
+ * and `Negate` call the `residue*` functions on their evaluated operands.
+ *
+ * Two classes with different moduli meet in ℤ/gcd(m, n), the largest ring
+ * that both reduce onto: `ResidueClass(2, 4) + ResidueClass(1, 6)` is
+ * `ResidueClass(1, 2)`. An integer, or a rational whose denominator is a unit,
+ * is read in the ring that the classes meet in.
+ */
+
+import type {
+  Expression,
+  IComputeEngine as ComputeEngine,
+  OperandDescriptor,
+} from '../global-types.js';
+import { asBigint, asRational } from './numerics.js';
+import { isFunction, isNumber } from './type-guards.js';
+import { gcd, modularInverse } from '../numerics/numeric-bigint.js';
+import { modPow } from '../numerics/primes.js';
+
+/** A class `k + nℤ`, with `0 ≤ k < n`. */
+export interface Residue {
+  readonly k: bigint;
+  readonly n: bigint;
+}
+
+const mod = (a: bigint, n: bigint): bigint => ((a % n) + n) % n;
+
+/** The class of the rational `num/den` in ℤ/nℤ, if `den` is a unit mod `n`. */
+export function residueOfRational(
+  num: bigint,
+  den: bigint,
+  n: bigint
+): Residue | undefined {
+  if (n < 1n) return undefined;
+  if (den === 1n) return { k: mod(num, n), n };
+  const inverse = modularInverse(den, n);
+  return inverse === null ? undefined : { k: mod(num * inverse, n), n };
+}
+
+/**
+ * The class of an integer or rational literal in ℤ/nℤ. `undefined` for
+ * anything else, a non-unit denominator included (`1/2` mod 4).
+ */
+export function residueOfLiteral(
+  expr: Expression,
+  n: bigint
+): Residue | undefined {
+  const integer = asBigint(expr);
+  if (integer !== null) return { k: mod(integer, n), n };
+  if (!isNumber(expr)) return undefined;
+  const q = asRational(expr);
+  if (q === undefined) return undefined;
+  return residueOfRational(BigInt(q[0]), BigInt(q[1]), n);
+}
+
+/** `ResidueClass(k, n)` with literal integer operands, read as a pair. */
+export function residueOf(expr: Expression | undefined): Residue | undefined {
+  if (!isFunction(expr, 'ResidueClass') || expr.nops !== 2) return undefined;
+  const k = asBigint(expr.op1);
+  const n = asBigint(expr.op2);
+  if (k === null || n === null || n < 1n || k < 0n || k >= n) return undefined;
+  return { k, n };
+}
+
+export function isResidueClass(expr: Expression | undefined): boolean {
+  return residueOf(expr) !== undefined;
+}
+
+/** An operand that is a `ResidueClass(…)` call, as a type handler sees it. */
+export function isResidueClassOperand(d: OperandDescriptor): boolean {
+  const s = d.structureOf?.();
+  return s?.kind === 'application' && s.head === 'ResidueClass';
+}
+
+export function residueExpression(ce: ComputeEngine, x: Residue): Expression {
+  return ce._fn('ResidueClass', [ce.number(x.k), ce.number(x.n)]);
+}
+
+/**
+ * `ResidueClass(k, n)` in canonical form: `k` reduced into [0, n). `undefined`
+ * when `n` is not an integer literal ≥ 1, or `k` is not an integer or a
+ * rational literal that is a unit mod `n`: the call stays as it is.
+ */
+export function residueClassOf(
+  ce: ComputeEngine,
+  k: Expression,
+  n: Expression
+): Expression | undefined {
+  const modulus = asBigint(n);
+  if (modulus === null || modulus < 1n) return undefined;
+  const r = residueOfLiteral(k, modulus);
+  return r === undefined ? undefined : residueExpression(ce, r);
+}
+
+//
+// Arithmetic on pairs
+//
+
+/** The same class read in ℤ/mℤ, for `m` dividing its modulus. */
+const reduce = (x: Residue, m: bigint): Residue => ({ k: mod(x.k, m), n: m });
+
+const meet = (x: Residue, y: Residue): [Residue, Residue] => {
+  const m = gcd(x.n, y.n);
+  return [reduce(x, m), reduce(y, m)];
+};
+
+const add = (x: Residue, y: Residue): Residue => {
+  const [a, b] = meet(x, y);
+  return { k: mod(a.k + b.k, a.n), n: a.n };
+};
+
+const multiply = (x: Residue, y: Residue): Residue => {
+  const [a, b] = meet(x, y);
+  return { k: mod(a.k * b.k, a.n), n: a.n };
+};
+
+const negate = (x: Residue): Residue => ({ k: mod(-x.k, x.n), n: x.n });
+
+/** x⁻¹, or `undefined` when gcd(k, n) ≠ 1. */
+const inverse = (x: Residue): Residue | undefined => {
+  const k = modularInverse(x.k, x.n);
+  return k === null ? undefined : { k, n: x.n };
+};
+
+/** xᵉ; a negative exponent inverts first, so it needs a unit. */
+const power = (x: Residue, e: bigint): Residue | undefined => {
+  const base = e < 0n ? inverse(x) : x;
+  if (base === undefined) return undefined;
+  return { k: modPow(base.k, e < 0n ? -e : e, x.n), n: x.n };
+};
+
+//
+// Operands
+//
+
+/**
+ * Every operand as a class of the ring that the classes among them meet in.
+ * `undefined` when an operand is neither a class nor an integer or rational
+ * literal that reads in that ring (a symbol, a float, `1/2` mod 4).
+ */
+function lift(ops: ReadonlyArray<Expression>): Residue[] | undefined {
+  const classes = ops.map(residueOf).filter((x) => x !== undefined);
+  if (classes.length === 0) return undefined;
+  const ring = classes.map((x) => x.n).reduce(gcd);
+  const values = ops.map((op) =>
+    isResidueClass(op) ? residueOf(op) : residueOfLiteral(op, ring)
+  );
+  return values.every((x) => x !== undefined) ? values : undefined;
+}
+
+function fold(
+  ce: ComputeEngine,
+  ops: ReadonlyArray<Expression>,
+  step: (x: Residue, y: Residue) => Residue | undefined
+): Expression | undefined {
+  const values = lift(ops);
+  if (values === undefined) return undefined;
+  let acc: Residue | undefined = values[0];
+  for (const next of values.slice(1))
+    acc = acc === undefined ? undefined : step(acc, next);
+  return acc === undefined ? undefined : residueExpression(ce, acc);
+}
+
+/** The sum of operands one of which is a class, or `undefined`. */
+export function residueAdd(
+  ce: ComputeEngine,
+  ops: ReadonlyArray<Expression>
+): Expression | undefined {
+  return fold(ce, ops, add);
+}
+
+export function residueMultiply(
+  ce: ComputeEngine,
+  ops: ReadonlyArray<Expression>
+): Expression | undefined {
+  return fold(ce, ops, multiply);
+}
+
+export function residueNegate(
+  ce: ComputeEngine,
+  x: Expression
+): Expression | undefined {
+  const r = residueOf(x);
+  return r === undefined ? undefined : residueExpression(ce, negate(r));
+}
+
+/** `num / den`: `undefined` when `den` is not a unit of the common ring. */
+export function residueDivide(
+  ce: ComputeEngine,
+  num: Expression,
+  den: Expression
+): Expression | undefined {
+  return fold(ce, [num, den], (x, y) => {
+    const [a, b] = meet(x, y);
+    const b1 = inverse(b);
+    return b1 === undefined ? undefined : multiply(a, b1);
+  });
+}
+
+/** `x ^ e` for an integer exponent: `undefined` for a negative one on a non-unit. */
+export function residuePower(
+  ce: ComputeEngine,
+  x: Expression,
+  e: Expression
+): Expression | undefined {
+  const r = residueOf(x);
+  const exponent = asBigint(e);
+  if (r === undefined || exponent === null) return undefined;
+  const p = power(r, exponent);
+  return p === undefined ? undefined : residueExpression(ce, p);
+}

src/compute-engine/boxed-expression/simplify.ts

@@ -8,6 +8,7 @@ import { sameSyntactic } from './compare.js';
 import { holdMap } from './hold.js';
 import { expToTrig } from './exp-to-trig.js';
 import { expand } from './expand.js';
+import { isResidueClass } from './residue-class.js';
 import {
   rebindEscaping,
   hasAssignedVariable,
@@ -113,7 +114,7 @@ function evaluateNumericSubexpressions(expr: Expression): Expression {
     !hasAssignedVariable(expr)
   ) {
     const evaluated = expr.evaluate();
-    if (isNumber(evaluated)) return evaluated;
+    if (isNumber(evaluated) || isResidueClass(evaluated)) return evaluated;
   }
 
   // Constant logarithms are folded to their exact value even when they are
@@ -921,7 +922,7 @@ function simplifyOperands(
         continue;
       }
       const evaluated = x.evaluate();
-      if (isNumber(evaluated)) {
+      if (isNumber(evaluated) || isResidueClass(evaluated)) {
         simplifiedOps.push(evaluated);
         continue;
       }
@@ -1008,6 +1009,22 @@ function simplifyExpression(
     expr = alt;
   }
 
+  // Arithmetic on residue classes folds to a class, as for numbers:
+  // `ResidueClass(5, 7) + ResidueClass(4, 7)` is `ResidueClass(2, 7)`. A
+  // class has no number literal, so the numeric fold of the operands does
+  // not reach it.
+  if (
+    isFunction(expr) &&
+    BASIC_ARITHMETIC.includes(expr.operator) &&
+    expr.ops.some(isResidueClass) &&
+    expr.unknowns.length === 0 &&
+    !hasAssignedVariable(expr)
+  ) {
+    const folded = expr.evaluate();
+    if (isResidueClass(folded))
+      return [...steps, { value: folded, because: 'residue class arithmetic' }];
+  }
+
   // A `NaN` or `Indeterminate` operand at a position whose NaN policy is
   // `propagate` decides the value, as it does for `evaluate()`:
   // `sin(NaN)` simplifies to `NaN` and `sin(Indeterminate)` to

src/compute-engine/latex-syntax/dictionary/definitions-complex.ts

@@ -1,7 +1,24 @@
-import { LatexDictionary, Parser } from '../types.js';
+import { LatexDictionary, Parser, Serializer } from '../types.js';
 import { MathJsonExpression } from '../../../math-json/types.js';
+import { isNumberObject, operand } from '../../../math-json/utils.js';
 import { singleArgSerializer } from './definitions-other.js';
 
+/**
+ * The digits of a non-negative integer literal (`7`, `{num: '123456…'}`), or
+ * `null` for anything else: a symbol, a sum, a negation, a decimal.
+ */
+function integerLiteralDigits(expr: MathJsonExpression | null): string | null {
+  if (typeof expr === 'number')
+    return Number.isSafeInteger(expr) && expr >= 0 ? String(expr) : null;
+  if (isNumberObject(expr)) {
+    const digits = String(expr.num);
+    return /^\d+$/.test(digits) ? digits : null;
+  }
+  return null;
+}
+
+const isPositive = (digits: string): boolean => /[1-9]/.test(digits);
+
 export const DEFINITIONS_COMPLEX: LatexDictionary = [
   {
     name: 'Real',
@@ -56,10 +73,39 @@ export const DEFINITIONS_COMPLEX: LatexDictionary = [
     latexTrigger: ['\\overline'],
     parse: (parser: Parser): MathJsonExpression => {
       const arg = parser.parseGroup();
-      return arg === null ? ['Conjugate'] : ['Conjugate', arg];
+      if (arg === null) return ['Conjugate'];
+
+      // `\overline{k}_{n}` with integer literals `k` and `n ≥ 1` is the
+      // residue class of `k` mod `n`. Anything else keeps the conjugate:
+      // `\overline{7}` and `\overline{z}_1` are not classes, and the
+      // subscript is left for the subscript parselet.
+      if (integerLiteralDigits(arg) !== null) {
+        const start = parser.index;
+        if (parser.match('_')) {
+          const modulus = parser.parseGroup() ?? parser.parseToken();
+          const digits = integerLiteralDigits(modulus);
+          if (digits !== null && isPositive(digits))
+            return ['ResidueClass', arg, modulus!];
+        }
+        parser.index = start;
+      }
+      return ['Conjugate', arg];
     },
     serialize: singleArgSerializer('\\overline'),
   },
+  // `ResidueClass(k, n)` is written `\overline{k}_{n}` when both operands are
+  // integer literals, the only form that reads back as a class; otherwise it
+  // is written as a function call.
+  {
+    name: 'ResidueClass',
+    serialize: (serializer: Serializer, expr: MathJsonExpression): string => {
+      const k = integerLiteralDigits(operand(expr, 1));
+      const n = integerLiteralDigits(operand(expr, 2));
+      if (k === null || n === null || !isPositive(n))
+        return serializer.serializeFunction(expr);
+      return `\\overline{${k}}_{${n}}`;
+    },
+  },
   // Function-style alias: `\operatorname{conj}(z)`, the spelling Desmos writes
   // for the complex conjugate. Without it `conj` lexed as an undeclared symbol:
   // `\operatorname{conj}(z)` was the product `conj·z`, and over a list argument