Added functionality and documentation for calculating roots of univariate polynomials #935

monien · 2021-04-11T21:06:26Z

Added functionality and documentation for calculating roots of univariate polynomials over the integers (fmpz_poly_t) over finite fields (fq_t), p-adic fields (padic_t) and q-adic fields (qadic_t).

…iate polynomials over the integers (fmpz_poly) over finite fields, p-adic fields and q-adic fields.

p-adic and q-adic roots.

fredrik-johansson · 2021-04-12T15:43:56Z

What is the rationale for making this a new module instead of putting root-finding methods in the fq, padic and qadic modules?

monien · 2021-04-12T15:58:39Z

Sage has the method roots attached to integer polynomials with input ring = GF, Zp .. Qq, where roots over Qq is not working correctly at the moment. Either way is fine with me.

…

On 12. Apr 2021, at 17:44, Fredrik Johansson ***@***.*** ***@***.***>> wrote: What is the rationale for making this a new module instead of putting root-finding methods in the fq, padic and qadic modules? — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#935 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAZTMAX2AAZTOV2IYUTQX5TTIMIM3ANCNFSM42YCE5TQ>.

monien · 2021-04-12T15:58:39Z

Sage has the method roots attached to integer polynomials with input ring = GF, Zp .. Qq, where roots over Qq is not working correctly at the moment. Either way is fine with me.

…

On 12. Apr 2021, at 17:44, Fredrik Johansson ***@***.*** ***@***.***>> wrote: What is the rationale for making this a new module instead of putting root-finding methods in the fq, padic and qadic modules? — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#935 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAZTMAX2AAZTOV2IYUTQX5TTIMIM3ANCNFSM42YCE5TQ>.

wbhart · 2021-04-12T16:50:20Z

These are finding padic, qadic or finite field roots of an integer polynomial (Z[x]). Are you sure you want it that way around @fredrik-johansson ?

fredrik-johansson · 2021-04-12T16:59:37Z

Well, someone who comes across this module might be disappointed that it doesn't do roots in R and C :-)

Does @tthsqe12 have an opinion on the interface?

tthsqe12 · 2021-04-12T18:28:31Z

I don't have much of an opinion on the interface. But I do have questions on the code. First of all, there is no test code yet, so we cannot be sure what the functions are supposed to do. The documentation is not specific on what roots it finds and to what precision, ect.
Here is what I have found so far: an example that words, a mundane crash, and weird root output (the root is listed as multiple and is only correct to precision 2, and the other root is missed):

#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>
#include "flint.h"
#include "fmpz.h"
#include "fmpz_poly.h"
#include "ulong_extras.h"
#include "padic.h"
#include "fmpz_poly_roots.h"

