Ask Your Question

Revision history [back]

click to hide/show revision 1
initial version

answered 2021-05-22 23:52:48 +0200

mwageringel gravatar image

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.

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.

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

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.

Update

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

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.

Update

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

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.

Update

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). In our case, we 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

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.

Update

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 ℚ.

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). In our case, we 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

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.

Update

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 ℚ.

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). In our case, we 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]
    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.
about | faq | help | privacy policy | terms of service
Powered by Askbot version 0.7.59