Zippel GCD algorithm

[Brent W. Baccala]

Hi -

I'm continuing the work I mentioned two or three weeks ago to improve FLINT's memory utilization.  I'm focusing right now on the Zippel GCD algorithm, and I'm having a lot of trouble understanding the code.  The only reference I've got is Zippel's 1979 paper, "Probabilistic Algorithms for Sparse Polynomials".

Are there any better references and/or documentation for this code?  The material in doc/source seems to document the public API without much description of the internals.

For example, I'm trying to figure the difference between the three functions nmod_mpolyu_gcd[p,s,m]_zippel which only differ in the one letter p, s, and m, and understand how they work together to implement the GCD algorithm.

Thanks.

    agape

    brent

[Bill Hart]

Hi Brent,

Thanks for your work on this. The guy who wrote this code will
hopefully reply, however it may take a few days to get an answer on
this. Hold tight, we aren't ignoring you.

Bill.

> --
>
> ---
> You received this message because you are subscribed to the Google Groups "flint-devel" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to flint-devel...@googlegroups.com.
> To view this discussion on the web, visit https://groups.google.com/d/msgid/flint-devel/4c748f9b-6ad6-4722-84e5-ec1675818914n%40googlegroups.com.

[da sh]

The old documentation is in mpoly/docs. I guess you will have to compile the latex file yourself. I have an update to this documentation, but it is not on trunk yet.

However, I do not like this part of the code and will be cleaning it up and simplifying it soon.

The only mpoly module with a clean implementation is fmpz_mod_mpoly.

[da sh]

Brent,

Could you join us for Flint Fridays?

[Fredrik Johansson]

I wonder if we should have a "Flint: algorithms and proofs" document similar to https://www.mpfr.org/algorithms.pdf, separate from the API documentation.

Granted, it would be years of work to document all the algorithms in Flint, but these notes for the multivariate polynomials could be a good place to start...

Fredrik

To view this discussion on the web, visit https://groups.google.com/d/msgid/flint-devel/aea3c902-75c4-4c7d-9c91-e1674c0ca4dcn%40googlegroups.com.

[fieker]

On Mon, May 17, 2021 at 02:45:08PM +0200, Fredrik Johansson wrote:
> I wonder if we should have a "Flint: algorithms and proofs" document
> similar to https://www.mpfr.org/algorithms.pdf, separate from the API
> documentation.
>
> Granted, it would be years of work to document all the algorithms in Flint,
> but these notes for the multivariate polynomials could be a good place to
> start...
>
> Fredrik

From the outside:
YEP!!!
plus comments in the code to explain stuff, not "just" the API...

enormous amout of work.

How does Calcium stand on this front?

> > <https://groups.google.com/d/msgid/flint-devel/aea3c902-75c4-4c7d-9c91-e1674c0ca4dcn%40googlegroups.com?utm_medium=email&utm_source=footer>
> > .

> >
>
> --
>
> ---
> You received this message because you are subscribed to the Google Groups "flint-devel" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to flint-devel...@googlegroups.com.

> To view this discussion on the web, visit https://groups.google.com/d/msgid/flint-devel/CAJdUXTL0GNx%2BtQaphC9_8JdKsXcdHKw6dKLtxDY0ZQj79xuchg%40mail.gmail.com.

[Fredrik Johansson]

On Mon, May 17, 2021 at 3:19 PM fieker <fie...@mathematik.uni-kl.de> wrote:

On Mon, May 17, 2021 at 02:45:08PM +0200, Fredrik Johansson wrote:
> I wonder if we should have a "Flint: algorithms and proofs" document
> similar to https://www.mpfr.org/algorithms.pdf, separate from the API
> documentation.
>
> Granted, it would be years of work to document all the algorithms in Flint,
> but these notes for the multivariate polynomials could be a good place to
> start...
>
> Fredrik

From the outside:
 YEP!!!
 plus comments in the code to explain stuff, not "just" the API...

enormous amout of work.

How does Calcium stand on this front?

Also an enormous amount of work, but something I have thought about. A large part of the problem is just to find the right "form" for such documentation...

Fredrik

[fieker]

I'd start with code-comments... side calulations in comments, links to
the source, translation table if notation is different.
A perfect, external document is good, however, code-comments are
easier...

Claus
>
> Fredrik

