Code coverage tests

This page documents the degree to which the PARI/GP source code is tested by our public test suite, distributed with the source distribution in directory src/test/. This is measured by the gcov utility; we then process gcov output using the lcov frond-end.

We test a few variants depending on Configure flags on the pari.math.u-bordeaux.fr machine (x86_64 architecture), and agregate them in the final report:

The target is 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

LCOV - code coverage report
Current view: top level - basemath - volcano.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.10.0 lcov report (development 21353-12523aa) Lines: 418 446 93.7 %
Date: 2017-11-24 06:20:58 Functions: 22 22 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2014  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : #include "pari.h"
      15             : #include "paripriv.h"
      16             : 
      17             : /* FIXME: Implement {ascend,descend}_volcano() in terms of the "new"
      18             :  * volcano traversal functions at the bottom of the file. */
      19             : 
      20             : /* Is j = 0 or 1728 (mod p)? */
      21             : INLINE int
      22      283883 : is_j_exceptional(ulong j, ulong p)
      23             : {
      24      283883 :   return j == 0 || j == 1728 % p;
      25             : }
      26             : 
      27             : 
      28             : INLINE long
      29       56499 : node_degree(GEN phi, long L, ulong j, ulong p, ulong pi)
      30             : {
      31       56499 :   pari_sp av = avma;
      32       56499 :   long n = Flx_nbroots(Flm_Fl_polmodular_evalx(phi, L, j, p, pi), p);
      33       56499 :   avma = av;
      34       56499 :   return n;
      35             : }
      36             : 
      37             : 
      38             : /*
      39             :  * Given an array path = [j0, j1] of length 2, return the polynomial
      40             :  *
      41             :  *   \Phi_L(X, j1) / (X - j0)
      42             :  *
      43             :  * where \Phi_L(X, Y) is the modular polynomial of level L.  An error
      44             :  * is raised if X - j0 does not divide \Phi_L(X, j1).
      45             :  */
      46             : INLINE GEN
      47      126723 : nhbr_polynomial(ulong path[], GEN phi, ulong p, ulong pi, long L)
      48             : {
      49      126723 :   pari_sp ltop = avma;
      50      126723 :   GEN modpol = Flm_Fl_polmodular_evalx(phi, L, path[0], p, pi);
      51             : 
      52             :   /* Note that, if the discriminant of End(path[0]) is less than L^2,
      53             :    * then it's possible for path[0] to appear among the roots of
      54             :    * nhbr_pol.  This case should have been obviated by earlier
      55             :    * choices. */
      56      126723 :   ulong rem = 0;
      57      126723 :   GEN nhbr_pol = Flx_div_by_X_x(modpol, path[-1], p, &rem);
      58      126723 :   if (rem)
      59           0 :     pari_err_BUG("nhbr_polynomial: invalid preceding j");
      60      126723 :   return gerepileupto(ltop, nhbr_pol);
      61             : }
      62             : 
      63             : 
      64             : /*
      65             :  * Assumes path is an array with space for at least max_len + 1
      66             :  * elements, whose first and second elements are the beginning of the
      67             :  * path.  I.e., the path starts
      68             :  *
      69             :  *   (path[0], path[1])
      70             :  *
      71             :  * If the result is less than max_len, then the last element of path
      72             :  * is definitely on the floor.  If the result equals max_len, then in
      73             :  * general it is unknown whether the last element of path is on the
      74             :  * floor or not.
      75             :  */
      76             : static long
      77      195339 : extend_path(
      78             :   ulong path[], GEN phi, ulong p, ulong pi, long L, long max_len)
      79             : {
      80      195339 :   pari_sp av = avma;
      81      195339 :   long d = 1;
      82      277319 :   for ( ; d < max_len; ++d) {
      83             :     ulong nhbr;
      84             :     GEN nhbr_pol;
      85             : 
      86       95611 :     nhbr_pol = nhbr_polynomial(path + d, phi, p, pi, L);
      87       95611 :     nhbr = Flx_oneroot(nhbr_pol, p);
      88       95611 :     avma = av;
      89             :     /* Flx_oneroot didn't find a root; so we must be on the floor. */
      90       95611 :     if (nhbr == p)
      91       13631 :       break;
      92       81980 :     path[d + 1] = nhbr;
      93             :   }
      94      195339 :   return d;
      95             : }
      96             : 
      97             : 
      98             : /*
      99             :  * This is Sutherland 2009 Algorithm Ascend (p12).
     100             :  */
     101             : ulong
     102      116820 : ascend_volcano(
     103             :   GEN phi, ulong j, ulong p, ulong pi, long level, long L,
     104             :   long depth, long steps)
     105             : {
     106      116820 :   pari_sp ltop = avma, av;
     107             :   /* path will never hold more than max_len elements, and max_len <
     108             :    * depth always. */
     109      116820 :   GEN path_g = cgetg(depth + 2, t_VECSMALL);
     110      116820 :   ulong *path = zv_to_ulongptr(path_g);
     111      116820 :   long max_len = depth - level;
     112      116820 :   int first_iter = 1;
     113             : 
     114      116820 :   if (steps <= 0 || max_len < 0)
     115           0 :     pari_err_BUG("ascend_volcano: bad params");
     116             : 
     117      116820 :   av = avma;
     118      381572 :   while (steps--) {
     119      147932 :     GEN nhbr_pol = first_iter
     120             :       ? Flm_Fl_polmodular_evalx(phi, L, j, p, pi)
     121      147932 :       : nhbr_polynomial(path + 1, phi, p, pi, L);
     122      147932 :     GEN nhbrs = Flx_roots(nhbr_pol, p);
     123      147932 :     long nhbrs_len = lg(nhbrs) - 1, i;
     124      147932 :     pari_sp btop = avma;
     125      147932 :     path[0] = j;
     126      147932 :     first_iter = 0;
     127             : 
     128      147932 :     j = nhbrs[nhbrs_len];
     129      177846 :     for (i = 1; i < nhbrs_len; ++i) {
     130       54823 :       ulong next_j = nhbrs[i], last_j;
     131             :       long len;
     132       54823 :       if (is_j_exceptional(next_j, p)) {
     133             :         /* According to Fouquet & Morain, Section 4.3, if j = 0 or
     134             :          * 1728, then it is necessarily on the surface.  So we just
     135             :          * return it. */
     136           0 :         if (steps) {
     137           0 :           pari_err_BUG("ascend_volcano: "
     138             :                        "Got to the top with more steps to go!");
     139             :         }
     140           0 :         j = next_j;
     141           0 :         break;
     142             :       }
     143       54823 :       path[1] = next_j;
     144       54823 :       len = extend_path(path, phi, p, pi, L, max_len);
     145       54823 :       last_j = path[len];
     146       54823 :       if (len == max_len
     147             :           /* Ended up on the surface */
     148       54823 :           && (is_j_exceptional(last_j, p)
     149       54823 :               || node_degree(phi, L, last_j, p, pi) > 1)) {
     150       24909 :         j = next_j;
     151       24909 :         break;
     152             :       }
     153       29914 :       avma = btop;
     154             :     }
     155      147932 :     path[1] = j; /* For nhbr_polynomial() at the top. */
     156             : 
     157      147932 :     ++max_len;
     158      147932 :     avma = av;
     159             :   }
     160      116820 :   avma = ltop;
     161      116820 :   return j;
     162             : }
     163             : 
     164             : 
     165             : static void
     166      174232 : random_distinct_neighbours_of(
     167             :   ulong *nhbr1, ulong *nhbr2,
     168             :   GEN phi, ulong j, ulong p, ulong pi, long L,
     169             :   long must_have_two_neighbours)
     170             : {
     171      174232 :   pari_sp ltop = avma;
     172      174232 :   GEN modpol = Flm_Fl_polmodular_evalx(phi, L, j, p, pi);
     173             :   ulong rem;
     174      174232 :   *nhbr1 = Flx_oneroot(modpol, p);
     175      174232 :   if (*nhbr1 == p) {
     176             :     /* Didn't even find one root! */
     177           0 :     char *err = stack_sprintf("random_distinct_neighbours_of: "
     178             :                               "No neighbours for j = %lu (mod %lu) "
     179             :                               "in %lu-volcano.", j, p, L);
     180           0 :     pari_err_BUG(err);
     181             :   }
     182      174232 :   modpol = Flx_div_by_X_x(modpol, *nhbr1, p, &rem);
     183      174232 :   *nhbr2 = Flx_oneroot(modpol, p);
     184      174232 :   if (must_have_two_neighbours && *nhbr2 == p) {
     185             :     /* Didn't find distinct root! */
     186           0 :     char *err = stack_sprintf("random_distinct_neighbours_of: "
     187             :                               "Only one neighbour for j = %lu (mod %lu) "
     188             :                               "in %lu-volcano.", j, p, L);
     189           0 :     pari_err_BUG(err);
     190             :   }
     191      174232 :   avma = ltop;
     192      174232 : }
     193             : 
     194             : 
     195             : /*
     196             :  * This is Sutherland 2009 Algorithm Descend (p12).
     197             :  */
     198             : ulong
     199        1586 : descend_volcano(
     200             :   GEN phi, ulong j, ulong p, ulong pi,
     201             :   long level, long L, long depth, long steps)
     202             : {
     203        1586 :   pari_sp ltop = avma;
     204             :   GEN path_g;
     205             :   ulong *path, res;
     206             :   long max_len;
     207             : 
     208        1586 :   if (steps <= 0 || level + steps > depth)
     209           0 :     pari_err_BUG("descend_volcano: bad params");
     210             : 
     211        1586 :   max_len = depth - level;
     212        1586 :   path_g = cgetg(max_len + 1 + 1, t_VECSMALL);
     213        1586 :   path = zv_to_ulongptr(path_g);
     214        1586 :   path[0] = j;
     215             :   /* level = 0 means we're on the volcano surface... */
     216        1586 :   if ( ! level) {
     217             :     /* Look for any path to the floor.  One of j's first three
     218             :      * neighbours must lead to the floor, since at most two neighbours
     219             :      * are on the surface. */
     220        1378 :     GEN nhbrs = Flx_roots(Flm_Fl_polmodular_evalx(phi, L, j, p, pi), p);
     221             :     long i;
     222        1488 :     for (i = 1; i <= 3; ++i) {
     223             :       long len;
     224        1488 :       path[1] = nhbrs[i];
     225        1488 :       len = extend_path(path, phi, p, pi, L, max_len);
     226             :       /* If nhbrs[i] took us to the floor: */
     227        1488 :       if (len < max_len || node_degree(phi, L, path[len], p, pi) == 1)
     228             :         break;
     229             :     }
     230             : 
     231        1378 :     if (i > 3) {
     232           0 :       pari_err_BUG("descend_volcano: "
     233             :                    "None of three neighbours lead to the floor");
     234             :     }
     235             :   } else {
     236             :     ulong nhbr1, nhbr2;
     237             :     long len;
     238         208 :     random_distinct_neighbours_of(&nhbr1, &nhbr2, phi, j, p, pi, L, 1);
     239         208 :     path[1] = nhbr1;
     240         208 :     len = extend_path(path, phi, p, pi, L, max_len);
     241             :     /* If last j isn't on the floor */
     242         208 :     if (len == max_len   /* Ended up on the surface. */
     243         208 :         && (is_j_exceptional(path[len], p)
     244         188 :             || node_degree(phi, L, path[len], p, pi) != 1)) {
     245             :       /* The other neighbour leads to the floor */
     246          76 :       path[1] = nhbr2;
     247          76 :       (void) extend_path(path, phi, p, pi, L, steps);
     248             :     }
     249             :   }
     250        1586 :   res = path[steps];
     251        1586 :   avma = ltop;
     252        1586 :   return res;
     253             : }
     254             : 
     255             : 
     256             : long
     257      174029 : j_level_in_volcano(
     258             :   GEN phi, ulong j, ulong p, ulong pi, long L, long depth)
     259             : {
     260      174029 :   pari_sp av = avma;
     261             :   GEN chunk;
     262             :   ulong *path1, *path2;
     263             :   long lvl;
     264             : 
     265      174029 :   if (depth == 0 || is_j_exceptional(j, p)) {
     266             :     /* According to Fouquet & Morain, Section 4.3, if j = 0 or 1728,
     267             :      * then it is necessarily on the surface.  Also, if the volcano
     268             :      * depth is zero then j necessarily has level 0. */
     269           5 :     return 0;
     270             :   }
     271             : 
     272      174024 :   chunk = new_chunk(2 * (depth + 1));
     273      174024 :   path1 = (ulong *) &chunk[0];
     274      174024 :   path2 = (ulong *) &chunk[depth + 1];
     275             : 
     276      174024 :   path1[0] = path2[0] = j;
     277             : 
     278      174024 :   random_distinct_neighbours_of(&path1[1], &path2[1],
     279             :                                 phi, j, p, pi, L, 0);
     280      174024 :   if (path2[1] == p) {
     281             :     /* Only found one neighbour, hence j is on the floor, hence
     282             :      * level == depth. */
     283      104652 :     lvl = depth;
     284             :   } else {
     285       69372 :     long path1_len = extend_path(path1, phi, p, pi, L, depth);
     286       69372 :     long path2_len = extend_path(path2, phi, p, pi, L, path1_len);
     287       69372 :     lvl = depth - path2_len;
     288             :   }
     289      174024 :   avma = av;
     290      174024 :   return lvl;
     291             : }
     292             : 
     293             : 
     294             : #define vecsmall_len(v) (lg(v) - 1)
     295             : 
     296             : INLINE GEN
     297    28144205 : Flx_remove_root(GEN f, ulong a, ulong p)
     298             : {
     299             :   ulong r;
     300    28144205 :   GEN g = Flx_div_by_X_x(f, a, p, &r);
     301    28153337 :   if (r) pari_err_BUG("Flx_remove_root");
     302    28163627 :   return g;
     303             : }
     304             : 
     305             : INLINE GEN
     306    21634489 : get_nbrs(GEN phi, long L, ulong J, const ulong *xJ, ulong p, ulong pi)
     307             : {
     308    21634489 :   pari_sp av = avma;
     309    21634489 :   GEN f = Flm_Fl_polmodular_evalx(phi, L, J, p, pi);
     310    21622741 :   if (xJ)
     311    21382698 :     f = Flx_remove_root(f, *xJ, p);
     312    21626561 :   return gerepileupto(av, Flx_roots(f, p));
     313             : }
     314             : 
     315             : /* Return a path of length n along the surface of an L-volcano of
     316             :  * height h starting from surface node j0.  Assumes (D|L) = 1 where D
     317             :  * is the discriminant of End(j0).
     318             :  *
     319             :  * Actually, if j0's endomorphism ring is a suborder, we return
     320             :  * the corresponding shorter path.
     321             :  *
     322             :  * W must hold space for n + h nodes.
     323             :  *
     324             :  * TODO: It could be worth having two versions of this function: one
     325             :  * that assumes J has the correct endomorphism ring (hence avoiding
     326             :  * several branches in the inner loop) and a second that doesn't
     327             :  * assume J has the correct endo ring and accordingly checks for
     328             :  * repetitions. */
     329             : static long
     330      175511 : surface_path(
     331             :   ulong W[],
     332             :   long n, GEN phi, long L, long h, ulong J, const ulong *nJ, ulong p, ulong pi)
     333             : {
     334      175511 :   pari_sp av = avma, bv;
     335             :   GEN T, v;
     336             :   long j, k, w, x;
     337             :   ulong W0;
     338             : 
     339      175511 :   W[0] = W0 = J;
     340      175511 :   if (n == 1)
     341           0 :     return 1;
     342             : 
     343      175511 :   T = cgetg(h + 2, t_VEC);
     344      175515 :   bv = avma;
     345             : 
     346      175515 :   v = get_nbrs(phi, L, J, nJ, p, pi);
     347             :   /* Insert known neighbour first */
     348      175513 :   if (nJ)
     349           0 :     v = gerepileupto(bv, vecsmall_prepend(v, *nJ));
     350      175514 :   gel(T, 1) = v;
     351      175514 :   k = vecsmall_len(v);
     352             : 
     353      175514 :   switch (k) {
     354             :   case 0:
     355             :     /* We must always have neighbours */
     356           0 :     pari_err_BUG("surface_path");
     357             :   case 1:
     358             :     /* If volcano is not flat, then we must have more than one
     359             :      * neighbour */
     360        5049 :     if (h) pari_err_BUG("surface_path");
     361        5049 :     W[1] = uel(v, 1);
     362        5049 :     avma = av;
     363             :     /* Check for bad endo ring */
     364        5049 :     if (W[1] == W[0])
     365          40 :       return 1;
     366        5009 :     return 2;
     367             :   case 2:
     368             :     /* If L=2 the only way we can have 2 neighbours is if we have a
     369             :      * double root.  This can only happen for |D| <= 16 (by Thm 2.2 of
     370             :      * Fouquet-Morain), and if it does we must have a 2-cycle.  This
     371             :      * case happens for D=-15. */
     372       20117 :     if (L == 2) {
     373             :       /* The double root is the neighbour on the surface, who will
     374             :        * have exactly one neighbour other than J; the other neighbour
     375             :        * of J has either 0 or 2 neighbours that are not J. */
     376           0 :       GEN u = get_nbrs(phi, L, uel(v, 1), &J, p, pi);
     377           0 :       long n = vecsmall_len(u) - !!vecsmall_isin(u, J);
     378           0 :       W[1] = n == 1 ? uel(v, 1) : uel(v, 2);
     379           0 :       avma = av;
     380           0 :       return 2;
     381             :     }
     382             : 
     383             :     /* Volcano is not flat, but we only found 2 neighbours for the
     384             :      * surface node J */
     385       20117 :     if (h) pari_err_BUG("surface_path");
     386             : 
     387       20262 :     W[1] = uel(v, 1);
     388             :     /* TODO: Can we use the other root -- uel(v,2) -- somehow? */
     389     4440925 :     for (w = 2; w < n; ++w) {
     390     4420878 :       v = get_nbrs(phi, L, W[w - 1], &W[w - 2], p, pi);
     391             :       /* A flat volcano must have exactly one non-previous
     392             :        * neighbour */
     393     4420279 :       if (vecsmall_len(v) != 1) pari_err_BUG("surface_path");
     394     4420733 :       W[w] = uel(v, 1);
     395     4420733 :       avma = av;
     396             :       /* Detect cycle in case J doesn't have the right endo ring. */
     397     4420733 :       if (W[w] == W0)
     398          70 :         return w;
     399             :     }
     400       20047 :     avma = av;
     401       20047 :     return n;
     402             :   }
     403             : 
     404             :   /* Can't have a flat volcano if k > 2. */
     405      150348 :   if (!h) pari_err_BUG("surface_path");
     406             : 
     407             :   /* At this point, each surface node must have L + 1 distinct
     408             :    * neighbours, 2 of which are on the surface. */
     409      150676 :   w = 1;
     410     5527094 :   for (x = 0; ; ++x) {
     411             :     /* Get next neighbour of last known surface node to attempt to
     412             :      * extend the path. */
     413     5527094 :     W[w] = umael(T, ((w - 1) % h) + 1, x + 1);
     414             : 
     415             :     /* Detect cycle in case the given J didn't have the right endo
     416             :      * ring. */
     417     5527094 :     if (W[w] == W0) {
     418          59 :       avma = av;
     419          59 :       return w;
     420             :     }
     421             : 
     422             :     /* If we have to test the last neighbour, we know it's on the
     423             :      * surface, and if we're done there's no need to extend. */
     424     5527035 :     if (x == k - 1 && w == n - 1) {
     425       97908 :       avma = av;
     426       97908 :       return n;
     427             :     }
     428             : 
     429             :     /* Walk forward until we hit the floor or finish. */
     430             :     /* NB: We don't keep the stack completely clean here; usage is
     431             :      * determined by the size of T, which will be in the order of Lh,
     432             :      * i.e. L roots for each level of the volcano of height h. */
     433     5429127 :     for (j = w; ; ) {
     434             :       long m;
     435             :       /* We must get 0 or L neighbours here. */
     436    16897985 :       v = get_nbrs(phi, L, W[j], &W[j - 1], p, pi);
     437    16892174 :       m = vecsmall_len(v);
     438    16892174 :       if ( ! m) {
     439             :         /* We hit the floor; save the neighbours of W[w-1] and throw
     440             :          * away the rest. */
     441     5374453 :         GEN nbrs = gel(T, ((w - 1) % h) + 1);
     442     5374453 :         gel(T, ((w - 1) % h) + 1) = gerepileupto(bv, nbrs);
     443     5376418 :         break;
     444             :       }
     445    11517721 :       if (m != L) pari_err_BUG("surface_path");
     446             : 
     447    11521240 :       gel(T, (j % h) + 1) = v;
     448             : 
     449    11521240 :       W[++j] = uel(v, 1);
     450             :       /* If we have our path by h nodes, we know W[w] is on the
     451             :        * surface */
     452    11521240 :       if (j == w + h) {
     453    10556904 :         ++w;
     454             :         /* Detect cycle in case the given J didn't have the right endo
     455             :          * ring. */
     456    10556904 :         if (W[w] == W0) {
     457       23707 :           avma = av;
     458       23707 :           return w;
     459             :         }
     460    10533197 :         x = 0;
     461    10533197 :         k = L;
     462             :       }
     463    11497533 :       if (w == n) {
     464       28675 :         avma = av;
     465       28675 :         return w;
     466             :       }
     467    11468858 :     }
     468     5376418 :   }
     469             :   return -1; /*LCOV_EXCL_LINE*/
     470             : }
     471             : 
     472             : long
     473      104970 : next_surface_nbr(
     474             :   ulong *nJ,
     475             :   GEN phi, long L, long h, ulong J, const ulong *pJ, ulong p, ulong pi)
     476             : {
     477      104970 :   pari_sp av = avma, bv;
     478             :   GEN S;
     479             :   ulong *P;
     480             :   long i, k;
     481             : 
     482      104970 :   S = get_nbrs(phi, L, J, pJ, p, pi);
     483      104971 :   k = vecsmall_len(S);
     484             :   /* If there is a double root and pJ is set, then k will be zero. */
     485      104971 :   if ( ! k) {
     486           0 :     avma = av;
     487           0 :     return 0;
     488             :   }
     489      104971 :   if (k == 1 || ( ! pJ && k == 2)) {
     490       87367 :     *nJ = uel(S, 1);
     491       87367 :     avma = av;
     492       87367 :     return 1;
     493             :   }
     494       17604 :   if ( ! h) pari_err_BUG("next_surface_nbr");
     495             : 
     496       17604 :   P = (ulong *) new_chunk(h + 1);
     497       17604 :   P[0] = J;
     498       17604 :   bv = avma;
     499       34902 :   for (i = 0; i < k; ++i) {
     500             :     long j;
     501       34902 :     P[1] = uel(S, i + 1);
     502       57238 :     for (j = 1; j <= h; ++j) {
     503       39634 :       GEN T = get_nbrs(phi, L, P[j], &P[j - 1], p, pi);
     504       39634 :       if ( ! vecsmall_len(T))
     505       17298 :         break;
     506       22336 :       P[j + 1] = uel(T, 1);
     507             :     }
     508       34902 :     if (j < h) pari_err_BUG("next_surface_nbr");
     509       34902 :     avma = bv;
     510       34902 :     if (j > h)
     511       17604 :       break;
     512             :   }
     513             :   /* TODO: We could save one get_nbrs call by iterating from i up to
     514             :    * k-1 and assume that the last (kth) nbr is the one we want. For
     515             :    * now we're careful and check that this last nbr really is on the
     516             :    * surface. */
     517       17604 :   if (i == k) pari_err_BUG("next_surf_nbr");
     518       17604 :   *nJ = uel(S, i + 1);
     519       17604 :   avma = av;
     520       17604 :   return 1;
     521             : }
     522             : 
     523             : /* If flag is set, requires a unique common neighor.  Otherwise it
     524             :  * will allow 2 common neighbors and return both. */
     525             : INLINE long
     526      170010 : common_nbr(
     527             :   ulong *nbr,
     528             :   ulong J1, GEN Phi1, long L1,
     529             :   ulong J2, GEN Phi2, long L2,
     530             :   ulong p, ulong pi, long flag)
     531             : {
     532      170010 :   pari_sp av = avma;
     533             :   GEN d, f, g, r;
     534             :   long rlen;
     535             : 
     536      170010 :   g = Flm_Fl_polmodular_evalx(Phi1, L1, J1, p, pi);
     537      170011 :   f = Flm_Fl_polmodular_evalx(Phi2, L2, J2, p, pi);
     538      170011 :   d = Flx_gcd(f, g, p);
     539      170008 :   if (degpol(d) == 1) {
     540      168895 :     *nbr = Flx_deg1_root(d, p);
     541      168898 :     avma = av;
     542      168898 :     return 1;
     543             :   }
     544        1113 :   if (flag || degpol(d) != 2) pari_err_BUG("common_neighbour");
     545        1113 :   r = Flx_roots(d, p);
     546        1113 :   rlen = vecsmall_len(r);
     547        1113 :   if ( ! rlen) pari_err_BUG("common_neighbour");
     548        1113 :   nbr[0] = uel(r, 1);
     549             :   /* rlen is 1 or 2 depending on whether the root is unique or not. */
     550        1113 :   nbr[1] = uel(r, rlen);
     551        1113 :   avma = av;
     552        1113 :   return 2;
     553             : }
     554             : 
     555             : /* Return gcd(Phi1(X,J1)/(X - J0), Phi2(X,J2)). Not stack clean. */
     556             : INLINE GEN
     557       42466 : common_nbr_pred_poly(
     558             :   ulong J1, GEN Phi1, long L1,
     559             :   ulong J2, GEN Phi2, long L2,
     560             :   ulong J0, ulong p, ulong pi)
     561             : {
     562             :   GEN f, g;
     563             : 
     564       42466 :   g = Flm_Fl_polmodular_evalx(Phi1, L1, J1, p, pi);
     565       42466 :   g = Flx_remove_root(g, J0, p);
     566       42466 :   f = Flm_Fl_polmodular_evalx(Phi2, L2, J2, p, pi);
     567       42466 :   return Flx_gcd(f, g, p);
     568             : }
     569             : 
     570             : /* Find common neighbour of J1 and J2, where J0 is an L1 predecessor
     571             :  * of J1.  Return 1 if successful, 0 if not. */
     572             : INLINE int
     573       41328 : common_nbr_pred(
     574             :   ulong *nbr,
     575             :   ulong J1, GEN Phi1, long L1,
     576             :   ulong J2, GEN Phi2, long L2,
     577             :   ulong J0, ulong p, ulong pi)
     578             : {
     579       41328 :   pari_sp av = avma;
     580             :   int res;
     581             :   GEN d;
     582             : 
     583       41328 :   d = common_nbr_pred_poly(J1, Phi1, L1, J2, Phi2, L2, J0, p, pi);
     584       41328 :   res = (degpol(d) == 1);
     585       41328 :   if (res)
     586       41328 :     *nbr = Flx_deg1_root(d, p);
     587       41328 :   avma = av;
     588       41328 :   return res;
     589             : }
     590             : 
     591             : INLINE long
     592        1138 : common_nbr_verify(
     593             :   ulong *nbr,
     594             :   ulong J1, GEN Phi1, long L1,
     595             :   ulong J2, GEN Phi2, long L2,
     596             :   ulong J0, ulong p, ulong pi)
     597             : {
     598        1138 :   pari_sp av = avma;
     599             :   GEN d;
     600             : 
     601        1138 :   d = common_nbr_pred_poly(J1, Phi1, L1, J2, Phi2, L2, J0, p, pi);
     602        1138 :   if ( ! degpol(d)) {
     603         765 :     avma = av;
     604         765 :     return 0;
     605             :   }
     606         373 :   if (degpol(d) > 1) pari_err_BUG("common_neighbour_verify");
     607         373 :   *nbr = Flx_deg1_root(d, p);
     608         373 :   avma = av;
     609         373 :   return 1;
     610             : }
     611             : 
     612             : INLINE ulong
     613         346 : Flm_Fl_polmodular_evalxy(
     614             :   GEN Phi, long L, ulong x, ulong y, ulong p, ulong pi)
     615             : {
     616         346 :   pari_sp av = avma;
     617         346 :   GEN f = Flm_Fl_polmodular_evalx(Phi, L, x, p, pi);
     618         346 :   ulong r = Flx_eval_pre(f, y, p, pi);
     619         346 :   avma = av;
     620         346 :   return r;
     621             : }
     622             : 
     623             : /* Find a common L1-neighbor of J1 and L2-neighbor of J2, given J0 an
     624             :  * L2-neighbor of J1 and an L1-neighbor of J2. Return 1 if successful,
     625             :  * 0 otherwise.  Will only fail if the initial J-invariant had the
     626             :  * wrong endo ring. */
     627             : INLINE int
     628       22607 : common_nbr_corner(
     629             :   ulong *nbr,
     630             :   ulong J1, GEN Phi1, long L1, long h1,
     631             :   ulong J2, GEN Phi2, long L2,
     632             :   ulong J0, ulong p, ulong pi)
     633             : {
     634             :   ulong nbrs[2];
     635             : 
     636       22607 :   if (common_nbr(nbrs, J1, Phi1, L1, J2, Phi2, L2, p, pi, 0) == 2
     637         225 :       && nbrs[0] != nbrs[1]) {
     638             :     ulong nJ1, nJ2;
     639             : 
     640         225 :     if ( ! next_surface_nbr(&nJ2, Phi1, L1, h1, J2, &J0, p, pi))
     641          13 :       return 0;  /* Error */
     642         225 :     if ( ! next_surface_nbr(&nJ1, Phi1, L1, h1, nbrs[0], &J1, p, pi))
     643           0 :       return 0;  /* Error */
     644             : 
     645         225 :     if (Flm_Fl_polmodular_evalxy(Phi2, L2, nJ1, nJ2, p, pi))
     646         104 :       nbrs[0] = nbrs[1];
     647             :     else {
     648         121 :       if ( ! next_surface_nbr(&nJ1, Phi1, L1, h1, nbrs[1], &J1, p, pi))
     649           0 :         return 0;  /* Error */
     650         121 :       if ( ! Flm_Fl_polmodular_evalxy(Phi2, L2, nJ1, nJ2, p, pi))
     651          13 :         return 0;  /* Error */
     652             :     }
     653             :   }
     654       22594 :   *nbr = nbrs[0];
     655       22594 :   return 1;
     656             : }
     657             : 
     658             : /*
     659             :  * Enumerates a surface L1-cycle using gcds with Phi_L2, where
     660             :  * alpha_L2=alpha_L1^e and |alpha_L1|=n (Here alpha_a denotes the
     661             :  * class represented by the pos def reduced primeform <a,b,c>).
     662             :  * Assumes n > e > 1 and that roots[0],...,roots[e-1] are already
     663             :  * present in W.
     664             : */
     665             : static long
     666       57116 : surface_gcd_cycle(
     667             :   ulong W[], ulong V[], long n,
     668             :   GEN Phi1, long L1, GEN Phi2, long L2, long e, ulong p, ulong pi)
     669             : {
     670       57116 :   pari_sp av = avma;
     671             :   long i1, i2, j1, j2;
     672             : 
     673       57116 :   i1 = j2 = 0;
     674       57116 :   i2 = j1 = e - 1;
     675             :   /* If W != V we assume V actually points to an L2-isogenous
     676             :    * parallel L1-path.  e should be 2 in this case */
     677       57116 :   if (W != V) {
     678       57116 :     i1 = j1 + 1;
     679       57116 :     i2 = n - 1;
     680             :   }
     681             :   do {
     682             :     ulong t0, t1, t2, h10, h11, h20, h21;
     683             :     long k;
     684             :     GEN f, g, h1, h2;
     685             : 
     686     3375059 :     f = Flm_Fl_polmodular_evalx(Phi2, L2, V[i1], p, pi);
     687     3366593 :     g = Flm_Fl_polmodular_evalx(Phi1, L1, W[j1], p, pi);
     688     3369268 :     g = Flx_remove_root(g, W[j1 - 1], p);
     689     3370553 :     h1 = Flx_gcd(f, g, p);
     690     3367933 :     if (degpol(h1) != 1) break; /* Error */
     691     3367696 :     h11 = Flx_coeff(h1, 1);
     692     3368565 :     h10 = Flx_coeff(h1, 0);
     693     3369626 :     avma = av;
     694             : 
     695     3369626 :     f = Flm_Fl_polmodular_evalx(Phi2, L2, V[i2], p, pi);
     696     3367015 :     g = Flm_Fl_polmodular_evalx(Phi1, L1, W[j2], p, pi);
     697     3369807 :     k = j2 + 1;
     698     3369807 :     if (k == n)
     699       51865 :       k = 0;
     700     3369807 :     g = Flx_remove_root(g, W[k], p);
     701     3370359 :     h2 = Flx_gcd(f, g, p);
     702     3368130 :     if (degpol(h2) != 1) break; /* Error */
     703     3367902 :     h21 = Flx_coeff(h2, 1);
     704     3369095 :     h20 = Flx_coeff(h2, 0);
     705     3369775 :     avma = av;
     706             : 
     707     3369775 :     ++i1;
     708     3369775 :     --i2;
     709     3369775 :     if (i2 < 0)
     710           0 :       i2 = n - 1;
     711     3369775 :     ++j1;
     712     3369775 :     --j2;
     713     3369775 :     if (j2 < 0)
     714       57115 :       j2 = n - 1;
     715             : 
     716     3369775 :     t0 = Fl_mul_pre(h11, h21, p, pi);
     717     3372958 :     t1 = Fl_inv(t0, p);
     718     3374198 :     t0 = Fl_mul_pre(t1, h21, p, pi);
     719     3375315 :     t2 = Fl_mul_pre(t0, h10, p, pi);
     720     3374820 :     W[j1] = Fl_neg(t2, p);
     721     3374187 :     t0 = Fl_mul_pre(t1, h11, p, pi);
     722     3375593 :     t2 = Fl_mul_pre(t0, h20, p, pi);
     723     3375583 :     W[j2] = Fl_neg(t2, p);
     724     3375059 :   } while (j2 > j1 + 1);
     725       57116 :   avma = av;
     726             :   /* Under normal circumstances, the loop exits when j2 == j1 + 1, in
     727             :    * which case we return n.  If we break early because of an error,
     728             :    * then j2 > j1 + 1, so (j2 - (j1 + 1)) > 0 (which is the number of
     729             :    * elements we haven't calculated yet), and we return n minus that
     730             :    * quantity. */
     731       57116 :   return n - j2 + (j1 + 1);
     732             : }
     733             : 
     734             : static long
     735        1036 : surface_gcd_path(
     736             :   ulong W[], ulong V[], long n,
     737             :   GEN Phi1, long L1, GEN Phi2, long L2, long e, ulong p, ulong pi)
     738             : {
     739        1036 :   pari_sp av = avma;
     740             :   long i, j;
     741             : 
     742        1036 :   i = 0;
     743        1036 :   j = e;
     744             :   /* If W != V then assume V actually points to a L2-isogenous
     745             :    * parallel L1-path.  e should be 2 in this case */
     746        1036 :   if (W != V)
     747        1036 :     i = j;
     748        4788 :   while (j < n) {
     749             :     GEN f, g, d;
     750             : 
     751        2716 :     f = Flm_Fl_polmodular_evalx(Phi2, L2, V[i], p, pi);
     752        2716 :     g = Flm_Fl_polmodular_evalx(Phi1, L1, W[j - 1], p, pi);
     753        2716 :     g = Flx_remove_root(g, W[j - 2], p);
     754        2716 :     d = Flx_gcd(f, g, p);
     755        2716 :     if (degpol(d) != 1) break; /* Error */
     756        2716 :     W[j] = Flx_deg1_root(d, p);
     757        2716 :     i++;
     758        2716 :     j++;
     759        2716 :     avma = av;
     760             :   }
     761        1036 :   avma = av;
     762        1036 :   return j;
     763             : }
     764             : 
     765             : /*
     766             :  * Given a path V of length n on an L1-volcano, and W[0] L2-isogenous
     767             :  * to V[0], extends the path W to length n on an L1-volcano, with W[i]
     768             :  * L2-isogenous to V[i].
     769             :  *
     770             :  * Uses gcds unless L2 is too large to make it helpful.  Always uses
     771             :  * GCD to get W[1] to ensure consistent orientation.
     772             :  *
     773             :  * Returns the new length of W.  This will of course almost always be
     774             :  * n, but could be lower if V was started with a J-invariant with bad
     775             :  * endomorphism ring.
     776             :  */
     777             : INLINE long
     778      147403 : surface_parallel_path(
     779             :   ulong W[], ulong V[], long n,
     780             :   GEN Phi1, long L1, GEN Phi2, long L2, ulong p, ulong pi, long cycle)
     781             : {
     782             :   ulong W2, nbrs[2];
     783             : 
     784      147403 :   if (common_nbr(nbrs, W[0], Phi1, L1, V[1], Phi2, L2, p, pi, 0) == 2
     785         888 :       && nbrs[0] != nbrs[1]) {
     786         790 :     if (n > 2) {
     787        1580 :       if ( ! common_nbr_verify(&W2, nbrs[0], Phi1, L1,
     788         790 :               V[2], Phi2, L2, W[0], p, pi)) {
     789             :         /* nbrs[1] must be the correct choice */
     790         442 :         nbrs[0] = nbrs[1];
     791         696 :       } else if (common_nbr_verify(&W2, nbrs[1], Phi1, L1,
     792         348 :               V[2], Phi2, L2, W[0], p, pi)) {
     793             :         /* Error: Both paths extend successfully */
     794          25 :         return 1;
     795             :       }
     796             :       /* Error: Two distinct choices with n = 2; ambiguous. */
     797           0 :     } else return 1;
     798             :   }
     799      147379 :   W[1] = nbrs[0];
     800      147379 :   if (n > 2) {
     801       58152 :     return cycle
     802             :       ? surface_gcd_cycle(W, V, n, Phi1, L1, Phi2, L2, 2, p, pi)
     803       58152 :       : surface_gcd_path(W, V, n, Phi1, L1, Phi2, L2, 2, p, pi);
     804             :   }
     805       89227 :   return n;
     806             : }
     807             : 
     808             : GEN
     809      175990 : enum_roots(ulong J0, norm_eqn_t ne, GEN fdb, classgp_pcp_t G)
     810             : {
     811             :   /* MAX_HEIGHT just needs to be bigger than the biggest value of
     812             :    * val_p(n) where p and n are ulongs. */
     813             :   enum { MAX_HEIGHT = BITS_IN_LONG };
     814      175990 :   pari_sp av, ltop = avma;
     815      175990 :   long s = !!G->L0;
     816      175990 :   long *n = G->n + s, *L = G->L + s, *o = G->o + s, k = G->k - s;
     817      175990 :   long i, t, vlen, *e, *h, *off, *poff, *M, N = G->enum_cnt;
     818      175990 :   ulong p = ne->p, pi = ne->pi, *roots;
     819             :   GEN vshape, vp, ve, roots_;
     820             :   GEN Phi;
     821             : 
     822      175990 :   if ( ! k) {
     823         480 :     roots_ = cgetg(2, t_VECSMALL);
     824         480 :     uel(roots_, 1) = J0;
     825         480 :     return roots_;
     826             :   }
     827             : 
     828      175510 :   roots_ = cgetg(N + MAX_HEIGHT, t_VECSMALL);
     829      175510 :   roots = zv_to_ulongptr(roots_);
     830      175510 :   av = avma;
     831             : 
     832             :   /* TODO: Shouldn't be factoring this every time. Store in *ne? */
     833      175510 :   vshape = factoru(ne->v);
     834      175510 :   vp = gel(vshape, 1);
     835      175510 :   ve = gel(vshape, 2);
     836             : 
     837      175510 :   vlen = vecsmall_len(vp);
     838      175510 :   Phi = new_chunk(k);
     839      175510 :   e = new_chunk(k);
     840      175511 :   off = new_chunk(k);
     841      175513 :   poff = new_chunk(k);
     842             :   /* TODO: Surely we can work these out ahead of time? */
     843             :   /* h[i] is the valuation of p[i] in v */
     844      175510 :   h = new_chunk(k);
     845      399020 :   for (i = 0; i < k; ++i) {
     846      223506 :     h[i] = 0;
     847      328256 :     for (t = 1; t <= vlen; ++t) {
     848      260679 :       if (vp[t] == L[i]) {
     849      155929 :         h[i] = uel(ve, t);
     850      155929 :         break;
     851             :       }
     852             :     }
     853      223506 :     e[i] = 0;
     854      223506 :     off[i] = 0;
     855      223506 :     gel(Phi, i) = polmodular_db_getp(fdb, L[i], p);
     856             :   }
     857             : 
     858      175514 :   M = new_chunk(k);
     859      223514 :   for (M[0] = 1, i = 1; i < k; ++i)
     860       48000 :     M[i] = M[i - 1] * n[i - 1];
     861             : 
     862      175514 :   t = surface_path(roots, n[0], gel(Phi, 0), L[0], h[0], J0, NULL, p, pi);
     863      175514 :   if (t < n[0]) {
     864             :     /* Error: J0 has bad endo ring */
     865         177 :     avma = ltop;
     866         177 :     return NULL;
     867             :   }
     868      175337 :   if (k == 1) {
     869      134011 :     avma = av;
     870      134011 :     setlg(roots_, t + 1);
     871      134011 :     return roots_;
     872             :   }
     873             : 
     874       41326 :   i = 1;
     875      230031 :   while (i < k) {
     876             :     long j, t0;
     877      147430 :     for (j = i + 1; j < k && ! e[j]; ++j);
     878      147430 :     if (j < k) {
     879       63935 :       if (e[i]) {
     880      330624 :         if (! common_nbr_pred(
     881      165312 :               &roots[t], roots[off[i]], gel(Phi, i), L[i],
     882      165312 :               roots[t - M[j]], gel(Phi, j), L[j], roots[poff[i]], p, pi)) {
     883           0 :           break; /* Error: J0 has bad endo ring */
     884             :         }
     885      203463 :       } else if ( ! common_nbr_corner(
     886      113035 :             &roots[t], roots[off[i]], gel(Phi, i), L[i], h[i],
     887       90428 :             roots[t - M[j]], gel(Phi, j), L[j], roots[poff[j]], p, pi)) {
     888          13 :         break; /* Error: J0 has bad endo ring */
     889             :       }
     890      536813 :     } else if ( ! next_surface_nbr(
     891      333980 :           &roots[t], gel(Phi, i), L[i], h[i],
     892      202833 :           roots[off[i]], e[i] ? &roots[poff[i]] : NULL, p, pi)) {
     893           0 :       break; /* Error: J0 has bad endo ring */
     894             :     }
     895             : 
     896      147418 :     if (roots[t] == roots[0]) break; /* Error: J0 has bad endo ring */
     897             : 
     898      147403 :     poff[i] = off[i];
     899      147403 :     off[i] = t;
     900      147403 :     ++e[i];
     901      170010 :     for (j = i - 1; j; --j) {
     902       22607 :       e[j] = 0;
     903       22607 :       off[j] = off[j + 1];
     904             :     }
     905             : 
     906      442209 :     t0 = surface_parallel_path(&roots[t], &roots[poff[i]], n[0],
     907      442209 :         gel(Phi, 0), L[0], gel(Phi, i), L[i], p, pi, n[0] == o[0]);
     908      147404 :     if (t0 < n[0]) break; /* Error: J0 has bad endo ring */
     909             : 
     910             :     /* TODO: Do I need to check if any of the new roots is a repeat in
     911             :      * the case where J0 has bad endo ring? */
     912      147379 :     t += n[0];
     913      147379 :     for (i = 1; i < k && e[i] == n[i] - 1; ++i);
     914             :   }
     915             :   /* Check if J0 had wrong endo ring */
     916       41328 :   if (t != N) { avma = ltop; return NULL; }
     917             : 
     918       41275 :   avma = av;
     919       41275 :   setlg(roots_, t + 1);
     920       41275 :   return roots_;
     921             : }

Generated by: LCOV version 1.11