1 | initial version |

Note that computing the lift of 1 with respect to your `C`

modulo a prime is fast – it takes only a few seconds. Doing this for many primes can allow to lift the result to ℚ. In this case, the result `F`

(see attachment) is huge, though, with coefficients as large as 4826 and 4818 digits for numerator and denominator.

```
sage: F, C = load("lift_one_result_F_C.sobj") # saved with Sage 9.3
sage: sum(f*c for f,c in zip(F, C)) # 1 minute
1
```

I am not sure if this solution would be unique. Maybe there is a smaller one.

To obtain the result with respect to `A1`

and `A2`

, you could try to lift `C1`

and `C2`

individually, but this might make the result even larger.

I will update this answer with more details and my ad-hoc implementation, tomorrow.

2 | No.2 Revision |

Note that computing the lift of 1 with respect to your `C`

modulo a prime is fast – it takes only a few seconds. Doing this for many primes can allow to lift the result to ℚ. In this case, the result `F`

~~(see ~~(~~see attachment~~ uploaded here (31MB), for now, as the file seems too large for an attachment) is huge, though, with coefficients as large as 4826 and 4818 digits for numerator and denominator.

```
sage: F, C = load("lift_one_result_F_C.sobj") # saved with Sage 9.3
sage: sum(f*c for f,c in zip(F, C)) # 1 minute
1
```

I am not sure if this solution would be unique. Maybe there is a smaller one.

To obtain the result with respect to `A1`

and `A2`

, you could try to lift `C1`

and `C2`

individually, but this might make the result even larger.

I will update this answer with more details and my ad-hoc implementation, tomorrow.

3 | No.3 Revision |

Note that computing the lift of 1 with respect to your `C`

modulo a prime is fast – it takes only a few seconds. Doing this for many primes can allow to lift the result to ℚ. In this case, the result `F`

(~~see attachment~~ uploaded here (31MB), for now, as the file seems too large for an attachment) is huge, though, with coefficients as large as 4826 and 4818 digits for numerator and denominator.

```
sage: F, C = load("lift_one_result_F_C.sobj") # saved with Sage 9.3
sage: sum(f*c for f,c in zip(F, C)) # 1 minute
1
```

I am not sure if this solution would be unique. Maybe there is a smaller one.

To obtain the result with respect to `A1`

and `A2`

, you could try to lift `C1`

and `C2`

individually, but this might make the result even larger.

When computing the lift of 1 mod p, most of the polynomials in the result have the same monomials as their support for varying p. This makes it very likely that the correct representative over ℚ has those same monomials as support and that the polynomials we computed mod p are the mod-p versions of the desired polynomials over ℚ.

This can fail for some primes which we need to ignore. In the function `likely_support`

, we use a simple ad-hoc heuristic to find the exponents that probably occur in the ℚ-result: if a monomial occurs in 50% of the mod-p polynomials (i.e. modulo 50% of the primes or more), we expect that it should occur in the ℚ-result as well. The 50% threshold does not matter much. We could have picked 10% or 90% instead, I suppose.

We only keep the primes p for which all the polynomials in the lift of 1 mod p have the expected supports. Out of the 2800 primes we started with, 2358 remain.

Next, we use the Chinese remainder theorem `CRT_list`

to lift the coefficients to ℤ modulo the product of the primes and then `rational_reconstruction`

to obtain rational coefficients, both implemented in `lift_poly_crt`

. For `rational_reconstruction`

to work, the product of the primes must have about twice as many digits as the numerators and denominators in the result over ℚ (which I do not know how to determine in advance).

I have used `cached_function.precompute`

in order to parallelize some computations because it is a bit less awkward to use than the `@parallel`

decorator, but also because it allows to increase the number of primes until rational reconstruction works, without redoing all the previous slow computations.

The slowest part of the computation is lifting 1 modulo the 2800 primes, so this benefits a lot from parallelization. Each prime takes about 10 seconds, so this would take about 7.5 hours on a single thread.

```
R1.<u0,u1,u2,v0,v1,v2,v3,v4,v5,v6>=PolynomialRing(QQ,10)
A1=[u0 + 7/5*u1 + 7/5*u2 - 4/125,
5*v0 + 5*v1 + 7*v3 + 7*v5 + 1/5,
25*v0^2 + 25*v1^2 + 35*v3^2 + 35*v5^2 - 4/5,
5*v0^3 + 5*v1^3 + 7*v3^3 + 7*v5^3 - v0^2 + 1/125,
5*v0*v1^2 + 5*v1*v2^2 + 7*v3*v4^2 + 7*v5*v6^2 + 1/125,
5*u0*v1 - v1^2 + 7*u1*v3 + 7*u2*v5 + 1/125,
5*v1 + 5*v2 + 7*v4 + 7*v6 + 1/5,
25*v0*v1 + 25*v1*v2 + 35*v3*v4 + 35*v5*v6 + 1/5,
5*v0^2*v1 + 5*v1^2*v2 + 7*v3^2*v4 + 7*v5^2*v6 - v1^2 + 1/125,
25*v1^2 + 25*v2^2 + 35*v4^2 + 35*v6^2 - 4/5,
5*v1^3 + 5*v2^3 + 7*v4^3 + 7*v6^3 - u0 + 1/125,
5*u0*v2 - v2^2 + 7*u1*v4 + 7*u2*v6 + 1/125]
I1=R1.ideal(A1)
C1=list(I1.groebner_basis())
R2.<y0,y1,y2,z0,z1,z2,z3,z4,z5,z6>=PolynomialRing(QQ,10)
A2=[y0 + 7/5*y1 + 7/5*y2 - 4/125,
5*z0 + 5*z1 + 7*z3 + 7*z5 + 1/5,
25*z0^2 + 25*z1^2 + 35*z3^2 + 35*z5^2 - 4/5,
5*z0^3 + 5*z1^3 + 7*z3^3 + 7*z5^3 - y0 + 1/125,
5*y0*z0 - z0^2 + 7*y1*z3 + 7*y2*z5 + 1/125,
5*z0*z1^2 + 5*z1*z2^2 + 7*z3*z4^2 + 7*z5*z6^2 - z1^2 + 1/125,
5*z1 + 5*z2 + 7*z4 + 7*z6 + 1/5,
25*z0*z1 + 25*z1*z2 + 35*z3*z4 + 35*z5*z6 + 1/5,
5*z0^2*z1 + 5*z1^2*z2 + 7*z3^2*z4 + 7*z5^2*z6 + 1/125,
5*y0*z1 - z1^2 + 7*y1*z4 + 7*y2*z6 + 1/125,
25*z1^2 + 25*z2^2 + 35*z4^2 + 35*z6^2 - 4/5,
5*z1^3 + 5*z2^3 + 7*z4^3 + 7*z6^3 - z2^2 + 1/125]
I2=R2.ideal(A2)
C2=list(I2.groebner_basis())
R.<u0,u1,u2,v0,v1,v2,v3,v4,v5,v6,y0,y1,y2,z0,z1,z2,z3,z4,z5,z6>=PolynomialRing(QQ,20)
C=C1+C2+[5*u0*z2 + 5*u1*z4 + 7*u2*z6 - u0 + 1/125] # the last one is the extra one
I=R.ideal(C)
########################################################################
@cached_function
def lift_one_modp(Ip):
Fp = Ip.ring()(1).lift(Ip.ideal())
print("completed p = %s" % Fp.ring().characteristic())
return Fp
@cached_function
def lift_poly_crt(fps):
ps = [fp.parent().characteristic() for fp in fps]
Ds = [fp.dict() for fp in fps]
supp = set([e for Di in Ds for e in Di])
return R({e: CRT_list([Di[e].lift() if e in Di else 0 for Di in Ds],
list(ps)).rational_reconstruction(prod(ps)) for e in supp})
def likely_support(fps):
supps = [set(fp.dict().keys()) for fp in fps]
supp_union = set.union(*supps)
supp = set()
half = len(supps) * 1/2 # exact amount does not matter
for e in supp_union:
k = 0
for s in supps:
if e in s:
k += 1
if k > half: # if exponent vector is in more than half the supports, we will
```~~update ~~probably want to keep it
supp.add(e)
break
return supp
ps = sorted(RecursivelyEnumeratedSet([next_prime(16000)], lambda p: [next_prime(p)])[:2800]) # many primes
Rps = [R.change_ring(GF(p)) for p in ps]
Ips = [I.change_ring(Rp).gens() for Rp in Rps]
CPUs = 16 # choose based on hardware
lift_one_modp.precompute(Ips, num_processes=CPUs)
assert all(type(r) != str for r in lift_one_modp.cache.values()) # otherwise, there was an exception
Fps = [lift_one_modp(Ip) for Ip in Ips]
supps = [likely_support(fps) for fps in zip(*Fps)]
# keep indices of primes for which all polynomials have the expected support
idxs = [j for j, Fp in enumerate(Fps) if all(set(fp.dict()) == supp for fp, supp in zip(Fp, supps))]
# restrict the results to those primes
ps1 = tuple(ps[j] for j in idxs)
Fps1 = [Fps[j] for j in idxs]
fpss = [tuple(fps) for fps in zip(*Fps1)]
modulus1 = prod(ps1)
floor(log(modulus1, 10)) # 10538 (must be large enough for rational reconstruction: this ~~answer with more details and my ad-hoc implementation, tomorrow. ~~needs about twice as many digits as the numerators or denominators in the result have)
lift_poly_crt.precompute([(fps,) for fps in fpss], num_processes=CPUs)
assert all(type(r) != str for r in lift_poly_crt.cache.values()) # otherwise, there was an exception (if so, remove it from cache by lift_poly_crt.cache.clear())
F = [lift_poly_crt(fps) for fps in fpss] # the final result

4 | No.4 Revision |

Note that computing the lift of 1 with respect to your `C`

modulo a prime is fast – it takes only a few seconds. Doing this for many primes can allow to lift the result to ℚ. In this case, the result `F`

(~~see attachment~~ uploaded here (31MB), for now, as the file seems too large for an attachment) is huge, though, with coefficients as large as 4826 and 4818 digits for numerator and denominator.

```
sage: F, C = load("lift_one_result_F_C.sobj") # saved with Sage 9.3
sage: sum(f*c for f,c in zip(F, C)) # 1 minute
1
```

I am not sure if this solution would be unique. Maybe there is a smaller one.

`A1`

and `A2`

, you could try to lift `C1`

and `C2`

individually, but this might make the result even larger.

When computing the lift of 1 mod p, most of the polynomials in the result have the same monomials as their support for varying p. This makes it very likely that the correct representative over ℚ has those same monomials as support and that the polynomials we computed mod p are the mod-p versions of the desired polynomials over ℚ.

This can fail for some primes which we need to ignore. In the function `likely_support`

, we use a simple ad-hoc heuristic to find the exponents that probably occur in the ℚ-result: if a monomial occurs in 50% of the mod-p polynomials (i.e. modulo 50% of the primes or more), we expect that it should occur in the ℚ-result as well. The 50% threshold does not matter much. We could have picked 10% or 90% instead, I suppose.

We only keep the primes p for which all the polynomials in the lift of 1 mod p have the expected supports. Out of the 2800 primes we started with, 2358 remain.

Next, we use the Chinese remainder theorem `CRT_list`

to lift the coefficients to ℤ modulo the product of the primes and then `rational_reconstruction`

to obtain rational coefficients, both implemented in `lift_poly_crt`

. For `rational_reconstruction`

to work, the product of the primes must have about twice as many digits as the numerators and denominators in the result over ℚ (which I do not know how to determine in advance).

I have used `cached_function.precompute`

in order to parallelize some computations because it is a bit less awkward to use than the `@parallel`

decorator, but also because it allows to increase the number of primes until rational reconstruction works, without redoing all the previous slow computations.

The slowest part of the computation is lifting 1 modulo the 2800 primes, so this benefits a lot from parallelization. Each prime takes about 10 seconds, so this would take about 7.5 hours on a single thread.

```
R1.<u0,u1,u2,v0,v1,v2,v3,v4,v5,v6>=PolynomialRing(QQ,10)
A1=[u0 + 7/5*u1 + 7/5*u2 - 4/125,
5*v0 + 5*v1 + 7*v3 + 7*v5 + 1/5,
25*v0^2 + 25*v1^2 + 35*v3^2 + 35*v5^2 - 4/5,
5*v0^3 + 5*v1^3 + 7*v3^3 + 7*v5^3 - v0^2 + 1/125,
5*v0*v1^2 + 5*v1*v2^2 + 7*v3*v4^2 + 7*v5*v6^2 + 1/125,
5*u0*v1 - v1^2 + 7*u1*v3 + 7*u2*v5 + 1/125,
5*v1 + 5*v2 + 7*v4 + 7*v6 + 1/5,
25*v0*v1 + 25*v1*v2 + 35*v3*v4 + 35*v5*v6 + 1/5,
5*v0^2*v1 + 5*v1^2*v2 + 7*v3^2*v4 + 7*v5^2*v6 - v1^2 + 1/125,
25*v1^2 + 25*v2^2 + 35*v4^2 + 35*v6^2 - 4/5,
5*v1^3 + 5*v2^3 + 7*v4^3 + 7*v6^3 - u0 + 1/125,
5*u0*v2 - v2^2 + 7*u1*v4 + 7*u2*v6 + 1/125]
I1=R1.ideal(A1)
C1=list(I1.groebner_basis())
R2.<y0,y1,y2,z0,z1,z2,z3,z4,z5,z6>=PolynomialRing(QQ,10)
A2=[y0 + 7/5*y1 + 7/5*y2 - 4/125,
5*z0 + 5*z1 + 7*z3 + 7*z5 + 1/5,
25*z0^2 + 25*z1^2 + 35*z3^2 + 35*z5^2 - 4/5,
5*z0^3 + 5*z1^3 + 7*z3^3 + 7*z5^3 - y0 + 1/125,
5*y0*z0 - z0^2 + 7*y1*z3 + 7*y2*z5 + 1/125,
5*z0*z1^2 + 5*z1*z2^2 + 7*z3*z4^2 + 7*z5*z6^2 - z1^2 + 1/125,
5*z1 + 5*z2 + 7*z4 + 7*z6 + 1/5,
25*z0*z1 + 25*z1*z2 + 35*z3*z4 + 35*z5*z6 + 1/5,
5*z0^2*z1 + 5*z1^2*z2 + 7*z3^2*z4 + 7*z5^2*z6 + 1/125,
5*y0*z1 - z1^2 + 7*y1*z4 + 7*y2*z6 + 1/125,
25*z1^2 + 25*z2^2 + 35*z4^2 + 35*z6^2 - 4/5,
5*z1^3 + 5*z2^3 + 7*z4^3 + 7*z6^3 - z2^2 + 1/125]
I2=R2.ideal(A2)
C2=list(I2.groebner_basis())
R.<u0,u1,u2,v0,v1,v2,v3,v4,v5,v6,y0,y1,y2,z0,z1,z2,z3,z4,z5,z6>=PolynomialRing(QQ,20)
C=C1+C2+[5*u0*z2 + 5*u1*z4 + 7*u2*z6 - u0 + 1/125] # the last one is the extra one
I=R.ideal(C)
########################################################################
@cached_function
def lift_one_modp(Ip):
Fp = Ip.ring()(1).lift(Ip.ideal())
print("completed p = %s" % Fp.ring().characteristic())
return Fp
@cached_function
def lift_poly_crt(fps):
ps = [fp.parent().characteristic() for fp in fps]
Ds = [fp.dict() for fp in fps]
supp = set([e for Di in Ds for e in Di])
return R({e: CRT_list([Di[e].lift() if e in Di else 0 for Di in Ds],
list(ps)).rational_reconstruction(prod(ps)) for e in supp})
def likely_support(fps):
supps = [set(fp.dict().keys()) for fp in fps]
supp_union = set.union(*supps)
supp = set()
half = len(supps) * 1/2 # exact amount does not matter
for e in supp_union:
k = 0
for s in supps:
if e in s:
k += 1
if k > half: # if exponent vector is in more than half the supports, we will probably want to keep it
supp.add(e)
break
return supp
ps = sorted(RecursivelyEnumeratedSet([next_prime(16000)], lambda p: [next_prime(p)])[:2800]) # many primes
Rps = [R.change_ring(GF(p)) for p in ps]
Ips = [I.change_ring(Rp).gens() for Rp in Rps]
CPUs = 16 # choose based on hardware
lift_one_modp.precompute(Ips, num_processes=CPUs)
assert all(type(r) != str for r in lift_one_modp.cache.values()) # otherwise, there was an exception
Fps = [lift_one_modp(Ip) for Ip in Ips]
supps = [likely_support(fps) for fps in zip(*Fps)]
# keep indices of primes for which all polynomials have the expected support
idxs = [j for j, Fp in enumerate(Fps) if all(set(fp.dict()) == supp for fp, supp in zip(Fp, supps))]
# restrict the results to those primes
ps1 = tuple(ps[j] for j in idxs)
Fps1 = [Fps[j] for j in idxs]
fpss = [tuple(fps) for fps in zip(*Fps1)]
modulus1 = prod(ps1)
floor(log(modulus1, 10)) # 10538 (must be large enough for rational reconstruction: this needs about twice as many digits as the numerators or denominators in the result have)
lift_poly_crt.precompute([(fps,) for fps in fpss], num_processes=CPUs)
assert all(type(r) != str for r in lift_poly_crt.cache.values()) # otherwise, there was an exception (if so, remove it from cache by lift_poly_crt.cache.clear())
F = [lift_poly_crt(fps) for fps in fpss] # the final result
```

5 | No.5 Revision |

Note that computing the lift of 1 with respect to your `C`

modulo a prime is fast – it takes only a few seconds. Doing this for many primes can allow to lift the result to ℚ. In this case, the result `F`

(~~see attachment~~ uploaded here (31MB), for now, as the file seems too large for an attachment) is huge, though, with coefficients as large as 4826 and 4818 digits for numerator and denominator.

```
sage: F, C = load("lift_one_result_F_C.sobj") # saved with Sage 9.3
sage: sum(f*c for f,c in zip(F, C)) # 1 minute
1
```

I am not sure if this solution would be unique. Maybe there is a smaller one.

`A1`

and `A2`

, you could try to lift `C1`

and `C2`

individually, but this might make the result even larger.

When computing the lift of 1 mod p, most of the polynomials in the result have the same monomials as their support for varying p. This makes it very likely that the correct representative over ℚ has those same monomials as support and that the polynomials we computed mod p are the mod-p versions of the desired polynomials over ℚ.

This can fail for some primes which we need to ignore. In the function `likely_support`

, we use a simple ad-hoc heuristic to find the exponents that probably occur in the ℚ-result: if a monomial occurs in 50% of the mod-p polynomials (i.e. modulo 50% of the primes or more), we expect that it should occur in the ℚ-result as well. The 50% threshold does not matter much. We could have picked 10% or 90% instead, I suppose.

We only keep the primes p for which all the polynomials in the lift of 1 mod p have the expected supports. Out of the 2800 primes we started with, 2358 remain.

Next, we use the Chinese remainder theorem `CRT_list`

to lift the coefficients to ℤ modulo the product of the primes and then `rational_reconstruction`

to obtain rational coefficients, both implemented in `lift_poly_crt`

. For `rational_reconstruction`

to work, the product of the primes must have about twice as many digits as the numerators and denominators in the result over ℚ (which I do not know how to determine in ~~advance).~~advance). In our case, obtain a modulus with 10538 digits.