int
main(void)
{
    int i, result;
    FLINT_TEST_INIT(state);

    flint_printf("roots....\n");
    fflush(stdout);

    {
        fmpz_poly_t a, b;
        fmpz_t p;
        padic_ctx_t Q5;
        fmpz_poly_roots_padic_t roots;

        fmpz_init_set_ui(p, 5);
        padic_ctx_init(Q5, p, 20, 20, PADIC_SERIES);
        fmpz_poly_init(a);

        flint_printf("------- roots works ---------\n");
        fmpz_poly_set_str(a, "3  3 1 1");
       
        flint_printf("a: ");
        fmpz_poly_print_pretty(a, "x");
        flint_printf("\n");

        fmpz_poly_roots_padic_init2(roots, 2);
        fmpz_poly_roots_padic(roots, a, Q5);
        fmpz_poly_roots_padic_print(roots, Q5);
        fmpz_poly_roots_padic_clear(roots);

        flint_printf("------- roots crashes ---------\n");
        fmpz_poly_set_str(a, "3  -30 1 1");
       
        flint_printf("a: ");
        fmpz_poly_print_pretty(a, "x");
        flint_printf("\n");

        fmpz_poly_roots_padic_init2(roots, 2);
        /* this crashes because an fq root is zero so its 'coeffs' cannot be dereferenced */
        /*fmpz_poly_roots_padic(roots, a, Q5);
        flint_printf("roots:\n");
        fmpz_poly_roots_padic_print(roots, Q5);*/
        fmpz_poly_roots_padic_clear(roots);

        flint_printf("------- roots gives wrong answer? ---------\n");
        fmpz_poly_set_str(a, "3  6 -7 1");
       
        flint_printf("a: ");
        fmpz_poly_print_pretty(a, "x");
        flint_printf("\n");

        fmpz_poly_roots_padic_init2(roots, 2);
        fmpz_poly_roots_padic(roots, a, Q5);
        flint_printf("roots:\n");
        fmpz_poly_roots_padic_print(roots, Q5);
        fmpz_poly_roots_padic_clear(roots);

        flint_printf("roots mod 5^20 reported by fmpz_mod_poly_roots:\n");

        {
            slong i;
            fmpz_mod_ctx_t mod55;
            fmpz_mod_poly_t f;
            fmpz_mod_poly_factor_t fac;
            fmpz_t t;
            fmpz_factor_t fac55;

            fmpz_init(t);
            fmpz_pow_ui(t, p, 20);
            fmpz_mod_ctx_init(mod55, t);
            fmpz_mod_poly_init(f, mod55);
            fmpz_mod_poly_factor_init(fac, mod55);
            fmpz_factor_init(fac55);
            fmpz_factor(fac55, fmpz_mod_ctx_modulus(mod55));

            fmpz_mod_poly_set_fmpz_poly(f, a, mod55);
            fmpz_mod_poly_roots_factored(fac, f, 1, fac55, mod55);

            for (i = 0; i < fac->num; i++)
            {
                fmpz_mod_neg(t, fac->poly[i].coeffs + 0, mod55);
                fmpz_print(t);
                flint_printf(" %wd\n", fac->exp[i]);
            }

            fmpz_factor_clear(fac55);
            fmpz_mod_poly_factor_clear(fac, mod55);
            fmpz_mod_poly_clear(f, mod55);
            fmpz_mod_ctx_clear(mod55);
            fmpz_clear(t);
        }

        fmpz_poly_clear(a);
        padic_ctx_clear(Q5);
        fmpz_clear(p);
    }

    FLINT_TEST_CLEANUP(state);
    
    flint_printf("PASS\n");
    return 0;
}

with output:

------- roots works ---------
a: x^2+x+3
3 + 1*5^1 + 1*5^2 + 3*5^3 + 4*5^4 + 4*5^5 + 3*5^7 + 1*5^8 + 1*5^9 + 4*5^11 + 2*5^13 + 1*5^15 + 1*5^17 + 4*5^18 + 1*5^19 1
1 + 3*5^1 + 3*5^2 + 1*5^3 + 4*5^6 + 1*5^7 + 3*5^8 + 3*5^9 + 4*5^10 + 4*5^12 + 2*5^13 + 4*5^14 + 3*5^15 + 4*5^16 + 3*5^17 + 3*5^19 1
------- roots crashes ---------
a: x^2+x-30
------- roots gives wrong answer? ---------
a: x^2-7*x+6
roots:
1 + 3*5^1 + 2*5^2 + 2*5^3 + 2*5^4 + 2*5^5 + 2*5^6 + 2*5^7 + 2*5^8 + 2*5^9 + 2*5^10 + 2*5^11 + 2*5^12 + 2*5^13 + 2*5^14 + 2*5^15 + 2*5^16 + 2*5^17 + 2*5^18 + 2*5^19 2
roots mod 5^20 reported by fmpz_mod_poly_roots:
1 1
19073486328126 1
38146972656251 1
57220458984376 1
76293945312501 1
6 1
19073486328131 1
38146972656256 1
57220458984381 1
76293945312506 1
PASS

tthsqe12 · 2021-04-12T18:33:02Z

maybe @thofma has an opinion on how this kind of stuff is supposed to work

tthsqe12 · 2021-04-12T19:43:52Z

So, the problem here for Qp is: the code finds roots mod p then attemps to lift them using newton's method on f or f', f'' ect if there are multiple roots mod p. The problem is this last case of multiple roots.

tthsqe12 · 2021-04-13T07:11:18Z