>
> --
>
> ---
> You received this message because you are subscribed to the Google Groups "flint-devel" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to flint-devel...@googlegroups.com.

> To view this discussion on the web, visit https://groups.google.com/d/msgid/flint-devel/CAJdUXTJDQz7y_1Hy6z%2BdgB%2B6fZcY_xK%2BOBFzXSJib4vVsps61g%40mail.gmail.com.

[Brent W. Baccala]

Dear Fredrik,

I'll join you for Flint Fridays!  Could you send me an access code, please?

It's 2 PM CEST (UTC +2), right?

    agape

    brent

--

---
You received this message because you are subscribed to a topic in the Google Groups "flint-devel" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/flint-devel/ICuIuVe6xno/unsubscribe.
To unsubscribe from this group and all its topics, send an email to flint-devel...@googlegroups.com.
To view this discussion on the web, visit https://groups.google.com/d/msgid/flint-devel/CAJdUXTJDQz7y_1Hy6z%2BdgB%2B6fZcY_xK%2BOBFzXSJib4vVsps61g%40mail.gmail.com.

[Brent W. Baccala]

Good evening -

Nice seeing you all this morning!

It was real helpful for me, because Dan, in particular, asked me some tough questions that I need to think about.

Definitely need to instrument this code to try and estimate how big these polynomials will be, to figure how practical my calculations are.

It'll take me a few days to put the Sage code together into something that makes sense to anybody but me.  You know how it goes :-).

In the meantime, I'm attaching a text file (at the end of this email) that contains my notes on the FLINT implementation of the ZIppel GCD algorithm.  I'm not sure just want to do with it other than share it.  It could ultimately work its way into an algorithms document.

    agape

    brent

Most common polynomial types:

fmpz_mpoly_t:        Z[a..d]