I have used `cached_function.precompute`

in order to parallelize some computations because it is a bit less awkward to use than the `@parallel`

decorator, but also because it allows to increase the number of primes until rational reconstruction works, without redoing all the previous slow computations.

The slowest part of the computation is lifting 1 modulo the 2800 primes, so this benefits a lot from parallelization. Each prime takes about 10 seconds, so this would take about 7.5 hours on a single thread.

```
R1.<u0,u1,u2,v0,v1,v2,v3,v4,v5,v6>=PolynomialRing(QQ,10)
A1=[u0 + 7/5*u1 + 7/5*u2 - 4/125,
5*v0 + 5*v1 + 7*v3 + 7*v5 + 1/5,
25*v0^2 + 25*v1^2 + 35*v3^2 + 35*v5^2 - 4/5,
5*v0^3 + 5*v1^3 + 7*v3^3 + 7*v5^3 - v0^2 + 1/125,
5*v0*v1^2 + 5*v1*v2^2 + 7*v3*v4^2 + 7*v5*v6^2 + 1/125,
5*u0*v1 - v1^2 + 7*u1*v3 + 7*u2*v5 + 1/125,
5*v1 + 5*v2 + 7*v4 + 7*v6 + 1/5,
25*v0*v1 + 25*v1*v2 + 35*v3*v4 + 35*v5*v6 + 1/5,
5*v0^2*v1 + 5*v1^2*v2 + 7*v3^2*v4 + 7*v5^2*v6 - v1^2 + 1/125,
25*v1^2 + 25*v2^2 + 35*v4^2 + 35*v6^2 - 4/5,
5*v1^3 + 5*v2^3 + 7*v4^3 + 7*v6^3 - u0 + 1/125,
5*u0*v2 - v2^2 + 7*u1*v4 + 7*u2*v6 + 1/125]
I1=R1.ideal(A1)
C1=list(I1.groebner_basis())
R2.<y0,y1,y2,z0,z1,z2,z3,z4,z5,z6>=PolynomialRing(QQ,10)
A2=[y0 + 7/5*y1 + 7/5*y2 - 4/125,
5*z0 + 5*z1 + 7*z3 + 7*z5 + 1/5,
25*z0^2 + 25*z1^2 + 35*z3^2 + 35*z5^2 - 4/5,
5*z0^3 + 5*z1^3 + 7*z3^3 + 7*z5^3 - y0 + 1/125,
5*y0*z0 - z0^2 + 7*y1*z3 + 7*y2*z5 + 1/125,
5*z0*z1^2 + 5*z1*z2^2 + 7*z3*z4^2 + 7*z5*z6^2 - z1^2 + 1/125,
5*z1 + 5*z2 + 7*z4 + 7*z6 + 1/5,
25*z0*z1 + 25*z1*z2 + 35*z3*z4 + 35*z5*z6 + 1/5,
5*z0^2*z1 + 5*z1^2*z2 + 7*z3^2*z4 + 7*z5^2*z6 + 1/125,
5*y0*z1 - z1^2 + 7*y1*z4 + 7*y2*z6 + 1/125,
25*z1^2 + 25*z2^2 + 35*z4^2 + 35*z6^2 - 4/5,
5*z1^3 + 5*z2^3 + 7*z4^3 + 7*z6^3 - z2^2 + 1/125]
I2=R2.ideal(A2)
C2=list(I2.groebner_basis())
R.<u0,u1,u2,v0,v1,v2,v3,v4,v5,v6,y0,y1,y2,z0,z1,z2,z3,z4,z5,z6>=PolynomialRing(QQ,20)
C=C1+C2+[5*u0*z2 + 5*u1*z4 + 7*u2*z6 - u0 + 1/125] # the last one is the extra one
I=R.ideal(C)
########################################################################
@cached_function
def lift_one_modp(Ip):
Fp = Ip.ring()(1).lift(Ip.ideal())
print("completed p = %s" % Fp.ring().characteristic())
return Fp
@cached_function
def lift_poly_crt(fps):
ps = [fp.parent().characteristic() for fp in fps]
Ds = [fp.dict() for fp in fps]
supp = set([e for Di in Ds for e in Di])
return R({e: CRT_list([Di[e].lift() if e in Di else 0 for Di in Ds],
list(ps)).rational_reconstruction(prod(ps)) for e in supp})
def likely_support(fps):
supps = [set(fp.dict().keys()) for fp in fps]
supp_union = set.union(*supps)
supp = set()
half = len(supps) * 1/2 # exact amount does not matter
for e in supp_union:
k = 0
for s in supps:
if e in s:
k += 1
if k > half: # if exponent vector is in more than half the supports, we will probably want to keep it
supp.add(e)
break
return supp
ps = sorted(RecursivelyEnumeratedSet([next_prime(16000)], lambda p: [next_prime(p)])[:2800]) # many primes
Rps = [R.change_ring(GF(p)) for p in ps]
Ips = [I.change_ring(Rp).gens() for Rp in Rps]
CPUs = 16 # choose based on hardware
lift_one_modp.precompute(Ips, num_processes=CPUs)
assert all(type(r) != str for r in lift_one_modp.cache.values()) # otherwise, there was an exception
Fps = [lift_one_modp(Ip) for Ip in Ips]
supps = [likely_support(fps) for fps in zip(*Fps)]
# keep indices of primes for which all polynomials have the expected support
idxs = [j for j, Fp in enumerate(Fps) if all(set(fp.dict()) == supp for fp, supp in zip(Fp, supps))]
# restrict the results to those primes
ps1 = tuple(ps[j] for j in idxs)
Fps1 = [Fps[j] for j in idxs]
fpss = [tuple(fps) for fps in zip(*Fps1)]
modulus1 = prod(ps1)
floor(log(modulus1, 10)) # 10538 (must be large enough for rational reconstruction: this needs about twice as many digits as the numerators or denominators in the result have)
lift_poly_crt.precompute([(fps,) for fps in fpss], num_processes=CPUs)
assert all(type(r) != str for r in lift_poly_crt.cache.values()) # otherwise, there was an exception (if so, remove it from cache by lift_poly_crt.cache.clear())
F = [lift_poly_crt(fps) for fps in fpss] # the final result
```