Im sure that @fredrik-johansson would have an opinion here, and I am going to give mine:
The roots function should behave as if it calculated the roots exactly (with multiplicities) and then truncated these exact roots to the precision of the padic_ctx_t. This means in particular for Q5:

The roots of (x-1)*(x-6) should reported as [[1 + O(5), 1], [1 + O(5), 1]] in precision 1 output and as [[1 + O(5^2), 1], [6 + O(5^2), 1]] in precision 2 output.
I am not sure how the padic_t handles relative precision, but there should be two simple roots output for x^2+x-30 in any precision.
5x-1 should have one simple root in any precision.
x^2-5 should have no roots (in any precision).
For the root of x (i.e. zero), the particular precision model used by the padic_t, of which I know nothing, is probably relevent.

This makes the function harder to code, but has the benefit of giving the function a well-defined mathematical prescription. This also means the output of the function is a priori reliable.

tthsqe12 · 2021-04-13T07:37:40Z

In order to accomplish this exact root stuff for f(x) in Z[x], we can assume:

f(0) != 0
f is squarefree over Qp by a squarefree factorization over Z
The leading coefficient of f is not divisible by p by a rescale of f.

At this point each exact root of f is identified by the following data:

An integer k > 0
An x0 in Z/(p^k Z) with f(x0) = 0 in Z/(p^k Z) and f'(x0) != 0 in Z/(p^floor((k+1)/2) Z)

Furthermore, it easy to see when (x0, k) and some other (x0', k') represent the same root, so we can calculate the roots mod p^k such that f' has at most half valuation at each of these roots mod p^k, and then delete duplicates.

monien · 2021-04-13T08:29:04Z

I will fix the embarrassing problems which of course did not occur in my generic test cases. I still would like to have a definite vote of you guys on either keeping the exact root stuff in one place (fmpz_poly) or distribute it over all fq, paid and qadic directories. Anybody a strong opinion? Cheers, H

…

On 13. Apr 2021, at 09:37, tthsqe12 ***@***.*** ***@***.***>> wrote: In order to accomplish this exact root stuff for f(x) in Z[x], we can assume: f(0) != 0 f is squarefree over Qp by a squarefree factorization over Z The leading coefficient of f is not divisible by p by a rescale of f. At this point each exact root of f is identified by the following data: An integer k0 > 0 An x0 in Z/p^kZ with f(x0) = 0 in Z/p^kZ and f'(x0) != 0 in Z/p^kZ Furthermore, it easy to see when (x0, k0) and some other (x0', k0') represent the same root, so we can calculate the roots mod p^k such that f' is nonzero at each of these roots mod p^k, and then delete duplicates. — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#935 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAZTMAR2LJGYN2LG2BG4EU3TIPYFHANCNFSM42YCE5TQ>.

tthsqe12 · 2021-04-13T10:24:34Z

Here is my opinion on the module thingy. In other modules we have a T_poly_t and and T_poly_factor_t and we put the root finding code in the factor module and return the roots ri and multiplicities ei as a "factorization" prod_i (x-ri)^ei. This is not great but was done to cut down on the number of flint types that need to be wrapped in foreign languages. The overhead of each x-ri is usually negligible compared to the root finding time.
At the moment we do not have qadic_poly nor do we have {p|q}adic_poly_factor_t.
I do not imagine that the factoring functions {p|q}adic_poly_factor will/can ever be implemented. Therefore, we are left with factoring fmpz_poly's over {p|q}adic. This root finding stuff does degree one factors. If we ever allow factors of arbitrary degree (which seems likely in the near future), there will definitely have to be a {p|q}adic_poly_factor_t to hold the answer. If the root finder and the hypothetical general factorer use different return types, it will be inconsistent with the other modules.

long summary: I don't care where the files are placed, but maybe the return type use for fmpz_poly -> {p|q}adic roots should be general enough to hold arbitrary degree factorizations to consistently support future extensions to the functions. Or, maybe roots in this case are a special enough case to warrant a special type... Arb has arb_fmpz_poly_complex_roots which takes squarefree exact input and just a raw pointer to the output. I don't think this is a great idea here because we don't know the number of roots in advance...