fmpz_mpolyu_t:       Z[a..d][x]      (x-exponents are long's)

fmpz_mpoly_univar_t: Z[a..d][x]      (x-exponents are fmpz's)

nmod_mpoly_t:        Fp[a..d]

nmod_mpolyu_t:       Fp[a..d][x]     (exps field are long's; coefficients are nmod_mpoly_t)

nmod_mpolyn_t:       Fp[v][a..d]     (exps field are mpoly-style; coefficients are nmod_poly_t)

nmod_mpolyun_t:      Fp[v][a..d][x]  (exps field are long's (the x-exponents); coefficients are nmod_mpolyn_t)

fmpz_mpoly_gcd(fmpz_mpoly_t G,

               const fmpz_mpoly_t A,

               const fmpz_mpoly_t B,

               const fmpz_mpoly_ctx_t ctx)

Accepts two multivariate polynomials over Z and computes their GCD.

Algorithm:

- if either polynomial is zero, GCD is the other one, maybe with sign flipped to ensure a positive leading coefficient

- if either polynomial has length 1, special case for monomial gcd

  - exponents are min exponents over both polynomials, cofficient is content of coefficients

- find min and max exponents for all variables in each polynomial

- if all variables have the same range (max-min) in both polynomials, check for monomial cofactors

  - coefficients must all have the same "cross ratio", i.e, (Ai / A0 == Bi / B0) for all i

  - exponents must all match a similar condition (x^a_i / x^a_0 == x^b_i / x^b_0) for all i

  - if coefficient and exponent conditions are met, then we have monomial cofactors

- if any min is greater than zero, factor out that variable^min

  (in principle; in fact we just record the min and use it later without rewriting the polynomial here)

- compute strides (GCD of each variable's exponents, after dividing out minimum exponents)

- if any stride > 1 and common between the two polynomials, reduce that variable's degree

  (in principle; in fact we just record the stride and use it later without rewriting the polynomial here)

- find which variables are in both polynomials

- if no variable in both, GCD is simply the GCD of the integer contents

- if only one variable is in both, and min==max for all other variables, special case for univariate GCD

- if a variable is in only one polynomial or the other, special case for missing variable

  - isolate the variable and compute GCD between all of its factors, and the other polynomial

  - multiply that GCD by variables raised to the minimum exponents we (in principle) factored out earlier

- all variables are now missing from both or present in both, and there are at least two present in both

- compute deflated degrees, which are (max-min)/stride,

  and the GCD's deflated upper degree bound, which is min of the two deflated degrees, for each variable

- if all GCD deflated upper degree bounds are zero, GCD is trivial (GCD of the integer contents, times x^min_deg)

- if all of A's min exponents are <= B's min exponents, and all of A's deflated degrees are == G's upper bounds,

  then check to see if B/A is exact and A is the GCD

- if all of B's min exponents are <= A's min exponents, and all of B's deflated degrees are == G's upper bounds,

  then check to see if A/B is exact and B is the GCD

- get fancy!

- if we've only got two variables, and Brown is usable, try Brown

- estimate time complexity of Brown and Berklekamp-Massey

- try Brown or Berklekamp-Massey, in order, depending on which is estimated to be faster

- fall back on Zippel

static int _try_zippel(

    fmpz_mpoly_t G,

    const fmpz_mpoly_t A,

    const fmpz_mpoly_t B,

    const mpoly_gcd_info_t I,

    const fmpz_mpoly_ctx_t ctx)

- the mpoly_gcd_info structure already contains a permutation of the variables that effectively

  selects one (or is it two?) variables that will be treated special (how were those variables picked?)

- deflate to fmpz_mpolyu: ZZ[a..d][x]

- for each polynomial, extract content: GCD of coefficients in ZZ[a..d], and divide it out

- compute the GCD using fmpz_mpolyu_gcdm_zippel

- inflate the answer back into our original ring and multiply by the content GCD

int fmpz_mpolyu_gcdm_zippel(

    fmpz_mpolyu_t G,

    fmpz_mpolyu_t Abar,

    fmpz_mpolyu_t Bbar,

    fmpz_mpolyu_t A,

    fmpz_mpolyu_t B,

    const fmpz_mpoly_ctx_t ctx,

    mpoly_zipinfo_t zinfo,

    flint_rand_t randstate)

At this point, the polynomials are in Z[a..d][x] (i.e, we've isolated a "privileged" variable)

- compute gamma (in ZZ), GCD of leading coefficients

- set GCD's degree bound (in x) to the minimum of the two polynomial's degrees

- find an upper bound on the coefficients of the GCD

    - Landau-Mignotte Bound for univariate polynomials (https://mathworld.wolfram.com/Landau-MignotteBound.html)

    - Lang and Mahler give bounds for multivariates (https://mathoverflow.net/questions/151556)

    - FLINT uses Landau-Mignotte combined with its standard Kronecker substitution

- pick a (outer) prime p

- reduce gamma mod p, and reduce both polynomials mod p

- try a different prime if either polynomial or gamma reduced to zero

- compute "skeleton" using nmod_mpolyu_gcdp_zippel (Zipeel's dense interpolation Algorithm D)

- try a different prime if the degree of the new GCD exceeds our degree bound

- update degree bound to be the degree of the computed GCD

- multiply GCD by gamma/lead(Gp) mod p, to make leading coeff of leading coeff = gamma (GCD of input poly's leading coeffs)

- if there's no x-dependency in the GCD (length==1 and exps[0]==0), return 1 as the result

- pick another (inner) prime

- reduce gamma and both polynomials, and loop if any of them reduced to zero

- compute nmod_mpoly_gcds_zippel using result of first calculation as Gform (Zippel's sparse interpolation Algorithm S)

- try inner prime again if scales_not_found, or eval_point_not_found, or eval_gcd_deg_too_high

- try outer prime again if form_main_degree_too_high or form_wrong or no_solution

- try inner prime again if leading coefficient of GCD is zero (how can this be?)

- multiple GCD by gamma/leadcoeff(Gp) mod p (inner prime)

- changed = fmpz_mpolyu_interp_mcrt_p (crt = Chinese Remainder Theorem)

- multiply modulus by inner prime

- if changed, loop to either outer prime (coeffbits > coeffbitbound) or inner prime otherwise

- if G=H/(content H) doesn't divide A and B, loop back to inner prime

- G is our answer

int nmod_mpolyu_gcdp_zippel(

    nmod_mpolyu_t G,

    nmod_mpolyu_t Abar,

    nmod_mpolyu_t Bbar,

    nmod_mpolyu_t A,

    nmod_mpolyu_t B,

    slong var,

    const nmod_mpoly_ctx_t ctx,

    mpoly_zipinfo_t zinfo,

    flint_rand_t randstate)

Zippel's dense interpolation Algorithm D

- special cases for univariate (either Euclidian or half-GCD, depending on length of longest poly) and bivariate

- pull out 'var' to get nmod_mpolyun_t's (An and Bn) in Fp[v][a..d][x]

- compute the GCD of the leading coefficients, producing a polynomial in Fp[v]

- bound the number of images required

- modulus/tempmod are in Fp[v]; set modulus to '1' and tempmod to 'v' (the coefficient variable) initially

- interpolant H is an mpolyun_t; set to zero initially

- pick an evaluation point (a single number for 'v')

- evaluate the input polys and the leading GCD at that point

- make sure it doesn't kill both lc(A) and lc(B), or A, or B

- recurse into self (nmod_mpolyu_gcdp_zippel) for the next variable to compute GCD w.r.t. other variables

- if the recursed GCD (Geval) is 1, then the result is 1 (return success)

- if degree modulus > 0, either try again if Geval's degree is greater than interpolants, or set modulus to 1 otherwise

- update interpolant (originally zero)

   - interpolant is an mpolyun_t

   - use nmod_mpolyun_interp_lift_sm_mpolyu if degree modulus == 0

     Geval (nmod_mpolyu_t Fp[a..d][x]) -> H (nmod_mpolyun_t Fp[v][a..d][x])

     operate on each x-coefficient independently

     set, ignoring previous value of interpolant; alpha and modulus both ignored

     map Fp -> Fp[v] as constants

   - use nmod_mpolyun_interp_crt_sm_mpolyu if degree modulus > 0

     H = H + modulus*(Geval - H(alpha))

     if this didn't change H, check to see if H divides A and B, if so, it's our answer, otherwise try another eval point

- multiple modulus by (v + (p - alpha))

- compute nmod_mpoly_gcds_zippel using result of first calculation as Gform (Zippel's sparse interpolation Algorithm S)

- fail if form_wrong or no_solution

- try another S-algorithm interpolation if scales_not_found, or eval_point_not_found, or eval_gcd_deg_too_high

- go back to D-algorithm if form_main_degree_too_high

nmod_gcds_ret_t nmod_mpolyu_gcds_zippel(nmod_mpolyu_t G,

               nmod_mpolyu_t A, nmod_mpolyu_t B,  nmod_mpolyu_t f, slong var,

               const nmod_mpoly_ctx_t ctx, flint_rand_t randstate, slong * degbound)

Zippel's sparse interpolation (S algorithm)

- (previous GCD is skeleton 'f')

- make sure it doesn't kill both lc(A) and lc(B), or A, or B

- if the skeleton has one term in x,

    - if it has multiple terms in its Fp[a..d] coefficient, return scales_not_found

    - set G = skeleton, then set its leading coefficient = 1

    - does G divide both A and B?  if so, return success; otherwise, return form_wrong

- based on length of skeleton and length of each Fp[a..d] coefficient, compute number of images needed

- for each remaining variable, pick random values

- for each image required, compute using (a,b,c,d), (a^2,b^2,c^2,d^2), (a^3,b^3,c^3,d^3), etc. as the coordinates

  this lets each monomial be computed by computing the power of its original value at (a,b,c,d)

- evaluate at those variables to get univariate polynomials and compute their GCD

- if that GCD's degree is greater than the skeleton's degree, try one more evaluation point, then return gcd_deg_too_high

- if that GCD's degree is less than the skeleton's degree, return form_main_degree_too_high (and update degbound)

- if any of the elements of the GCD correspond to a missing exponent in the skeleton, return form_wrong

- assemble the coefficients of the univariate GCD into the columns of a matrix (W), but only the rows

     corresponding to actual x-exponents in the skeleton are used

- build a matrix for each non-zero exponent in the skeleton

     a row for each image (coordinate)

     as many columns as coefficients for that exponent in the skeleton

- convert them to reduced row Echelon form

- if any of them are rank-deficient, pick a new random coordinate

- concatenate them column-wise into one big matrix

- if it has no solution (nullity = 0), trial form must be wrong.  return gcds_form_wrong

- if it's underdetermined (nullity > 1), try once to pick a new evaluation point, then return scales_not_found

- solve it using vandsolve (another advantage of picking coordinates as powers - it's a Vandemonde matrix)

- the solution gives the coefficients for the skeleton