6 | No.6 Revision |

`C`

modulo a prime is fast – it takes only a few seconds. Doing this for many primes can allow to lift the result to ℚ. In this case, the result `F`

(~~see attachment~~ uploaded here (31MB), for now, as the file seems too large for an attachment) is huge, though, with coefficients as large as 4826 and 4818 digits for numerator and denominator.

```
sage: F, C = load("lift_one_result_F_C.sobj") # saved with Sage 9.3
sage: sum(f*c for f,c in zip(F, C)) # 1 minute
1
```

I am not sure if this solution would be unique. Maybe there is a smaller one.

`A1`

and `A2`

, you could try to lift `C1`

and `C2`

individually, but this might make the result even larger.

`likely_support`

, we use a simple ad-hoc heuristic to find the exponents that probably occur in the ℚ-result: if a monomial occurs in 50% of the mod-p polynomials (i.e. modulo 50% of the primes or more), we expect that it should occur in the ℚ-result as well. The 50% threshold does not matter much. We could have picked 10% or 90% instead, I suppose.

Next, we use the Chinese remainder theorem `CRT_list`

to lift the coefficients to ℤ modulo the product of the primes and then `rational_reconstruction`

to obtain rational coefficients, both implemented in `lift_poly_crt`

. For `rational_reconstruction`