fredrik-johansson · 2021-04-13T11:02:52Z

Life will maybe become easier if we introduce non-underscore vector types with length management (fmpz_vec_t, etc). I've been doing this in Calcium and so far it "feels" right.

tthsqe12 · 2021-04-13T11:44:16Z

Ok, but do you have any issues with my prescription for the behavior of the function? There can be controversy here...

fredrik-johansson · 2021-04-13T12:33:52Z

No, I agree with everything you said.

monien · 2021-04-24T05:01:30Z

Adding roots functionality to fmpz_poly_factor will not work.

One would need to include "fq.h" and "quadic.h" in "fmpz_poly_factor.h" which introduces circular dependencies in the include tree.

tthsqe12 · 2021-04-24T09:25:38Z

I have no good ideas on the interface. Maybe @fredrik-johansson or @wbhart can chime in.

wbhart · 2021-05-03T13:11:58Z

The way we normally handle the circular reference issue is to put the signatures in both files, but as a comment in the one that can't support the actual signatures. Usually these are placed at the very bottom of the .h file, with a note in the comment explaining where the real signatures are.

This only affects the signature. The implementation remains in the module with the comments rather than the one with the real signatures. See charpoly and minpoly for examples of where this is already done, see for example nmod_mat.h/nmod_poly.h.

MGessinger · 2021-08-02T14:18:12Z

If I may, I believe the algorithm used to compute padic roots is far from optimal. The roots as computed over the finite field Fp are exactly the "first" (i.e. least significant) digits of the padic expansion. The Euler iteration does not respect that digit, so after just one iteration we lose that piece of info. It also doesn't "refine" roots that agree in low precision, but differ in high precision.

Instead, the following method computes the "digits" one at a time:
Suppose we call the polynomial q, and the prime p (i.e. we want to find roots of q over Qp).

Compute the roots with multiplicity (ri, mi) of q over Z/pZ.
For each pair (ri, mi) define qi(x) = q(p*x+ri)/p^mi.
Return to step 1 with qi instead of q
This iteration will produce a series of digits {ri}, which correspond exactly to the digits of the root r. Up to a small caveat (that I will come back to in a moment), this algorithm converges exponentially (w.r.t. the padic metric).

To further justify why this method might be of interest, allow me to address the specific properties that @tthsqe12 listed above:

Unfortunately, this method doesn't distinguish roots that only differ after a certain precision. Those roots are reported as duplicate, even though they aren't. Maybe a refinement of this method can allow that?
The two simple roots are reported as [(0,1),(1,1)] in precision 1 and [(5,1),(6,1)] in any higher precision.
This is the caveat mentioned earlier: It doesn't do roots in Qp, it does roots in Zp. I am sure there is some easy rescaling of x that can bring all the roots into Zp (probably something like "divide x by powers of p until the leading coefficient is coprime to p"), and then transform them back after calculating the expansion.
This case would need to be checked explicitly, but that can be done with relative ease: As used in the definition of qi, if ri is a root of qi of multiplicity (up to O(p)), then q(p*x+ri) is divisible by p^mi. Otherwise, this is an "erroneous" root, and can be ignored. Such is the case for q(x) = x^2-5, which reports a double root at ri=0, but q(5x) = 25x^2-5 is not divisible by 25.
The constant polynomial x (and in fact all monomials x^n) return only the digit zero with multiplicity 1 (or n in general), so set however many digits you want.

I am aware that memory managing this will be a nightmare, but there are some optimizations that can be made right off the bat.
If (ri, mi) is the root used to define qi, then the modular projection of qi will be of degree mi, which will in general be less than the degree of q.
Also, if at some point qi == q (defined over Z, not Z/pZ), then the root ri is indeed of multiplicity mi, and its expansion becomes periodic from this point on.

tthsqe12 · 2021-08-09T20:58:00Z

Indeed, this is a intricate function to write. There are many possible approaches, and I do not expect anything that works to come in at fewer than 1k LOC.

fieker · 2023-10-09T16:15:47Z

fmpz_poly_roots/padic_roots.c

