diff --git a/lib/src/Main.qs b/lib/src/Main.qs index 4b92cd5..ca98b79 100644 --- a/lib/src/Main.qs +++ b/lib/src/Main.qs @@ -1,13 +1,13 @@ -import QuantumArithmetic.Yuan2022.Divide; +import QuantumArithmetic.CG20192.ModExpWindow; import TestUtils.*; import QuantumArithmetic.Utils; // For debugging, run with Ctrl+F5, operation Main() : Unit { let n = 6; - let m = 2; - let a = 37; - let b = 1; - let ans = Test_Divide_Unequal(n, a, m, b, Divide); + let a = 2L; + let x = 42L; + let N = 11L; + let ans = TestModExp(n, a, x, N, ModExpWindow(_, _, _, _, 2, 2)); Message($"ans={ans}"); } \ No newline at end of file diff --git a/lib/src/QuantumArithmetic/CG20192.qs b/lib/src/QuantumArithmetic/CG20192.qs index b9be4e5..4eeac65 100644 --- a/lib/src/QuantumArithmetic/CG20192.qs +++ b/lib/src/QuantumArithmetic/CG20192.qs @@ -18,6 +18,13 @@ import QuantumArithmetic.WindowedArithmeticUtils.Util.ForceMeasureResetBigInt; import QuantumArithmetic.WindowedArithmeticUtils.Util.BitLength; import QuantumArithmetic.WindowedArithmeticUtils.Xor.XorEqualConst; import QuantumArithmetic.WindowedArithmeticUtils.MulAdd_Window.PlusEqualConstTimesLEWindowed; +import QuantumArithmetic.Utils; +import Std.Arrays; +import Std.Convert; +import Std.Math; +import Std.TableLookup.*; +import Std.Arithmetic.RippleCarryCGIncByLE; +import QuantumArithmetic.LYY2021.ModAdd; operation Multiply (nx : Int, ny : Int, result_t : BigInt, classical_factor_x : BigInt, quantum_factor_y: BigInt) : BigInt { @@ -52,4 +59,108 @@ operation MultiplyWindow (nx : Int, ny : Int, result_t : BigInt, return result; } -export Multiply, MultiplyWindow; \ No newline at end of file +internal function Skip2Data( + generator: BigInt, period: BigInt, num_exponents: Int, num_bits: Int +) : Bool[][] { + mutable total = 1L; + let num_entries = 1 <<< num_exponents; + mutable table = [[false, size = num_bits], size = num_entries]; + + for k2 in 0..num_entries-1 { + // Convert total to bits (little-endian) + set table w/= k2 <- Convert.BigIntAsBoolArray(total, num_bits); + set total *= generator; + set total %= period; + } + + return table; +} + +/// Computes ans=(base^exponent)%modulus. +/// ans must be prepared in zero state. +/// base must be coprime with modulus. +/// Doesn't change exponent. +/// Fig. 7 in the paper +operation ModExpWindow(exponent : Qubit[], ans : Qubit[], base : BigInt, modulus : BigInt, + expWindowLen : Int, mulWindowLen : Int +) : Unit is Adj + Ctl { + let n1 = Length(exponent); + let n2 = Length(ans); + + let expWindows = Arrays.Chunks(expWindowLen, exponent); + + // skip the first two expWindows with a direct lookup + // based on [Gidney 2025](https://arxiv.org/abs/2505.15917) + let skipNum = + if (expWindowLen * 2 < n1) { + expWindowLen * 2 + } else { + n1 + }; + let data = Skip2Data(base, modulus, skipNum, n2); + use output = Qubit[n2]; + within { + Select(data, expWindows[0] + expWindows[1], output); + } apply { + Utils.ParallelCNOT(output, ans); + } + + for i in 2..Length(expWindows)-1 { + let adjustedBase = Math.ExpModL(base, 1L <<< (i * expWindowLen), modulus); + if (i % 2 == 1) { + AddExpModWindowed(adjustedBase, modulus, 1, mulWindowLen, expWindows[i], output, ans); + AddExpModWindowed(Math.InverseModL(adjustedBase, modulus), modulus, -1, mulWindowLen, expWindows[i], ans, output); + } else{ + AddExpModWindowed(adjustedBase, modulus, 1, mulWindowLen, expWindows[i], ans, output); + AddExpModWindowed(Math.InverseModL(adjustedBase, modulus), modulus, -1, mulWindowLen, expWindows[i], output, ans); + } + } + if (Length(expWindows) % 2 == 1) { + // Handle the case where there are more than 2 exponent windows + Utils.ParallelSWAP(ans, output); + } +} + +internal function ModExpData(factor : BigInt, expLength : Int, mulLength : Int, base : BigInt, mod : BigInt, sign : Int, numBits : Int) : Bool[][] { + mutable data = [[false, size = numBits], size = 2^(expLength + mulLength)]; + for b in 0..2^mulLength - 1 { + for a in 0..2^expLength - 1 { + let idx = b * 2^expLength + a; + let value = Math.ModulusL(factor * Convert.IntAsBigInt(b) * Convert.IntAsBigInt(sign) * (base^a), mod); + set data w/= idx <- Convert.BigIntAsBoolArray(value, numBits); + } + } + + data +} + +/// Computes zs += ys * (base ^ xs) % mod (for small registers xs and ys) +/// based on Fig. 5 in [Gidney 2019](https://arxiv.org/abs/1905.07682) +internal operation AddExpModWindowed( + base : BigInt, + mod : BigInt, + sign : Int, + mulWindowLen : Int, + xs : Qubit[], + ys : Qubit[], + zs : Qubit[] +) : Unit is Adj + Ctl { + // split factor into parts + let factorWindows = Arrays.Chunks(mulWindowLen, ys); + + for i in 0..Length(factorWindows)-1 { + // compute data for table lookup + let factorValue = Math.ExpModL(2L, Convert.IntAsBigInt(i * mulWindowLen), mod); + let data = ModExpData(factorValue, Length(xs), Length(factorWindows[i]), base, mod, sign, Length(zs)); + + use output = Qubit[Length(data[0])]; + + within { + Select(data, xs + factorWindows[i], output); + } apply { + ModAdd(output, zs, mod); + } + } +} + +export Multiply, MultiplyWindow, ModExpWindow; \ No newline at end of file diff --git a/test/CG20192_test.py b/test/CG20192_test.py index 1d76752..9318266 100644 --- a/test/CG20192_test.py +++ b/test/CG20192_test.py @@ -2,6 +2,7 @@ from qsharp import init, eval import random +import test_utils @pytest.fixture(scope="session", autouse=True) def setup(): @@ -14,3 +15,13 @@ def test_Multiply(n: int): x, y, t = random.randint(0, 2**n - 1), random.randint(0, 2**n - 1), random.randint(0, 2**n - 1) ans = eval(f"{op}({n},{n},{t}L,{x}L,{y}L)") assert ans == x * y + t + + +@pytest.mark.parametrize("n", [3, 4, 8, 16, 32, 64, 80]) +def test_ModExp(n: int): + op = "QuantumArithmetic.CG20192.ModExpWindow(_,_,_,_,2,2)" + N = 1+2*random.randint(1, 2**(n-1)-1) + a = test_utils.random_coprime(N) + x = random.randint(0, 2**n-1) + ans = eval(f"TestUtils.TestModExp({n},{a}L,{x}L,{N}L,{op})") + assert ans == test_utils.pow_mod(a, x, N) \ No newline at end of file