to work, the product of the primes must have about twice as many digits as the numerators and denominators in the result over ℚ (which I do not know how to determine in advance). In our case, we obtain a modulus with 10538 digits.

`cached_function.precompute`

in order to parallelize some computations because it is a bit less awkward to use than the `@parallel`

decorator, but also because it allows to increase the number of primes until rational reconstruction works, without redoing all the previous slow computations.

```
R1.<u0,u1,u2,v0,v1,v2,v3,v4,v5,v6>=PolynomialRing(QQ,10)
A1=[u0 + 7/5*u1 + 7/5*u2 - 4/125,
5*v0 + 5*v1 + 7*v3 + 7*v5 + 1/5,
25*v0^2 + 25*v1^2 + 35*v3^2 + 35*v5^2 - 4/5,
5*v0^3 + 5*v1^3 + 7*v3^3 + 7*v5^3 - v0^2 + 1/125,
5*v0*v1^2 + 5*v1*v2^2 + 7*v3*v4^2 + 7*v5*v6^2 + 1/125,
5*u0*v1 - v1^2 + 7*u1*v3 + 7*u2*v5 + 1/125,
5*v1 + 5*v2 + 7*v4 + 7*v6 + 1/5,
25*v0*v1 + 25*v1*v2 + 35*v3*v4 + 35*v5*v6 + 1/5,
5*v0^2*v1 + 5*v1^2*v2 + 7*v3^2*v4 + 7*v5^2*v6 - v1^2 + 1/125,
25*v1^2 + 25*v2^2 + 35*v4^2 + 35*v6^2 - 4/5,
5*v1^3 + 5*v2^3 + 7*v4^3 + 7*v6^3 - u0 + 1/125,
5*u0*v2 - v2^2 + 7*u1*v4 + 7*u2*v6 + 1/125]
I1=R1.ideal(A1)
C1=list(I1.groebner_basis())
R2.<y0,y1,y2,z0,z1,z2,z3,z4,z5,z6>=PolynomialRing(QQ,10)
A2=[y0 + 7/5*y1 + 7/5*y2 - 4/125,
5*z0 + 5*z1 + 7*z3 + 7*z5 + 1/5,
25*z0^2 + 25*z1^2 + 35*z3^2 + 35*z5^2 - 4/5,
5*z0^3 + 5*z1^3 + 7*z3^3 + 7*z5^3 - y0 + 1/125,
5*y0*z0 - z0^2 + 7*y1*z3 + 7*y2*z5 + 1/125,
5*z0*z1^2 + 5*z1*z2^2 + 7*z3*z4^2 + 7*z5*z6^2 - z1^2 + 1/125,
5*z1 + 5*z2 + 7*z4 + 7*z6 + 1/5,
25*z0*z1 + 25*z1*z2 + 35*z3*z4 + 35*z5*z6 + 1/5,
5*z0^2*z1 + 5*z1^2*z2 + 7*z3^2*z4 + 7*z5^2*z6 + 1/125,
5*y0*z1 - z1^2 + 7*y1*z4 + 7*y2*z6 + 1/125,
25*z1^2 + 25*z2^2 + 35*z4^2 + 35*z6^2 - 4/5,
5*z1^3 + 5*z2^3 + 7*z4^3 + 7*z6^3 - z2^2 + 1/125]
I2=R2.ideal(A2)
C2=list(I2.groebner_basis())
R.<u0,u1,u2,v0,v1,v2,v3,v4,v5,v6,y0,y1,y2,z0,z1,z2,z3,z4,z5,z6>=PolynomialRing(QQ,20)
C=C1+C2+[5*u0*z2 + 5*u1*z4 + 7*u2*z6 - u0 + 1/125] # the last one is the extra one
I=R.ideal(C)
########################################################################
@cached_function
def lift_one_modp(Ip):
Fp = Ip.ring()(1).lift(Ip.ideal())
print("completed p = %s" % Fp.ring().characteristic())
return Fp
@cached_function
def lift_poly_crt(fps):
ps = [fp.parent().characteristic() for fp in fps]
Ds = [fp.dict() for fp in fps]
supp = set([e for Di in Ds for e in Di])
return R({e: CRT_list([Di[e].lift() if e in Di else 0 for Di in Ds],
list(ps)).rational_reconstruction(prod(ps)) for e in supp})
def likely_support(fps):
supps = [set(fp.dict().keys()) for fp in fps]
supp_union = set.union(*supps)
supp = set()
half = len(supps) * 1/2 # exact amount does not matter
for e in supp_union:
k = 0
for s in supps:
if e in s:
k += 1
if k > half: # if exponent vector is in more than half the supports, we will probably want to keep it
supp.add(e)
break
return supp
ps = sorted(RecursivelyEnumeratedSet([next_prime(16000)], lambda p: [next_prime(p)])[:2800]) # many primes
Rps = [R.change_ring(GF(p)) for p in ps]
Ips = [I.change_ring(Rp).gens() for Rp in Rps]
CPUs = 16 # choose based on hardware
lift_one_modp.precompute(Ips, num_processes=CPUs)
assert all(type(r) != str for r in lift_one_modp.cache.values()) # otherwise, there was an exception
Fps = [lift_one_modp(Ip) for Ip in Ips]
supps = [likely_support(fps) for fps in zip(*Fps)]
# keep indices of primes for which all polynomials have the expected support
idxs = [j for j, Fp in enumerate(Fps) if all(set(fp.dict()) == supp for fp, supp in zip(Fp, supps))]
# restrict the results to those primes
ps1 = tuple(ps[j] for j in idxs)
Fps1 = [Fps[j] for j in idxs]
fpss = [tuple(fps) for fps in zip(*Fps1)]
modulus1 = prod(ps1)
floor(log(modulus1, 10)) # 10538 (must be large enough for rational reconstruction: this needs about twice as many digits as the numerators or denominators in the result have)
lift_poly_crt.precompute([(fps,) for fps in fpss], num_processes=CPUs)
assert all(type(r) != str for r in lift_poly_crt.cache.values()) # otherwise, there was an exception (if so, remove it from cache by lift_poly_crt.cache.clear())
F = [lift_poly_crt(fps) for fps in fpss] # the final result
```

