Merge pull request #113 from sumiya11/better-modular-computation

Add tests for rational reconstruction
sumiya11 · Jan 30, 2024 · fce221f · fce221f
2 parents 4d24a65 + 7a56315
commit fce221f
Show file tree

Hide file tree

Showing 3 changed files with 117 additions and 89 deletions.
diff --git a/src/reconstruction/crt.jl b/src/reconstruction/crt.jl
@@ -308,40 +308,51 @@ function crt_vec_full!(
     tables_ff::Vector{Vector{Vector{T}}},
     moduli::Vector{T}
 ) where {T <: Integer}
-    indices = [(j, i) for j in 1:length(table_zz) for i in 1:length(table_zz[j])]
-    crt_vec_partial!(table_zz, modulo, tables_ff, moduli, indices)
-end
+    # indices = [(j, i) for j in 1:length(table_zz) for i in 1:length(table_zz[j])]
+    # crt_vec_partial!(table_zz, modulo, tables_ff, moduli, indices)
+    @assert isbitstype(T)
+    @assert length(moduli) > 0
 
-###
-# A couple of inline tests
-
-begin
-    table_zz = [[BigInt(0)], [BigInt(0), BigInt(0)]]
-    modulo = BigInt(0)
-    tables_ff = [[UInt64[2], UInt64[1, 3]], [UInt64[3], UInt64[4, 7]]]
-    moduli = UInt64[7, 11]
-    indices = [(1, 1), (2, 1)]
-    Groebner.crt_vec_partial!(table_zz, modulo, tables_ff, moduli, indices)
-    @assert table_zz[1][1] == BigInt(58)
-    @assert table_zz[2][1] == BigInt(15)
-    Groebner.crt_vec_full!(table_zz, modulo, tables_ff, moduli)
-    @assert table_zz[2][2] == BigInt(73)
-
-    ###
-    table_zz = [[BigInt(0)]]
-    modulo = BigInt(0)
-    tables_ff = [[UInt64[2]]]
-    moduli = UInt64[7]
-    indices = [(1, 1)]
-    Groebner.crt_vec_partial!(table_zz, modulo, tables_ff, moduli, indices)
-    @assert table_zz[1][1] == BigInt(2)
-
-    ###
-    table_zz = [[BigInt(0)]]
-    modulo = BigInt(0)
-    tables_ff = [[Int32[2]], [Int32[3]]]
-    moduli = Int32[11, 13]
-    indices = [(1, 1)]
-    Groebner.crt_vec_partial!(table_zz, modulo, tables_ff, moduli, indices)
-    @assert table_zz[1][1] == BigInt(68)
+    n = length(moduli)
+    # Base case
+    if n == 1
+        table_ff = tables_ff[1]
+        Base.GMP.MPZ.set_ui!(modulo, UInt64(moduli[1]))
+        @inbounds for i in 1:length(table_zz)
+            @assert length(table_zz[i]) == length(table_ff[i])
+            for j in 1:length(table_zz[i])
+                rem_ij = UInt64(table_ff[i][j])
+                @assert 0 <= rem_ij < moduli[1]
+                table_zz[i][j] = rem_ij
+            end
+        end
+        return nothing
+    end
+
+    # n > 1
+    # Precompute CRT multipliers
+    buf = BigInt()
+    n1, n2 = BigInt(), BigInt()
+    mults = Vector{BigInt}(undef, n)
+    for i in 1:length(mults)
+        mults[i] = BigInt(0)
+    end
+    crt_precompute!(modulo, n1, n2, mults, map(UInt64, moduli))
+
+    rems = Vector{UInt64}(undef, n)
+    @inbounds for i in 1:length(table_zz)
+        for j in 1:length(table_zz[i])
+            for t in 1:n
+                @assert length(table_zz[i]) == length(tables_ff[t][i])
+                @assert 0 <= tables_ff[t][i][j] < moduli[t]
+                rems[t] = UInt64(tables_ff[t][i][j])
+            end
+
+            crt!(modulo, buf, n1, n2, rems, mults)
+
+            Base.GMP.MPZ.set!(table_zz[i][j], buf)
+        end
+    end
+
+    nothing
 end
diff --git a/src/reconstruction/ratrec.jl b/src/reconstruction/ratrec.jl
@@ -204,61 +204,24 @@ function ratrec_vec_full!(
     table_zz::Vector{Vector{BigInt}},
     modulo::BigInt
 )
-    indices = [(j, i) for j in 1:length(table_zz) for i in 1:length(table_zz[j])]
-    ratrec_vec_partial!(table_qq, table_zz, modulo, indices)
-end
+    # indices = [(j, i) for j in 1:length(table_zz) for i in 1:length(table_zz[j])]
+    # ratrec_vec_partial!(table_qq, table_zz, modulo, indices)
+    @assert length(table_qq) == length(table_zz)
+    @assert modulo > 1
 