+	  padic_add (y0, y0, tmp, ctx);
+	}
+      /* Newton step: x -> x - poly / poly' */
+      padic_inv (y1, y1, ctx);


I'd suggest a double iteration to compute rt and inv(f'(rt)) at the same time. Otherwise, this inv repeats the previous inve all the time.

Also, in case on many roots, a lifting of the factorisation is much faster than lifting of individual roots.

The precision should be dynamically increased using a lifting chain...

I guess most of this applies to the qadic case as well

https://arxiv.org/pdf/1610.06837.pdf
algo 13, steps 1-3

P. Panayi, Computation of Leopoldt’s p-adic regulator, Dissertation, University of East Anglia, 1995.
contains a "certified" p-adic root finder also in the hard cases

monien added 3 commits April 11, 2021 23:01

Added functionality and documentation for calculating roots of univar…

ce8fab1

…iate polynomials over the integers (fmpz_poly) over finite fields, p-adic fields and q-adic fields.

fmpz_poly_roots: Added code to deal with multiplicity correctly for

98535b5

p-adic and q-adic roots.

fmpz_poly_roots: fixed typo in documentation.

056cddb

fieker reviewed Oct 9, 2023

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Added functionality and documentation for calculating roots of univariate polynomials #935

Added functionality and documentation for calculating roots of univariate polynomials #935

monien commented Apr 11, 2021

fredrik-johansson commented Apr 12, 2021

monien commented Apr 12, 2021 via email

monien commented Apr 12, 2021 via email

wbhart commented Apr 12, 2021

fredrik-johansson commented Apr 12, 2021

tthsqe12 commented Apr 12, 2021

tthsqe12 commented Apr 12, 2021

tthsqe12 commented Apr 12, 2021

tthsqe12 commented Apr 13, 2021 •

edited

Loading

tthsqe12 commented Apr 13, 2021 •

edited

Loading

monien commented Apr 13, 2021 via email

tthsqe12 commented Apr 13, 2021 •

edited

Loading

fredrik-johansson commented Apr 13, 2021

tthsqe12 commented Apr 13, 2021

fredrik-johansson commented Apr 13, 2021

monien commented Apr 24, 2021

tthsqe12 commented Apr 24, 2021

wbhart commented May 3, 2021 •

edited

Loading

MGessinger commented Aug 2, 2021 •

edited

Loading

tthsqe12 commented Aug 9, 2021

fieker Oct 9, 2023

fieker Oct 9, 2023

Added functionality and documentation for calculating roots of univariate polynomials #935

Are you sure you want to change the base?

Added functionality and documentation for calculating roots of univariate polynomials #935

Conversation

monien commented Apr 11, 2021

fredrik-johansson commented Apr 12, 2021

monien commented Apr 12, 2021 via email

monien commented Apr 12, 2021 via email

wbhart commented Apr 12, 2021

fredrik-johansson commented Apr 12, 2021

tthsqe12 commented Apr 12, 2021

tthsqe12 commented Apr 12, 2021

tthsqe12 commented Apr 12, 2021

tthsqe12 commented Apr 13, 2021 • edited Loading

tthsqe12 commented Apr 13, 2021 • edited Loading

monien commented Apr 13, 2021 via email

tthsqe12 commented Apr 13, 2021 • edited Loading

fredrik-johansson commented Apr 13, 2021

tthsqe12 commented Apr 13, 2021

fredrik-johansson commented Apr 13, 2021

monien commented Apr 24, 2021

tthsqe12 commented Apr 24, 2021

wbhart commented May 3, 2021 • edited Loading

MGessinger commented Aug 2, 2021 • edited Loading

tthsqe12 commented Aug 9, 2021

fieker Oct 9, 2023

Choose a reason for hiding this comment

fieker Oct 9, 2023

Choose a reason for hiding this comment

tthsqe12 commented Apr 13, 2021 •

edited

Loading

tthsqe12 commented Apr 13, 2021 •

edited

Loading

tthsqe12 commented Apr 13, 2021 •

edited

Loading

wbhart commented May 3, 2021 •

edited

Loading

MGessinger commented Aug 2, 2021 •

edited

Loading