7 | No.7 Revision |

`C`

modulo a prime is fast – it takes only a few seconds. Doing this for many primes can allow to lift the result to ℚ. In this case, the result `F`

(~~see attachment~~ uploaded here (31MB), for now, as the file seems too large for an attachment) is huge, though, with coefficients as large as 4826 and 4818 digits for numerator and denominator.

```
sage: F, C = load("lift_one_result_F_C.sobj") # saved with Sage 9.3
sage: sum(f*c for f,c in zip(F, C)) # 1 minute
1
```

I am not sure if this solution would be unique. Maybe there is a smaller one.

`A1`

and `A2`

, you could try to lift `C1`

and `C2`

individually, but this might make the result even larger.

When computing the lift of 1 mod p, most of the polynomials in the result have the same monomials ~~as ~~for their ~~support for varying p. ~~support, if p varies. This makes it very likely that the correct representative over ℚ has those same monomials as support and that the polynomials we computed mod p are the mod-p versions of the desired polynomials over ℚ.

`likely_support`

, we use a simple ad-hoc heuristic to find the exponents that probably occur in the ℚ-result: if a monomial occurs in 50% of the mod-p polynomials (i.e. modulo 50% of the primes or more), we expect that it should occur in the ℚ-result as well. The 50% threshold does not matter much. We could have picked 10% or 90% instead, I suppose.

