thofma
CRT for prime power moduli from Hensel lifting
For my application, I need CRT over ℤ[x]/(p^k ⋅ f). Now, this turns out to be decently easy to do with existing functionality in Hecke, but I had to touch a bunch of what seem to be internal functions. I'm wondering if this functionality, could be exposed through some nicer interface (another method of crt_env perhaps?). Please also let me know if this functionality is actually already available through some nice interface that I'm missing.
My code is below. If this doesn't already exist in this form and my implementation seems reasonable, I'm happy to submit a PR, but given my lack of overview of this complete library, I figured I'd check in before attempting that.
using Nemo, Hecke
const p = 2
const k = 8
const pk = 2^8
const ℤ = ZZ
const ℤx, x = PolynomialRing(ℤ, "x")
const ℤp = ResidueRing(ℤ, p)
const ℤxp = PolynomialRing(ℤp, "x")[1]
const ℤpk = ResidueRing(ℤ, pk)
const ℤxpk = PolynomialRing(ℤpk, "x")[1]
cycpoly = cyclotomic(31, x)
using Hecke: HenselCtx, fmpz_poly_raw
ctx = HenselCtx(cycpoly, Hecke.fmpz(p))
# TODO: Why 2k here? The flint docs say that the xgcd factors are 2^l for some l,
# but don't elaborate what l is here. Do we need to use a different method to guarantee that # l>=k?
Hecke.start_lift(ctx, 2*k)
function _crt(p, idx, (v, w, link))
(l, r) = link[idx], link[idx+1]
left = l < 0 ? p[-l] : _crt(p, l+1, (v,w,link))
right = r < 0 ? p[-r] : _crt(p, r+1, (v,w,link))
left * v[idx+1] * w[idx+1] + right * v[idx] * w[idx]
end
poly_array(ctx, a) = [(f = ℤx(); ccall((:fmpz_poly_set, :libflint), Nothing,
(Ref{fmpz_poly}, Ref{fmpz_poly_raw}),
f, a+(i-1)*sizeof(fmpz_poly_raw)); f) for i = 1:(2ctx.r-2)]
function crt(p, h::HenselCtx)
v = map(parent(p[1]), poly_array(ctx, ctx.v))
w = map(parent(p[1]), poly_array(ctx, ctx.w))
link = [unsafe_load(ctx.link, i) for i = 1:(2ctx.r-2)]
_crt(p, length(link)-1, (v,w,link))
end
function crt_inv(p, h::HenselCtx)
factors = map(parent(p), poly_array(ctx, ctx.v)[1:ctx.r])
map(f->rem(p, f), factors)
end
crtpoly = crt(map(ℤxpk, [1:6;]), ctx)
crt_inv(crtpoly, ctx)
Just for clarification. Is the ring in your first sentence equal to Z[x]/(p^k * f) or Z[x]/(p^k, f)?
The latter. Polynomials with coefficients in Z/p^kZ modded out by some polynomial f (usually cyclotomic). This is the standard setting for many lattice cryptography algorithms.
Could you elaborate on what you are trying to do with the polynomials (more abstractly, i.e., more mathemtically)? It might be that we already have some functionality for it.
Sure, I'm currently in particular interested in homomorphic encryption schemes where the plain text space is that particular ring. Homomorphic encryption basically, means that I have some protocol for which Decrypt(f(Encrypt(p1), Encrypt(p2))) == f(p1, p2) for a suitable set of interesting functions f. Now, for schemes of practical interest, we can generally have p1,p2 ∈ ℤ[x]/(n, ϕm) for arbitrary choices of n and where ϕm is the m-th cyclotomic polynomial for some m. For security there needs to be a lower bound on m (usually ~1000-4000). Now, imagine you want to emulate computation on UInt8 data. One way to do this is to set n=256, embed your two UInt8s into the constant terms of the polynomials p1 and p2, and voila you can do your computation. However, as you might imagine, this is fairly wasteful as you're doing the whole polynomial multiplication but only caring about the constant term. A smarter thing to do is to choose some m for which ϕm factors into i factors over p (2 in the case of n=256) and use the polynomial chinese remainder theorem to identify ℤ[x]/(p^k, ϕm) ≃ ℤ[x]/(p^k, f₁) × ⋯ × ℤ[x]/(p^k, fᵢ). That way you can perform i scalar operations for every operation you're doing on encrypted data. I was using Hecke to construct this isomorphism (as shown in the code above), but as far as I understand this is a fairly standard construction, so I thought Hecke might either already have it (and I didn't find it) or it should, in which case I was hoping for some guidance and discussion on the interface.
For the constructive CRT, the crt functions should "just work". Here they don't, since the modulus (p^k) is not prime. But we have a gcdx function, which works also for non-prime modulus. In the following example, I just overwrite the gcdx function to make it work.
julia> Zx, x = FlintZZ["x"]
(Univariate Polynomial Ring in x over Integer Ring, x)
julia> f = cyclotomic(31, x)
x^30+x^29+x^28+x^27+x^26+x^25+x^24+x^23+x^22+x^21+x^20+x^19+x^18+x^17+x^16+x^15+x^14+x^13+x^12+x^11+x^10+x^9+x^8+x^7+x^6+x^5+x^4+x^3+x^2+x+1
julia> H = Hecke.factor_mod_pk_init(f, 2)
factorisation of x^30+x^29+x^28+x^27+x^26+x^25+x^24+x^23+x^22+x^21+x^20+x^19+x^18+x^17+x^16+x^15+x^14+x^13+x^12+x^11+x^10+x^9+x^8+x^7+x^6+x^5+x^4+x^3+x^2+x+1 modulo 2^0
julia> fa = Hecke.factor_mod_pk(H, 8);
julia> factors = collect(keys(ans))
6-element Array{fmpz_poly,1}:
x^5+62*x^4-95*x^3+92*x^2-8*x-1
x^5+18*x^4-85*x^3-35*x^2+67*x-1
x^5+8*x^4-92*x^3+95*x^2-62*x-1
x^5+77*x^4+59*x^3+80*x^2+97*x-1
x^5-97*x^4-80*x^3-59*x^2-77*x-1
x^5-67*x^4+35*x^3+85*x^2-18*x-1
julia> R = ResidueRing(ZZ, 2^8)
Integers modulo 256
julia> Rt, t = R["t"]
(Univariate Polynomial Ring in t over Integers modulo 256, t)
julia> crt([t, t + 1, t + 2], [Rt(factors[1]), Rt(factors[2]), Rt(factors[3])])
ERROR: Modulus not prime in gcdx
Stacktrace:
[1] error(::String) at ./error.jl:33
[2] crt(::nmod_poly, ::nmod_poly, ::nmod_poly, ::nmod_poly) at /home/thofmann/.julia/dev/Nemo/src/flint/nmod_poly.jl:501
julia> Base.gcdx(x::nmod_poly, y::nmod_poly) = Hecke.rresx(x, y)
julia> crt([t + 1, t], [Rt(factors[1]), Rt(factors[2])])
74*t^9+57*t^8+119*t^7+235*t^6+148*t^5+84*t^4+196*t^3+149*t^2+133*t+8
julia> mod(ans, Rt(factors[1]))
t+1
We will make a few adjustments in the CRT code to accommodate for non-prime modulus.
Just a note here that there's a similar assumption in crt_env though there's another branch currently guarded by if false that does work with the above gcdx definition.