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 to exceed 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.16.2 lcov report (development 29115-f22e516b23) Lines: 340 345 98.6 %
Date: 2024-03-18 08:03:28 Functions: 28 28 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; either version 2 of the License, or (at your option) any later
       8             : version. It is distributed in the hope that it will be useful, but WITHOUT
       9             : ANY WARRANTY WHATSOEVER.
      10             : 
      11             : Check the License for details. You should have received a copy of it, along
      12             : with the package; see the file 'COPYING'. If not, write to the Free Software
      13             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14             : 
      15             : #include "pari.h"
      16             : #include "paripriv.h"
      17             : 
      18             : static GEN
      19      213934 : pcp_get_L(GEN G) { return gmael(G,1,1)+1; }
      20             : static GEN
      21      213934 : pcp_get_n(GEN G) { return gmael(G,1,2)+1; }
      22             : static GEN
      23      213935 : pcp_get_o(GEN G) { return gmael(G,1,3)+1; }
      24             : static long
      25      213934 : pcp_get_L0(GEN G) { return mael(G,2,1); }
      26             : static long
      27      213935 : pcp_get_k(GEN G) { return mael(G,2,2); }
      28             : static long
      29      213935 : pcp_get_enum_cnt(GEN G) { return mael(G,2,3); }
      30             : 
      31             : /* FIXME: Implement {ascend,descend}_volcano() in terms of the "new"
      32             :  * volcano traversal functions at the bottom of the file. */
      33             : 
      34             : /* Is j = 0 or 1728 (mod p)? */
      35             : INLINE int
      36      356101 : is_j_exceptional(ulong j, ulong p) { return j == 0 || j == 1728 % p; }
      37             : 
      38             : INLINE long
      39       79929 : node_degree(GEN phi, long L, ulong j, ulong p, ulong pi)
      40             : {
      41       79929 :   pari_sp av = avma;
      42       79929 :   long n = Flx_nbroots(Flm_Fl_polmodular_evalx(phi, L, j, p, pi), p);
      43       79921 :   return gc_long(av, n);
      44             : }
      45             : 
      46             : /* Given an array path = [j0, j1] of length 2, return the polynomial
      47             :  *
      48             :  *   \Phi_L(X, j1) / (X - j0)
      49             :  *
      50             :  * where \Phi_L(X, Y) is the modular polynomial of level L.  An error
      51             :  * is raised if X - j0 does not divide \Phi_L(X, j1) */
      52             : INLINE GEN
      53      140142 : nhbr_polynomial(ulong path[], GEN phi, ulong p, ulong pi, long L)
      54             : {
      55      140142 :   GEN modpol = Flm_Fl_polmodular_evalx(phi, L, path[0], p, pi);
      56             :   ulong rem;
      57      140146 :   GEN nhbr_pol = Flx_div_by_X_x(modpol, path[-1], p, &rem);
      58             :   /* If disc End(path[0]) <= L^2, it's possible for path[0] to appear among the
      59             :    * roots of nhbr_pol. This should have been obviated by earlier choices */
      60      140132 :   if (rem) pari_err_BUG("nhbr_polynomial: invalid preceding j");
      61      140132 :   return nhbr_pol;
      62             : }
      63             : 
      64             : /* Path is an array with space for at least max_len+1 * elements, whose first
      65             :  * and second elements are the beginning of the path.  I.e., the path starts
      66             :  *   (path[0], path[1])
      67             :  * If the result is less than max_len, then the last element of path is on the
      68             :  * floor.  If the result equals max_len, then it is unknown whether the last
      69             :  * element of path is on the floor or not */
      70             : static long
      71      274343 : extend_path(ulong path[], GEN phi, ulong p, ulong pi, long L, long max_len)
      72             : {
      73      274343 :   pari_sp av = avma;
      74      274343 :   long d = 1;
      75      353288 :   for ( ; d < max_len; d++)
      76             :   {
      77      101764 :     GEN nhbr_pol = nhbr_polynomial(path + d, phi, p, pi, L);
      78      101757 :     ulong nhbr = Flx_oneroot_pre(nhbr_pol, p, pi);
      79      101762 :     set_avma(av);
      80      101764 :     if (nhbr == p) break; /* no root: we are on the floor. */
      81       78945 :     path[d+1] = nhbr;
      82             :   }
      83      274343 :   return d;
      84             : }
      85             : 
      86             : /* This is Sutherland 2009 Algorithm Ascend (p12) */
      87             : ulong
      88      125151 : ascend_volcano(GEN phi, ulong j, ulong p, ulong pi, long level, long L,
      89             :   long depth, long steps)
      90             : {
      91      125151 :   pari_sp ltop = avma, av;
      92             :   /* path will never hold more than max_len < depth elements */
      93      125151 :   GEN path_g = cgetg(depth + 2, t_VECSMALL);
      94      125151 :   ulong *path = zv_to_ulongptr(path_g);
      95      125151 :   long max_len = depth - level;
      96      125151 :   int first_iter = 1;
      97             : 
      98      125151 :   if (steps <= 0 || max_len < 0) pari_err_BUG("ascend_volcano: bad params");
      99      125151 :   av = avma;
     100      288680 :   while (steps--)
     101             :   {
     102      125151 :     GEN nhbr_pol = first_iter? Flm_Fl_polmodular_evalx(phi, L, j, p, pi)
     103      163533 :                              : nhbr_polynomial(path+1, phi, p, pi, L);
     104      163536 :     GEN nhbrs = Flx_roots_pre(nhbr_pol, p, pi);
     105      163530 :     long nhbrs_len = lg(nhbrs)-1, i;
     106      163530 :     pari_sp btop = avma;
     107      163530 :     path[0] = j;
     108      163530 :     first_iter = 0;
     109             : 
     110      163530 :     j = nhbrs[nhbrs_len];
     111      206046 :     for (i = 1; i < nhbrs_len; i++)
     112             :     {
     113       76775 :       ulong next_j = nhbrs[i], last_j;
     114             :       long len;
     115       76775 :       if (is_j_exceptional(next_j, p))
     116             :       {
     117             :         /* Fouquet & Morain, Section 4.3, if j = 0 or 1728, then it is on the
     118             :          * surface.  So we just return it. */
     119          40 :         if (steps)
     120           0 :           pari_err_BUG("ascend_volcano: Got to the top with more steps to go!");
     121          40 :         j = next_j; break;
     122             :       }
     123       76735 :       path[1] = next_j;
     124       76735 :       len = extend_path(path, phi, p, pi, L, max_len);
     125       76735 :       last_j = path[len];
     126       76735 :       if (len == max_len
     127             :           /* Ended up on the surface */
     128       76735 :           && (is_j_exceptional(last_j, p)
     129       76734 :               || node_degree(phi, L, last_j, p, pi) > 1)) { j = next_j; break; }
     130       42513 :       set_avma(btop);
     131             :     }
     132      163526 :     path[1] = j; /* For nhbr_polynomial() at the top. */
     133             : 
     134      163526 :     max_len++; set_avma(av);
     135             :   }
     136      125147 :   return gc_long(ltop, j);
     137             : }
     138             : 
     139             : static void
     140      202593 : random_distinct_neighbours_of(ulong *nhbr1, ulong *nhbr2,
     141             :   GEN phi, ulong j, ulong p, ulong pi, long L, long must_have_two_neighbours)
     142             : {
     143      202593 :   pari_sp av = avma;
     144      202593 :   GEN modpol = Flm_Fl_polmodular_evalx(phi, L, j, p, pi);
     145             :   ulong rem;
     146      202596 :   *nhbr1 = Flx_oneroot_pre(modpol, p, pi);
     147      202598 :   if (*nhbr1 == p) pari_err_BUG("random_distinct_neighbours_of [no neighbour]");
     148      202598 :   modpol = Flx_div_by_X_x(modpol, *nhbr1, p, &rem);
     149      202585 :   *nhbr2 = Flx_oneroot_pre(modpol, p, pi);
     150      202595 :   if (must_have_two_neighbours && *nhbr2 == p)
     151           0 :     pari_err_BUG("random_distinct_neighbours_of [single neighbour]");
     152      202595 :   set_avma(av);
     153      202594 : }
     154             : 
     155             : /*
     156             :  * This is Sutherland 2009 Algorithm Descend (p12).
     157             :  */
     158             : ulong
     159        2931 : descend_volcano(GEN phi, ulong j, ulong p, ulong pi,
     160             :   long level, long L, long depth, long steps)
     161             : {
     162        2931 :   pari_sp ltop = avma;
     163             :   GEN path_g;
     164             :   ulong *path;
     165             :   long max_len;
     166             : 
     167        2931 :   if (steps <= 0 || level + steps > depth) pari_err_BUG("descend_volcano");
     168        2931 :   max_len = depth - level;
     169        2931 :   path_g = cgetg(max_len + 1 + 1, t_VECSMALL);
     170        2931 :   path = zv_to_ulongptr(path_g);
     171        2931 :   path[0] = j;
     172             :   /* level = 0 means we're on the volcano surface... */
     173        2931 :   if (!level)
     174             :   {
     175             :     /* Look for any path to the floor. One of j's first three neighbours leads
     176             :      * to the floor, since at most two neighbours are on the surface. */
     177        2647 :     GEN nhbrs = Flx_roots_pre(Flm_Fl_polmodular_evalx(phi, L, j, p, pi), p, pi);
     178             :     long i;
     179        2950 :     for (i = 1; i <= 3; i++)
     180             :     {
     181             :       long len;
     182        2950 :       path[1] = nhbrs[i];
     183        2950 :       len = extend_path(path, phi, p, pi, L, max_len);
     184             :       /* If nhbrs[i] took us to the floor: */
     185        2950 :       if (len < max_len || node_degree(phi, L, path[len], p, pi) == 1) break;
     186             :     }
     187        2647 :     if (i > 3) pari_err_BUG("descend_volcano [2]");
     188             :   }
     189             :   else
     190             :   {
     191             :     ulong nhbr1, nhbr2;
     192             :     long len;
     193         284 :     random_distinct_neighbours_of(&nhbr1, &nhbr2, phi, j, p, pi, L, 1);
     194         284 :     path[1] = nhbr1;
     195         284 :     len = extend_path(path, phi, p, pi, L, max_len);
     196             :     /* If last j isn't on the floor */
     197         284 :     if (len == max_len   /* Ended up on the surface. */
     198         284 :         && (is_j_exceptional(path[len], p)
     199         244 :             || node_degree(phi, L, path[len], p, pi) != 1)) {
     200             :       /* The other neighbour leads to the floor */
     201         120 :       path[1] = nhbr2;
     202         120 :       (void) extend_path(path, phi, p, pi, L, steps);
     203             :     }
     204             :   }
     205        2931 :   return gc_ulong(ltop, path[steps]);
     206             : }
     207             : 
     208             : long
     209      202315 : j_level_in_volcano(
     210             :   GEN phi, ulong j, ulong p, ulong pi, long L, long depth)
     211             : {
     212      202315 :   pari_sp av = avma;
     213             :   GEN chunk;
     214             :   ulong *path1, *path2;
     215             :   long lvl;
     216             : 
     217             :   /* Fouquet & Morain, Section 4.3, if j = 0 or 1728 then it is on the
     218             :    * surface.  Also, if the volcano depth is zero then j has level 0 */
     219      202315 :   if (depth == 0 || is_j_exceptional(j, p)) return 0;
     220             : 
     221      202309 :   chunk = new_chunk(2 * (depth + 1));
     222      202309 :   path1 = (ulong *) &chunk[0];
     223      202309 :   path2 = (ulong *) &chunk[depth + 1];
     224      202309 :   path1[0] = path2[0] = j;
     225      202309 :   random_distinct_neighbours_of(&path1[1], &path2[1], phi, j, p, pi, L, 0);
     226      202315 :   if (path2[1] == p)
     227      105186 :     lvl = depth; /* Only one neighbour => j is on the floor => level = depth */
     228             :   else
     229             :   {
     230       97129 :     long path1_len = extend_path(path1, phi, p, pi, L, depth);
     231       97134 :     long path2_len = extend_path(path2, phi, p, pi, L, path1_len);
     232       97134 :     lvl = depth - path2_len;
     233             :   }
     234      202320 :   return gc_long(av, lvl);
     235             : }
     236             : 
     237             : INLINE GEN
     238    32116701 : Flx_remove_root(GEN f, ulong a, ulong p)
     239             : {
     240             :   ulong r;
     241    32116701 :   GEN g = Flx_div_by_X_x(f, a, p, &r);
     242    31952001 :   if (r) pari_err_BUG("Flx_remove_root");
     243    31952349 :   return g;
     244             : }
     245             : 
     246             : INLINE GEN
     247    24212665 : get_nbrs(GEN phi, long L, ulong J, const ulong *xJ, ulong p, ulong pi)
     248             : {
     249    24212665 :   pari_sp av = avma;
     250    24212665 :   GEN f = Flm_Fl_polmodular_evalx(phi, L, J, p, pi);
     251    24220596 :   if (xJ) f = Flx_remove_root(f, *xJ, p);
     252    24123704 :   return gerepileupto(av, Flx_roots_pre(f, p, pi));
     253             : }
     254             : 
     255             : /* Return a path of length n along the surface of an L-volcano of height h
     256             :  * starting from surface node j0.  Assumes (D|L) = 1 where D = disc End(j0).
     257             :  *
     258             :  * Actually, if j0's endomorphism ring is a suborder, we return the
     259             :  * corresponding shorter path. W must hold space for n + h nodes.
     260             :  *
     261             :  * TODO: have two versions of this function: one that assumes J has the correct
     262             :  * endomorphism ring (hence avoiding several branches in the inner loop) and a
     263             :  * second that does not and accordingly checks for repetitions */
     264             : static long
     265      211582 : surface_path(
     266             :   ulong W[],
     267             :   long n, GEN phi, long L, long h, ulong J, const ulong *nJ, ulong p, ulong pi)
     268             : {
     269      211582 :   pari_sp av = avma, bv;
     270             :   GEN T, v;
     271             :   long j, k, w, x;
     272             :   ulong W0;
     273             : 
     274      211582 :   W[0] = W0 = J;
     275      211582 :   if (n == 1) return 1;
     276             : 
     277      211582 :   T = cgetg(h+2, t_VEC); bv = avma;
     278      211580 :   v = get_nbrs(phi, L, J, nJ, p, pi);
     279             :   /* Insert known neighbour first */
     280      211588 :   if (nJ) v = gerepileupto(bv, vecsmall_prepend(v, *nJ));
     281      211588 :   gel(T,1) = v; k = lg(v)-1;
     282             : 
     283      211588 :   switch (k) {
     284           0 :   case 0: pari_err_BUG("surface_path"); /* We must always have neighbours */
     285        8772 :   case 1:
     286             :     /* If volcano is not flat, then we must have more than one neighbour */
     287        8772 :     if (h) pari_err_BUG("surface_path");
     288        8772 :     W[1] = uel(v, 1);
     289        8772 :     set_avma(av);
     290             :     /* Check for bad endo ring */
     291        8772 :     if (W[1] == W[0]) return 1;
     292        8591 :     return 2;
     293       25981 :   case 2:
     294             :     /* If L=2 the only way we can have 2 neighbours is if we have a double root
     295             :      * which can only happen for |D| <= 16 (Thm 2.2 of Fouquet-Morain)
     296             :      * and if it does we must have a 2-cycle. Happens for D=-15. */
     297       25981 :     if (L == 2)
     298             :     { /* The double root is the neighbour on the surface, with exactly one
     299             :        * neighbour other than J; the other neighbour of J has either 0 or 2
     300             :        * neighbours that are not J */
     301          84 :       GEN u = get_nbrs(phi, L, uel(v, 1), &J, p, pi);
     302          84 :       long n = lg(u)-1 - !!vecsmall_isin(u, J);
     303          84 :       W[1] = n == 1 ? uel(v,1) : uel(v,2);
     304          84 :       return gc_long(av, 2);
     305             :     }
     306             :     /* Volcano is not flat but found only 2 neighbours for the surface node J */
     307       25897 :     if (h) pari_err_BUG("surface_path");
     308             : 
     309       25897 :     W[1] = uel(v,1); /* TODO: Can we use the other root uel(v,2) somehow? */
     310     4451789 :     for (w = 2; w < n; w++)
     311             :     {
     312     4426450 :       v = get_nbrs(phi, L, W[w-1], &W[w-2], p, pi);
     313             :       /* A flat volcano must have exactly one non-previous neighbour */
     314     4426477 :       if (lg(v) != 2) pari_err_BUG("surface_path");
     315     4426477 :       W[w] = uel(v, 1);
     316             :       /* Detect cycle in case J doesn't have the right endo ring. */
     317     4426477 :       set_avma(av); if (W[w] == W0) return w;
     318             :     }
     319       25339 :     return gc_long(av, n);
     320             :   }
     321      176835 :   if (!h) pari_err_BUG("surface_path"); /* Can't have a flat volcano if k > 2 */
     322             : 
     323             :   /* At this point, each surface node has L+1 distinct neighbours, 2 of which
     324             :    * are on the surface */
     325      176835 :   w = 1;
     326     6373999 :   for (x = 0;; x++)
     327             :   {
     328             :     /* Get next neighbour of last known surface node to attempt to
     329             :      * extend the path. */
     330     6373999 :     W[w] = umael(T, ((w-1) % h) + 1, x + 1);
     331             : 
     332             :     /* Detect cycle in case the given J didn't have the right endo ring */
     333     6373999 :     if (W[w] == W0) return gc_long(av,w);
     334             : 
     335             :     /* If we have to test the last neighbour, we know it's on the
     336             :      * surface, and if we're done there's no need to extend. */
     337     6373967 :     if (x == k-1 && w == n-1) return gc_long(av,n);
     338             : 
     339             :     /* Walk forward until we hit the floor or finish. */
     340             :     /* NB: We don't keep the stack clean here; usage is in the order of Lh,
     341             :      * i.e. L roots for each level of the volcano of height h. */
     342     6257871 :     for (j = w;;)
     343    13135158 :     {
     344             :       long m;
     345             :       /* We must get 0 or L neighbours here. */
     346    19393029 :       v = get_nbrs(phi, L, W[j], &W[j-1], p, pi);
     347    19342198 :       m = lg(v)-1;
     348    19342198 :       if (!m) {
     349             :         /* We hit the floor: save the neighbours of W[w-1] and dump the rest */
     350     6197355 :         GEN nbrs = gel(T, ((w-1) % h) + 1);
     351     6197355 :         gel(T, ((w-1) % h) + 1) = gerepileupto(bv, nbrs);
     352     6197164 :         break;
     353             :       }
     354    13144843 :       if (m != L) pari_err_BUG("surface_path");
     355             : 
     356    13195868 :       gel(T, (j % h) + 1) = v;
     357             : 
     358    13195868 :       W[++j] = uel(v, 1);
     359             :       /* If we have our path by h nodes, we know W[w] is on the surface */
     360    13195868 :       if (j == w + h) {
     361    12209339 :         ++w;
     362             :         /* Detect cycle in case the given J didn't have the right endo ring */
     363    12209339 :         if (W[w] == W0) return gc_long(av,w);
     364    12182011 :         x = 0; k = L;
     365             :       }
     366    13168540 :       if (w == n) return gc_long(av,w);
     367             :     }
     368             :   }
     369             : }
     370             : 
     371             : long
     372      144919 : next_surface_nbr(
     373             :   ulong *nJ,
     374             :   GEN phi, long L, long h, ulong J, const ulong *pJ, ulong p, ulong pi)
     375             : {
     376      144919 :   pari_sp av = avma, bv;
     377             :   GEN S;
     378             :   ulong *P;
     379             :   long i, k;
     380             : 
     381      144919 :   S = get_nbrs(phi, L, J, pJ, p, pi); k = lg(S)-1;
     382             :   /* If there is a double root and pJ is set, then k will be zero. */
     383      144916 :   if (!k) return gc_long(av,0);
     384      144916 :   if (k == 1 || ( ! pJ && k == 2)) { *nJ = uel(S, 1); return gc_long(av,1); }
     385       24099 :   if (!h) pari_err_BUG("next_surface_nbr");
     386             : 
     387       24099 :   P = (ulong *) new_chunk(h + 1); bv = avma;
     388       24100 :   P[0] = J;
     389       51403 :   for (i = 0; i < k; i++)
     390             :   {
     391             :     long j;
     392       51403 :     P[1] = uel(S, i + 1);
     393       81029 :     for (j = 1; j <= h; j++)
     394             :     {
     395       56935 :       GEN T = get_nbrs(phi, L, P[j], &P[j - 1], p, pi);
     396       56932 :       if (lg(T) == 1) break;
     397       29626 :       P[j + 1] = uel(T, 1);
     398             :     }
     399       51400 :     if (j < h) pari_err_BUG("next_surface_nbr");
     400       51400 :     set_avma(bv); if (j > h) break;
     401             :   }
     402             :   /* TODO: We could save one get_nbrs call by iterating from i up to k-1 and
     403             :    * assume that the last (kth) nbr is the one we want. For now we're careful
     404             :    * and check that this last nbr really is on the surface */
     405       24098 :   if (i == k) pari_err_BUG("next_surf_nbr");
     406       24098 :   *nJ = uel(S, i+1); return gc_long(av,1);
     407             : }
     408             : 
     409             : /* Return the number of distinct neighbours (1 or 2) */
     410             : INLINE long
     411      236994 : common_nbr(ulong *nbr,
     412             :   ulong J1, GEN Phi1, long L1,
     413             :   ulong J2, GEN Phi2, long L2, ulong p, ulong pi)
     414             : {
     415      236994 :   pari_sp av = avma;
     416             :   GEN d, f, g, r;
     417             :   long rlen;
     418             : 
     419      236994 :   g = Flm_Fl_polmodular_evalx(Phi1, L1, J1, p, pi);
     420      237042 :   f = Flm_Fl_polmodular_evalx(Phi2, L2, J2, p, pi);
     421      237038 :   d = Flx_gcd(f, g, p);
     422      236841 :   if (degpol(d) == 1) { *nbr = Flx_deg1_root(d, p); return gc_long(av,1); }
     423        1237 :   if (degpol(d) != 2) pari_err_BUG("common_neighbour");
     424        1237 :   r = Flx_roots_pre(d, p, pi);
     425        1237 :   rlen = lg(r)-1;
     426        1237 :   if (!rlen) pari_err_BUG("common_neighbour");
     427             :   /* rlen is 1 or 2 depending on whether the root is unique or not. */
     428        1237 :   nbr[0] = uel(r, 1);
     429        1237 :   nbr[1] = uel(r, rlen); return gc_long(av,rlen);
     430             : }
     431             : 
     432             : /* Return gcd(Phi1(X,J1)/(X - J0), Phi2(X,J2)). Not stack clean. */
     433             : INLINE GEN
     434       44093 : common_nbr_pred_poly(
     435             :   ulong J1, GEN Phi1, long L1,
     436             :   ulong J2, GEN Phi2, long L2, ulong J0, ulong p, ulong pi)
     437             : {
     438             :   GEN f, g;
     439             : 
     440       44093 :   g = Flm_Fl_polmodular_evalx(Phi1, L1, J1, p, pi);
     441       44100 :   g = Flx_remove_root(g, J0, p);
     442       44098 :   f = Flm_Fl_polmodular_evalx(Phi2, L2, J2, p, pi);
     443       44103 :   return Flx_gcd(f, g, p);
     444             : }
     445             : 
     446             : /* Find common neighbour of J1 and J2, where J0 is an L1 predecessor of J1.
     447             :  * Return 1 if successful, 0 if not. */
     448             : INLINE int
     449       43018 : common_nbr_pred(ulong *nbr,
     450             :   ulong J1, GEN Phi1, long L1,
     451             :   ulong J2, GEN Phi2, long L2, ulong J0, ulong p, ulong pi)
     452             : {
     453       43018 :   pari_sp av = avma;
     454       43018 :   GEN d = common_nbr_pred_poly(J1, Phi1, L1, J2, Phi2, L2, J0, p, pi);
     455       43009 :   int res = (degpol(d) == 1);
     456       43007 :   if (res) *nbr = Flx_deg1_root(d, p);
     457       43023 :   return gc_bool(av,res);
     458             : }
     459             : 
     460             : INLINE long
     461        1074 : common_nbr_verify(ulong *nbr,
     462             :   ulong J1, GEN Phi1, long L1,
     463             :   ulong J2, GEN Phi2, long L2, ulong J0, ulong p, ulong pi)
     464             : {
     465        1074 :   pari_sp av = avma;
     466        1074 :   GEN d = common_nbr_pred_poly(J1, Phi1, L1, J2, Phi2, L2, J0, p, pi);
     467             : 
     468        1074 :   if (!degpol(d)) return gc_long(av,0);
     469         408 :   if (degpol(d) > 1) pari_err_BUG("common_neighbour_verify");
     470         408 :   *nbr = Flx_deg1_root(d, p);
     471         408 :   return gc_long(av,1);
     472             : }
     473             : 
     474             : INLINE ulong
     475         477 : Flm_Fl_polmodular_evalxy(GEN Phi, long L, ulong x, ulong y, ulong p, ulong pi)
     476             : {
     477         477 :   pari_sp av = avma;
     478         477 :   GEN f = Flm_Fl_polmodular_evalx(Phi, L, x, p, pi);
     479         477 :   return gc_ulong(av, Flx_eval_pre(f, y, p, pi));
     480             : }
     481             : 
     482             : /* Find a common L1-neighbor of J1 and L2-neighbor of J2, given J0 an
     483             :  * L2-neighbor of J1 and an L1-neighbor of J2. Return 1 if successful, 0
     484             :  * otherwise. Will only fail if initial J-invariant had the wrong endo ring */
     485             : INLINE int
     486       37284 : common_nbr_corner(ulong *nbr,
     487             :   ulong J1, GEN Phi1, long L1, long h1,
     488             :   ulong J2, GEN Phi2, long L2, ulong J0, ulong p, ulong pi)
     489             : {
     490             :   ulong nbrs[2];
     491       37284 :   if (common_nbr(nbrs, J1,Phi1,L1, J2,Phi2,L2, p, pi) == 2)
     492             :   {
     493             :     ulong nJ1, nJ2;
     494         640 :     if (!next_surface_nbr(&nJ2, Phi1, L1, h1, J2, &J0, p, pi) ||
     495         354 :         !next_surface_nbr(&nJ1, Phi1, L1, h1, nbrs[0], &J1, p, pi)) return 0;
     496             : 
     497         320 :     if (Flm_Fl_polmodular_evalxy(Phi2, L2, nJ1, nJ2, p, pi))
     498         163 :       nbrs[0] = nbrs[1];
     499         314 :     else if (!next_surface_nbr(&nJ1, Phi1, L1, h1, nbrs[1], &J1, p, pi) ||
     500         191 :              !Flm_Fl_polmodular_evalxy(Phi2, L2, nJ1, nJ2, p, pi)) return 0;
     501             :   }
     502       37251 :   *nbr = nbrs[0]; return 1;
     503             : }
     504             : 
     505             : /* Enumerate a surface L1-cycle using gcds with Phi_L2, where c_L2=c_L1^e and
     506             :  * |c_L1|=n, where c_a is the class of the pos def reduced primeform <a,b,c>.
     507             :  * Assumes n > e > 1 and roots[0],...,roots[e-1] are already present in W */
     508             : static long
     509       92778 : surface_gcd_cycle(
     510             :   ulong W[], ulong V[], long n,
     511             :   GEN Phi1, long L1, GEN Phi2, long L2, long e, ulong p, ulong pi)
     512             : {
     513       92778 :   pari_sp av = avma;
     514             :   long i1, i2, j1, j2;
     515             : 
     516       92778 :   i1 = j2 = 0;
     517       92778 :   i2 = j1 = e - 1;
     518             :   /* If W != V we assume V actually points to an L2-isogenous parallel L1-path.
     519             :    * e should be 2 in this case */
     520       92778 :   if (W != V) { i1 = j1+1; i2 = n-1; }
     521             :   do {
     522             :     ulong t0, t1, t2, h10, h11, h20, h21;
     523             :     long k;
     524             :     GEN f, g, h1, h2;
     525             : 
     526     4111568 :     f = Flm_Fl_polmodular_evalx(Phi2, L2, V[i1], p, pi);
     527     4106262 :     g = Flm_Fl_polmodular_evalx(Phi1, L1, W[j1], p, pi);
     528     4105583 :     g = Flx_remove_root(g, W[j1 - 1], p);
     529     4092994 :     h1 = Flx_gcd(f, g, p);
     530     4082996 :     if (degpol(h1) != 1) break; /* Error */
     531     4084071 :     h11 = Flx_coeff(h1, 1);
     532     4085371 :     h10 = Flx_coeff(h1, 0); set_avma(av);
     533             : 
     534     4086076 :     f = Flm_Fl_polmodular_evalx(Phi2, L2, V[i2], p, pi);
     535     4106310 :     g = Flm_Fl_polmodular_evalx(Phi1, L1, W[j2], p, pi);
     536     4106087 :     k = j2 + 1;
     537     4106087 :     if (k == n) k = 0;
     538     4106087 :     g = Flx_remove_root(g, W[k], p);
     539     4094486 :     h2 = Flx_gcd(f, g, p);
     540     4084535 :     if (degpol(h2) != 1) break; /* Error */
     541     4086219 :     h21 = Flx_coeff(h2, 1);
     542     4086939 :     h20 = Flx_coeff(h2, 0); set_avma(av);
     543             : 
     544     4088503 :     i1++; i2--; if (i2 < 0) i2 = n-1;
     545     4088503 :     j1++; j2--; if (j2 < 0) j2 = n-1;
     546             : 
     547     4088503 :     t0 = Fl_mul_pre(h11, h21, p, pi);
     548     4110515 :     t1 = Fl_inv(t0, p);
     549     4110212 :     t0 = Fl_mul_pre(t1, h21, p, pi);
     550     4110238 :     t2 = Fl_mul_pre(t0, h10, p, pi);
     551     4110461 :     W[j1] = Fl_neg(t2, p);
     552     4109934 :     t0 = Fl_mul_pre(t1, h11, p, pi);
     553     4112206 :     t2 = Fl_mul_pre(t0, h20, p, pi);
     554     4112034 :     W[j2] = Fl_neg(t2, p);
     555     4111460 :   } while (j2 > j1 + 1);
     556             :   /* Usually the loop exits when j2 = j1 + 1, in which case we return n.
     557             :    * If we break early because of an error, then (j2 - (j1+1)) > 0 is the
     558             :    * number of elements we haven't calculated yet, and we return n minus that
     559             :    * quantity */
     560       92670 :   return gc_long(av, n - j2 + j1 + 1);
     561             : }
     562             : 
     563             : static long
     564        1212 : surface_gcd_path(
     565             :   ulong W[], ulong V[], long n,
     566             :   GEN Phi1, long L1, GEN Phi2, long L2, long e, ulong p, ulong pi)
     567             : {
     568        1212 :   pari_sp av = avma;
     569             :   long i, j;
     570             : 
     571        1212 :   i = 0; j = e;
     572             :   /* If W != V then assume V actually points to a L2-isogenous
     573             :    * parallel L1-path.  e should be 2 in this case */
     574        1212 :   if (W != V) i = j;
     575        4945 :   while (j < n)
     576             :   {
     577             :     GEN f, g, d;
     578             : 
     579        3733 :     f = Flm_Fl_polmodular_evalx(Phi2, L2, V[i], p, pi);
     580        3734 :     g = Flm_Fl_polmodular_evalx(Phi1, L1, W[j - 1], p, pi);
     581        3734 :     g = Flx_remove_root(g, W[j - 2], p);
     582        3733 :     d = Flx_gcd(f, g, p);
     583        3731 :     if (degpol(d) != 1) break; /* Error */
     584        3731 :     W[j] = Flx_deg1_root(d, p);
     585        3732 :     i++; j++; set_avma(av);
     586             :   }
     587        1212 :   return gc_long(av, j);
     588             : }
     589             : 
     590             : /* Given a path V of length n on an L1-volcano, and W[0] L2-isogenous to V[0],
     591             :  * extends the path W to length n on an L1-volcano, with W[i] L2-isogenous
     592             :  * to V[i]. Uses gcds unless L2 is too large to make it helpful. Always uses
     593             :  * GCD to get W[1] to ensure consistent orientation.
     594             :  *
     595             :  * Returns the new length of W. This will almost always be n, but could be
     596             :  * lower if V was started with a J-invariant with bad endomorphism ring */
     597             : INLINE long
     598      199730 : surface_parallel_path(
     599             :   ulong W[], ulong V[], long n,
     600             :   GEN Phi1, long L1, GEN Phi2, long L2, ulong p, ulong pi, long cycle)
     601             : {
     602             :   ulong W2, nbrs[2];
     603      199730 :   if (common_nbr(nbrs, W[0], Phi1, L1, V[1], Phi2, L2, p, pi) == 2)
     604             :   {
     605         721 :     if (n <= 2) return 1; /* Error: Two choices with n = 2; ambiguous */
     606         721 :     if (!common_nbr_verify(&W2,nbrs[0], Phi1,L1,V[2], Phi2,L2,W[0], p,pi))
     607         368 :       nbrs[0] = nbrs[1]; /* nbrs[1] must be the correct choice */
     608         353 :     else if (common_nbr_verify(&W2,nbrs[1], Phi1,L1,V[2], Phi2,L2,W[0], p,pi))
     609          55 :       return 1; /* Error: Both paths extend successfully */
     610             :   }
     611      199673 :   W[1] = nbrs[0];
     612      199673 :   if (n <= 2) return n;
     613       92778 :   return cycle? surface_gcd_cycle(W, V, n, Phi1, L1, Phi2, L2, 2, p, pi)
     614      186768 :               : surface_gcd_path (W, V, n, Phi1, L1, Phi2, L2, 2, p, pi);
     615             : }
     616             : 
     617             : GEN
     618      213935 : enum_roots(ulong J0, norm_eqn_t ne, GEN fdb, GEN G, GEN vshape)
     619             : { /* MAX_HEIGHT >= max_{p,n} val_p(n) where p and n are ulongs */
     620             :   enum { MAX_HEIGHT = BITS_IN_LONG };
     621      213935 :   pari_sp av, ltop = avma;
     622      213935 :   long s = !!pcp_get_L0(G);
     623      213933 :   long *n = pcp_get_n(G)+s, *L = pcp_get_L(G)+s, *o = pcp_get_o(G)+s, k = pcp_get_k(G)-s;
     624      213935 :   long i, t, vlen, *e, *h, *off, *poff, *M, N = pcp_get_enum_cnt(G);
     625      213933 :   ulong p = ne->p, pi = ne->pi, *roots;
     626             :   GEN Phi, vp, ve, roots_;
     627             : 
     628      213933 :   if (!k) return mkvecsmall(J0);
     629             : 
     630      211584 :   roots_ = cgetg(N + MAX_HEIGHT, t_VECSMALL);
     631      211583 :   roots = zv_to_ulongptr(roots_);
     632      211583 :   av = avma;
     633             : 
     634      211583 :   if (!vshape) vshape = factoru(ne->v);
     635      211583 :   vp = gel(vshape, 1); vlen = lg(vp)-1;
     636      211583 :   ve = gel(vshape, 2);
     637             : 
     638      211583 :   Phi = new_chunk(k);
     639      211583 :   e = new_chunk(k);
     640      211584 :   off = new_chunk(k);
     641      211584 :   poff = new_chunk(k);
     642             :   /* TODO: Surely we can work these out ahead of time? */
     643             :   /* h[i] is the valuation of p[i] in v */
     644      211584 :   h = new_chunk(k);
     645      489605 :   for (i = 0; i < k; ++i) {
     646      278023 :     h[i] = 0;
     647      413177 :     for (t = 1; t <= vlen; ++t)
     648      320998 :       if (vp[t] == L[i]) { h[i] = uel(ve, t); break; }
     649      278023 :     e[i] = 0;
     650      278023 :     off[i] = 0;
     651      278023 :     gel(Phi, i) = polmodular_db_getp(fdb, L[i], p);
     652             :   }
     653             : 
     654      211582 :   t = surface_path(roots, n[0], gel(Phi, 0), L[0], h[0], J0, NULL, p, pi);
     655      211585 :   if (t < n[0]) return gc_NULL(ltop); /* J0 has bad endo ring */
     656      210872 :   if (k == 1) { setlg(roots_, t + 1); return gc_const(av,roots_); }
     657             : 
     658       53837 :   M = new_chunk(k);
     659      119093 :   for (M[0] = 1, i = 1; i < k; ++i) M[i] = M[i-1] * n[i-1];
     660       53837 :   i = 1;
     661      253505 :   while (i < k) {
     662             :     long j, t0;
     663      226501 :     for (j = i + 1; j < k && ! e[j]; ++j);
     664      199819 :     if (j < k) {
     665       80298 :       if (e[i]) {
     666       43023 :         if (! common_nbr_pred(
     667       43018 :               &roots[t], roots[off[i]], gel(Phi,i), L[i],
     668       43018 :               roots[t - M[j]], gel(Phi, j), L[j], roots[poff[i]], p, pi)) {
     669           0 :           break; /* J0 has bad endo ring */
     670             :         }
     671       37285 :       } else if ( ! common_nbr_corner(
     672       37280 :             &roots[t], roots[off[i]], gel(Phi,i), L[i], h[i],
     673       37280 :             roots[t - M[j]], gel(Phi, j), L[j], roots[poff[j]], p, pi)) {
     674          34 :         break; /* J0 has bad endo ring */
     675             :       }
     676      184662 :     } else if ( ! next_surface_nbr(
     677      119521 :           &roots[t], gel(Phi,i), L[i], h[i],
     678      184658 :           roots[off[i]], e[i] ? &roots[poff[i]] : NULL, p, pi))
     679           0 :       break; /* J0 has bad endo ring */
     680      199799 :     if (roots[t] == roots[0]) break; /* J0 has bad endo ring */
     681             : 
     682      199742 :     poff[i] = off[i];
     683      199742 :     off[i] = t;
     684      199742 :     e[i]++;
     685      237031 :     for (j = i-1; j; --j) { e[j] = 0; off[j] = off[j+1]; }
     686             : 
     687      199742 :     t0 = surface_parallel_path(&roots[t], &roots[poff[i]], n[0],
     688      199742 :         gel(Phi, 0), L[0], gel(Phi, i), L[i], p, pi, n[0] == o[0]);
     689      199723 :     if (t0 < n[0]) break; /* J0 has bad endo ring */
     690             : 
     691             :     /* TODO: Do I need to check if any of the new roots is a repeat in
     692             :      * the case where J0 has bad endo ring? */
     693      199668 :     t += n[0];
     694      301908 :     for (i = 1; i < k && e[i] == n[i]-1; i++);
     695             :   }
     696       53832 :   if (t != N) return gc_NULL(ltop); /* J0 has wrong endo ring */
     697       53686 :   setlg(roots_, t + 1); return gc_const(av,roots_);
     698             : }

Generated by: LCOV version 1.14