Next, we use the Chinese remainder theorem `CRT_list`

to lift the coefficients to ℤ modulo the product of the primes and then `rational_reconstruction`

to obtain rational coefficients, both implemented in `lift_poly_crt`

. For `rational_reconstruction`

to work, the product of the primes must have about twice as many digits as the numerators and denominators in the result over ℚ (which I do not know how to determine in advance). In our case, we obtain a modulus with 10538 digits.

`cached_function.precompute`

in order to parallelize some computations because it is a bit less awkward to use than the `@parallel`

decorator, but also because it allows to increase the number of primes until rational reconstruction works, without redoing all the previous slow computations.

```
R1.<u0,u1,u2,v0,v1,v2,v3,v4,v5,v6>=PolynomialRing(QQ,10)
A1=[u0 + 7/5*u1 + 7/5*u2 - 4/125,
5*v0 + 5*v1 + 7*v3 + 7*v5 + 1/5,
25*v0^2 + 25*v1^2 + 35*v3^2 + 35*v5^2 - 4/5,
5*v0^3 + 5*v1^3 + 7*v3^3 + 7*v5^3 - v0^2 + 1/125,
5*v0*v1^2 + 5*v1*v2^2 + 7*v3*v4^2 + 7*v5*v6^2 + 1/125,
5*u0*v1 - v1^2 + 7*u1*v3 + 7*u2*v5 + 1/125,
5*v1 + 5*v2 + 7*v4 + 7*v6 + 1/5,
25*v0*v1 + 25*v1*v2 + 35*v3*v4 + 35*v5*v6 + 1/5,
5*v0^2*v1 + 5*v1^2*v2 + 7*v3^2*v4 + 7*v5^2*v6 - v1^2 + 1/125,
25*v1^2 + 25*v2^2 + 35*v4^2 + 35*v6^2 - 4/5,
5*v1^3 + 5*v2^3 + 7*v4^3 + 7*v6^3 - u0 + 1/125,
5*u0*v2 - v2^2 + 7*u1*v4 + 7*u2*v6 + 1/125]
I1=R1.ideal(A1)
C1=list(I1.groebner_basis())
R2.<y0,y1,y2,z0,z1,z2,z3,z4,z5,z6>=PolynomialRing(QQ,10)
A2=[y0 + 7/5*y1 + 7/5*y2 - 4/125,
5*z0 + 5*z1 + 7*z3 + 7*z5 + 1/5,
25*z0^2 + 25*z1^2 + 35*z3^2 + 35*z5^2 - 4/5,
5*z0^3 + 5*z1^3 + 7*z3^3 + 7*z5^3 - y0 + 1/125,
5*y0*z0 - z0^2 + 7*y1*z3 + 7*y2*z5 + 1/125,
5*z0*z1^2 + 5*z1*z2^2 + 7*z3*z4^2 + 7*z5*z6^2 - z1^2 + 1/125,
5*z1 + 5*z2 + 7*z4 + 7*z6 + 1/5,
25*z0*z1 + 25*z1*z2 + 35*z3*z4 + 35*z5*z6 + 1/5,
5*z0^2*z1 + 5*z1^2*z2 + 7*z3^2*z4 + 7*z5^2*z6 + 1/125,
5*y0*z1 - z1^2 + 7*y1*z4 + 7*y2*z6 + 1/125,
25*z1^2 + 25*z2^2 + 35*z4^2 + 35*z6^2 - 4/5,
5*z1^3 + 5*z2^3 + 7*z4^3 + 7*z6^3 - z2^2 + 1/125]
I2=R2.ideal(A2)
C2=list(I2.groebner_basis())
R.<u0,u1,u2,v0,v1,v2,v3,v4,v5,v6,y0,y1,y2,z0,z1,z2,z3,z4,z5,z6>=PolynomialRing(QQ,20)
C=C1+C2+[5*u0*z2 + 5*u1*z4 + 7*u2*z6 - u0 + 1/125] # the last one is the extra one
I=R.ideal(C)
########################################################################
@cached_function
def lift_one_modp(Ip):
Fp = Ip.ring()(1).lift(Ip.ideal())
print("completed p = %s" % Fp.ring().characteristic())
return Fp
@cached_function
def lift_poly_crt(fps):
ps = [fp.parent().characteristic() for fp in fps]
Ds = [fp.dict() for fp in fps]
supp = set([e for Di in Ds for e in Di])
return R({e: CRT_list([Di[e].lift() if e in Di else 0 for Di in Ds],
list(ps)).rational_reconstruction(prod(ps)) for e in supp})
def likely_support(fps):
supps = [set(fp.dict().keys()) for fp in fps]
supp_union = set.union(*supps)
supp = set()
half = len(supps) * 1/2 # exact amount does not matter
for e in supp_union:
k = 0
for s in supps:
if e in s:
k += 1
if k > half: # if exponent vector is in more than half the supports, we will probably want to keep it
supp.add(e)
break
return supp
ps = sorted(RecursivelyEnumeratedSet([next_prime(16000)], lambda p: [next_prime(p)])[:2800]) # many primes
Rps = [R.change_ring(GF(p)) for p in ps]
Ips = [I.change_ring(Rp).gens() for Rp in Rps]
CPUs = 16 # choose based on hardware
lift_one_modp.precompute(Ips, num_processes=CPUs)
assert all(type(r) != str for r in lift_one_modp.cache.values()) # otherwise, there was an exception
Fps = [lift_one_modp(Ip) for Ip in Ips]
supps = [likely_support(fps) for fps in zip(*Fps)]
# keep indices of primes for which all polynomials have the expected support
idxs = [j for j, Fp in enumerate(Fps) if all(set(fp.dict()) == supp for fp, supp in zip(Fp, supps))]
# restrict the results to those primes
ps1 = tuple(ps[j] for j in idxs)
Fps1 = [Fps[j] for j in idxs]
fpss = [tuple(fps) for fps in zip(*Fps1)]
modulus1 = prod(ps1)
floor(log(modulus1, 10)) # 10538 (must be large enough for rational reconstruction: this needs about twice as many digits as the numerators or denominators in the result have)
lift_poly_crt.precompute([(fps,) for fps in fpss], num_processes=CPUs)
assert all(type(r) != str for r in lift_poly_crt.cache.values()) # otherwise, there was an exception (if so, remove it from cache by lift_poly_crt.cache.clear())
F = [lift_poly_crt(fps) for fps in fpss] # the final result
```