-###
-# A couple of inline tests
-
-begin
-    table_qq = [[BigInt(0) // BigInt(1)], [BigInt(0) // BigInt(1), BigInt(0) // BigInt(1)]]
-    table_zz = [[BigInt(58)], [BigInt(15), BigInt(73)]]
-    modulo = BigInt(77)
-    indices = [(1, 1), (2, 1)]
-    success = Groebner.ratrec_vec_partial!(table_qq, table_zz, modulo, indices)
-    @assert success
-    @assert table_qq[1][1] == BigInt(1) // BigInt(4)
-    @assert table_qq[2][1] == BigInt(-2) // BigInt(5)
-    success = Groebner.ratrec_vec_full!(table_qq, table_zz, modulo)
-    @assert success
-    @assert table_qq[2][2] == BigInt(-4) // BigInt(1)
-
-    # Impossible reconstruction
-    table_qq = [[BigInt(0) // BigInt(1)]]
-    table_zz = [[BigInt(643465418)]]
-    modulo = BigInt(2^31 - 1)
-    indices = [(1, 1)]
-    success = Groebner.ratrec_vec_partial!(table_qq, table_zz, modulo, indices)
-    @assert !success
-
-    ### A complete example. ###
-    table_zz = [[BigInt(0)], [BigInt(0), BigInt(0)]]
-    table_qq = [[BigInt(0) // BigInt(1)], [BigInt(0) // BigInt(1), BigInt(0) // BigInt(1)]]
-    modulo = BigInt(0)
-    # Collected remainders
-    tables_ff = [[UInt64[2], UInt64[1, 3]], [UInt64[3], UInt64[4, 7]]]
-    moduli = UInt64[7, 11]
-    # Try partial reconstruction first
-    indices = [(1, 1), (2, 1)]
-    Groebner.crt_vec_partial!(table_zz, modulo, tables_ff, moduli, indices)
-    success = Groebner.ratrec_vec_partial!(table_qq, table_zz, modulo, indices)
-    if success
-        # Partial reconstrction succeeded;
-        # Try full reconstruction
-        Groebner.crt_vec_full!(table_zz, modulo, tables_ff, moduli)
-        success = Groebner.ratrec_vec_full!(table_qq, table_zz, modulo)
-        if success
-            # Full reconstruction succeeded
-            @assert table_qq[1][1] == BigInt(1) // BigInt(4)
-            @assert table_qq[2][1] == BigInt(-2) // BigInt(5)
-            @assert table_qq[2][2] == BigInt(-4) // BigInt(1)
-        else
-            # Partial reconstruction succeeded, but full reconstrction failed;
-            # Go for another iteration ...
+    modulo_nemo = Nemo.ZZRingElem(modulo)
+    @inbounds for i in 1:length(table_zz)
+        @assert length(table_zz[i]) == length(table_qq[i])
+        for j in 1:length(table_zz[i])
+            rem_nemo = Nemo.ZZRingElem(table_zz[i][j])
+            @assert rem_nemo >= 0
+
+            success, (num, den) = ratrec_nemo(rem_nemo, modulo_nemo)
+            table_qq[i][j] = Base.unsafe_rational(num, den)
+
+            !success && return false
         end
-    else
-        # Partial reconstrction failed;
-        # Go for another iteration ...
     end
+
+    true
 end
diff --git a/test/reconstruction/ratrec.jl b/test/reconstruction/ratrec.jl
@@ -43,3 +43,57 @@ import Primes
         end
     end
 end
+
+function table_compute_res(table, m, T)
+    [
+        map(c -> T(mod(mod(numerator(c), m) * invmod(denominator(c), m), m)), table[j]) for
+        j in 1:length(table)
+    ]
+end
+
+@testset "vec ratrec + crt" begin
+    table_qq = [[BigInt(0) // BigInt(1)], [BigInt(0) // BigInt(1), BigInt(0) // BigInt(1)]]
+    table_zz = [[BigInt(58)], [BigInt(15), BigInt(73)]]
+    modulo = BigInt(77)
+    success = Groebner.ratrec_vec_full!(table_qq, table_zz, modulo)
+    @assert success
+    @assert table_qq[1][1] == BigInt(1) // BigInt(4)
+    @assert table_qq[2][1] == BigInt(-2) // BigInt(5)
+    @test table_qq[2][2] == BigInt(-4) // BigInt(1)
+
+    # Impossible reconstruction
+    table_qq = [[BigInt(0) // BigInt(1)]]
+    table_zz = [[BigInt(643465418)]]
+    modulo = BigInt(2^31 - 1)
+    indices = [(1, 1)]
+    success = Groebner.ratrec_vec_full!(table_qq, table_zz, modulo)
+    @test !success
+
+    boot = 10
+    for _ in 1:boot
+        k = rand(1:100)
+        moduli = map(UInt32, Primes.nextprimes(rand((2^20):(2^28)), k))
+        m = rand(1:10)
+        P0 = prod(map(BigInt, moduli))
+        P1 = floor(BigInt, sqrt(P0 / 2) - 1)
+        N, D = P1, P1
+        _s = rand(1:10, m)
+        # Generate rational numbers
+        table_ans = [rand((-N):N, _s[_m]) .// rand(1:D, _s[_m]) for _m in 1:m]
+        # Compute k tables of residuals
+        tables_res = [table_compute_res(table_ans, moduli[i], UInt32) for i in 1:k]
+        # Run CRT
+        table_zz =
+            [[BigInt(0) for _ in 1:length(table_ans[i])] for i in 1:length(table_ans)]
+        Groebner.crt_vec_full!(table_zz, modulo, tables_res, moduli)
+        @test modulo == prod(BigInt, moduli)
+        @test table_zz == table_compute_res(table_ans, modulo, BigInt)
+        # Run rat. rec.
+        table_qq = [
+            [Rational{BigInt}(0) for _ in 1:length(table_ans[i])] for
+            i in 1:length(table_ans)
+        ]
+        success = Groebner.ratrec_vec_full!(table_qq, table_zz, modulo)
+        @test table_ans == table_qq && success
+    end
+end