diff --git a/src/algorithms/MontgomeryMultiplication.test.ts b/src/algorithms/MontgomeryMultiplication.test.ts new file mode 100644 index 0000000..2468208 --- /dev/null +++ b/src/algorithms/MontgomeryMultiplication.test.ts @@ -0,0 +1,41 @@ +import { bigIntToUint32Array, uint32ArrayToBigInt } from "../gpu/utils"; +import { FIELD_MODULUS } from "../params/BLS12_377Constants"; +import { Montgomery, Montgomery_K, Montgomery_MultiPrecisionREDC } from "./MontgomeryMultiplication"; + +describe('Montgomery', () => { + it.each([ + [BigInt(8), BigInt(7), BigInt(56)], + [FIELD_MODULUS, BigInt(0), BigInt(0)], + [FIELD_MODULUS, BigInt(2), BigInt(0)], + [FIELD_MODULUS, FIELD_MODULUS, BigInt(0)], + [FIELD_MODULUS + BigInt(2), BigInt(2), BigInt(4)], + ])('properly does modulus multiplication', (a: bigint, b: bigint, expected: bigint) => { + const result = Montgomery(a, b); + expect(result).toEqual(expected); + }); + + it.each([ + [BigInt(8), BigInt(7), BigInt(56)], + [FIELD_MODULUS, BigInt(0), BigInt(0)], + [FIELD_MODULUS, BigInt(2), BigInt(0)], + [FIELD_MODULUS, FIELD_MODULUS, BigInt(0)], + [FIELD_MODULUS + BigInt(2), BigInt(2), BigInt(4)], + ])('properly does modulus multiplication', (a: bigint, b: bigint, expected: bigint) => { + const result = Montgomery_K(a, b); + expect(result).toEqual(expected); + }); + + it.each([ + [BigInt(8), BigInt(7), BigInt(56)], + [FIELD_MODULUS, BigInt(0), BigInt(0)], + [FIELD_MODULUS, BigInt(2), BigInt(0)], + [FIELD_MODULUS, FIELD_MODULUS, BigInt(0)], + [FIELD_MODULUS + BigInt(2), BigInt(2), BigInt(4)], + ])('properly does modulus multiplication', (a: bigint, b: bigint, expected: bigint) => { + const mul = a * b; + const t = bigIntToUint32Array(mul, false); + const result = Montgomery_MultiPrecisionREDC(t); + const res = uint32ArrayToBigInt(result, false); + expect(res).toEqual(expected); + }); +}); \ No newline at end of file diff --git a/src/algorithms/MontgomeryMultiplication.ts b/src/algorithms/MontgomeryMultiplication.ts new file mode 100644 index 0000000..f943296 --- /dev/null +++ b/src/algorithms/MontgomeryMultiplication.ts @@ -0,0 +1,303 @@ +import { uint32ArrayToBigInt } from "../gpu/utils"; +import { + FIELD_MODULUS as N, + MONTGOMERY_RADIX as R, + MONTGOMERY_RINVERSE as R_PRIME, + FIELD_MODULUS_NPRIME as N_PRIME +} from "../params/BLS12_377Constants"; + +import { Kochanski } from "./Kochanski"; + +const Rsq = BigInt('508595941311779472113692600146818027278633330499214071737745792929336755579'); // R^2 mod N. Used to transform out of Mont form. + +/** + * Computes a x b mod AleoFieldModulus using Montgomery multiplication. + * Note that converting to and from Montgomery form are the most expensive operations as they require modulo AleoFieldModulus. + * + * @param a A bigint + * @param b Another bigint + * @returns a * b mod AleoFieldModulus + */ +export const Montgomery = (a: bigint, b: bigint) => { + // Convert to montgomery form. Can rely on Montgomery_redc for this utilizing the fact that mont_redc(x * (R^2 mod N)) = x * R mod N. + // let a_mont = (a * R) % N; + // let b_mont = (b * R) % N; + let a_mont = Montgomery_redc(a * Rsq); // Note that if a * Rsq is a 512 bit number. u256 wrapping makes this incorrect. + let b_mont = Montgomery_redc(b * Rsq); + + // Do montgomery multiplication: + let t = a_mont * b_mont; + let u = Montgomery_redc(t); + + // Convert out back out of montgomery form. Takes one more montgomery reduction. + return Montgomery_redc(u); +}; + +export const Montgomery_K = (a: bigint, b: bigint) => { + // convert to montgomery form. Can rely on Montgomery_redc for this utilizing the fact that mont_redc(x * (R^2 mod N)) = x * R mod N. + let a_mont = Kochanski(a, R, N); + let b_mont = Kochanski(b, R, N); + + // Do montgomery multiplication: + let t = a_mont * b_mont; + let u = Montgomery_redc(t); + + // Convert out back out of montgomery form. Also slow. Could use Kochanski for this step. + return Montgomery_redc(u); +}; + +/** + * https://en.wikipedia.org/wiki/Montgomery_modular_multiplication + * FUTURE: See Montgomery arithmetic on multiprecision integers section for a reduction utilizing limbs. + * Note that this would need a rechoosing of R from R = 2^253 to R = 2^256. + * + * REDC is a way to efficiently do tR' mod N, where R' is the inverse of R mod N and t is the product of 2 montgomery numbers. + * + * From Wikipedia: + * function REDC is + input: Integers R and N with gcd(R, N) = 1, + Integer N′ in [0, R − 1] such that NN′ ≡ −1 mod R, + Integer T in the range [0, RN − 1]. + output: Integer S in the range [0, N − 1] such that S ≡ TR−1 mod N + + m ← ((T mod R)N′) mod R + t ← (T + mN) / R + if t ≥ N then + return t − N + else + return t + end if +end function + */ + +const Montgomery_redc = (t: bigint) => { + const m = (t * N_PRIME) % R; // Note than this mod R is fast as you can right shift by 256 bits to get the remainder + const newT = (t + (m * N)) >> BigInt(256); // divide by R + if (newT >= N) { + return newT - N; + } + else { + return newT; + } +} + +/** + * Given 2 number a and b, returns the two numbers x and y such that ax + by = gcd(a, b). (Bezout's identity) + * + * @param a + * @param b + * @returns [gcd, x, y] + */ +function extendedGCD(a: bigint, b: bigint): [bigint, bigint, bigint] { + if (a === BigInt(0)) { + return [b, BigInt(0), BigInt(1)]; + } else { + const [g, x, y] = extendedGCD(b % a, a); + return [g, y - (b / a) * x, x]; + } +} + +// Used for finding R' such that RR' mod N = 1 +function modInv(a: bigint, m: bigint): bigint { + const [g, x, _] = extendedGCD(a, m); + if (g !== BigInt(1)) { + throw new Error('Modular inverse does not exist'); + } else { + return (x % m + m) % m; // Ensure the result is positive + } +} + +const FIELD_MODULUS_AS_U32_ARRAY: Uint32Array = new Uint32Array([ + 313222494, + 2586617174, + 1622428958, + 1547153409, + 1504343806, + 3489660929, + 168919040, + 1 +]); + +/** + * Wikipedia version. Can't get working :/ + * + * function MultiPrecisionREDC is + Input: Integer N with gcd(B, N) = 1, stored as an array of p words, + Integer R = Br, --thus, r = logB R + Integer N′ in [0, B − 1] such that NN′ ≡ −1 (mod B), + Integer T in the range 0 ≤ T < RN, stored as an array of r + p words. + + Output: Integer S in [0, N − 1] such that TR−1 ≡ S (mod N), stored as an array of p words. + + Set T[r + p] = 0 (extra carry word) + for 0 ≤ i < r do + --loop1- Make T divisible by Bi+1 + + c ← 0 + m ← T[i] ⋅ N′ mod B + for 0 ≤ j < p do + --loop2- Add the low word of m ⋅ N[j] and the carry from earlier, and find the new carry + + x ← T[i + j] + m ⋅ N[j] + c + T[i + j] ← x mod B + c ← ⌊x / B⌋ + end for + for p ≤ j ≤ r + p − i do + --loop3- Continue carrying + + x ← T[i + j] + c + T[i + j] ← x mod B + c ← ⌊x / B⌋ + end for + end for + + for 0 ≤ i < p do + S[i] ← T[i + r] + end for + + if S ≥ N then + return S − N + else + return S + end if +end function + */ +export const Montgomery_MultiPrecisionREDC = (t: Uint32Array): Uint32Array => { + console.log('t array: ', t); + console.log('t: ', uint32ArrayToBigInt(t, false)); + const r = 8; // We use 8 limbs for our u256s + const p = FIELD_MODULUS_AS_U32_ARRAY.length; + const B = 2**32; + const nBPrime = 4294967295; // N*nBPrime mod B = -1 + const NlittleEndian = FIELD_MODULUS_AS_U32_ARRAY.slice().reverse(); + let s = new Uint32Array(p).fill(0); + // let newT = new Uint32Array(p + r).fill(0); + t[r + p] = 0; // Extra carry word + console.log('t array: ', t); + for (let i = 0; i < r; i++) { + let c = 0; + let m = (t[i] * nBPrime) % B; // m = T[i] * N' mod B + for (let j = 0; j < p; j++) { + let x = BigInt(t[i + j]) + BigInt(m * NlittleEndian[j]) + BigInt(c); + t[i + j] = Number(x % BigInt(B)); + c = Math.floor(Number(x / BigInt(B))); + } + for (let j = p; j <= r + p - i; j++) { + let x = BigInt(t[i + j]) + BigInt(c); + t[i + j] = Number(x % BigInt(B)); + c = Math.floor(Number(x / BigInt(B))); + } + } + + for (let i = 0; i < p; i++) { + s[i] = t[i + r]; + } + console.log('t array: ', t); + console.log('s array: ', s); + console.log('S:', uint32ArrayToBigInt(s, false)); + + let compare = compareUint32Arrays(s, NlittleEndian); + if (compare || compare == 0) { + return subtractUint32Arrays(s, NlittleEndian); + } + else { + return s; + } +} + +/** + * Algorithm 14.36 from Handbook of Applied Cryptography + * + * + * INPUT: integers m = (mn−1 ··· m1m0)b, x = (xn−1 ··· x1x0)b, y = (yn−1 ··· y1y0)b + * with 0 ≤ x, y < m, R = bn with gcd(m, b)=1, and m0 = −m−1 mod b. + * OUTPUT: xyR−1 mod m. + * 1. A←0. (Notation: A = (anan−1 ··· a1a0)b.) + * 2. For i from 0 to (n − 1) do the following: + * 2.1 ui←(a0 + xiy0)m0 mod b. + * 2.2 A←(A + xiy + uim)/b. + * 3. If A ≥ m the + */ +export const MontgomeryMultiplication = (x: Uint32Array, y: Uint32Array): Uint32Array => { + const m = FIELD_MODULUS_AS_U32_ARRAY; + const m0 = 4294967295; // m * m0 mod b = -1 + let A = new Uint32Array(1).fill(0); +} + +function addUint32Arrays(a: Uint32Array, b: Uint32Array, isBigEndian = false): Uint32Array { + const maxLength = Math.max(a.length, b.length); + const result = new Uint32Array(maxLength); + let carry = 0; + + const loopStart = isBigEndian ? maxLength - 1 : 0; + const loopEnd = isBigEndian ? -1 : maxLength; + const loopStep = isBigEndian ? -1 : 1; + + for (let i = loopStart; isBigEndian ? i > loopEnd : i < loopEnd; i += loopStep) { + const sum = (a[i] || 0) + (b[i] || 0) + carry; + result[i] = sum >>> 0; + carry = sum < a[i] || sum < b[i] ? 1 : 0; + } + + if (carry > 0) { + const extendedResult = new Uint32Array(maxLength + 1); + extendedResult.set(result); + if (isBigEndian) { + extendedResult[maxLength] = carry; + } else { + extendedResult[0] = carry; + } + return extendedResult; + } + + return result; +} + +function subtractUint32Arrays(a: Uint32Array, b: Uint32Array, isBigEndian = false): Uint32Array { + const maxLength = Math.max(a.length, b.length); + const result = new Uint32Array(maxLength); + let borrow = 0; + + const loopStart = isBigEndian ? maxLength - 1 : 0; + const loopEnd = isBigEndian ? -1 : maxLength; + const loopStep = isBigEndian ? -1 : 1; + + for (let i = loopStart; isBigEndian ? i > loopEnd : i < loopEnd; i += loopStep) { + const diff = (a[i] || 0) - (b[i] || 0) - borrow; + result[i] = diff >>> 0; + borrow = diff < 0 ? 1 : 0; + } + + let lastIndex = isBigEndian ? maxLength - 1 : 0; + while (lastIndex >= 0 && result[lastIndex] === 0) { + lastIndex -= loopStep; + } + + return result.subarray(isBigEndian ? lastIndex + 1 : 0, isBigEndian ? undefined : lastIndex + 1); +} + +// Compares two Uint32Arrays of arbitrary length. +// Returns -1 if a < b, 1 if a > b, and 0 if a == b. +// Takes into account the specified endianness. +function compareUint32Arrays(a: Uint32Array, b: Uint32Array, isBigEndian = false): number { + const maxLength = Math.max(a.length, b.length); + + const loopStart = isBigEndian ? 0 : maxLength - 1; + const loopEnd = isBigEndian ? maxLength : -1; + const loopStep = isBigEndian ? 1 : -1; + + for (let i = loopStart; isBigEndian ? i < loopEnd : i > loopEnd; i += loopStep) { + const aValue = a[i] || 0; + const bValue = b[i] || 0; + + if (aValue < bValue) { + return 1; + } + if (aValue > bValue) { + return -1; + } + } + + return 0; +} + diff --git a/src/algorithms/U32.ts b/src/algorithms/U32.ts new file mode 100644 index 0000000..f046682 --- /dev/null +++ b/src/algorithms/U32.ts @@ -0,0 +1,190 @@ +export class U32 { + private value: number; + + constructor(input: number | BigInt) { + if (typeof input === 'number') { + // Truncate the decimal part + input = Math.trunc(input); + // Wrap around to fit into 32 bits + this.value = input >>> 0; + } else if (typeof input === 'bigint') { + // Wrap around to fit into 32 bits + this.value = Number(input % BigInt(2 ** 32)); + } else { + throw new Error('Invalid type for U32'); + } + } + + // Getter to retrieve the value + get(): number { + return this.value; + } + + // Setter to update the value + set(newValue: number | BigInt): void { + if (typeof newValue === 'number') { + // Truncate the decimal part + newValue = Math.trunc(newValue); + // Wrap around to fit into 32 bits + this.value = newValue >>> 0; + } else if (typeof newValue === 'bigint') { + // Wrap around to fit into 32 bits + this.value = Number(newValue % BigInt(2 ** 32)); + } else { + throw new Error('Invalid type for U32'); + } + } + + // Addition + add(other: U32): U32 { + return new U32(this.get() + other.get()); + } + + // Subtraction + sub(other: U32): U32 { + return new U32(this.get() - other.get()); + } + + // Multiplication + mul(other: U32): U32 { + return new U32(this.get() * other.get()); + } + + // Division + div(other: U32): U32 { + if (other.get() === 0) { + throw new Error('Division by zero'); + } + return new U32(Math.floor(this.get() / other.get())); + } +} + +export function calculateMultCarry(a: number, b: number): number { + // Break a and b into 16-bit halves + const aLow = a & 0xFFFF; + const aHigh = a >>> 16; + const bLow = b & 0xFFFF; + const bHigh = b >>> 16; + + // Multiply the low 16 bits of a and b + const lowProduct = aLow * bLow; + + // Multiply the high 16 bits of a and b + const highProduct = aHigh * bHigh; + + // Calculate the carry + const carry = highProduct + ((lowProduct >>> 16) + aLow * bHigh + aHigh * bLow >>> 16); + + return carry; +} + +export function calculateAdditionCarry(a: number, b: number): number { + const sum = a + b; + + // Check for carry + const carry = (sum < a || sum < b) ? 1 : 0; + + return carry; +} + +type U32Array = U32[]; + +// Multiplies two U32Arrays of any length. +// Returns a U32Array whose length is the sum of the lengths of the input arrays. +function multiplyArbitraryLength(a: U32Array, b: U32Array): U32Array { + const resultLength = a.length + b.length; + const result: U32Array = new Array(resultLength).fill(0); + + for (let i = 0; i < a.length; i++) { + let carry = new U32(0); + for (let j = 0; j < b.length; j++) { + const temp = a[i].mul(b[j]).add(result[i+j]).add(carry); + result[i + j] = temp; + carry = temp.div(new U32(0x100000000)); // Extract carry + } + result[i + b.length] =result[i + b.length].add(carry); + } + + return result; +} + +export function compareU32Arrays(a: U32[], b: U32[], bigEndian: boolean = false): number { + // Make sure both arrays have the same length by padding with zeros if needed + const maxLength = Math.max(a.length, b.length); + const aPadded = Array(maxLength - a.length).fill(0).concat(a); + const bPadded = Array(maxLength - b.length).fill(0).concat(b); + + const iterStart = bigEndian ? 0 : maxLength - 1; + const iterEnd = bigEndian ? maxLength - 1 : 0; + + // Compare arrays element by element, starting from the most significant + for (let i of rangeBetween(iterStart, iterEnd)) { + if (aPadded[i] > bPadded[i]) { + return 1; + } else if (aPadded[i] < bPadded[i]) { + return -1; + } + } + + // If we reach this point, the arrays are equal + return 0; +} + +// Adds two U32Arrays of arbitrary length. +// Returns a new U32Array representing the sum. +export function addU32Arrays(a: U32[], b: U32[], bigEndian = false): U32Array { + const maxLength = Math.max(a.length, b.length); + const result: U32Array = new Array(maxLength).fill(0); + let carry = 0; + + const iterStart = bigEndian ? maxLength - 1 : 0; + const iterEnd = bigEndian ? 0 : maxLength - 1; + + for (let i of rangeBetween(iterStart, iterEnd)) { + const sum = (a[i].get() || 0) + (b[i].get() || 0) + carry; + result[i] = new U32(sum >>> 0); // Truncate to 32 bits + carry = sum < a[i].get() || sum < b[i].get() ? 1 : 0; // Check for carry + } + + if (carry > 0) { + result.push(new U32(carry)); + } + + return result; +} + +// Subtracts two U32Arrays of arbitrary length. +// Returns a new U32Array representing the difference. +export function subtractU32Arrays(a: U32Array, b: U32Array, bigEndian = false): U32Array { + const maxLength = Math.max(a.length, b.length); + const result: U32Array = new Array(maxLength).fill(0); + let borrow = 0; + + const iterStart = bigEndian ? maxLength - 1 : 0; + const iterEnd = bigEndian ? 0 : maxLength - 1; + + for (let i of rangeBetween(iterStart, iterEnd)) { + const diff = (a[i].get() || 0) - (b[i].get() || 0) - borrow; + result[i] = new U32(diff >>> 0); // Truncate to 32 bits + borrow = diff < 0 ? 1 : 0; // Check for borrow + } + + // Remove leading zeros from the result + while (result.length > 1 && result[result.length - 1].get() === 0) { + result.pop(); + } + + return result; +} + +// Inclusive on both ends +function rangeBetween(a: number, b: number): number[] { + const result: number[] = []; + const step = a <= b ? 1 : -1; + + for (let i = a; step > 0 ? i <= b : i >= b; i += step) { + result.push(i); + } + + return result; +} \ No newline at end of file diff --git a/src/gpu/utils.ts b/src/gpu/utils.ts index 6fcc603..98715fe 100644 --- a/src/gpu/utils.ts +++ b/src/gpu/utils.ts @@ -44,6 +44,40 @@ export const u32ArrayToBigInts = (u32Array: Uint32Array): bigint[] => { return bigInts; } +// Not u256 specific +export function bigIntToUint32Array(value: bigint, isBigEndian: boolean): Uint32Array { + // Calculate the required size of the Uint32Array + const size = Math.ceil(Number(value.toString(2).length) / 32); + const result = new Uint32Array(size); + + const loopStart = isBigEndian ? 0 : size - 1; + const loopEnd = isBigEndian ? size : -1; + const loopStep = isBigEndian ? 1 : -1; + + let tempValue = value; + for (let i = loopStart; isBigEndian ? i < loopEnd : i > loopEnd; i += loopStep) { + result[i] = Number(tempValue & BigInt(0xFFFFFFFF)); + tempValue >>= BigInt(32); + } + + return result; +} + +// Not u256 specific +export function uint32ArrayToBigInt(arr: Uint32Array, isBigEndian: boolean): BigInt { + let result = BigInt(0); + + const loopStart = isBigEndian ? 0 : arr.length - 1; + const loopEnd = isBigEndian ? arr.length : -1; + const loopStep = isBigEndian ? 1 : -1; + + for (let i = loopStart; isBigEndian ? i < loopEnd : i > loopEnd; i += loopStep) { + result = (result << BigInt(32)) | BigInt(arr[i]); + } + + return result; +} + export const generateRandomFields = (inputSize: number): bigint[] => { const randomBigInts = []; for (let i = 0; i < inputSize; i++) { diff --git a/src/params/BLS12_377Constants.ts b/src/params/BLS12_377Constants.ts index 5bbde5e..e9e2847 100644 --- a/src/params/BLS12_377Constants.ts +++ b/src/params/BLS12_377Constants.ts @@ -1,5 +1,8 @@ // base field modulus https://docs.rs/ark-ed-on-bls12-377/latest/ark_ed_on_bls12_377/ -export const FIELD_MODULUS = BigInt('8444461749428370424248824938781546531375899335154063827935233455917409239041'); +export const FIELD_MODULUS = BigInt('8444461749428370424248824938781546531375899335154063827935233455917409239041'); // N -> Chosen by aleo +export const FIELD_MODULUS_NPRIME = BigInt('47752251086953357377073236701509605140872345086634869599321669320666611974143'); // N' -> NN' mod R = -1. RR' - NN' = 1 by extend Euclidean algorithm +export const MONTGOMERY_RADIX = BigInt('115792089237316195423570985008687907853269984665640564039457584007913129639936'); // R -> 2 ^ 256 +export const MONTGOMERY_RINVERSE = BigInt('3482466379256973933331601287759811764685972354380176549708408303012390300674'); // R' -> RR' mod N = 1 export const EDWARDS_A = BigInt('8444461749428370424248824938781546531375899335154063827935233455917409239040'); export const EDWARDS_D = BigInt('3021'); export const SUBGROUP_CHARACTERISTIC = BigInt('2111115437357092606062206234695386632838870926408408195193685246394721360383'); \ No newline at end of file