8 | No.8 Revision |

`C`

modulo a prime is fast – it takes only a few seconds. Doing this for many primes can allow to lift the result to ℚ. In this case, the result `F`

(~~see attachment~~ uploaded here (31MB), for now, as the file seems too large for an attachment) is huge, though, with coefficients as large as 4826 and 4818 digits for numerator and denominator.

```
sage: F, C = load("lift_one_result_F_C.sobj") # saved with Sage 9.3
sage: sum(f*c for f,c in zip(F, C)) # 1 minute
1
```

I am not sure if this solution would be unique. Maybe there is a smaller one.

`A1`

and `A2`

, you could try to lift `C1`

and `C2`

individually, but this might make the result even larger.

When computing the lift of 1 mod p, most of the polynomials in the result have the same monomials for their support, if p varies. This makes it very likely that the correct representative over ℚ has those same monomials as support and that the polynomials we computed mod p are the mod-p versions of the desired polynomials over ℚ.

`likely_support`

, we use a simple ad-hoc heuristic to find the exponents that probably occur in the ℚ-result: if a monomial occurs in 50% of the mod-p polynomials (i.e. modulo 50% of the primes or more), we expect that it should occur in the ℚ-result as well. The 50% threshold does not matter much. We could have picked 10% or 90% instead, I suppose.

Next, we use the Chinese remainder theorem `CRT_list`

to lift the coefficients to ℤ modulo the product of the primes and then `rational_reconstruction`

to obtain rational coefficients, both implemented in `lift_poly_crt`

. For `rational_reconstruction`

to work, the product of the primes must have about twice as many digits as the numerators and denominators in the result over ℚ (which I do not know how to determine in advance). In our case, we obtain a modulus with 10538 digits.

`cached_function.precompute`

in order to parallelize some computations because it is a bit less awkward to use than the `@parallel`

