From 7a56315af359b006e374e49a0924d7953e644490 Mon Sep 17 00:00:00 2001 From: Alexander Demin Date: Wed, 31 Jan 2024 02:39:06 +0300 Subject: [PATCH] add tests for rational reconstruction --- src/reconstruction/crt.jl | 81 ++++++++++++++++++++--------------- src/reconstruction/ratrec.jl | 71 ++++++++---------------------- test/reconstruction/ratrec.jl | 54 +++++++++++++++++++++++ 3 files changed, 117 insertions(+), 89 deletions(-) diff --git a/src/reconstruction/crt.jl b/src/reconstruction/crt.jl index ee9c0b33..e2952561 100644 --- 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 index 5ab526a4..e8607c24 100644 --- 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 index 7f038d61..e7f08531 100644 --- 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