Add windowed modexp from Gidney 2019

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}");
 }

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;
+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;

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)