decorator, but also because it allows to increase the number of primes until rational reconstruction works, without redoing all the previous slow computations.

```
R1.<u0,u1,u2,v0,v1,v2,v3,v4,v5,v6>=PolynomialRing(QQ,10)
A1=[u0 + 7/5*u1 + 7/5*u2 - 4/125,
5*v0 + 5*v1 + 7*v3 + 7*v5 + 1/5,
25*v0^2 + 25*v1^2 + 35*v3^2 + 35*v5^2 - 4/5,
5*v0^3 + 5*v1^3 + 7*v3^3 + 7*v5^3 - v0^2 + 1/125,
5*v0*v1^2 + 5*v1*v2^2 + 7*v3*v4^2 + 7*v5*v6^2 + 1/125,
5*u0*v1 - v1^2 + 7*u1*v3 + 7*u2*v5 + 1/125,
5*v1 + 5*v2 + 7*v4 + 7*v6 + 1/5,
25*v0*v1 + 25*v1*v2 + 35*v3*v4 + 35*v5*v6 + 1/5,
5*v0^2*v1 + 5*v1^2*v2 + 7*v3^2*v4 + 7*v5^2*v6 - v1^2 + 1/125,
25*v1^2 + 25*v2^2 + 35*v4^2 + 35*v6^2 - 4/5,
5*v1^3 + 5*v2^3 + 7*v4^3 + 7*v6^3 - u0 + 1/125,
5*u0*v2 - v2^2 + 7*u1*v4 + 7*u2*v6 + 1/125]
I1=R1.ideal(A1)
C1=list(I1.groebner_basis())
R2.<y0,y1,y2,z0,z1,z2,z3,z4,z5,z6>=PolynomialRing(QQ,10)
A2=[y0 + 7/5*y1 + 7/5*y2 - 4/125,
5*z0 + 5*z1 + 7*z3 + 7*z5 + 1/5,
25*z0^2 + 25*z1^2 + 35*z3^2 + 35*z5^2 - 4/5,
5*z0^3 + 5*z1^3 + 7*z3^3 + 7*z5^3 - y0 + 1/125,
5*y0*z0 - z0^2 + 7*y1*z3 + 7*y2*z5 + 1/125,
5*z0*z1^2 + 5*z1*z2^2 + 7*z3*z4^2 + 7*z5*z6^2 - z1^2 + 1/125,
5*z1 + 5*z2 + 7*z4 + 7*z6 + 1/5,
25*z0*z1 + 25*z1*z2 + 35*z3*z4 + 35*z5*z6 + 1/5,
5*z0^2*z1 + 5*z1^2*z2 + 7*z3^2*z4 + 7*z5^2*z6 + 1/125,
5*y0*z1 - z1^2 + 7*y1*z4 + 7*y2*z6 + 1/125,
25*z1^2 + 25*z2^2 + 35*z4^2 + 35*z6^2 - 4/5,
5*z1^3 + 5*z2^3 + 7*z4^3 + 7*z6^3 - z2^2 + 1/125]
I2=R2.ideal(A2)
C2=list(I2.groebner_basis())
R.<u0,u1,u2,v0,v1,v2,v3,v4,v5,v6,y0,y1,y2,z0,z1,z2,z3,z4,z5,z6>=PolynomialRing(QQ,20)
C=C1+C2+[5*u0*z2 + 5*u1*z4 + 7*u2*z6 - u0 + 1/125] # the last one is the extra one
I=R.ideal(C)
########################################################################
@cached_function
def lift_one_modp(Ip):
Fp = Ip.ring()(1).lift(Ip.ideal())
print("completed p = %s" % Fp.ring().characteristic())
return Fp
@cached_function
def lift_poly_crt(fps):
ps = [fp.parent().characteristic() for fp in fps]
modulus = prod(ps)
Ds = [fp.dict() for fp in fps]
supp = set([e for Di in Ds for e in Di])
return R({e: CRT_list([Di[e].lift() if e in Di else 0 for Di in Ds],
```~~list(ps)).rational_reconstruction(prod(ps)) ~~list(ps)).rational_reconstruction(modulus) for e in supp})
def likely_support(fps):
supps = [set(fp.dict().keys()) for fp in fps]
supp_union = set.union(*supps)
supp = set()
half = len(supps) * 1/2 # exact amount does not matter
for e in supp_union:
k = 0
for s in supps:
if e in s:
k += 1
if k > half: # if exponent vector is in more than half the supports, we will probably want to keep it
supp.add(e)
break
return supp
ps = sorted(RecursivelyEnumeratedSet([next_prime(16000)], lambda p: [next_prime(p)])[:2800]) # many primes
Rps = [R.change_ring(GF(p)) for p in ps]
Ips = [I.change_ring(Rp).gens() for Rp in Rps]
CPUs = 16 # choose based on hardware
lift_one_modp.precompute(Ips, num_processes=CPUs)
assert all(type(r) != str for r in lift_one_modp.cache.values()) # otherwise, there was an exception
Fps = [lift_one_modp(Ip) for Ip in Ips]
supps = [likely_support(fps) for fps in zip(*Fps)]
# keep indices of primes for which all polynomials have the expected support
idxs = [j for j, Fp in enumerate(Fps) if all(set(fp.dict()) == supp for fp, supp in zip(Fp, supps))]
# restrict the results to those primes
ps1 = tuple(ps[j] for j in idxs)
Fps1 = [Fps[j] for j in idxs]
fpss = [tuple(fps) for fps in zip(*Fps1)]
modulus1 = prod(ps1)
floor(log(modulus1, 10)) # 10538 (must be large enough for rational reconstruction: this needs about twice as many digits as the numerators or denominators in the result have)
lift_poly_crt.precompute([(fps,) for fps in fpss], num_processes=CPUs)
assert all(type(r) != str for r in lift_poly_crt.cache.values()) # otherwise, there was an exception (if so, remove it from cache by lift_poly_crt.cache.clear())
F = [lift_poly_crt(fps) for fps in fpss] # the final result

Copyright Sage, 2010. Some rights reserved under creative commons license. Content on this site is licensed under a Creative Commons Attribution Share Alike 3.0 license.