Blob Blame History Raw
--- config/Makefile.SH
+++ config/Makefile.SH
@@ -791,6 +791,9 @@ for dir in basemath modules language gp
     depend="$ANAL_H $src/language/opcode.h"
     cflags="$cflags \$(DLCFLAGS)"
     ;;
+  ratpoints)
+    cflags="$cflags \$(DLCFLAGS) \$(KERNELCFLAGS)"
+    ;;
   *)
     cflags="$cflags \$(DLCFLAGS)"
     ;;
--- src/funclist
+++ src/funclist
@@ -226,6 +226,7 @@
 2702502306 1174 ../functions/elliptic_curves/ellperiods
 2075345511 2181 ../functions/elliptic_curves/ellpointtoz
 1928877975 183 ../functions/elliptic_curves/ellpow
+1259705161 842 ../functions/elliptic_curves/ellratpoints
 2904045120 763 ../functions/elliptic_curves/ellrootno
 3056570066 2988 ../functions/elliptic_curves/ellsea
 1274315519 1818 ../functions/elliptic_curves/ellsearch
@@ -243,6 +244,7 @@
 3073123427 6008 ../functions/elliptic_curves/genus2red
 4179269820 731 ../functions/elliptic_curves/hyperellcharpoly
 1622426669 1038 ../functions/elliptic_curves/hyperellpadicfrobenius
+2980508906 1348 ../functions/elliptic_curves/hyperellratpoints
 1668730040 123 ../functions/gp2c/DEBUGLEVEL
 673657102 528 ../functions/gp2c/clone
 4282475994 385 ../functions/gp2c/copy
--- src/functions/elliptic_curves/ellratpoints
+++ src/functions/elliptic_curves/ellratpoints
@@ -0,0 +1,21 @@
+Function: ellratpoints
+Section: elliptic_curves
+C-Name: ellratpoints
+Prototype: GGD0,L,
+Help: ellratpoints(E,h,{flag=0}): E being an integral model of an
+ elliptic curve, return a vector containing the affine rational points on the curve
+ of naive height less than h.
+ If fl=1, stop as soon as a point is found.
+Doc: $E$ being an integral model of elliptic curve , return a vector
+ containing the affine rational points on the curve of naive height less than
+ $h$. If $\fl=1$, stop as soon as a point is found; return either an empty
+ vector or a vector containing a single point.
+ See \kbd{hyperellratpoints} for how $h$ can be specified.
+ \bprog
+ ? E=ellinit([-25,1]);
+ ? ellratpoints(E,10)
+ %2 = [[-5,1],[-5,-1],[-3,7],[-3,-7],[-1,5],[-1,-5],
+       [0,1],[0,-1],[5,1],[5,-1],[7,13],[7,-13]]
+ ? ellratpoints(E,10,1)
+ %3 = [[-5,1]]
+ @eprog
--- src/functions/elliptic_curves/hyperellratpoints
+++ src/functions/elliptic_curves/hyperellratpoints
@@ -0,0 +1,29 @@
+Function: hyperellratpoints
+Section: elliptic_curves
+C-Name: hyperellratpoints
+Prototype: GGD0,L,
+Help: hyperellratpoints(X,h,{flag=0}): X being a non-singular hyperelliptic
+ curve given by an integral model, return a vector containing the affine rational points on the curve
+ of naive height less than h.
+ If fl=1, stop as soon as a point is found.
+ X can be given either by a squarefree polynomial P such that
+ X:y^2=P(x) or by a vector [P,Q] such that X:y^2+Q(x)y=P(x) and Q^2+4P is
+ squarefree.
+Doc: $X$ being a non-singular hyperelliptic curve given by an integral model,
+ return a vector containing the affine rational points on the curve of naive height less than $h$.
+ If $\fl=1$, stop as soon as a point is found; return either an empty vector or a vector containing a single point.
+
+ $X$ is given either by a squarefree polynomial $P$ such that $X: y^2=P(x)$
+ or by a vector $[P,Q]$ such that $X: y^2+Q(x)\*y=P(x)$ and $Q^2+4\*P$ is
+ squarefree.
+
+ \noindent The parameter $h$ can be
+
+ \item an integer $H$: find the points $[n/d,y]$ whose abscissas $x = n/d$ have
+ naive height (= $\max(|n|, d)$) less than $H$;
+
+ \item a vector $[N,D]$ with $D\leq N$: find the points $[n/d,y]$ with
+ $|n| \leq N$, $d \leq D$.
+
+ \item a vector $[N,[D_1,D_2]]$ with $D_1<D_2\leq N$  find the points
+ $[n/d,y]$ with $|n| \leq N$ and $D_1 \leq d \leq D_2$.
--- src/headers/paridecl.h
+++ src/headers/paridecl.h
@@ -4071,6 +4071,11 @@ ulong   random_Fl(ulong n);
 long    random_bits(long k);
 void    setrand(GEN seed);
 
+/* ratpoints.c */
+
+GEN     ellratpoints(GEN E, GEN h, long flag);
+GEN     hyperellratpoints(GEN P, GEN h, long flag);
+
 /* rootpol.c */
 
 GEN     QX_complex_roots(GEN p, long l);
@@ -4535,6 +4540,7 @@ INLINE long   smodis(GEN x, long y);
 INLINE long   smodss(long x, long y);
 INLINE void   stackdummy(pari_sp av, pari_sp ltop);
 INLINE char  *stack_malloc(size_t N);
+INLINE char  *stack_malloc_align(size_t N, long k);
 INLINE char  *stack_calloc(size_t N);
 INLINE GEN    stoi(long x);
 INLINE GEN    stor(long x, long prec);
--- src/kernel/none/level1.h
+++ src/kernel/none/level1.h
@@ -128,6 +128,14 @@ stack_malloc(size_t N)
 }
 
 INLINE char *
+stack_malloc_align(size_t N, long k)
+{
+  ulong d = ((ulong)avma) % k;
+  if (d) (void)new_chunk(d/sizeof(long));
+  return (char*) new_chunk(nchar2nlong(N));
+}
+
+INLINE char *
 stack_calloc(size_t N)
 {
   char *p = stack_malloc(N);
--- src/language/init.h
+++ src/language/init.h
@@ -372,6 +372,7 @@ entree functions_basic[]={
 {"ellperiods",0,(void*)ellperiods,5,"GD0,L,p","ellperiods(w, {flag = 0}): w describes a complex period lattice ([w1,w2] or an ellinit structure). Returns normalized periods [W1,W2] generating the same lattice such that tau := W1/W2 satisfies Im(tau) > 0 and lies in the standard fundamental domain for SL2. If flag is 1, the return value is [[W1,W2], [eta1,eta2]], where eta1, eta2 are the quasi-periods attached to [W1,W2], satisfying eta1 W2 - eta2 W1 = 2 I Pi."},
 {"ellpointtoz",0,(void*)zell,5,"GGp","ellpointtoz(E,P): lattice point z corresponding to the point P on the elliptic curve E."},
 {"ellpow",0,(void*)ellmul,5,"GGG","ellpow(E,z,n): deprecated alias for ellmul."},
+{"ellratpoints",0,(void*)ellratpoints,12,"GGD0,L,","ellratpoints(E,h,{flag=0}): E being an integral model of an elliptic curve, return a vector containing the affine rational points on the curve of naive height less than h. If fl=1, stop as soon as a point is found."},
 {"ellrootno",0,(void*)ellrootno,5,"lGDG","ellrootno(E,{p}): root number for the L-function of the elliptic curve E/Q at a prime p (including 0, for the infinite place); global root number if p is omitted."},
 {"ellsea",0,(void*)ellsea,5,"GD0,U,","ellsea(E,{tors=0}): computes the order of the group E(Fq) for the elliptic curve E, defined over a finite field, using SEA algorithm, with early abort for curves with non prime orders."},
 {"ellsearch",0,(void*)ellsearch,5,"G","ellsearch(N): returns all curves in the elldata database matching constraint N:  given name (N = \"11a1\" or [11,0,1]), given isogeny class (N = \"11a\" or [11,0]), or given conductor (N = 11, \"11\", or [11])."},
@@ -457,6 +458,7 @@ entree functions_basic[]={
 {"hilbert",0,(void*)hilbert,4,"lGGDG","hilbert(x,y,{p}): Hilbert symbol at p of x,y."},
 {"hyperellcharpoly",0,(void*)hyperellcharpoly,5,"G","hyperellcharpoly(X): X being a non-singular hyperelliptic curve defined over a finite field, return the characteristic polynomial of the Frobenius automorphism.  X can be given either by a squarefree polynomial P such that X:y^2=P(x) or by a vector [P,Q] such that X:y^2+Q(x)*y=P(x) and Q^2+4P is squarefree."},
 {"hyperellpadicfrobenius",0,(void*)hyperellpadicfrobenius,5,"GUL","hyperellpadicfrobenius(Q,p,n): Q being a  rational polynomial of degree d and X being the curve defined by y^2=Q(x), return the matrix of the Frobenius at p>=d in the standard basis of H^1_dR(X) to absolute p-adic precision p^n."},
+{"hyperellratpoints",0,(void*)hyperellratpoints,12,"GGD0,L,","hyperellratpoints(X,h,{flag=0}): X being a non-singular hyperelliptic curve given by an integral model, return a vector containing the affine rational points on the curve of naive height less than h. If fl=1, stop as soon as a point is found. X can be given either by a squarefree polynomial P such that X:y^2=P(x) or by a vector [P,Q] such that X:y^2+Q(x)y=P(x) and Q^2+4P is squarefree."},
 {"hyperu",0,(void*)hyperu,3,"GGGp","hyperu(a,b,x): U-confluent hypergeometric function."},
 {"idealadd",0,(void*)idealadd,8,"GGG","idealadd(nf,x,y): sum of two ideals x and y in the number field defined by nf."},
 {"idealaddtoone",0,(void*)idealaddtoone0,8,"GGDG","idealaddtoone(nf,x,{y}): if y is omitted, when the sum of the ideals in the number field K defined by nf and given in the vector x is equal to Z_K, gives a vector of elements of the corresponding ideals who sum to 1. Otherwise, x and y are ideals, and if they sum up to 1, find one element in each of them such that the sum is 1."},
--- src/modules/ratpoints.c
+++ src/modules/ratpoints.c
@@ -0,0 +1,1749 @@
+/* Copyright (C) 2017  The PARI group.
+
+This file is part of the PARI/GP package.
+
+PARI/GP is free software; you can redistribute it and/or modify it under the
+terms of the GNU General Public License as published by the Free Software
+Foundation. It is distributed in the hope that it will be useful, but WITHOUT
+ANY WARRANTY WHATSOEVER.
+
+Check the License for details. You should have received a copy of it, along
+with the package; see the file 'COPYING'. If not, write to the Free Software
+Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
+
+/* This file is based on ratpoints-2.1.3 by Michael Stoll, see
+ * http://www.mathe2.uni-bayreuth.de/stoll/programs/
+ * Original copyright / license: */
+/***********************************************************************
+ * ratpoints-2.1.3                                                     *
+ * Copyright (C) 2008, 2009  Michael Stoll                             *
+ *  - A program to find rational points on hyperelliptic curves        *
+ *                                                                     *
+ * This program is free software: you can redistribute it and/or       *
+ * modify it under the terms of the GNU General Public License         *
+ * as published by the Free Software Foundation, either version 2 of   *
+ * the License, or (at your option) any later version.                 *
+ ***********************************************************************/
+
+#include "pari.h"
+#include "paripriv.h"
+
+#define pel(a,b)  gel((a),(b)+2)
+
+#define RATPOINTS_ARRAY_SIZE 256           /* Array size in longs */
+#define RATPOINTS_DEFAULT_SP1 9            /* Default value for sp1 */
+#define RATPOINTS_DEFAULT_SP2 16           /* Default value for sp2 */
+#define RATPOINTS_DEFAULT_NUM_PRIMES 28    /* Default value for num_primes */
+#define RATPOINTS_DEFAULT_BIT_PRIMES 7     /* Default value for bit_primes */
+#define RATPOINTS_DEFAULT_MAX_FORBIDDEN 30 /* Default value for max_forbidden */
+
+typedef struct {double low; double up;} ratpoints_interval;
+
+/* Define the flag bits for the flags component: */
+#define RATPOINTS_NO_REVERSE      0x0004UL
+
+#define RATPOINTS_FLAGS_INPUT_MASK (RATPOINTS_NO_REVERSE)
+
+/* Flags bits for internal purposes */
+#define RATPOINTS_REVERSED        0x0100UL
+#define RATPOINTS_CHECK_DENOM     0x0200UL
+#define RATPOINTS_USE_SQUARES     0x0400UL
+#define RATPOINTS_USE_SQUARES1    0x0800UL
+
+#define LONG_MASK (~(-(1UL<<TWOPOTBITS_IN_LONG)))
+
+#define CEIL(a,b) (((a) <= 0) ? -(-(a) / (b)) : 1 + ((a)-1) / (b))
+
+/* Some Macros for working with SSE registers */
+#ifdef HAS_SSE2
+#include <emmintrin.h>
+
+#define AND(a,b) ((ratpoints_bit_array)__builtin_ia32_andps((__v4sf)(a), (__v4sf)(b)))
+#define EXT0(a) ((ulong)__builtin_ia32_vec_ext_v2di((__v2di)(a), 0))
+#define EXT1(a) ((ulong)__builtin_ia32_vec_ext_v2di((__v2di)(a), 1))
+#define TEST(a) (EXT0(a) || EXT1(a))
+#define RBA(a) ((__v2di){((long) a), ((long) a)})
+
+/* Use SSE 128 bit registers for the bit arrays */
+typedef __v2di ratpoints_bit_array;
+
+#define RBA_LENGTH (128)
+#define RBA_SHIFT (7)
+#define RBA_ALIGN  (sizeof(ratpoints_bit_array))
+#define RBA_MASK (~(-(1UL<<RBA_SHIFT)))
+
+#else
+/* Use ulong for the bit arrays */
+typedef ulong ratpoints_bit_array;
+#define AND(a,b) ((a)&(b))
+#define EXT0(a) (a)
+#define TEST(a) (a)
+#define RBA(a) (a)
+
+#define RBA_LENGTH BITS_IN_LONG
+#define RBA_SHIFT TWOPOTBITS_IN_LONG
+#define RBA_ALIGN  (sizeof(long))
+#define RBA_MASK LONG_MASK
+
+#endif
+
+typedef struct { long p; long offset; ratpoints_bit_array *ptr; } sieve_spec;
+
+typedef enum { num_all, num_even, num_odd, num_none } bit_selection;
+
+typedef struct {
+  long p; int *is_f_square;
+  const long *inverses;
+  long offset; ratpoints_bit_array** sieve;
+} ratpoints_sieve_entry;
+
+typedef struct { long p;
+                 ulong *start;
+                 ulong *end;
+                 ulong *curr; }
+               forbidden_entry;
+
+typedef struct {
+  GEN cof, listprime;
+  ratpoints_interval *domain;
+  long height, b_low, b_high, sp1, sp2, array_size;
+  long num_inter, num_primes, bit_primes, max_forbidden;
+  ulong flags;
+/* from here: private data */
+  GEN bc;
+  ratpoints_sieve_entry *se_buffer;
+  ratpoints_sieve_entry *se_next;
+  ratpoints_bit_array *ba_buffer;
+  ratpoints_bit_array *ba_next;
+  int *int_buffer, *int_next;
+  forbidden_entry *forb_ba;
+  long *forbidden;
+  GEN inverses, offsets, den_info, divisors;
+  ulong **sieves0;
+} ratpoints_args;
+
+#ifdef HAS_SSE2
+#define CODE_INIT_SIEVE_COPY for (a = 0; a < p; a++) si[a+p] = si[a];
+#else
+#define CODE_INIT_SIEVE_COPY
+#endif
+
+static ratpoints_bit_array *
+#if defined(__GNUC__)
+__attribute__((noinline))
+#endif
+sieve_init1(long p, ratpoints_sieve_entry *se1, long b1, ratpoints_args *args1)
+{
+  ratpoints_sieve_entry *se = se1;
+  ratpoints_args *args = args1;
+  register int *isfs = se->is_f_square;
+  register long b = b1;
+  long lmp = BITS_IN_LONG % p;
+  long ldp = BITS_IN_LONG / p;
+  long p1 = (ldp + 1) * p;
+  long diff_shift = p1 & LONG_MASK;
+  long diff = BITS_IN_LONG - diff_shift;
+  register ulong help0;
+  register long a;
+  register long d = se->inverses[b];
+  register long ab = 0; /* a/b mod p */
+  register ulong test = 1UL;
+  register ulong he0 = 0UL;
+  for (a = 0; a < p; a++)
+  {
+    if (isfs[ab]) he0 |= test;
+    ab += d;
+    if (ab >= p) ab -= p;
+    test <<= 1;
+  }
+  help0 = he0;
+  {
+    register ulong help1;
+     /* repeat bit pattern floor(BITS_IN_LONG/p) times */
+    register ulong pattern = help0;
+    register long i;
+    /* the p * (floor(BITS_IN_LONG/p) + 1) - BITS_IN_LONG
+            = p - (BITS_IN_LONG mod p)
+       upper bits into help[b][1] :
+       shift away the  BITS_IN_LONG mod p  lower bits */
+    help1 = pattern >> lmp;
+    for (i = p; i < BITS_IN_LONG; i <<= 1)
+      help0 |= help0 << i;
+    { /* fill the bit pattern from help0/help1 into sieve[b][].
+          sieve[b][a0] has the same semantics as help0/help1,
+          but here, a0 runs from 0 to p-1 and all bits are filled. */
+      register long a;
+      ulong *si = (ulong *)args->ba_next;
+
+      args->ba_next += p;
+      /* copy the first chunk into sieve[b][] */
+      si[0] = help0;
+      /* now keep repeating the bit pattern,
+         rotating it in help0/help1 */
+      for (a = 1 ; a < p; a++)
+      {
+        register ulong temp = help0 >> diff;
+        help0 = help1 | (help0 << diff_shift);
+        si[a] = help0;
+        help1 = temp;
+      }
+      CODE_INIT_SIEVE_COPY
+      /* set sieve array */
+      se->sieve[b] = (ratpoints_bit_array *)si;
+      return (ratpoints_bit_array *)si;
+    }
+  }
+}
+
+/* This is for p > BITS_IN_LONG */
+static ratpoints_bit_array *
+#if defined(__GNUC__)
+__attribute__((noinline))
+#endif
+sieve_init2(long p, ratpoints_sieve_entry *se1, long b1, ratpoints_args *args1)
+{
+  ratpoints_sieve_entry *se = se1;
+  ratpoints_args *args = args1;
+  register int *isfs = se->is_f_square;
+  register long b = b1;
+  /* long ldp = 0;  = BITS_IN_LONG / p */
+  /* long p1 = p; = (ldp + 1) * p; */
+  long wp = p >> TWOPOTBITS_IN_LONG;
+  long diff_shift = p & LONG_MASK;
+  long diff = BITS_IN_LONG - diff_shift;
+  ulong help[(p>>TWOPOTBITS_IN_LONG) + 2];
+
+  /* initialize help */
+  {
+    register ulong *he = &help[0];
+    register ulong *he1 = &he[(p>>TWOPOTBITS_IN_LONG) + 2];
+    while (he1 != he) { he1--; *he1 = 0UL; }
+  }
+  { register ulong work = 0UL;
+    register long a;
+    register long ab = 0; /* a/b mod p */
+    register long d = se->inverses[b];
+    register long n = 0;
+    register ulong test = 1UL;
+    for (a = 0; a < p; )
+    {
+      if (isfs[ab]) work |= test;
+      ab += d;
+      if (ab >= p) ab -= p;
+      test <<= 1;
+      a++;
+      if ((a & LONG_MASK) == 0)
+      { help[n] = work; n++; work = 0UL; test = 1UL; }
+    }
+    help[n] = work;
+  }
+
+  { /* fill the bit pattern from help[] into sieve[b][].
+       sieve[b][a0] has the same semantics as help[b][a0],
+       but here, a0 runs from 0 to p-1 and all bits are filled. */
+    register ulong *si = (ulong *)args->ba_next;
+    register long a1;
+    register long a;
+
+    args->ba_next += p;
+    /* copy the first chunk from help[] into sieve[num][b][] */
+    for (a = 0; a < wp; a++) si[a] = help[a];
+    /* now keep repeating the bit pattern, rotating it in help */
+    for (a1 = a ; a < p; a++)
+    {
+      register long t = (a1 == wp) ? 0 : a1+1;
+      help[a1] |= help[t]<<diff_shift;
+      si[a] = help[a1];
+      a1 = t;
+      help[a1] >>= diff;
+    }
+     CODE_INIT_SIEVE_COPY
+    /* set sieve array */
+    se->sieve[b] = (ratpoints_bit_array *)si;
+    return (ratpoints_bit_array *)si;
+  }
+}
+
+static GEN
+gen_squares(GEN listprime)
+{
+  long nbprime = lg(listprime)-1;
+  GEN sq = cgetg(nbprime+1, t_VEC);
+  long n;
+  for (n = 1; n <= nbprime; n++)
+  {
+    ulong i, p = uel(listprime,n);
+    GEN w = zero_zv(p), work = w+1;
+    work[0] = 1;
+    /* record non-zero squares mod p, p odd */
+    for (i = 1; i < p; i += 2) work[(i*i) % p] = 1;
+    gel(sq, n) = w;
+  }
+  return sq;
+}
+
+static GEN
+gen_offsets(GEN listprime)
+{
+  long n;
+  long nbprime = lg(listprime)-1;
+  GEN of = cgetg(nbprime+1, t_VEC);
+  for (n = 1; n <= nbprime; n++)
+  {
+    ulong p = uel(listprime,n);
+    uel(of, n) = Fl_inv((2*RBA_LENGTH)%p, p);
+  }
+  return of;
+}
+
+static GEN
+gen_inverses(GEN listprime)
+{
+  long n;
+  long nbprime = lg(listprime)-1;
+  GEN iv = cgetg(nbprime+1, t_VEC);
+  for (n = 1; n <= nbprime; n++)
+  {
+    ulong i, p = uel(listprime,n);
+    GEN w = cgetg(p, t_VECSMALL);
+    for (i = 1; i < p; i++)
+      uel(w, i) = Fl_inv(i, p);
+    gel(iv, n) = w;
+  }
+  return iv;
+}
+
+static ulong **
+gen_sieves0(GEN listprime)
+{
+  long n;
+  long nbprime = lg(listprime)-1;
+  ulong ** si = (ulong**) new_chunk(nbprime+1);
+  for (n = 1; n <= nbprime; n++)
+  {
+    ulong i, p = uel(listprime,n);
+    ulong *w = (ulong *) stack_malloc_align(2*p*sizeof(ulong), RBA_ALIGN);
+    for (i = 0; i < p; i++) uel(w,i) = ~0UL;
+    for (i = 0; i < BITS_IN_LONG; i++)
+      uel(w,(p*i)>>TWOPOTBITS_IN_LONG) &= ~(1UL<<((p*i) & LONG_MASK));
+    for (i = 0; i < p; i++) uel(w,i+p) = uel(w,i);
+    si[n] = w;
+  }
+  return si;
+}
+
+static void
+gen_sieve(ratpoints_args *args)
+{
+  GEN listprimes = args->listprime;
+  args->offsets = gen_offsets(listprimes);
+  args->inverses = gen_inverses(listprimes);
+  args->sieves0 = gen_sieves0(listprimes);
+}
+
+static GEN
+ZX_positive_region(GEN P, long h, long bitprec)
+{
+  long prec = nbits2prec(bitprec);
+  GEN it = mkvec2(stoi(-h),stoi(h));
+  GEN R = ZX_Uspensky(P, it, 1, bitprec);
+  long nR = lg(R)-1;
+  long s = signe(poleval(P,gel(it,1)));
+  long i=1, j;
+  GEN iv, st, en;
+  if (s<0 && nR==0) return NULL;
+  iv = cgetg(((nR+1+(s>=0))>>1)+1, t_VEC);
+  if (s>=0) st = itor(gel(it,1),prec);
+  else    { st = gel(R,i); i++; }
+  for (j=1; i<nR; j++)
+  {
+    gel(iv, j) = mkvec2(st, gel(R,i));
+    st = gel(R,i+1);
+    i+=2;
+  }
+  if (i==nR) en = gel(R,i); else en = itor(gel(it,2),prec);
+  gel(iv,j) = mkvec2(st, en);
+  return iv;
+}
+
+static long
+posint(ratpoints_interval *ivlist, GEN P, long h)
+{
+  GEN R = ZX_positive_region(P, h, 53);
+  const double eps = 1e-5;
+  long nR, i;
+
+  if (!R) return 0;
+  nR = lg(R)-1;
+  i = 1;
+  for (i=1; i<=nR; i++)
+  {
+    ivlist[i-1].low = rtodbl(gmael(R,i,1))-eps;
+    ivlist[i-1].up  = rtodbl(gmael(R,i,2))+eps;
+  }
+  return nR;
+}
+
+static long
+ratpoints_compute_sturm(ratpoints_args *args)
+{
+  ratpoints_interval *ivlist = args->domain;
+  args->num_inter = posint(ivlist, args->cof, (long) ivlist[0].up);
+  return args->num_inter;
+}
+
+/**************************************************************************
+ * Try to avoid divisions                                                 *
+ **************************************************************************/
+INLINE long
+mod(long a, long b)
+{
+  long b1 = b << 4; /* b1 = 16*b */
+
+  if (a < -b1) { a %= b; if (a < 0) { a += b; } return a ; }
+  if (a < 0) { a += b1; }
+  else { if (a >= b1) { return a % b; } }
+  b1 >>= 1; /* b1 = 8*b */
+  if (a >= b1) { a -= b1; }
+  b1 >>= 1; /* b1 = 4*b */
+  if (a >= b1) { a -= b1; }
+  b1 >>= 1; /* b1 = 2*b */
+  if (a >= b1) { a -= b1; }
+  if (a >= b) { a -= b; }
+  return a;
+}
+
+static void
+set_bc(long b, ratpoints_args *args)
+{
+  GEN w0 = gen_1;
+  GEN c = args->cof, bc;
+  long k, degree = degpol(c);
+  bc = cgetg(degree+2, t_POL);
+  for (k = degree-1; k >= 0; k--)
+  {
+    w0 = muliu(w0, b);
+    gel(bc,k+2) = mulii(gel(c,k+2), w0);
+  }
+  args->bc = bc;
+}
+
+/**************************************************************************
+ * Check a `survivor' of the sieve if it really gives a point.            *
+ **************************************************************************/
+
+static long
+_ratpoints_check_point(long a, long b, ratpoints_args *args, int *quit,
+                 int process(long, long, GEN, void*, int*), void *info)
+{
+  pari_sp av = avma;
+  GEN w0, w2, c = args->cof, bc = args->bc;
+  long k, degree = degpol(c);
+  int reverse = args->flags & RATPOINTS_REVERSED;
+
+  /* Compute F(a, b), where F is the homogenized version of f
+     of smallest possible even degree  */
+  w2 = gel(c, degree+2);
+  for (k = degree-1; k >= 0; k--)
+  {
+    w2 = mulis(w2, a);
+    w2 = addii(w2, gel(bc,k+2));
+  }
+  if (odd(degree)) w2 = muliu(w2, b);
+  /* check if f(x,z) is a square; if so, process the point(s) */
+  if (signe(w2) >= 0 && Z_issquareall(w2, &w0))
+  {
+    if (reverse)
+    {
+      if (a >= 0) (void)process(b, a, w0, info, quit);
+      else        (void)process(-b, -a, w0, info, quit);
+    }
+    else (void)process(a, b, w0, info, quit);
+    if (!*quit && signe(w0) != 0)
+    {
+      GEN nw0 = negi(w0);
+      if (reverse)
+      {
+        if (a >= 0) (void)process(b, a, nw0, info, quit);
+        else        (void)process(-b, -a, nw0, info, quit);
+      }
+      else (void)process(a, b, nw0, info, quit);
+    }
+    return 1;
+  }
+  avma = av;
+  return 0;
+}
+
+/**************************************************************************
+ * The inner loop of the sieving procedure                                *
+ **************************************************************************/
+static long
+_ratpoints_sift0(long b, long w_low, long w_high,
+           ratpoints_args *args, bit_selection which_bits,
+           ratpoints_bit_array *survivors, sieve_spec *sieves, int *quit,
+           int process(long, long, GEN, void*, int*), void *info)
+{
+  long range = w_high - w_low;
+  long sp1 = args->sp1;
+  long sp2 = args->sp2;
+  long i, n, nb = 0, absb = labs(b);
+  ratpoints_bit_array *surv0;
+
+  /* now do the sieving (fast!) */
+  for (n = 0; n < sp1; n++)
+  {
+    ratpoints_bit_array *sieve_n = sieves[n].ptr;
+    register long p = sieves[n].p;
+    long r = mod(-w_low-sieves[n].offset, p);
+    register ratpoints_bit_array *surv = survivors;
+
+    if (w_high < w_low + r)
+    { /* if we get here, r > 0, since w_high >= w_low always */
+      register ratpoints_bit_array *siv1 = &sieve_n[p-r];
+      register ratpoints_bit_array *siv0 = siv1 + range;
+
+      while (siv1 != siv0) { *surv = AND(*surv, *siv1++); surv++; }
+    }
+    else
+    {
+      register ratpoints_bit_array *siv1 = &sieve_n[p-r];
+      register ratpoints_bit_array *surv_end = &survivors[range - p];
+      register long i;
+      for (i = r; i; i--) { *surv = AND(*surv, *siv1++); surv++; }
+      siv1 -= p;
+      while (surv <= surv_end)
+      {
+        for (i = p; i; i--) { *surv = AND(*surv, *siv1++); surv++; }
+        siv1 -= p;
+      }
+      surv_end += p;
+      while (surv < surv_end) { *surv = AND(*surv, *siv1++); surv++; }
+    }
+  } /* for n */
+  /* 2nd phase of the sieve: test each surviving bit array with more primes */
+  surv0 = &survivors[0];
+  for (i = w_low; i < w_high; i++)
+  {
+    register ratpoints_bit_array nums = *surv0++;
+    sieve_spec *ssp = &sieves[sp1];
+    register long n;
+
+    for (n = sp2-sp1; n && TEST(nums); n--)
+    {
+      register long p = ssp->p;
+      nums = AND(nums, ssp->ptr[mod(i + ssp->offset, p)]);
+      ssp++;
+    }
+
+    /* Check the survivors of the sieve if they really give points */
+    if (TEST(nums))
+    {
+      long a0, a, d;
+#ifdef HAS_SSE2
+      long da = (which_bits == num_all) ? RBA_LENGTH/2 : RBA_LENGTH;
+#endif
+           /* a will be the numerator corresponding to the selected bit */
+      if (which_bits == num_all)
+      {
+        d = 1; a0 = i << RBA_SHIFT;
+      }
+      else
+      {
+        d = 2; a0 = i << (RBA_SHIFT+1);
+        if (which_bits == num_odd) a0++;
+      }
+      {
+#ifdef HAS_SSE2
+        ulong nums0 = EXT0(nums);
+        ulong nums1 = EXT1(nums);
+#else
+        ulong nums0 = nums;
+#endif
+
+        for (a = a0; nums0; a += d, nums0 >>= 1)
+        { /* test one bit */
+          if (odd(nums0) && ugcd(labs(a), absb)==1)
+          {
+            if (!args->bc) set_bc(b, args);
+            nb += _ratpoints_check_point(a, b, args, quit, process, info);
+            if (*quit) return nb;
+          }
+        }
+#ifdef HAS_SSE2
+        for (a = a0 + da; nums1; a += d, nums1 >>= 1)
+        { /* test one bit */
+          if (odd(nums1) && ugcd(labs(a), absb)==1)
+          {
+            if (!args->bc) set_bc(b, args);
+            nb += _ratpoints_check_point(a, b, args, quit, process, info);
+            if (*quit) return nb;
+          }
+        }
+#endif
+      }
+    }
+  }
+  return nb;
+}
+
+#define MAX_DIVISORS 512
+ /* Maximal length of array for squarefree divisors of leading coefficient */
+
+typedef struct { double r; ratpoints_sieve_entry *ssp; } entry;
+
+static const int squares16[16] = {1,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0};
+ /* Says if a is a square mod 16, for a = 0..15 */
+
+/**************************************************************************
+ * Initialization and cleanup of ratpoints_args structure                 *
+ **************************************************************************/
+
+static ratpoints_sieve_entry*
+alloc_sieve(long nbprime, long maxprime)
+{
+  long i;
+  ratpoints_sieve_entry * s = (ratpoints_sieve_entry*)
+                        stack_malloc(nbprime*sizeof(ratpoints_sieve_entry));
+  for (i=0; i<nbprime; i++)
+    s[i].sieve = (ratpoints_bit_array**) new_chunk(maxprime);
+  return s;
+}
+
+/* NOTE: args->degree must be set */
+static void
+find_points_init(ratpoints_args *args, long bit_primes)
+{
+  long need = 0;
+  long n, nbprime,maxprime;
+  args->listprime = primes_interval_zv(3, 1<<bit_primes);
+  nbprime = lg(args->listprime)-1;
+  maxprime = args->listprime[nbprime];
+
+  /* allocate space for se_buffer */
+  args->se_buffer = alloc_sieve(nbprime, maxprime);
+  args->se_next = args->se_buffer;
+  for (n = 1; n <= nbprime; n++)
+  {
+    ulong p = args->listprime[n];
+    need += p*p;
+  }
+  args->ba_buffer = (ratpoints_bit_array*)
+     stack_malloc_align(need*sizeof(ratpoints_bit_array),RBA_ALIGN);
+  args->ba_next = args->ba_buffer;
+
+  /* allocate space for int_buffer */
+  args->int_buffer = (int *) stack_malloc(nbprime*(maxprime+1)*sizeof(int));
+  args->int_next = args->int_buffer;
+
+  args->forb_ba   = (forbidden_entry*)
+    stack_malloc((nbprime + 1)*sizeof(forbidden_entry));
+  args->forbidden = new_chunk(nbprime + 1);
+  gen_sieve(args);
+  return;
+}
+
+/* f = leading coeff; b = b1*b2, b1 maximal with (b1, 2*f) = 1;
+ * return Jacobi symbol (f, b1) */
+INLINE int
+rpjacobi(long b, GEN lcf)
+{
+  ulong f;
+  b >>= vals(b);
+  f = umodiu(lcf, b);
+  return krouu(f, u_ppo(b,f));
+}
+
+/************************************************************************
+ * Set up information on possible denominators                          *
+ * when polynomial is of odd degree with leading coefficient != +-1     *
+ ************************************************************************/
+
+static void
+setup_us1(ratpoints_args *args, GEN w0)
+{
+  GEN divisors;
+  GEN F = Z_factor_limit(w0, 1000), P = gel(F,1), E = gel(F,2), S;
+  long i, l = lg(P);
+  if (cmpiu(gel(P,l-1), 1000) > 0)
+    return;
+  P = ZV_to_zv(P); E = ZV_to_zv(E);
+  divisors  = cgetg(1+(1<<(l-1)), t_VECSMALL);
+  /* factorization is complete, set up array of squarefree divisors */
+  divisors[1] = 1;
+  for (i = 1; i < l; i++)
+  { /* multiply all divisors known so far by next prime */
+    long k, n = 1<<(i-1);
+    for (k=0; k<n; k++)
+      uel(divisors,1+n+k) = uel(divisors,1+k) * P[i];
+  }
+  S = cgetg(l, t_VECSMALL);
+  /* set slopes in den_info */
+  for (i = 1; i < l; i++)
+  { /* compute min{n : (d-k)*n > v_p(f_d) - v_p(f_k), k = 0,...,d-1} */
+    GEN c = args->cof;
+    long p = P[i], v = E[i];
+    long k, n = 1, d = degpol(c);
+
+    for (k = d - 1; k >= 0; k--)
+    {
+      long t = 1 + v - Z_lval(gel(c,k+2), p);
+      long m = CEIL(t, (d - k));
+
+      if (m > n) n = m;
+    }
+    S[i] = n;
+  }
+  args->divisors = divisors;
+  args->flags |= RATPOINTS_USE_SQUARES1;
+  args->den_info = mkvec3(P, E, S);
+}
+
+/************************************************************************
+ * Consider 2-adic information                                          *
+ ************************************************************************/
+
+static bit_selection
+get_2adic_info(ratpoints_args *args, ulong *den_bits,
+               ratpoints_bit_array *num_bits)
+{
+  GEN c = args->cof;
+  long degree = degpol(c);
+  int is_f_square16[24];
+  long *cmp = new_chunk(degree+1);
+  long npe = 0, npo = 0;
+  bit_selection result;
+  long n, a, b;
+
+  /* compute coefficients mod 16 */
+  for (n = 0; n <= degree; n++) cmp[n] = Mod16(gel(c,n+2));
+  for (a = 0 ; a < 16; a++)
+  {
+    ulong s = cmp[degree];
+    long n;
+    for (n = degree - 1 ; n >= 0 ; n--) s = s*a + cmp[n];
+    s &= 0xf;
+    if ((is_f_square16[a] = squares16[s])) { if (odd(a)) npo++; else npe++; }
+  }
+
+  /* even denominators:
+     is_f_square16[16+k] says if f((2k+1)/2) is a square, k = 0..3
+     is_f_square16[20+k] says if f((2k+1)/4) is a square, k = 0,1
+     is_f_square16[22]   says if f(odd/8) is a square
+     is_f_square16[23]   says if f(odd/2^n), n >= 4, can be a square */
+  {
+    long np1 = 0, np2 = 0, np3 = 0, np4 = 0;
+
+    if (odd(degree))
+    {
+      long a, cf = 4*cmp[degree-1];
+
+      if (degree >= 2) cf += 8*cmp[degree-2];
+      for (a = 0; a < 4; a++)
+      { /* Compute  2 c[d] k^d + 4 c[d-1] k^(d-1) + 8 c[d-2] k^(d-2), k = 2a+1.
+           Note that k^d = k mod 8, k^(d-1) = 1 mod 8. */
+        long k = 2*a+1;
+        long s = (2*k*cmp[degree] + cf) & 0xf;
+        if ((is_f_square16[16+a] = squares16[s])) np1++;
+      }
+      if ((is_f_square16[20] = squares16[(4*cmp[degree])  & 0xf])) np2++;
+      if ((is_f_square16[21] = squares16[(12*cmp[degree]) & 0xf])) np2++;
+      if ((is_f_square16[22] = squares16[(8*cmp[degree])  & 0xf])) np3++;
+      is_f_square16[23] = 1; np4++;
+    }
+    else
+    {
+      long a, cf = (degree >= 2) ? 4*cmp[degree-2] : 0;
+
+      if (degree >= 3) cf += 8*cmp[degree-3];
+      for (a = 0; a < 4; a++)
+      { /* compute c[d] k^d + 2 c[d-1] k^(d-1) + ... + 8 c[d-3] k^(d-3),
+           k = 2a+1. Note that k^d = k^2 mod 16, k^(d-1) = k mod 8. */
+        long k = 2*a+1;
+        long s = ((cmp[degree]*k + 2*cmp[degree-1])*k + cf) & 0xf;
+        if ((is_f_square16[16+a] = squares16[s])) np1++;
+      }
+      if ((is_f_square16[20] = squares16[(cmp[degree]+4*cmp[degree-1])  & 0xf]))
+        np2++;
+      if ((is_f_square16[21] = squares16[(cmp[degree]+12*cmp[degree-1]) & 0xf]))
+        np2++;
+      if ((is_f_square16[22] = squares16[(cmp[degree]+8*cmp[degree-1])  & 0xf]))
+        np3++;
+      if ((is_f_square16[23] = squares16[cmp[degree]])) np4++;
+    }
+
+    /* set den_bits */
+    { ulong db = 0;
+      long i;
+
+      if (npe + npo > 0) db |= 0xaaaaUL; /* odd denominators */
+      if (np1 > 0)       db |= 0x4444UL; /* v_2(den) = 1 */
+      if (np2 > 0)       db |= 0x1010UL; /* v_2(den) = 2 */
+      if (np3 > 0)       db |= 0x0100UL; /* v_2(den) = 3 */
+      if (np4 > 0)       db |= 0x0001UL; /* v_2(den) >= 4 */
+      if (db == 0) { *den_bits = 0UL; return num_none; }
+
+      for (i = 16; i < BITS_IN_LONG; i <<= 1) db |= db << i;
+      *den_bits = db;
+    }
+    result = (npe == 0) ? (npo == 0) ? num_none : num_odd
+                        : (npo == 0) ? num_even : num_all;
+  }
+
+  /* set up num_bits[16] */
+
+  /* odd denominators */
+  switch(result)
+  {
+    case num_all:
+      for (b = 1; b < 16; b += 2)
+      {
+        ulong work = 0, bit = 1;
+        long i, invb = b; /* inverse of b mod 16 */
+        if (b & 2) invb ^= 8;
+        if (b & 4) invb ^= 8;
+        for (i = 0; i < 16; i++)
+        {
+          if (is_f_square16[(invb*i) & 0xf]) work |= bit;
+          bit <<= 1;
+        }
+        /* now repeat the 16 bits */
+        for (i = 16; i < BITS_IN_LONG; i <<= 1) work |= work << i;
+        num_bits[b] = RBA(work);
+      }
+      break;
+
+    case num_odd:
+      for (b = 1; b < 16; b += 2)
+      {
+        ulong work = 0, bit = 1;
+        long i, invb = b; /* inverse of b mod 16 */
+        if (b & 2) invb ^= 8;
+        if (b & 4) invb ^= 8;
+        for (i = 1; i < 16; i += 2)
+        {
+          if (is_f_square16[(invb*i) & 0xf]) work |= bit;
+          bit <<= 1;
+        }
+        /* now repeat the 8 bits */
+        for (i = 8; i < BITS_IN_LONG; i <<= 1) { work |= work << i; }
+        num_bits[b] = RBA(work);
+      }
+      break;
+
+    case num_even:
+      for (b = 1; b < 16; b += 2)
+      {
+        ulong work = 0, bit = 1;
+        long i, invb = b; /* inverse of b mod 16 */
+        if (b & 2) invb ^= 8;
+        if (b & 4) invb ^= 8;
+        for (i = 0; i < 16; i += 2)
+        {
+          if (is_f_square16[(invb*i) & 0xf]) work |= bit;
+          bit <<= 1;
+        }
+        /* now repeat the 8 bits */
+        for (i = 8; i < BITS_IN_LONG; i <<= 1) work |= work << i;
+        num_bits[b] = RBA(work);
+      }
+      break;
+
+    case num_none:
+      for (b = 1; b < 16; b += 2) num_bits[b] = RBA(0UL);
+      break;
+  }
+
+  /* v_2(den) = 1 : only odd numerators */
+  for (b = 1; b < 8; b += 2)
+  {
+    ulong work = 0, bit = 1;
+    long i;
+    for (i = 1; i < 16; i += 2)
+    {
+      if (is_f_square16[16 + (((b*i)>>1) & 0x3)]) work |= bit;
+      bit <<= 1;
+    }
+    /* now repeat the 8 bits */
+    for (i = 8; i < BITS_IN_LONG; i <<= 1) work |= work << i;
+    num_bits[2*b] = RBA(work);
+  }
+
+  /* v_2(den) = 2 : only odd numerators */
+  for (b = 1; b < 4; b += 2)
+  {
+    ulong work = 0, bit = 1;
+    long i;
+    for (i = 1; i < 8; i += 2)
+    {
+      if (is_f_square16[20 + (((b*i)>>1) & 0x1)]) work |= bit;
+      bit <<= 1;
+    }
+    /* now repeat the 4 bits */
+    for (i = 4; i < BITS_IN_LONG; i <<= 1) work |= work << i;
+    num_bits[4*b] = RBA(work);
+  }
+
+  /* v_2(den) = 3, >= 4 : only odd numerators */
+  num_bits[8] = (is_f_square16[22]) ? RBA(~(0UL)) : RBA(0UL);
+  num_bits[0] = (is_f_square16[23]) ? RBA(~(0UL)) : RBA(0UL);
+
+  return result;
+}
+
+/**************************************************************************
+ * This is a comparison function needed for sorting in order to determine *
+ * the `best' primes for sieving.                                         *
+ **************************************************************************/
+
+static int
+compare_entries(const void *a, const void *b)
+{
+  double diff = ((entry *)a)->r - ((entry *)b)->r;
+  return (diff > 0) ? 1 : (diff < 0) ? -1 : 0;
+}
+
+/************************************************************************
+ * Collect the sieving information                                      *
+ ************************************************************************/
+
+static long
+sieving_info(ratpoints_args *args,
+             ratpoints_sieve_entry **sieve_list)
+{
+  GEN c = args->cof;
+  GEN inverses = args->inverses, squares;
+  GEN offsets = args->offsets;
+  ulong ** sieves0 = args->sieves0;
+  long degree = degpol(c);
+  long fba = 0, fdc = 0;
+  long pn, pnp = 0;
+  long n;
+  long nbprime = lg(args->listprime)-1;
+  entry *prec = (entry*) stack_malloc(nbprime*sizeof(entry));
+    /* This array is used for sorting in order to
+       determine the `best' sieving primes. */
+
+  forbidden_entry *forb_ba = args->forb_ba;
+  long *forbidden = args->forbidden;
+  ulong bound = (1UL)<<(BITS_IN_LONG - args->bit_primes);
+  pari_sp av = avma;
+  squares = gen_squares(args->listprime);
+
+  /* initialize sieve in se_buffer */
+  for (pn = 1; pn <= args->num_primes; pn++)
+  {
+    long coeffs_mod_p[degree+1]; /* The coefficients of f reduced modulo p */
+    ulong a, p = args->listprime[pn], np;
+    long n;
+    int *is_f_square = args->int_next;
+
+    args->int_next += p + 1; /* need space for (p+1) int's */
+
+    /* compute coefficients mod p */
+    for (n = 0; n <= degree; n++) coeffs_mod_p[n] = umodiu(pel(c,n), p);
+
+    np = umael(squares,pn,coeffs_mod_p[0]+1);
+    is_f_square[0] = np;
+    for (a = 1 ; a < p; a++)
+    {
+      ulong s = coeffs_mod_p[degree];
+      if ((degree+1)*args->bit_primes <= BITS_IN_LONG)
+      {
+        for (n = degree - 1 ; n >= 0 ; n--) s = s*a + coeffs_mod_p[n];
+        /* here, s < p^(degree+1) <= max. long */
+        s %= p;
+      }
+      else
+      {
+        for (n = degree - 1 ; n >= 0 ; n--)
+        {
+          s = s*a + coeffs_mod_p[n];
+          if (s+1 >= bound) s %= p;
+        }
+        s %= p;
+      }
+      if ((is_f_square[a] = mael(squares,pn,s+1))) np++;
+    }
+    is_f_square[p] = odd(degree) || mael(squares,pn,coeffs_mod_p[degree]+1);
+
+    /* check if there are no solutions mod p */
+    if (np == 0 && !is_f_square[p]) { avma = av; return p; }
+
+    /* Fill arrays with info for p */
+    if (np < p)
+    { /* only when there is some information */
+      ulong i;
+      ratpoints_sieve_entry *se = args->se_next;
+      double r = is_f_square[p] ? ((double)(np*(p-1) + p))/((double)(p*p))
+                                : (double)np/(double)p;
+      prec[pnp].r = r;
+      args->se_next ++;
+      se->p = p;
+      se->is_f_square = is_f_square;
+      se->inverses = gel(inverses,pn);
+      se->offset = offsets[pn];
+      se->sieve[0] = (ratpoints_bit_array *)sieves0[pn];
+      for (i = 1; i < p; i++) se->sieve[i] = NULL;
+      prec[pnp].ssp = se;
+      pnp++;
+    }
+
+    if ((args->flags & RATPOINTS_CHECK_DENOM)
+         && fba + fdc < args->max_forbidden
+         && !is_f_square[p])
+    { /* record forbidden divisors of the denominator */
+      if (coeffs_mod_p[degree] == 0)
+      { /* leading coeff. divisible by p */
+        GEN r;
+        long v = Z_lvalrem(pel(c,degree), p, &r);
+
+        if (odd(v) || !mael(squares,pn, umodiu(r,p)+1))
+        { /* Can only get something when valuation is odd
+             or when valuation is even and lcf is not a p-adic square.
+             Compute smallest n such that if v(den) >= n, the leading
+             term determines the valuation. Then we must have v(den) < n. */
+          long k, n = 1;
+          for (k = degree-1; k >= 0; k--)
+          {
+            if (coeffs_mod_p[k] == 0)
+            {
+              long t = 1 + v - Z_lval(pel(c,k), p);
+              long m = CEIL(t, (degree-k));
+              if (m > n) n = m;
+            }
+          }
+          if (n == 1)
+          {
+            forb_ba[fba].p     = p;
+            forb_ba[fba].start = sieves0[pn];
+            forb_ba[fba].end   = sieves0[pn]+p;
+            forb_ba[fba].curr  = forb_ba[fba].start;
+            fba++;
+          }
+          else
+          {
+            forbidden[fdc] = upowuu(p, n);
+            fdc++;
+          }
+        }
+      }
+      else /* leading coefficient is a non-square mod p */
+      { /* denominator divisible by p is excluded */
+        forb_ba[fba].p     = p;
+        forb_ba[fba].start = sieves0[pn];
+        forb_ba[fba].end   = sieves0[pn]+p;
+        forb_ba[fba].curr  = forb_ba[fba].start;
+        fba++;
+      }
+    }
+  } /* end for pn */
+
+  /* update sp2 and sp1 if necessary */
+  if (args->sp2 > pnp)       args->sp2 = pnp;
+  if (args->sp1 > args->sp2) args->sp1 = args->sp2;
+
+  /* sort the array to get at the best primes */
+  qsort(prec, pnp, sizeof(entry), compare_entries);
+
+  /* put the sorted entries into sieve_list */
+  for (n = 0; n < args->sp2; n++) sieve_list[n] = prec[n].ssp;
+
+  /* terminate array of forbidden divisors */
+  if (args->flags & RATPOINTS_CHECK_DENOM)
+  {
+    long n;
+
+    for (n = args->num_primes+1;
+        fba + fdc < args->max_forbidden && n <= nbprime; n++)
+    {
+      ulong p = args->listprime[n];
+
+      if (p*p > (ulong) args->b_high) break;
+      if (kroiu(pel(c,degree), p) == -1)
+      {
+        forb_ba[fba].p     = p;
+        forb_ba[fba].start = sieves0[n];
+        forb_ba[fba].end   = sieves0[n]+p;
+        forb_ba[fba].curr  = forb_ba[fba].start;
+        fba++;
+      }
+    }
+    forb_ba[fba].p = 0; /* terminating zero */
+    forbidden[fdc] = 0; /* terminating zero */
+    args->max_forbidden = fba + fdc; /* note actual number */
+  }
+
+  if (fba + fdc == 0) args->flags &= ~RATPOINTS_CHECK_DENOM;
+  avma = av; return 0;
+}
+
+/**************************************************************************
+ * The sieving procedure itself                                           *
+ **************************************************************************/
+static void
+sift(long b, ratpoints_bit_array *survivors, ratpoints_args *args,
+     bit_selection which_bits, ratpoints_bit_array bits16,
+     ratpoints_sieve_entry **sieve_list, long *bp_list, int *quit,
+     int process(long, long, GEN, void*, int*), void *info)
+{
+  pari_sp av = avma;
+  sieve_spec ssp[args->sp2];
+  int do_setup = 1;
+  long k, height = args->height, nb;
+
+  if (odd(b) == 0) which_bits = num_odd; /* even denominator */
+
+  /* Note that b is new */
+  args->bc = NULL;
+  nb = 0;
+
+  for (k = 0; k < args->num_inter; k++)
+  {
+    ratpoints_interval inter = args->domain[k];
+    long low, high;
+
+    /* Determine relevant interval [low, high] of numerators. */
+    if (b*inter.low <= -height)
+      low = -height;
+    else
+    {
+      if (b*inter.low > height) break;
+      low = ceil(b*inter.low);
+    }
+    if (b*inter.up >= height)
+      high = height;
+    else
+    {
+      if (b*inter.up < -height) continue;
+      high = floor(b*inter.up);
+    }
+
+    if (do_setup)
+    { /* set up the sieve information */
+      long n;
+
+      do_setup = 0; /* only do it once for every b */
+      for (n = 0; n < args->sp2; n++)
+      {
+        ratpoints_sieve_entry *se = sieve_list[n];
+        long p = se->p;
+        long bp = bp_list[n];
+        ratpoints_bit_array *sptr;
+
+        if (which_bits != num_all) /* divide by 2 mod p */
+          bp = odd(bp) ? (bp+p) >> 1 : bp >> 1;
+        sptr = se->sieve[bp];
+
+        ssp[n].p = p;
+        ssp[n].offset = (which_bits == num_odd) ? se->offset : 0;
+
+        /* copy if already initialized, else initialize */
+        ssp[n].ptr = sptr ? sptr : (p<BITS_IN_LONG?sieve_init1(p, se, bp, args)
+                                                  :sieve_init2(p, se, bp, args));
+      }
+    }
+
+    switch(which_bits)
+    {
+      case num_all: break;
+      case num_none: break;
+      case num_odd: low >>= 1; high--; high >>= 1; break;
+      case num_even: low++; low >>= 1; high >>= 1; break;
+    }
+
+    /* now turn the bit interval into [low, high[ */
+    high++;
+
+    if (low < high)
+    {
+      long w_low, w_high, w_low0, w_high0, range = args->array_size;
+
+      /* Now the range of longwords (= bit_arrays) */
+      w_low = low >> RBA_SHIFT;
+      w_high = (high + (long)(RBA_LENGTH-1)) >> RBA_SHIFT;
+      w_low0 = w_low;
+      w_high0 = w_low0 + range;
+      for ( ; w_low0 < w_high; w_low0 = w_high0, w_high0 += range)
+      {
+        long i;
+        if (w_high0 > w_high) { w_high0 = w_high; range = w_high0 - w_low0; }
+        /* initialise the bits */
+        for (i = range; i; i--) survivors[i-1] = bits16;
+        /* boundary words */
+        if (w_low0 == w_low)
+        {
+          long sh = low - RBA_LENGTH * w_low;
+          ulong *survl = (ulong *)survivors;
+
+#ifdef HAS_SSE2
+          if (sh >= BITS_IN_LONG)
+          {
+            survl[0] = 0UL;
+            survl[1] &= (~0UL)<<(sh - BITS_IN_LONG);
+          }
+          else
+            survl[0] &= ~(0UL)<<sh;
+#else
+          survl[0] &= ~(0UL)<<sh;
+#endif
+        }
+        if (w_high0 == w_high)
+        {
+          long sh = RBA_LENGTH * w_high - high;
+          ulong *survl = (ulong *)&survivors[range-1];
+
+#ifdef HAS_SSE2
+          if (sh >= BITS_IN_LONG)
+          {
+            survl[0] &= ~(0UL)>>(sh - BITS_IN_LONG);
+            survl[1] = 0UL;
+          }
+          else
+            survl[1] &= ~(0UL)>>sh;
+#else
+          survl[0] &= ~(0UL)>>sh;
+#endif
+        }
+        nb += _ratpoints_sift0(b, w_low0, w_high0, args, which_bits,
+                         survivors, &ssp[0], quit, process, info);
+        if (*quit) return;
+      }
+    }
+  }
+  if (nb==0) avma = av;
+}
+
+/**************************************************************************
+ * Find points by looping over the denominators and sieving numerators    *
+ **************************************************************************/
+static void
+find_points_work(ratpoints_args *args,
+                 int process(long, long, GEN, void*, int*), void *info)
+{
+  int quit = 0;
+  GEN c = args->cof;
+  long degree = degpol(c);
+  long nbprime = lg(args->listprime)-1;
+  long height = args->height;
+
+  int point_at_infty = 0; /* indicates if there are points at infinity */
+  int lcfsq = Z_issquare(pel(c,degree));
+
+  forbidden_entry *forb_ba = args->forb_ba;
+  long *forbidden = args->forbidden;
+    /* The forbidden divisors, a zero-terminated array.
+       Used when degree is even and leading coefficient is not a square */
+
+    /* These are used when degree is odd and leading coeff. is not +-1 */
+
+  ratpoints_sieve_entry **sieve_list = (ratpoints_sieve_entry**)
+     stack_malloc(nbprime*sizeof(ratpoints_sieve_entry*));
+  bit_selection which_bits = num_all;
+  ulong den_bits;
+  ratpoints_bit_array num_bits[16];
+
+  args->flags &= RATPOINTS_FLAGS_INPUT_MASK;
+  args->flags |= RATPOINTS_CHECK_DENOM;
+
+  /* initialize memory management */
+  args->se_next = args->se_buffer;
+  args->ba_next = args->ba_buffer;
+  args->int_next = args->int_buffer;
+
+  /* Some sanity checks */
+  args->num_inter = 0;
+
+  if (args->num_primes > nbprime) args->num_primes = nbprime;
+  if (args->sp2 > args->num_primes) args->sp2 = args->num_primes;
+  if (args->sp1 > args->sp2)        args->sp1 = args->sp2;
+
+  if (args->b_low < 1)  args->b_low = 1;
+  if (args->b_high < 1) args->b_high = height;
+  if (args->b_high > height) args->b_high = height;
+  if (args->max_forbidden < 0)
+    args->max_forbidden = RATPOINTS_DEFAULT_MAX_FORBIDDEN;
+  if (args->max_forbidden > nbprime)
+    args->max_forbidden = nbprime;
+  if (args->array_size <= 0) args->array_size = RATPOINTS_ARRAY_SIZE;
+  {
+    long s = 2*CEIL(height, BITS_IN_LONG);
+    if (args->array_size > s) args->array_size = s;
+  }
+  /* Don't reverse if intervals are specified or limits for the denominator
+     are given */
+  if (args->num_inter > 0 || args->b_low > 1 || args->b_high < height)
+    args->flags |= RATPOINTS_NO_REVERSE;
+
+  /* Check if reversal of polynomial might be better:
+   * case 1: degree is even, but trailing coefficient is zero
+   * case 2: degree is even, leading coefficient is a square, but
+   *         trailing coefficient is not
+   * case 3: degree is odd, |leading coefficient| > 1,
+   *         trailing coefficient is zero, |coeff. of x| = 1 */
+  if (!(args->flags & RATPOINTS_NO_REVERSE))
+  {
+    if (!odd(degree))
+    {
+      if (signe(pel(c,0)) == 0)
+      { /* case 1 */
+        long n;
+        args->flags |= RATPOINTS_REVERSED;
+        for (n = 0; n < degree>>1; n++) swap(pel(c,n), pel(c,degree-n));
+        degree--;
+        setlg(c,degree+3);
+      }
+      else if (lcfsq && !Z_issquare(pel(c,0)))
+      { /* case 2 */
+        long n;
+        args->flags |= RATPOINTS_REVERSED;
+        for (n = 0; n < degree>>1; n++) swap(pel(c,n), pel(c,degree-n));
+        lcfsq = 0;
+      }
+    }
+    else
+    { /* odd degree, case 3*/
+      if (!is_pm1(pel(c,degree)) && !signe(pel(c,0)) && is_pm1(pel(c,1)))
+      {
+        long n;
+        args->flags |= RATPOINTS_REVERSED;
+        for (n = 1; n < degree>>1; n++) swap(pel(c,n),pel(c,degree+1-n));
+      }
+    }
+  }
+
+  /* Deal with the intervals */
+  if (args->num_inter == 0)
+  { /* default interval (effectively ]-oo,oo[) if none is given */
+    args->domain = (ratpoints_interval*) stack_malloc(2*degree*sizeof(ratpoints_interval));
+    args->domain[0].low = -height; args->domain[0].up = height;
+    args->num_inter = 1;
+  }
+
+  ratpoints_compute_sturm(args);
+
+  /* Point(s) at infinity? */
+  if (odd(degree) || lcfsq)
+  {
+    args->flags &= ~RATPOINTS_CHECK_DENOM;
+    point_at_infty = 1;
+  }
+
+  /* Can use only squares as denoms if degree is odd and poly is +-monic */
+  if (odd(degree))
+  {
+    GEN w1 = pel(c,degree);
+    if (is_pm1(w1))
+      args->flags |= RATPOINTS_USE_SQUARES;
+    else /* set up information on divisors of leading coefficient */
+      setup_us1(args, absi(w1));
+  }
+
+  /* deal with f mod powers of 2 */
+  which_bits = get_2adic_info(args, &den_bits, &num_bits[0]);
+  /* which_bits says whether to consider even and/or odd numerators
+   * when the denominator is odd.
+   *
+   * Bit k in den_bits is 0 if b congruent to k mod BITS_IN_LONG need
+   * not be considered as a denominator.
+   *
+   * Bit k in num_bits[b] is 0 is numerators congruent to
+   *  k (which_bits = den_all)
+   *  2k (which_bits = den_even)
+   *  2k+1 (which_bits = den_odd)
+   * need not be considered for denominators congruent to b mod 16. */
+
+  /* set up the sieve data structure */
+  if (sieving_info(args, sieve_list)) return;
+
+  /* deal with point(s) at infinity */
+  if (point_at_infty)
+  {
+    long a = 1, b = 0;
+
+    if (args->flags & RATPOINTS_REVERSED) { a = 0; b = 1; }
+    if (odd(degree))
+      (void)process(a, b, gen_0, info, &quit);
+    else
+    {
+      GEN w0 = sqrti(pel(c,degree));
+      (void)process(a, b, w0, info, &quit);
+      (void)process(a, b, negi(w0), info, &quit);
+    }
+    if (quit) return;
+  }
+  /* now do the sieving */
+  {
+    ratpoints_bit_array *survivors = (ratpoints_bit_array *)
+      stack_malloc_align((args->array_size)*sizeof(ratpoints_bit_array), RBA_ALIGN);
+    if (args->flags & (RATPOINTS_USE_SQUARES | RATPOINTS_USE_SQUARES1))
+    {
+      if (args->flags & RATPOINTS_USE_SQUARES)
+      /* need only take squares as denoms */
+      {
+        long b, bb;
+        long bp_list[args->sp2];
+        long last_b = args->b_low;
+        long n;
+        for (n = 0; n < args->sp2; n++)
+          bp_list[n] = mod(args->b_low, sieve_list[n]->p);
+
+        for (b = 1; bb = b*b, bb <= args->b_high; b++)
+          if (bb >= args->b_low)
+          {
+            ratpoints_bit_array bits = num_bits[bb & 0xf];
+            if (TEST(bits))
+            {
+              long n;
+              long d = bb - last_b;
+
+              /* fill bp_list */
+              for (n = 0; n < args->sp2; n++)
+                bp_list[n] = mod(bp_list[n] + d, sieve_list[n]->p);
+              last_b = bb;
+
+              sift(bb, survivors, args, which_bits, bits,
+                   sieve_list, &bp_list[0], &quit, process, info);
+              if (quit) break;
+            }
+          }
+      }
+      else /* args->flags & RATPOINTS_USE_SQUARES1 */
+      {
+        GEN den_info = args->den_info;
+        GEN divisors = args->divisors;
+        long ld = lg(divisors);
+        long k;
+        long b, bb;
+        long bp_list[args->sp2];
+
+        for (k = 1; k < ld; k++)
+        {
+          long d = divisors[k];
+          long last_b = d;
+          long n;
+          if (d > args->b_high) continue;
+          for (n = 0; n < args->sp2; n++)
+            bp_list[n] = mod(d, sieve_list[n]->p);
+
+          for (b = 1; bb = d*b*b, bb <= args->b_high; b++)
+          {
+            if (bb >= args->b_low)
+            {
+              int flag = 1;
+              ratpoints_bit_array bits = num_bits[bb & 0xf];
+
+              if (EXT0(bits))
+              {
+                long i, n, l = lg(gel(den_info,1));
+                long d = bb - last_b;
+
+                /* fill bp_list */
+                for (n = 0; n < args->sp2; n++)
+                  bp_list[n] = mod(bp_list[n] + d, sieve_list[n]->p);
+                last_b = bb;
+
+                for(i = 1; i < l; i++)
+                {
+                  long v = z_lval(bb, mael(den_info,1,i));
+                  if((v >= mael(den_info,3,i))
+                      && odd(v + (mael(den_info,2,i)))) { flag = 0; break; }
+                }
+                if (flag)
+                {
+                  sift(bb, survivors, args, which_bits, bits,
+                       sieve_list, &bp_list[0], &quit, process, info);
+                  if (quit) break;
+                }
+              }
+            }
+          }
+          if (quit) break;
+        }
+      }
+    }
+    else
+    {
+      if (args->flags & RATPOINTS_CHECK_DENOM)
+      {
+        long *forb;
+        long b;
+        long bp_list[args->sp2];
+        long last_b = args->b_low;
+        ulong b_bits;
+        long n;
+        for (n = 0; n < args->sp2; n++)
+          bp_list[n] = mod(args->b_low, sieve_list[n]->p);
+        {
+          forbidden_entry *fba = &forb_ba[0];
+          long b_low = args->b_low;
+          long w_low = (b_low-1) >> TWOPOTBITS_IN_LONG;
+
+          b_bits = den_bits;
+          while (fba->p)
+          {
+            fba->curr = fba->start + mod(w_low, fba->p);
+            b_bits &= *(fba->curr);
+            fba++;
+          }
+          b_bits >>= (b_low-1) & LONG_MASK;
+        }
+
+        for (b = args->b_low; b <= args->b_high; b++)
+        {
+          ratpoints_bit_array bits = num_bits[b & 0xf];
+
+          if ((b & LONG_MASK) == 0)
+          { /* next b_bits */
+            forbidden_entry *fba = &forb_ba[0];
+
+            b_bits = den_bits;
+            while (fba->p)
+            {
+              fba->curr++;
+              if (fba->curr == fba->end) fba->curr = fba->start;
+              b_bits &= *(fba->curr);
+              fba++;
+            }
+          }
+          else
+            b_bits >>= 1;
+
+          if (odd(b_bits) && EXT0(bits))
+          { /* check if denominator is excluded */
+            for (forb = &forbidden[0] ; *forb && (b % (*forb)); forb++)
+              continue;
+            if (*forb == 0 && rpjacobi(b, pel(c,degree)) == 1)
+            {
+              long n, d = b - last_b;
+
+              /* fill bp_list */
+              for (n = 0; n < args->sp2; n++)
+              {
+                long bp = bp_list[n] + d;
+                long p = sieve_list[n]->p;
+
+                while (bp >= p) bp -= p;
+                bp_list[n] = bp;
+              }
+              last_b = b;
+
+              sift(b, survivors, args, which_bits, bits,
+                   sieve_list, &bp_list[0], &quit, process, info);
+              if (quit) break;
+            }
+          }
+        }
+      } /* if (args->flags & RATPOINTS_CHECK_DENOM) */
+      else
+      {
+        long b, n;
+        long bp_list[args->sp2];
+        long last_b = args->b_low;
+        for (n = 0; n < args->sp2; n++)
+          bp_list[n] = mod(args->b_low, sieve_list[n]->p);
+        for (b = args->b_low; b <= args->b_high; b++)
+        {
+          ratpoints_bit_array bits = num_bits[b & 0xf];
+          if (EXT0(bits))
+          {
+            long n;
+            long d = b - last_b;
+
+            /* fill bp_list */
+            for (n = 0; n < args->sp2; n++)
+            {
+              long bp = bp_list[n] + d;
+              long p = sieve_list[n]->p;
+
+              while (bp >= p) bp -= p;
+              bp_list[n] = bp;
+            }
+            last_b = b;
+
+            sift(b, survivors, args, which_bits, bits,
+                 sieve_list, &bp_list[0], &quit, process, info);
+            if (quit) break;
+          }
+        }
+      }
+    }
+  }
+}
+
+static GEN
+vec_append_grow(GEN z, long i, GEN x)
+{
+  long n = lg(z)-1;
+  if (i > n)
+  {
+    n <<= 1;
+    z = vec_lengthen(z,n);
+  }
+  gel(z,i) = x;
+  return z;
+}
+
+struct points
+{
+  GEN z;
+  long i, f;
+};
+
+static int
+process(long a, long b, GEN y, void *info0, int *quit)
+{
+  struct points *p = (struct points *) info0;
+  if (b==0 && (p->f&2L)) return 0;
+  *quit = (p->f&1);
+  p->z = vec_append_grow(p->z, ++p->i, mkvec3(stoi(a),y,stoi(b)));
+  return 1;
+}
+
+static GEN
+ZX_hyperellratpoints(GEN P, GEN h, long flag)
+{
+  pari_sp av = avma;
+  ratpoints_args args;
+  struct points data;
+  ulong flags = 0;
+
+  if (!ZX_is_squarefree(P))
+    pari_err_DOMAIN("hyperellratpoints","issquarefree(pol)","=",gen_0, P);
+  if (typ(h)==t_INT && signe(h)>0)
+  {
+    long H = itos(h);
+    args.height        = H;
+    args.b_low         = 1;
+    args.b_high        = H;
+  } else if (typ(h)==t_VEC && lg(h)==3)
+  {
+    args.height        = gtos(gel(h,1));
+    if (typ(gel(h,2))==t_INT)
+    {
+      args.b_low         = 1;
+      args.b_high        = itos(gel(h,2));
+    } else if (typ(gel(h,2))==t_VEC && lg(gel(h,2))==3)
+    {
+      args.b_low         = gtos(gmael(h,2,1));
+      args.b_high        = gtos(gmael(h,2,2));
+    } else pari_err_TYPE("hyperellratpoints",h);
+  } else pari_err_TYPE("hyperellratpoints",h);
+
+  find_points_init(&args, RATPOINTS_DEFAULT_BIT_PRIMES);
+
+  args.cof           = shallowcopy(P);
+  args.num_inter     = 0;
+  args.sp1           = RATPOINTS_DEFAULT_SP1;
+  args.sp2           = RATPOINTS_DEFAULT_SP2;
+  args.array_size    = RATPOINTS_ARRAY_SIZE;
+  args.num_primes    = RATPOINTS_DEFAULT_NUM_PRIMES;
+  args.bit_primes    = RATPOINTS_DEFAULT_BIT_PRIMES;
+  args.max_forbidden = RATPOINTS_DEFAULT_MAX_FORBIDDEN;
+  args.flags         = flags;
+  data.z = cgetg(17,t_VEC);
+  data.i = 0;
+  data.f = flag;
+  find_points_work(&args, process, (void *)&data);
+
+  setlg(data.z, data.i+1);
+  return gerepilecopy(av, data.z);
+}
+
+static GEN
+ZX_homogenous_evalpow(GEN Q, GEN A, GEN B)
+{
+  pari_sp av = avma;
+  long i, d = degpol(Q);
+  GEN s = gel(Q, d+2);
+  for (i = d-1; i >= 0; i--)
+    s = addii(mulii(s, A), mulii(gel(B,d+1-i), gel(Q,i+2)));
+  return gerepileuptoint(av, s);
+}
+
+static GEN
+to_ZX(GEN a, long v) { return typ(a)==t_INT? scalarpol(a,v): a; }
+
+static int
+hyperell_check(GEN PQ, GEN *P, GEN *Q)
+{
+  *P = *Q = NULL;
+  if (typ(PQ) == t_POL)
+  {
+    if (!RgX_is_ZX(PQ)) return 0;
+    *P = PQ;
+  }
+  else
+  {
+    long v = gvar(PQ);
+    if (v == NO_VARIABLE || typ(PQ) != t_VEC || lg(PQ) != 3) return 0;
+    *P = to_ZX(gel(PQ,1), v); if (!RgX_is_ZX(*P)) return 0;
+    *Q = to_ZX(gel(PQ,2), v); if (!RgX_is_ZX(*Q)) return 0;
+    if (!signe(*Q)) *Q = NULL;
+  }
+  return 1;
+}
+
+GEN
+hyperellratpoints(GEN PQ, GEN h, long flag)
+{
+  pari_sp av = avma;
+  GEN P, Q, H, L;
+  long i, l, dy, dQ;
+
+  if (flag<0 || flag>1) pari_err_FLAG("ellratpoints");
+  if (!hyperell_check(PQ, &P, &Q)) pari_err_TYPE("hyperellratpoints",PQ);
+  if (!Q)
+  {
+    L = ZX_hyperellratpoints(P, h, flag|2L);
+    dy = (degpol(P)+1)>>1;
+    l = lg(L);
+    for (i = 1; i < l; i++)
+    {
+      GEN Li = gel(L,i), x = gel(Li,1), y = gel(Li,2), z = gel(Li,3);
+      gel(L,i) = mkvec2(gdiv(x,z), gdiv(y, powiu(z, dy)));
+    }
+    return gerepilecopy(av, L);
+  }
+  H = ZX_add(ZX_shifti(P,2), ZX_sqr(Q));
+  dy = (degpol(H)+1)>>1; dQ = degpol(Q);
+  L = ZX_hyperellratpoints(H, h, flag|2L);
+  l = lg(L);
+  for (i = 1; i < l; i++)
+  {
+    GEN Li = gel(L,i), x = gel(Li,1), y = gel(Li,2), z = gel(Li,3);
+    GEN B = gpowers(z, dQ);
+    GEN Qx = gdiv(ZX_homogenous_evalpow(Q, x, B), gel(B, dQ+1));
+    gel(L,i) = mkvec2(gdiv(x,z), gmul2n(gsub(gdiv(y,powiu(z, dy)),Qx),-1));
+  }
+  return gerepilecopy(av, L);
+}
+
+GEN
+ellratpoints(GEN E, GEN h, long flag)
+{
+  pari_sp av = avma;
+  GEN L, a1, a3;
+  long i, l;
+  checkell_Q(E);
+  if (flag<0 || flag>1) pari_err_FLAG("ellratpoints");
+  if (!RgV_is_ZV(vecslice(E,1,5))) pari_err_TYPE("ellratpoints",E);
+  a1 = ell_get_a1(E);
+  a3 = ell_get_a3(E);
+  L = ZX_hyperellratpoints(ec_bmodel(E), h, flag|2L);
+  l = lg(L);
+  for (i = 1; i < l; i++)
+  {
+    GEN P, Li = gel(L,i), x = gel(Li,1), y = gel(Li,2), z = gel(Li,3);
+    if (!signe(z))
+      P = ellinf();
+    else
+    {
+      GEN z2 = sqri(z);
+      y = subii(y, addii(mulii(a1, mulii(x,z)), mulii(a3,z2)));
+      P = mkvec2(gdiv(x,z), gdiv(y,shifti(z2,1)));
+    }
+    gel(L,i) = P;
+  }
+  return gerepilecopy(av, L);
+}
--- src/test/32/ellratpoints
+++ src/test/32/ellratpoints
@@ -0,0 +1,29 @@
+201
+173
+203
+217
+201
+209
+143
+183
+137
+127
+[[-5, 1]]
+[]
+[[-3, 0], [-3, -1], [-2, 3], [-2, -4], [-1, 3], [-1, -4], [0, 2], [0, -3], [
+1, 0], [1, -1], [2, 0], [2, -1], [3, 3], [3, -4], [4, 6], [4, -7], [8, 21], 
+[8, -22], [11, 35], [11, -36], [14, 51], [14, -52], [21, 95], [21, -96], [37
+, 224], [37, -225], [52, 374], [52, -375], [93, 896], [93, -897], [342, 6324
+], [342, -6325], [406, 8180], [406, -8181], [816, 23309], [816, -23310]]
+[[-26/9, 28/27], [-26/9, -55/27], [4/9, 35/27], [4/9, -62/27], [7/9, 17/27],
+ [7/9, -44/27], [25/9, 64/27], [25/9, -91/27], [31/9, 116/27], [31/9, -143/2
+7]]
+[[-4, 831896], [-4, -831896], [-1, 3128], [-1, -3128], [0, 15740], [0, -1574
+0], [4, 167392], [4, -167392], [5, 399050], [5, -399050], [6, 886754], [6, -
+886754]]
+470
+8
+  ***   at top-level: E=ellinit([0,0,1,-7,6]);ellratpoints(E,[10^5
+  ***                                         ^--------------------
+  *** ellratpoints: incorrect type in hyperellratpoints (t_VEC).
+Total time spent: 72
--- src/test/32/ratpoints
+++ src/test/32/ratpoints
@@ -0,0 +1,1031 @@
+1:[]
+2:[[0, 2, 1], [0, -2, 1]]
+3:[[0, 2, 1], [0, -2, 1], [-2, 346, 5], [-2, -346, 5]]
+4:[[0, 0, 1]]
+5:[]
+6:[[0, 0, 1], [-1, 0, 1]]
+7:[[-3, 6, 2], [-3, -6, 2]]
+8:[]
+9:[]
+10:[]
+11:[]
+12:[[3, 31, 2], [3, -31, 2], [2, 51, 3], [2, -51, 3]]
+13:[]
+14:[[0, 0, 1], [-1, 2, 1], [-1, -2, 1]]
+15:[]
+16:[]
+17:[]
+18:[[1, 0, 1]]
+19:[[-1, 1, 1], [-1, -1, 1], [1, 20, 2], [1, -20, 2]]
+20:[[-1, 3, 1], [-1, -3, 1]]
+21:[[0, 2, 1], [0, -2, 1]]
+22:[[-2, 7, 1], [-2, -7, 1], [0, 1, 1], [0, -1, 1]]
+23:[[-1, 2, 1], [-1, -2, 1]]
+24:[[1, 11, 2], [1, -11, 2]]
+25:[]
+26:[[1, 1, 1], [1, -1, 1]]
+27:[]
+28:[]
+29:[[-3, 38, 1], [-3, -38, 1], [-1, 2, 1], [-1, -2, 1], [1, 2, 1], [1, -2, 1
+], [7, 602, 1], [7, -602, 1], [1, 513, 6], [1, -513, 6]]
+30:[]
+31:[[0, 2, 1], [0, -2, 1]]
+32:[[1, 1, 1], [1, -1, 1]]
+33:[]
+34:[[0, 0, 1], [-1, 3, 2], [-1, -3, 2]]
+35:[[-1, 3, 1], [-1, -3, 1]]
+36:[]
+37:[]
+38:[[0, 0, 1], [-1, 3, 1], [-1, -3, 1]]
+39:[[0, 2, 1], [0, -2, 1]]
+40:[]
+41:[[0, 1, 1], [0, -1, 1]]
+42:[]
+43:[[0, 1, 1], [0, -1, 1]]
+44:[[1, 2, 1], [1, -2, 1]]
+45:[[0, 3, 1], [0, -3, 1], [1, 23, 2], [1, -23, 2]]
+46:[]
+47:[]
+48:[[0, 2, 1], [0, -2, 1]]
+49:[[-1, 4, 1], [-1, -4, 1]]
+50:[[1, 3, 1], [1, -3, 1]]
+51:[]
+52:[[1, 4, 1], [1, -4, 1]]
+53:[[-1, 2, 1], [-1, -2, 1], [-1, 4, 2], [-1, -4, 2]]
+54:[[-1, 1, 1], [-1, -1, 1]]
+55:[[0, 0, 1]]
+56:[[1, 3, 1], [1, -3, 1]]
+57:[[-1, 4, 1], [-1, -4, 1], [0, 3, 1], [0, -3, 1], [2, 107, 3], [2, -107, 3
+], [38, 1939745, 83], [38, -1939745, 83]]
+58:[[0, 1, 1], [0, -1, 1], [1, 3, 2], [1, -3, 2], [-13, 6499, 9], [-13, -649
+9, 9], [-49, 198463, 12], [-49, -198463, 12]]
+59:[[1, 0, 1]]
+60:[]
+61:[]
+62:[[0, 3, 1], [0, -3, 1]]
+63:[]
+64:[]
+65:[[0, 3, 1], [0, -3, 1]]
+66:[[0, 0, 1]]
+67:[]
+68:[]
+69:[]
+70:[[1, 3, 1], [1, -3, 1]]
+71:[[1, 0, 1]]
+72:[]
+73:[]
+74:[[0, 0, 1]]
+75:[[0, 3, 1], [0, -3, 1]]
+76:[]
+77:[[0, 3, 1], [0, -3, 1]]
+78:[[-1, 4, 1], [-1, -4, 1]]
+79:[]
+80:[]
+81:[]
+82:[]
+83:[]
+84:[[-1, 2, 1], [-1, -2, 1]]
+85:[]
+86:[]
+87:[]
+88:[]
+89:[]
+90:[]
+91:[]
+92:[[1, 1, 1], [1, -1, 1]]
+93:[]
+94:[[-1, 0, 1]]
+95:[[0, 3, 1], [0, -3, 1], [4, 23, 1], [4, -23, 1]]
+96:[[0, 1, 1], [0, -1, 1]]
+97:[[0, 2, 1], [0, -2, 1]]
+98:[[-1, 2, 1], [-1, -2, 1], [1, 2, 1], [1, -2, 1]]
+99:[]
+100:[[2, 3461, 11], [2, -3461, 11]]
+101:[[-1, 1, 1], [-1, -1, 1]]
+102:[[0, 1, 1], [0, -1, 1]]
+103:[]
+104:[]
+105:[]
+106:[]
+107:[]
+108:[[-1, 5, 1], [-1, -5, 1]]
+109:[[-1, 1, 1], [-1, -1, 1]]
+110:[[-2, 8, 1], [-2, -8, 1], [-1, 1, 1], [-1, -1, 1], [2, 16, 1], [2, -16, 
+1]]
+111:[[0, 3, 1], [0, -3, 1]]
+112:[]
+113:[[1, 144, 4], [1, -144, 4]]
+114:[]
+115:[[3, 82, 2], [3, -82, 2]]
+116:[]
+117:[[2, 4, 1], [2, -4, 1]]
+118:[[-1, 2, 1], [-1, -2, 1]]
+119:[]
+120:[]
+121:[[0, 0, 1], [1, 2, 1], [1, -2, 1], [-2, 16, 1], [-2, -16, 1]]
+122:[[0, 2, 1], [0, -2, 1], [3, 23, 1], [3, -23, 1]]
+123:[]
+124:[]
+125:[[1, 3, 1], [1, -3, 1]]
+126:[]
+127:[[-1, 1, 1], [-1, -1, 1]]
+128:[]
+129:[]
+130:[]
+131:[]
+132:[[0, 3, 1], [0, -3, 1]]
+133:[]
+134:[[0, 1, 1], [0, -1, 1]]
+135:[]
+136:[]
+137:[[0, 2, 1], [0, -2, 1]]
+138:[]
+139:[]
+140:[[-1, 3, 1], [-1, -3, 1], [2, 3, 1], [2, -3, 1], [3, 6, 2], [3, -6, 2]]
+141:[[0, 3, 1], [0, -3, 1], [1, 6, 1], [1, -6, 1], [-4, 419, 7], [-4, -419, 
+7]]
+142:[[1, 2, 1], [1, -2, 1], [7, 2996, 10], [7, -2996, 10]]
+143:[]
+144:[[1, 2, 1], [1, -2, 1]]
+145:[[0, 2, 1], [0, -2, 1]]
+146:[]
+147:[[3, 49, 2], [3, -49, 2]]
+148:[[-1, 18, 2], [-1, -18, 2]]
+149:[]
+150:[]
+151:[[0, 0, 1]]
+152:[[0, 3, 1], [0, -3, 1]]
+153:[]
+154:[[0, 2, 1], [0, -2, 1]]
+155:[[-1, 1, 1], [-1, -1, 1]]
+156:[[-2, 22, 1], [-2, -22, 1]]
+157:[]
+158:[[0, 1, 1], [0, -1, 1]]
+159:[[1, 1, 1], [1, -1, 1]]
+160:[]
+161:[]
+162:[]
+163:[]
+164:[]
+165:[]
+166:[[-1, 1, 1], [-1, -1, 1], [0, 3, 1], [0, -3, 1]]
+167:[]
+168:[[3, 1, 1], [3, -1, 1]]
+169:[]
+170:[]
+171:[[1, 4, 1], [1, -4, 1]]
+172:[[-1, 1, 1], [-1, -1, 1]]
+173:[[0, 3, 1], [0, -3, 1], [1, 0, 1]]
+174:[]
+175:[]
+176:[[1, 5, 1], [1, -5, 1]]
+177:[[10, 3648, 7], [10, -3648, 7]]
+178:[[-20, 21952, 1], [-20, -21952, 1], [0, 2, 1], [0, -2, 1]]
+179:[[-1, 2, 1], [-1, -2, 1]]
+180:[]
+181:[[1, 1, 1], [1, -1, 1], [-27, 1795891, 89], [-27, -1795891, 89]]
+182:[[-1, 5, 2], [-1, -5, 2]]
+183:[]
+184:[]
+185:[]
+186:[]
+187:[]
+188:[]
+189:[[-1, 9, 2], [-1, -9, 2]]
+190:[]
+191:[]
+192:[]
+193:[[1, 3, 1], [1, -3, 1]]
+194:[[-1, 1, 2], [-1, -1, 2], [1, 1, 1], [1, -1, 1]]
+195:[]
+196:[[-1, 1, 1], [-1, -1, 1]]
+197:[]
+198:[]
+199:[[0, 2, 1], [0, -2, 1]]
+200:[]
+201:[[7, 357, 4], [7, -357, 4]]
+202:[[0, 1, 1], [0, -1, 1]]
+203:[]
+204:[]
+205:[]
+206:[[1, 3, 1], [1, -3, 1]]
+207:[]
+208:[[0, 0, 1]]
+209:[[-2, 9, 1], [-2, -9, 1]]
+210:[[1, 0, 1]]
+211:[[-1, 7, 2], [-1, -7, 2]]
+212:[[0, 2, 1], [0, -2, 1]]
+213:[[0, 1, 1], [0, -1, 1]]
+214:[[0, 2, 1], [0, -2, 1], [2, 12, 3], [2, -12, 3]]
+215:[]
+216:[]
+217:[[-2, 3, 1], [-2, -3, 1], [-1, 0, 1], [0, 1, 1], [0, -1, 1], [-6, 133, 5
+], [-6, -133, 5]]
+218:[]
+219:[[-1, 75, 4], [-1, -75, 4]]
+220:[]
+221:[[-1, 3, 1], [-1, -3, 1]]
+222:[]
+223:[[1, 3, 1], [1, -3, 1]]
+224:[]
+225:[]
+226:[]
+227:[]
+228:[[0, 3, 1], [0, -3, 1], [2, 21, 1], [2, -21, 1]]
+229:[]
+230:[]
+231:[[0, 1, 1], [0, -1, 1], [1, 5, 2], [1, -5, 2]]
+232:[]
+233:[[0, 3, 1], [0, -3, 1]]
+234:[]
+235:[]
+236:[[0, 2, 1], [0, -2, 1]]
+237:[[-1, 1, 1], [-1, -1, 1], [1, 101, 4], [1, -101, 4]]
+238:[[-1, 2, 1], [-1, -2, 1], [1, 13, 2], [1, -13, 2]]
+239:[[-3, 58, 1], [-3, -58, 1]]
+240:[[1, 4, 1], [1, -4, 1], [3, 1408, 8], [3, -1408, 8]]
+241:[]
+242:[[0, 1, 1], [0, -1, 1]]
+243:[[1, 2, 1], [1, -2, 1]]
+244:[]
+245:[]
+246:[[-2, 9, 1], [-2, -9, 1]]
+247:[[1, 4, 1], [1, -4, 1]]
+248:[[-1, 3, 1], [-1, -3, 1]]
+249:[]
+250:[]
+251:[]
+252:[]
+253:[]
+254:[]
+255:[[-1, 1, 1], [-1, -1, 1]]
+256:[]
+257:[]
+258:[]
+259:[[0, 1, 1], [0, -1, 1]]
+260:[]
+261:[]
+262:[]
+263:[[0, 0, 1]]
+264:[[0, 1, 1], [0, -1, 1]]
+265:[]
+266:[]
+267:[[1, 1, 1], [1, -1, 1]]
+268:[]
+269:[]
+270:[]
+271:[]
+272:[]
+273:[[0, 2, 1], [0, -2, 1]]
+274:[[-4, 114, 1], [-4, -114, 1]]
+275:[[0, 3, 1], [0, -3, 1]]
+276:[]
+277:[]
+278:[]
+279:[]
+280:[]
+281:[]
+282:[[0, 2, 1], [0, -2, 1]]
+283:[[-1, 3, 1], [-1, -3, 1]]
+284:[]
+285:[]
+286:[[0, 3, 1], [0, -3, 1]]
+287:[[20, 25315, 19], [20, -25315, 19]]
+288:[]
+289:[[-1, 1, 1], [-1, -1, 1]]
+290:[]
+291:[[0, 3, 1], [0, -3, 1], [3, 291, 5], [3, -291, 5]]
+292:[]
+293:[]
+294:[]
+295:[]
+296:[]
+297:[]
+298:[]
+299:[]
+300:[]
+301:[]
+302:[[0, 3, 1], [0, -3, 1]]
+303:[]
+304:[[-2, 37, 3], [-2, -37, 3]]
+305:[[-1, 2, 1], [-1, -2, 1], [1, 2, 1], [1, -2, 1]]
+306:[]
+307:[]
+308:[]
+309:[]
+310:[[0, 3, 1], [0, -3, 1]]
+311:[]
+312:[[0, 3, 1], [0, -3, 1]]
+313:[]
+314:[]
+315:[]
+316:[]
+317:[[0, 2, 1], [0, -2, 1], [1, 0, 1]]
+318:[]
+319:[[0, 3, 1], [0, -3, 1], [-704, 40680297, 933], [-704, -40680297, 933]]
+320:[[0, 0, 1]]
+321:[[-1, 1, 1], [-1, -1, 1], [1, 5, 1], [1, -5, 1]]
+322:[]
+323:[[36, 131678, 1], [36, -131678, 1]]
+324:[[0, 2, 1], [0, -2, 1]]
+325:[]
+326:[]
+327:[]
+328:[[-1, 1, 1], [-1, -1, 1], [1, 11, 2], [1, -11, 2], [1, 503, 6], [1, -503
+, 6], [2, 13, 1], [2, -13, 1], [3, 113, 5], [3, -113, 5]]
+329:[[0, 0, 1]]
+330:[[1, 3, 1], [1, -3, 1]]
+331:[]
+332:[[1, 3, 1], [1, -3, 1]]
+333:[]
+334:[]
+335:[]
+336:[[-3, 10, 1], [-3, -10, 1], [0, 2, 1], [0, -2, 1], [-1, 46, 3], [-1, -46
+, 3]]
+337:[]
+338:[[0, 0, 1]]
+339:[]
+340:[[0, 1, 1], [0, -1, 1], [-1, 16, 2], [-1, -16, 2]]
+341:[[93, 594431, 28], [93, -594431, 28]]
+342:[[0, 1, 1], [0, -1, 1]]
+343:[[-1, 2, 1], [-1, -2, 1]]
+344:[]
+345:[[1, 0, 1]]
+346:[[0, 2, 1], [0, -2, 1]]
+347:[[0, 0, 1]]
+348:[]
+349:[[0, 2, 1], [0, -2, 1], [1, 5, 1], [1, -5, 1]]
+350:[]
+351:[[-1, 0, 1], [0, 2, 1], [0, -2, 1], [1, 0, 1]]
+352:[[0, 2, 1], [0, -2, 1], [1, 0, 1]]
+353:[]
+354:[]
+355:[]
+356:[[0, 3, 1], [0, -3, 1]]
+357:[[0, 2, 1], [0, -2, 1]]
+358:[[0, 0, 1], [-1, 2, 1], [-1, -2, 1], [1, 11, 2], [1, -11, 2]]
+359:[]
+360:[]
+361:[]
+362:[]
+363:[]
+364:[]
+365:[[1, 3, 1], [1, -3, 1]]
+366:[]
+367:[]
+368:[[0, 0, 1]]
+369:[]
+370:[[0, 3, 1], [0, -3, 1]]
+371:[[-1, 0, 1], [-2, 0, 7]]
+372:[]
+373:[]
+374:[]
+375:[[-1, 0, 1]]
+376:[]
+377:[]
+378:[]
+379:[[0, 0, 1]]
+380:[]
+381:[[0, 2, 1], [0, -2, 1]]
+382:[]
+383:[]
+384:[[-2, 92, 3], [-2, -92, 3]]
+385:[]
+386:[[1, 0, 1]]
+387:[[0, 0, 1]]
+388:[]
+389:[]
+390:[[-1, 4, 1], [-1, -4, 1]]
+391:[]
+392:[[-4, 190, 1], [-4, -190, 1]]
+393:[]
+394:[[-1, 0, 1]]
+395:[]
+396:[[-1, 3, 1], [-1, -3, 1]]
+397:[]
+398:[]
+399:[[-1, 2, 1], [-1, -2, 1]]
+400:[[1, 2, 1], [1, -2, 1]]
+401:[[-5, 159, 3], [-5, -159, 3]]
+402:[]
+403:[[-1, 3, 1], [-1, -3, 1]]
+404:[[-2, 11, 1], [-2, -11, 1]]
+405:[]
+406:[]
+407:[]
+408:[[2, 6, 1], [2, -6, 1]]
+409:[]
+410:[]
+411:[[-2, 13, 1], [-2, -13, 1]]
+412:[]
+413:[]
+414:[[-5, 95, 4], [-5, -95, 4]]
+415:[]
+416:[]
+417:[[-1, 6, 1], [-1, -6, 1]]
+418:[[2, 1, 1], [2, -1, 1]]
+419:[]
+420:[[1, 2, 2], [1, -2, 2]]
+421:[]
+422:[]
+423:[]
+424:[[1, 1, 1], [1, -1, 1]]
+425:[]
+426:[]
+427:[[-1, 2, 1], [-1, -2, 1], [0, 1, 1], [0, -1, 1]]
+428:[]
+429:[[-2, 9, 1], [-2, -9, 1], [0, 3, 1], [0, -3, 1]]
+430:[]
+431:[]
+432:[[1, 5, 1], [1, -5, 1]]
+433:[[4, 43, 5], [4, -43, 5]]
+434:[]
+435:[[1, 0, 1], [-3, 42, 4], [-3, -42, 4]]
+436:[]
+437:[]
+438:[[1, 3, 1], [1, -3, 1]]
+439:[[1, 0, 1]]
+440:[[1, 2, 1], [1, -2, 1]]
+441:[[5, 1735, 8], [5, -1735, 8], [-11, 7265, 16], [-11, -7265, 16]]
+442:[]
+443:[[0, 3, 1], [0, -3, 1]]
+444:[[-2, 9, 1], [-2, -9, 1]]
+445:[[1, 2, 1], [1, -2, 1]]
+446:[[-1, 2, 1], [-1, -2, 1]]
+447:[]
+448:[]
+449:[[0, 3, 1], [0, -3, 1]]
+450:[]
+451:[]
+452:[]
+453:[]
+454:[[-1, 2, 1], [-1, -2, 1]]
+455:[]
+456:[]
+457:[]
+458:[]
+459:[[-1, 0, 1]]
+460:[]
+461:[]
+462:[[-1, 3, 1], [-1, -3, 1]]
+463:[]
+464:[]
+465:[[0, 1, 1], [0, -1, 1]]
+466:[[2, 11, 1], [2, -11, 1]]
+467:[[0, 0, 1], [1, 10, 2], [1, -10, 2]]
+468:[[-1, 460, 6], [-1, -460, 6]]
+469:[[-2, 20, 1], [-2, -20, 1]]
+470:[[-1, 3, 1], [-1, -3, 1], [-6, 368, 5], [-6, -368, 5]]
+471:[[1, 3, 1], [1, -3, 1]]
+472:[]
+473:[]
+474:[[0, 2, 1], [0, -2, 1]]
+475:[[-1, 3, 1], [-1, -3, 1], [11, 2535, 1], [11, -2535, 1]]
+476:[]
+477:[[1, 19, 2], [1, -19, 2]]
+478:[]
+479:[[-1, 3, 1], [-1, -3, 1], [0, 2, 1], [0, -2, 1]]
+480:[]
+481:[]
+482:[]
+483:[]
+484:[[1, 2, 1], [1, -2, 1], [-1, 11, 2], [-1, -11, 2]]
+485:[]
+486:[]
+487:[[1, 3, 1], [1, -3, 1]]
+488:[]
+489:[[-2, 4, 1], [-2, -4, 1]]
+490:[[-1, 27, 2], [-1, -27, 2]]
+491:[]
+492:[]
+493:[]
+494:[[0, 0, 1]]
+495:[[0, 1, 1], [0, -1, 1]]
+496:[]
+497:[[-1, 3, 1], [-1, -3, 1], [0, 1, 1], [0, -1, 1], [3, 74, 2], [3, -74, 2]
+, [-8, 1591, 7], [-8, -1591, 7]]
+498:[]
+499:[]
+500:[]
+501:[[-1, 5, 1], [-1, -5, 1]]
+502:[[-1, 4, 1], [-1, -4, 1]]
+503:[[0, 1, 1], [0, -1, 1]]
+504:[]
+505:[[2, 76, 3], [2, -76, 3], [7, 7780, 15], [7, -7780, 15]]
+506:[[0, 1, 1], [0, -1, 1], [1, 1, 1], [1, -1, 1], [3, 10, 2], [3, -10, 2]]
+507:[[0, 1, 1], [0, -1, 1], [-4, 95, 3], [-4, -95, 3]]
+508:[[1, 1, 1], [1, -1, 1], [-1, 55, 3], [-1, -55, 3], [2, 64, 3], [2, -64, 
+3]]
+509:[[3, 57, 2], [3, -57, 2]]
+510:[[2, 14, 1], [2, -14, 1]]
+511:[[-1, 4, 1], [-1, -4, 1]]
+512:[]
+513:[[-3, 5, 1], [-3, -5, 1]]
+514:[[0, 3, 1], [0, -3, 1]]
+515:[]
+516:[[0, 3, 1], [0, -3, 1]]
+517:[[-1, 5, 1], [-1, -5, 1], [1, 1, 1], [1, -1, 1], [1, 22, 2], [1, -22, 2]
+]
+518:[]
+519:[]
+520:[[-1, 0, 1], [0, 1, 1], [0, -1, 1]]
+521:[[-7, 521, 5], [-7, -521, 5]]
+522:[]
+523:[]
+524:[]
+525:[]
+526:[[-2, 20, 1], [-2, -20, 1], [1, 2, 1], [1, -2, 1], [22, 31638, 15], [22,
+ -31638, 15]]
+527:[]
+528:[]
+529:[[1, 3, 1], [1, -3, 1]]
+530:[]
+531:[]
+532:[[-1, 1, 1], [-1, -1, 1]]
+533:[]
+534:[]
+535:[]
+536:[]
+537:[]
+538:[[1, 4, 1], [1, -4, 1], [1, 23, 2], [1, -23, 2]]
+539:[]
+540:[]
+541:[[0, 3, 1], [0, -3, 1], [-5, 257, 4], [-5, -257, 4]]
+542:[]
+543:[[-1, 3, 1], [-1, -3, 1]]
+544:[]
+545:[[-1, 0, 1], [0, 3, 1], [0, -3, 1]]
+546:[[-1, 0, 1]]
+547:[[0, 0, 1]]
+548:[[0, 3, 1], [0, -3, 1], [-5, 348, 6], [-5, -348, 6]]
+549:[]
+550:[]
+551:[]
+552:[[-1, 1, 1], [-1, -1, 1]]
+553:[]
+554:[]
+555:[]
+556:[]
+557:[[0, 3, 1], [0, -3, 1]]
+558:[]
+559:[[0, 3, 1], [0, -3, 1]]
+560:[]
+561:[[-1, 7, 2], [-1, -7, 2]]
+562:[]
+563:[]
+564:[]
+565:[]
+566:[[3, 3, 2], [3, -3, 2]]
+567:[]
+568:[]
+569:[]
+570:[]
+571:[]
+572:[]
+573:[]
+574:[[0, 3, 1], [0, -3, 1]]
+575:[]
+576:[[1, 4, 1], [1, -4, 1]]
+577:[]
+578:[]
+579:[]
+580:[[1, 0, 1]]
+581:[]
+582:[[0, 1, 1], [0, -1, 1]]
+583:[]
+584:[[-1, 4, 1], [-1, -4, 1]]
+585:[[0, 1, 1], [0, -1, 1]]
+586:[[0, 1, 1], [0, -1, 1]]
+587:[]
+588:[]
+589:[[4, 100, 1], [4, -100, 1]]
+590:[[1, 0, 1], [-313, 27928221051, 2496], [-313, -27928221051, 2496]]
+591:[]
+592:[[0, 1, 1], [0, -1, 1]]
+593:[]
+594:[]
+595:[]
+596:[]
+597:[]
+598:[[0, 2, 1], [0, -2, 1]]
+599:[]
+600:[[1, 3, 1], [1, -3, 1]]
+601:[]
+602:[[-1, 0, 1]]
+603:[]
+604:[[-1, 4, 1], [-1, -4, 1]]
+605:[[0, 3, 1], [0, -3, 1]]
+606:[[1, 1, 1], [1, -1, 1]]
+607:[]
+608:[[-3, 81, 1], [-3, -81, 1]]
+609:[[-47, 3295394, 104], [-47, -3295394, 104]]
+610:[[0, 2, 1], [0, -2, 1]]
+611:[[-1, 5, 1], [-1, -5, 1]]
+612:[[1, 1, 1], [1, -1, 1]]
+613:[[-1, 5, 1], [-1, -5, 1]]
+614:[[1, 0, 1]]
+615:[]
+616:[[-1, 53, 4], [-1, -53, 4]]
+617:[]
+618:[[-2, 16, 1], [-2, -16, 1], [-1, 0, 1], [0, 2, 1], [0, -2, 1], [5, 198, 
+1], [5, -198, 1], [1, 27, 2], [1, -27, 2], [7, 5094, 11], [7, -5094, 11]]
+619:[[0, 2, 1], [0, -2, 1], [3, 20, 4], [3, -20, 4]]
+620:[]
+621:[[0, 0, 1]]
+622:[]
+623:[]
+624:[[1, 4, 1], [1, -4, 1]]
+625:[]
+626:[]
+627:[]
+628:[[0, 1, 1], [0, -1, 1], [-57, 543145, 112], [-57, -543145, 112]]
+629:[[0, 3, 1], [0, -3, 1]]
+630:[[0, 3, 1], [0, -3, 1]]
+631:[]
+632:[[-1, 2, 1], [-1, -2, 1], [1, 0, 1]]
+633:[]
+634:[]
+635:[]
+636:[]
+637:[[0, 3, 1], [0, -3, 1]]
+638:[[-1, 3, 1], [-1, -3, 1], [0, 3, 1], [0, -3, 1]]
+639:[[-1, 3, 1], [-1, -3, 1]]
+640:[]
+641:[]
+642:[[-4, 170, 1], [-4, -170, 1], [0, 2, 1], [0, -2, 1]]
+643:[[-8, 604, 1], [-8, -604, 1], [-2, 22, 1], [-2, -22, 1], [-1, 5, 1], [-1
+, -5, 1], [0, 0, 1], [1, 1, 1], [1, -1, 1], [-1, 381, 9], [-1, -381, 9], [1,
+ 4, 2], [1, -4, 2], [-9, 3228, 8], [-9, -3228, 8], [-1, 95, 5], [-1, -95, 5]
+, [9, 1200, 10], [9, -1200, 10]]
+644:[[0, 1, 1], [0, -1, 1], [1, 3, 1], [1, -3, 1]]
+645:[]
+646:[]
+647:[[1, 1, 1], [1, -1, 1]]
+648:[[0, 1, 1], [0, -1, 1]]
+649:[[0, 1, 1], [0, -1, 1]]
+650:[]
+651:[[1, 0, 1], [-4, 123, 5], [-4, -123, 5]]
+652:[]
+653:[[0, 3, 1], [0, -3, 1]]
+654:[]
+655:[[0, 2, 1], [0, -2, 1]]
+656:[]
+657:[]
+658:[[0, 1, 1], [0, -1, 1]]
+659:[]
+660:[]
+661:[[0, 0, 1], [-1, 2, 1], [-1, -2, 1], [1, 8, 2], [1, -8, 2]]
+662:[[0, 2, 1], [0, -2, 1], [1, 2, 1], [1, -2, 1]]
+663:[[-1, 5, 1], [-1, -5, 1], [0, 3, 1], [0, -3, 1], [2, 17, 1], [2, -17, 1]
+]
+664:[[-1, 3, 1], [-1, -3, 1]]
+665:[[-3, 112, 2], [-3, -112, 2]]
+666:[]
+667:[]
+668:[]
+669:[[-1, 2, 1], [-1, -2, 1]]
+670:[[1, 2, 1], [1, -2, 1]]
+671:[]
+672:[]
+673:[[-1, 1, 1], [-1, -1, 1], [-1, 7, 2], [-1, -7, 2]]
+674:[]
+675:[]
+676:[[0, 1, 1], [0, -1, 1]]
+677:[]
+678:[[-1, 1, 1], [-1, -1, 1]]
+679:[]
+680:[]
+681:[]
+682:[[0, 0, 1]]
+683:[[0, 3, 1], [0, -3, 1]]
+684:[]
+685:[]
+686:[[1, 2, 1], [1, -2, 1]]
+687:[]
+688:[]
+689:[[0, 0, 1]]
+690:[[2, 5, 1], [2, -5, 1], [3, 50, 1], [3, -50, 1]]
+691:[[0, 3, 1], [0, -3, 1]]
+692:[[-1, 25, 3], [-1, -25, 3]]
+693:[]
+694:[]
+695:[]
+696:[]
+697:[]
+698:[]
+699:[]
+700:[[-1, 4, 1], [-1, -4, 1]]
+701:[]
+702:[[-1, 4, 1], [-1, -4, 1]]
+703:[[0, 2, 1], [0, -2, 1], [1, 1, 1], [1, -1, 1], [63, 338641, 71], [63, -3
+38641, 71], [168, 8223850, 205], [168, -8223850, 205]]
+704:[[1, 2, 1], [1, -2, 1]]
+705:[[4, 10, 3], [4, -10, 3]]
+706:[]
+707:[]
+708:[]
+709:[[0, 1, 1], [0, -1, 1]]
+710:[[-1, 2, 1], [-1, -2, 1]]
+711:[[0, 2, 1], [0, -2, 1]]
+712:[[-1, 5, 1], [-1, -5, 1], [-2, 24, 1], [-2, -24, 1]]
+713:[]
+714:[[-4, 16, 1], [-4, -16, 1], [0, 2, 1], [0, -2, 1], [1, 187, 5], [1, -187
+, 5]]
+715:[]
+716:[[-1, 3, 1], [-1, -3, 1]]
+717:[]
+718:[[-1, 0, 1], [0, 2, 1], [0, -2, 1]]
+719:[[-1, 2, 1], [-1, -2, 1]]
+720:[]
+721:[]
+722:[[-1, 3, 1], [-1, -3, 1]]
+723:[]
+724:[[1, 16, 2], [1, -16, 2], [-1, 71, 3], [-1, -71, 3], [-5, 820, 6], [-5, 
+-820, 6], [71, 117115, 51], [71, -117115, 51]]
+725:[]
+726:[[1, 0, 1]]
+727:[]
+728:[]
+729:[]
+730:[]
+731:[]
+732:[[1, 5, 2], [1, -5, 2], [5, 209, 8], [5, -209, 8]]
+733:[]
+734:[]
+735:[[-1, 1, 1], [-1, -1, 1], [1, 1, 1], [1, -1, 1], [4, 1, 1], [4, -1, 1]]
+736:[]
+737:[]
+738:[[0, 1, 1], [0, -1, 1]]
+739:[[1, 0, 1]]
+740:[]
+741:[[0, 2, 1], [0, -2, 1], [1, 0, 1], [10, 2328, 1], [10, -2328, 1], [-9, 7
+67, 4], [-9, -767, 4]]
+742:[]
+743:[[1, 4, 1], [1, -4, 1]]
+744:[[-1, 1, 1], [-1, -1, 1]]
+745:[[-2, 10, 1], [-2, -10, 1], [1, 2, 1], [1, -2, 1]]
+746:[]
+747:[[3, 67, 1], [3, -67, 1], [-1, 4, 2], [-1, -4, 2]]
+748:[]
+749:[]
+750:[[-2, 12, 1], [-2, -12, 1]]
+751:[]
+752:[]
+753:[]
+754:[[0, 1, 1], [0, -1, 1]]
+755:[]
+756:[]
+757:[]
+758:[]
+759:[[0, 2, 1], [0, -2, 1]]
+760:[[5, 2086, 13], [5, -2086, 13]]
+761:[[0, 3, 1], [0, -3, 1]]
+762:[]
+763:[]
+764:[[0, 1, 1], [0, -1, 1]]
+765:[[-2, 23, 3], [-2, -23, 3]]
+766:[]
+767:[]
+768:[[1, 4, 1], [1, -4, 1]]
+769:[[-1, 537, 6], [-1, -537, 6], [41, 636063, 51], [41, -636063, 51]]
+770:[]
+771:[]
+772:[[-1, 1, 1], [-1, -1, 1]]
+773:[]
+774:[]
+775:[]
+776:[[0, 3, 1], [0, -3, 1], [1, 25, 2], [1, -25, 2], [-1, 209, 4], [-1, -209
+, 4]]
+777:[[0, 2, 1], [0, -2, 1]]
+778:[[2, 10, 1], [2, -10, 1]]
+779:[[0, 1, 1], [0, -1, 1], [3, 241, 5], [3, -241, 5]]
+780:[]
+781:[]
+782:[]
+783:[]
+784:[[0, 0, 1], [-3, 102, 2], [-3, -102, 2], [12, 34950, 23], [12, -34950, 2
+3]]
+785:[]
+786:[[1, 2, 1], [1, -2, 1]]
+787:[[-1, 3, 1], [-1, -3, 1]]
+788:[[0, 2, 1], [0, -2, 1]]
+789:[]
+790:[]
+791:[]
+792:[[-1, 3, 1], [-1, -3, 1], [0, 1, 1], [0, -1, 1]]
+793:[]
+794:[]
+795:[[0, 0, 1], [1, 0, 2]]
+796:[[-2, 6, 1], [-2, -6, 1], [-1, 0, 1]]
+797:[[0, 3, 1], [0, -3, 1]]
+798:[[-1, 1, 1], [-1, -1, 1]]
+799:[]
+800:[[-1, 1, 1], [-1, -1, 1], [0, 3, 1], [0, -3, 1]]
+801:[[-2, 17, 1], [-2, -17, 1], [5, 368, 1], [5, -368, 1]]
+802:[]
+803:[]
+804:[[0, 1, 1], [0, -1, 1], [1, 0, 1]]
+805:[[-1, 0, 1]]
+806:[]
+807:[[0, 2, 1], [0, -2, 1], [3, 2357, 10], [3, -2357, 10]]
+808:[]
+809:[]
+810:[]
+811:[[-1, 0, 1], [0, 3, 1], [0, -3, 1]]
+812:[]
+813:[]
+814:[]
+815:[]
+816:[]
+817:[]
+818:[]
+819:[]
+820:[]
+821:[[1, 3, 1], [1, -3, 1]]
+822:[]
+823:[[0, 3, 1], [0, -3, 1]]
+824:[]
+825:[[-1, 6, 1], [-1, -6, 1], [0, 3, 1], [0, -3, 1]]
+826:[[-1, 70, 3], [-1, -70, 3]]
+827:[]
+828:[[1, 1, 1], [1, -1, 1]]
+829:[]
+830:[[0, 2, 1], [0, -2, 1]]
+831:[[-1, 0, 1], [4, 190, 1], [4, -190, 1]]
+832:[]
+833:[[-2, 1, 1], [-2, -1, 1]]
+834:[[-3, 27, 2], [-3, -27, 2], [-1, 9, 2], [-1, -9, 2]]
+835:[[-1, 4, 1], [-1, -4, 1]]
+836:[]
+837:[[0, 2, 1], [0, -2, 1]]
+838:[]
+839:[]
+840:[[1, 12, 2], [1, -12, 2], [-1, 180, 4], [-1, -180, 4], [4, 9039, 17], [4
+, -9039, 17]]
+841:[[0, 2, 1], [0, -2, 1], [-1, 25, 2], [-1, -25, 2]]
+842:[]
+843:[[-1, 4, 1], [-1, -4, 1]]
+844:[]
+845:[]
+846:[[-1, 4, 1], [-1, -4, 1]]
+847:[[3, 265, 5], [3, -265, 5]]
+848:[]
+849:[]
+850:[]
+851:[]
+852:[]
+853:[]
+854:[[-1, 3, 1], [-1, -3, 1]]
+855:[]
+856:[[1, 3, 1], [1, -3, 1]]
+857:[]
+858:[[1, 11, 2], [1, -11, 2]]
+859:[]
+860:[]
+861:[]
+862:[]
+863:[[2, 12, 1], [2, -12, 1]]
+864:[[0, 2, 1], [0, -2, 1]]
+865:[]
+866:[[0, 2, 1], [0, -2, 1]]
+867:[]
+868:[]
+869:[]
+870:[]
+871:[]
+872:[]
+873:[]
+874:[[1, 2, 1], [1, -2, 1]]
+875:[[13, 6807, 9], [13, -6807, 9]]
+876:[]
+877:[[-1, 1, 1], [-1, -1, 1]]
+878:[[0, 2, 1], [0, -2, 1], [1, 2, 1], [1, -2, 1], [4, 154, 1], [4, -154, 1]
+]
+879:[[1, 24, 2], [1, -24, 2]]
+880:[]
+881:[]
+882:[[1, 3, 1], [1, -3, 1]]
+883:[[-1, 1, 1], [-1, -1, 1]]
+884:[]
+885:[[-3, 80, 2], [-3, -80, 2]]
+886:[]
+887:[[1, 0, 1]]
+888:[]
+889:[[-1, 4, 1], [-1, -4, 1]]
+890:[]
+891:[[0, 3, 1], [0, -3, 1], [92, 14946963, 209], [92, -14946963, 209]]
+892:[]
+893:[]
+894:[]
+895:[[-2, 27, 1], [-2, -27, 1]]
+896:[]
+897:[[0, 1, 1], [0, -1, 1], [-1, 12, 2], [-1, -12, 2], [3, 56, 2], [3, -56, 
+2]]
+898:[[0, 1, 1], [0, -1, 1]]
+899:[]
+900:[[0, 3, 1], [0, -3, 1]]
+901:[[0, 2, 1], [0, -2, 1]]
+902:[[1, 0, 1]]
+903:[]
+904:[]
+905:[[1, 0, 1]]
+906:[]
+907:[]
+908:[]
+909:[[-1, 3, 1], [-1, -3, 1]]
+910:[[-1, 2, 1], [-1, -2, 1]]
+911:[]
+912:[]
+913:[[-1, 1, 1], [-1, -1, 1]]
+914:[]
+915:[]
+916:[]
+917:[[-1, 3, 1], [-1, -3, 1], [0, 1, 1], [0, -1, 1]]
+918:[[-1, 3, 1], [-1, -3, 1], [-1, 23, 2], [-1, -23, 2], [14, 4236, 13], [14
+, -4236, 13]]
+919:[]
+920:[[0, 3, 1], [0, -3, 1]]
+921:[]
+922:[]
+923:[]
+924:[]
+925:[]
+926:[[-1, 1, 1], [-1, -1, 1], [1, 1, 1], [1, -1, 1]]
+927:[]
+928:[]
+929:[[0, 2, 1], [0, -2, 1], [1, 0, 1], [7, 1990, 11], [7, -1990, 11]]
+930:[]
+931:[[0, 0, 1]]
+932:[[-1, 14, 2], [-1, -14, 2]]
+933:[]
+934:[]
+935:[[1, 1, 1], [1, -1, 1]]
+936:[[0, 0, 1]]
+937:[]
+938:[]
+939:[]
+940:[]
+941:[]
+942:[[0, 1, 1], [0, -1, 1]]
+943:[]
+944:[[0, 2, 1], [0, -2, 1], [1, 2, 2], [1, -2, 2]]
+945:[]
+946:[[1, 1, 1], [1, -1, 1]]
+947:[]
+948:[[-2, 8, 1], [-2, -8, 1], [-1, 4, 1], [-1, -4, 1], [0, 2, 1], [0, -2, 1]
+, [3, 38, 1], [3, -38, 1], [6, 200, 1], [6, -200, 1], [127, 3151792, 49], [1
+27, -3151792, 49], [1, 16, 2], [1, -16, 2], [325, 2360754224, 1058], [325, -
+2360754224, 1058]]
+949:[[3, 81, 4], [3, -81, 4]]
+950:[]
+951:[]
+952:[[1, 2, 1], [1, -2, 1]]
+953:[]
+954:[[-1, 4, 1], [-1, -4, 1]]
+955:[]
+956:[[1, 2, 2], [1, -2, 2]]
+957:[]
+958:[]
+959:[[1, 19, 2], [1, -19, 2]]
+960:[]
+961:[]
+962:[]
+963:[[-3, 39, 1], [-3, -39, 1], [-2, 9, 1], [-2, -9, 1]]
+964:[]
+965:[[2, 3, 1], [2, -3, 1]]
+966:[[1, 1, 1], [1, -1, 1]]
+967:[]
+968:[]
+969:[]
+970:[[0, 1, 1], [0, -1, 1], [-3, 287, 5], [-3, -287, 5]]
+971:[]
+972:[[1, 2, 1], [1, -2, 1]]
+973:[]
+974:[]
+975:[]
+976:[[-1, 3, 1], [-1, -3, 1], [1, 1, 1], [1, -1, 1]]
+977:[]
+978:[[0, 2, 1], [0, -2, 1]]
+979:[[0, 3, 1], [0, -3, 1], [1, 1, 1], [1, -1, 1]]
+980:[[0, 2, 1], [0, -2, 1]]
+981:[[-1, 1, 1], [-1, -1, 1]]
+982:[]
+983:[[1, 1, 1], [1, -1, 1], [-1, 23, 2], [-1, -23, 2]]
+984:[[0, 0, 1], [-1, 0, 1], [3, 0, 2]]
+985:[]
+986:[[0, 3, 1], [0, -3, 1]]
+987:[]
+988:[]
+989:[]
+990:[]
+991:[[-1, 3, 1], [-1, -3, 1]]
+992:[]
+993:[[1, 3, 1], [1, -3, 1]]
+994:[]
+995:[[1, 1, 1], [1, -1, 1]]
+996:[[0, 0, 1]]
+997:[]
+998:[[0, 2, 1], [0, -2, 1]]
+999:[]
+1000:[[0, 2, 1], [0, -2, 1]]
+Total time spent: 4527
--- src/test/in/ellratpoints
+++ src/test/in/ellratpoints
@@ -0,0 +1,43 @@
+W=
+{[
+[0, 0, 1, -79, 342],
+[1, 0, 0, -22, 219],
+[0, 0, 1, -247, 1476],
+[1, -1, 0, -415, 3481],
+[0, 0, 0, -532, 4420],
+[1, 1, 0, -2582, 48720],
+[0, 0, 1, -7077, 235516],
+[1, -1, 0, -2326, 43456],
+[1, -1, 0, -16249, 799549],
+[1, -1, 1, -63147, 6081915] ];
+}
+check(E)=
+{
+  E=ellinit(E);
+  L=ellratpoints(E,1000);
+  if(#L!=#Set(L),error([E,L]));
+  for(i=1,#L,if(!ellisoncurve(E,L[i]),error([E,L[i]])));
+  #L+1;
+}
+for(i=1,#W,print(check(W[i])))
+
+E=ellinit([-25,1]);ellratpoints(E,10,1)
+E=ellinit([-25,2]);ellratpoints(E,10,1)
+E=ellinit([0,0,1,-7,6]);ellratpoints(E,[10^5,1])
+E=ellinit([0,0,1,-7,6]);ellratpoints(E,[10^5,[5,10]])
+
+checkhyp(P,Q)=
+{
+  L=hyperellratpoints([P,Q],10000);
+  if(#L!=#Set(L),error([P,Q,L]));
+  for(i=1,#L,
+    my([x,y]=L[i]);
+    if(y^2+y*subst(Q,'x,x)!=subst(P,'x,x),error([P,Q,L[i]])));
+  #L;
+}
+P=82342800 *x^6 - 470135160 *x^5 + 52485681 *x^4 + 2396040466 *x^3 + 567207969 *x^2 - 985905640 *x + 247747600;
+hyperellratpoints(P,[10,1])
+checkhyp(P,0)
+checkhyp(-x^6+x^3+x+1,2*x^3)
+
+E=ellinit([0,0,1,-7,6]);ellratpoints(E,[10^5,[5]])
--- src/test/in/ratpoints
+++ src/test/in/ratpoints
@@ -0,0 +1,1014 @@
+affine2proj(P)=
+{
+  my([X,Y]=P);
+  my(x=numerator(X),z=denominator(X),y=Y*z^3);
+  [x,y,z];
+}
+{
+  V=[
+ [-10,-8,1,-2,0,-5,-4],
+ [4,-6,-4,-7,-8,5,-3],
+ [4,-8,2,-3,-1,0,-6],
+ [0,5,6,2,-7,4,-8],
+ [-5,10,-9,-6,7,-1,3],
+ [0,-8,-6,-9,-8,-6,-9],
+ [-6,5,7,-10,-1,1,-2],
+ [-2,9,-3,-3,-4,-2,-3],
+ [7,-7,-2,5,-1,8,7],
+ [-6,10,-6,9,-3,-7,-2],
+ [-3,-8,10,6,3,7,-7],
+ [-5,5,8,7,-2,2,-3],
+ [-8,5,9,-6,-2,-6,-1],
+ [0,6,0,-9,6,2,-3],
+ [2,1,-5,9,-9,1,-2],
+ [-5,8,-4,-5,10,-10,-5],
+ [-2,3,1,4,-5,7,3],
+ [-8,-8,4,8,7,3,-6],
+ [6,6,-5,-8,-5,-5,-2],
+ [6,-10,-6,-6,-7,-8,-8],
+ [4,10,-9,-1,6,10,3],
+ [1,2,7,3,-7,-9,-2],
+ [-7,-5,-2,0,-3,-9,2],
+ [3,1,-7,-3,4,7,3],
+ [-3,7,6,-3,-1,7,-8],
+ [-8,-3,0,-3,4,7,4],
+ [10,10,5,-7,0,3,-2],
+ [-7,-3,10,-4,-3,-7,3],
+ [7,-8,-2,7,-4,1,3],
+ [3,-3,-9,-10,5,-2,-1],
+ [4,9,-7,8,-3,1,10],
+ [-4,-9,-7,-3,9,7,8],
+ [-2,5,7,4,0,6,4],
+ [0,1,1,-1,4,4,9],
+ [8,4,-4,1,-3,-10,3],
+ [-10,-10,-3,9,7,0,-4],
+ [-7,-1,7,-9,9,-5,-7],
+ [0,5,1,-4,10,-2,-3],
+ [4,6,-6,1,-7,-8,-9],
+ [2,-7,-4,3,-2,-8,2],
+ [1,8,-8,-5,-2,9,-4],
+ [-6,9,-10,2,9,5,2],
+ [1,-4,-6,-9,-1,-3,-4],
+ [-2,4,3,1,-2,-10,10],
+ [9,5,-10,1,-10,-10,5],
+ [-1,6,-3,4,-6,-6,-2],
+ [3,4,-6,-6,-9,5,-1],
+ [4,-8,-2,6,6,3,8],
+ [2,0,3,0,10,4,5],
+ [-10,3,0,1,5,8,2],
+ [8,6,-2,5,-3,3,-9],
+ [10,-10,-3,-4,8,8,7],
+ [6,10,-4,-3,-9,-10,8],
+ [-5,-1,-3,-4,5,10,9],
+ [0,8,10,10,-8,8,-2],
+ [8,-4,-2,-10,-1,9,9],
+ [9,2,5,1,4,9,10],
+ [1,-2,6,-7,-4,-8,1],
+ [-6,8,-2,0,7,-5,-2],
+ [-4,-5,-7,8,-4,-4,-1],
+ [6,9,-9,3,-8,10,10],
+ [9,-5,5,2,2,1,1],
+ [-7,-2,1,8,2,-8,-1],
+ [5,-6,7,-7,6,-6,9],
+ [9,10,4,1,-1,7,-7],
+ [0,1,10,7,-9,-8,9],
+ [-10,9,5,-9,4,-3,6],
+ [-7,-6,5,9,-7,-10,-7],
+ [-8,-5,6,-7,-1,-4,-2],
+ [-4,7,6,0,3,1,-4],
+ [2,4,5,7,-10,-9,1],
+ [-3,9,-5,9,2,-8,2],
+ [7,2,10,-1,-4,-10,-2],
+ [0,-4,-9,-4,-4,-6,-10],
+ [9,1,-4,-4,4,7,6],
+ [-2,-1,6,8,4,-8,-10],
+ [9,-8,-5,9,2,-5,8],
+ [5,0,10,-10,2,5,-6],
+ [6,0,9,-7,-8,0,-2],
+ [6,3,-6,3,8,10,-9],
+ [-7,-2,-6,-1,10,-2,-5],
+ [-7,3,-5,9,1,-8,6],
+ [-5,7,-8,10,1,-9,1],
+ [6,-6,-7,0,6,3,-4],
+ [-2,-2,-10,-7,-5,-3,0],
+ [-6,-6,-3,-4,-9,-10,6],
+ [7,-6,0,10,-1,1,-8],
+ [10,9,-9,2,1,9,-3],
+ [-6,-4,-9,5,-10,10,8],
+ [-10,-8,-10,8,-1,-9,9],
+ [3,-3,-9,-4,-6,-8,-4],
+ [-2,8,9,-6,5,-7,-6],
+ [6,9,-9,3,7,10,-5],
+ [6,10,2,-5,0,-6,-9],
+ [9,-6,6,-9,-8,7,-1],
+ [1,0,-9,7,-10,-7,-10],
+ [4,-3,3,5,-4,-5,2],
+ [-9,-3,8,-6,2,9,3],
+ [10,-4,1,4,10,1,-5],
+ [5,10,0,-8,-9,6,3],
+ [-9,5,-4,-5,4,-4,6],
+ [1,-2,5,-5,-9,3,2],
+ [-8,-3,6,1,9,0,10],
+ [-3,-5,-3,-7,5,-4,-8],
+ [-7,8,7,1,-3,10,-4],
+ [8,0,8,9,9,-7,0],
+ [-5,-8,-3,10,-9,4,-8],
+ [3,-9,10,3,4,-9,-7],
+ [6,-8,-1,10,2,-5,-9],
+ [8,4,-6,-1,3,3,2],
+ [9,5,9,4,5,3,-9],
+ [-3,4,-5,9,-7,8,-7],
+ [3,9,-5,8,0,1,-4],
+ [-1,3,10,0,-3,2,1],
+ [7,-10,0,8,-10,3,10],
+ [7,-6,-4,2,0,-4,0],
+ [-8,4,-8,6,10,-9,2],
+ [-1,-10,-5,-2,1,9,6],
+ [-10,-4,2,9,4,7,-3],
+ [5,1,4,2,-7,4,-4],
+ [0,4,2,-10,-3,5,6],
+ [4,4,3,6,4,3,-1],
+ [-10,10,-10,8,-8,-6,-9],
+ [7,7,10,-9,-5,-1,-3],
+ [-7,10,4,-3,1,1,3],
+ [-10,5,-8,4,1,4,-4],
+ [-3,1,-1,8,10,4,8],
+ [2,10,3,-7,-1,-3,-6],
+ [-7,8,0,7,10,2,-6],
+ [-3,-4,-5,-2,1,-1,-5],
+ [-2,10,-4,8,-7,-4,-2],
+ [9,0,10,8,-4,4,-7],
+ [3,9,-5,0,1,4,7],
+ [1,10,4,9,1,-10,5],
+ [7,-6,-10,10,-8,2,-7],
+ [7,-3,7,10,-7,-6,6],
+ [4,4,2,4,10,-1,3],
+ [-4,-6,10,10,7,-5,-9],
+ [6,-8,1,-9,2,8,-7],
+ [3,-5,-6,-5,7,3,-2],
+ [9,8,-7,8,7,3,8],
+ [7,4,-4,6,4,-7,-6],
+ [-3,10,-1,1,4,-9,0],
+ [-3,5,-2,-6,9,-6,7],
+ [4,7,3,9,4,2,-7],
+ [-10,-3,1,10,5,-9,5],
+ [-5,5,-9,7,10,2,-3],
+ [7,1,0,6,-5,9,-6],
+ [-5,1,-7,-9,-1,9,-4],
+ [7,-6,-9,5,-2,-3,-6],
+ [0,6,9,-5,8,-5,10],
+ [9,10,-2,9,9,-3,5],
+ [-4,7,8,0,-6,-1,-1],
+ [4,5,7,-3,5,-9,4],
+ [-5,-7,8,5,-8,-5,-1],
+ [-10,5,-10,0,-2,2,10],
+ [7,7,5,-9,-10,2,1],
+ [1,-10,-8,1,3,4,-8],
+ [-8,4,10,1,-10,2,2],
+ [7,2,-4,-2,-10,-3,3],
+ [3,-5,-6,-2,-5,2,9],
+ [-3,-6,-9,-7,-3,3,-8],
+ [-6,-4,0,-1,9,-3,-8],
+ [-7,9,9,-7,4,0,-1],
+ [2,4,4,1,-9,10,-10],
+ [9,-1,-8,6,8,1,-2],
+ [10,-10,5,10,9,-10,-3],
+ [10,0,-10,-6,6,-7,2],
+ [10,0,-1,3,10,2,-3],
+ [-10,10,-2,7,5,-9,10],
+ [-2,-8,7,-6,6,9,10],
+ [6,5,2,-1,3,2,-4],
+ [9,5,-10,-6,-1,3,0],
+ [-2,-10,5,4,8,-6,-5],
+ [8,-4,-3,2,4,-5,-5],
+ [2,4,10,-4,8,-2,7],
+ [-4,10,8,9,1,-2,8],
+ [4,5,6,4,-8,9,8],
+ [-3,-6,4,1,6,2,-6],
+ [-8,8,-9,9,-7,-1,-3],
+ [7,-2,-9,10,1,4,-10],
+ [-1,-5,-5,-7,-9,4,-3],
+ [-3,8,2,0,-6,9,9],
+ [-4,5,0,-10,-4,-6,-8],
+ [-10,-6,-7,-6,-9,-2,4],
+ [-2,-5,-1,0,8,3,-10],
+ [-4,-3,6,9,-5,1,-8],
+ [2,0,6,8,-8,-2,-9],
+ [2,4,2,-7,-2,-4,-7],
+ [2,-7,-4,7,-2,9,7],
+ [3,-7,-3,0,1,1,-4],
+ [-4,-7,1,7,1,-4,-3],
+ [8,7,-4,-3,0,-2,3],
+ [-3,-7,1,8,3,-2,1],
+ [3,0,-7,2,-6,5,-7],
+ [-4,-2,-6,-3,-6,-2,10],
+ [5,10,1,10,10,-6,-4],
+ [-4,6,5,-10,2,6,-3],
+ [4,-5,2,3,-10,3,8],
+ [-5,4,-10,4,6,0,-2],
+ [-7,-3,-9,7,-3,2,1],
+ [1,-9,-9,2,8,8,7],
+ [-7,1,-7,-7,1,-6,-10],
+ [-9,1,-2,-5,3,5,-3],
+ [-3,-9,-3,9,6,-9,4],
+ [-2,9,-9,5,5,-7,8],
+ [5,-10,8,-8,-3,-8,-10],
+ [0,1,-7,7,6,3,-8],
+ [-1,-9,4,-8,-5,-6,-2],
+ [-6,5,1,8,-8,3,-3],
+ [6,8,-4,1,3,9,-1],
+ [4,-6,7,10,-9,-7,-8],
+ [1,7,-3,3,0,-10,-1],
+ [4,-6,-2,3,-1,7,-6],
+ [5,-8,-2,-4,-6,10,-5],
+ [10,-3,2,-1,-6,-8,-1],
+ [1,2,-5,-10,-7,-4,-1],
+ [7,-8,4,0,5,4,1],
+ [3,6,-2,2,6,-6,1],
+ [8,-4,6,1,-9,8,-8],
+ [-4,-6,8,5,-5,-4,5],
+ [-7,-3,8,-6,1,-6,-4],
+ [8,-4,-10,-3,10,4,4],
+ [-4,4,3,-7,1,-3,-8],
+ [6,-5,-4,-5,-7,-5,9],
+ [-10,6,-9,-8,2,2,2],
+ [-7,1,10,-9,10,-2,-4],
+ [9,0,-8,8,-3,10,2],
+ [3,-4,2,8,10,2,-1],
+ [2,8,-8,6,-6,7,8],
+ [1,5,-8,-4,-8,-7,7],
+ [-3,4,-5,-6,-10,-1,7],
+ [9,-6,-3,3,-6,-1,-1],
+ [-6,4,-10,6,-10,-4,-1],
+ [-4,-7,-4,4,-2,-8,7],
+ [4,8,-2,-2,-3,-7,-5],
+ [5,-10,0,2,-9,-4,-7],
+ [8,-7,-5,-6,0,7,-5],
+ [7,-9,4,10,-7,1,6],
+ [10,-8,1,8,0,-3,8],
+ [3,3,2,-2,-6,-2,-10],
+ [1,0,-8,5,6,-1,-8],
+ [-5,6,0,10,7,-4,-10],
+ [10,-5,2,5,6,8,7],
+ [7,2,7,-6,-8,-2,2],
+ [3,7,3,8,-3,-4,1],
+ [8,-2,9,0,-5,9,-3],
+ [6,0,5,-6,-9,-3,-2],
+ [-5,-8,-10,2,1,-2,-1],
+ [3,1,-9,5,-10,8,-8],
+ [6,-7,-4,-8,-2,-1,-4],
+ [-1,0,0,5,-2,-9,2],
+ [-7,10,-2,-2,-7,-7,-2],
+ [6,4,-3,-3,-1,-9,5],
+ [10,6,5,4,-1,-1,-4],
+ [-8,-4,-2,7,3,9,1],
+ [5,0,3,-6,0,10,-2],
+ [5,-10,7,7,5,2,5],
+ [1,3,4,-10,3,2,-8],
+ [-2,-1,4,-10,-10,-3,-7],
+ [-9,1,6,1,4,0,9],
+ [6,9,10,-4,2,6,10],
+ [0,2,-3,-1,-4,6,-1],
+ [1,-4,-7,-2,3,1,6],
+ [-1,-4,0,-2,-3,4,8],
+ [-8,4,-9,-4,-5,5,-6],
+ [-1,4,6,-9,6,-5,0],
+ [-6,-1,-3,5,-6,4,-7],
+ [6,-10,10,2,3,5,1],
+ [-10,-8,-1,7,-8,-2,5],
+ [3,-2,3,8,4,7,10],
+ [7,-4,8,5,-1,3,5],
+ [4,10,7,0,5,-3,5],
+ [-8,-3,-8,7,-3,-10,1],
+ [9,7,2,3,-7,-3,-3],
+ [-9,-7,-2,9,-2,-4,4],
+ [-4,-7,3,9,5,4,-5],
+ [3,3,4,10,-6,2,5],
+ [-9,2,9,0,-9,2,3],
+ [-6,-7,-1,-1,-7,-10,1],
+ [-6,8,9,-10,-3,-4,8],
+ [4,9,7,-4,4,1,10],
+ [6,6,4,-5,-2,-8,-6],
+ [-8,5,10,-2,1,1,3],
+ [-3,-2,2,-3,0,8,-4],
+ [9,2,6,8,9,8,-1],
+ [5,9,5,-10,10,3,-8],
+ [-4,1,-9,-4,-2,0,-9],
+ [10,7,4,3,-4,-7,-6],
+ [-5,9,-6,-9,-6,-6,4],
+ [9,0,-6,3,-9,-8,-6],
+ [-4,0,4,-5,8,-10,9],
+ [3,-7,5,5,9,-6,-10],
+ [-4,-5,2,9,5,-3,-10],
+ [5,6,0,2,-6,10,0],
+ [-7,-8,-1,10,-10,-3,-8],
+ [-9,6,2,-2,3,7,-9],
+ [-7,-3,-10,-6,-6,9,1],
+ [3,-10,9,2,10,7,9],
+ [-5,9,-9,9,-7,1,-5],
+ [6,-2,3,-3,3,-8,-8],
+ [9,-7,1,-9,-9,-6,-5],
+ [-3,-2,-5,1,0,0,-8],
+ [3,3,-6,-7,3,-8,-2],
+ [10,-9,-8,1,-4,8,6],
+ [-8,3,5,3,-4,-8,0],
+ [-3,3,-9,-7,5,-5,-8],
+ [8,10,-1,9,7,10,-3],
+ [-2,1,10,9,-2,-6,1],
+ [9,4,-6,9,-7,2,4],
+ [7,-8,-6,2,-6,2,-2],
+ [9,-9,0,8,4,0,-9],
+ [-4,-6,-3,6,9,9,9],
+ [-3,-6,-2,-7,6,0,10],
+ [5,-4,-7,0,-3,8,-8],
+ [7,1,5,-2,-7,-7,-7],
+ [4,-3,5,-6,7,-1,-6],
+ [6,-2,-4,4,-9,7,8],
+ [9,4,-3,8,1,7,3],
+ [0,-2,-8,-9,-2,4,10],
+ [-3,6,9,8,-3,-2,10],
+ [-8,-10,7,-8,6,6,-5],
+ [-8,-5,0,9,-9,-1,8],
+ [4,-1,10,1,5,-8,1],
+ [-4,5,4,0,-4,-9,-1],
+ [-9,8,-6,-8,2,-7,-6],
+ [-7,-8,8,9,5,6,9],
+ [7,-9,-3,2,-5,6,1],
+ [0,1,-1,0,-1,-5,0],
+ [3,-9,2,0,5,6,2],
+ [8,-7,-3,3,-8,-5,-4],
+ [10,8,-9,-7,3,2,2],
+ [10,5,3,-2,-6,9,-5],
+ [-10,7,5,-1,-6,3,-6],
+ [-10,-9,-1,-4,-10,-1,-1],
+ [4,7,10,-4,-1,3,1],
+ [-8,-7,0,-6,0,-6,4],
+ [0,2,-4,10,-7,2,9],
+ [8,-8,-1,9,5,-8,7],
+ [1,-10,-8,1,0,-5,-2],
+ [10,9,-10,6,8,9,-3],
+ [1,-3,-8,4,-7,-10,0],
+ [6,-7,-1,4,-1,-2,-5],
+ [8,6,-4,2,-2,-7,-6],
+ [10,-2,7,-10,-6,9,-8],
+ [4,2,9,6,-2,1,4],
+ [0,9,-4,-5,3,-4,-6],
+ [2,4,-4,-10,-9,-1,-1],
+ [4,6,10,2,2,-5,6],
+ [-9,3,-1,-9,-7,3,-10],
+ [4,9,-6,-5,3,-4,-1],
+ [4,9,-10,-7,-8,9,3],
+ [-2,8,0,8,5,5,-7],
+ [-7,9,1,-3,2,9,-8],
+ [-4,-1,-10,6,8,-6,-2],
+ [9,7,4,0,2,8,-4],
+ [4,9,-7,10,10,-10,-4],
+ [0,6,0,-8,-1,-2,1],
+ [-4,3,-5,1,-2,-9,-2],
+ [-7,-5,-8,-5,7,-5,8],
+ [8,1,-10,-8,-9,-10,6],
+ [-8,-3,-2,4,-5,2,-2],
+ [-5,-1,9,0,-6,-8,-5],
+ [-2,-7,-1,-3,7,-9,7],
+ [2,4,2,0,6,2,-7],
+ [7,-3,-1,7,-5,4,-2],
+ [-9,-8,-1,-2,10,-7,-4],
+ [0,-10,10,7,-8,-10,3],
+ [2,8,2,-4,4,0,-7],
+ [9,6,7,-8,10,-8,-10],
+ [-2,-9,-7,0,-2,-9,-7],
+ [-3,8,1,3,7,1,-2],
+ [-9,-2,-2,1,4,0,4],
+ [-9,4,-5,2,-9,7,2],
+ [-2,-9,5,2,-1,4,-5],
+ [2,2,-7,-10,5,-7,6],
+ [5,-7,2,-4,0,-1,-1],
+ [8,-1,2,10,10,-5,4],
+ [0,9,-2,-4,2,-5,-1],
+ [-9,4,-6,1,-7,-7,-3],
+ [4,-1,5,-4,3,-7,2],
+ [-3,9,-5,6,5,2,0],
+ [6,7,-4,9,-9,9,-4],
+ [8,-6,-5,-7,-3,2,7],
+ [-6,-8,1,2,8,5,-7],
+ [5,9,-7,2,-6,-8,5],
+ [0,-2,-7,-5,-9,-7,-7],
+ [-6,5,10,6,0,6,0],
+ [10,4,-6,0,-5,-8,4],
+ [-1,-3,-9,-7,1,-9,6],
+ [-3,-8,-4,-1,1,-7,5],
+ [-4,-6,7,10,7,6,10],
+ [8,7,-1,-8,-10,1,-9],
+ [3,7,1,8,1,-4,6],
+ [-7,1,-3,9,0,10,-7],
+ [-7,-6,8,-7,-6,2,3],
+ [3,0,4,-3,-5,-8,-2],
+ [6,1,2,-3,1,7,9],
+ [3,-6,5,9,-1,5,5],
+ [6,4,-3,6,0,-5,-4],
+ [-1,-4,1,-9,3,3,0],
+ [-2,-7,-8,-2,0,-2,-5],
+ [8,0,-3,1,-4,-3,6],
+ [-1,-7,-9,10,-4,-5,2],
+ [6,-4,-2,-1,-4,4,-2],
+ [3,-2,-7,-1,-9,2,-7],
+ [-9,-8,-8,-10,-7,6,-10],
+ [10,9,-4,1,5,-4,1],
+ [5,10,-9,-3,-3,-6,-10],
+ [10,6,-10,-7,3,8,-2],
+ [-7,10,5,6,-8,-1,5],
+ [-8,-8,-10,-6,6,-8,-8],
+ [-3,-8,9,-1,3,1,4],
+ [-10,-8,-9,8,-1,-10,1],
+ [3,-8,5,1,3,-10,-1],
+ [-8,5,1,9,-9,0,-3],
+ [8,-9,10,-10,-1,-10,-10],
+ [-1,5,-10,-4,-8,-10,8],
+ [-6,9,-8,10,-10,-3,4],
+ [3,-10,2,9,4,1,10],
+ [8,6,8,4,1,4,-9],
+ [3,-10,8,-5,-2,-10,8],
+ [10,-5,-3,7,-1,6,-2],
+ [6,-10,-4,-9,0,8,10],
+ [2,10,-6,6,10,10,-1],
+ [7,-8,4,-10,3,9,3],
+ [1,2,-3,-10,3,10,5],
+ [-5,-8,-6,9,-1,9,-6],
+ [9,-10,9,-6,8,3,-1],
+ [-7,1,-8,0,-10,-4,4],
+ [-10,1,0,3,4,3,10],
+ [10,-3,10,2,7,-4,3],
+ [-3,9,-6,-3,7,0,-6],
+ [-9,-3,0,1,4,1,-4],
+ [6,2,-4,6,0,-6,-4],
+ [7,-5,5,4,-9,6,-6],
+ [-10,-10,4,-2,7,-2,-7],
+ [8,-8,2,-6,10,2,1],
+ [5,-10,-2,1,-10,8,8],
+ [8,6,3,-9,-5,8,-7],
+ [5,2,8,8,0,6,-7],
+ [-4,3,-10,4,-3,4,-1],
+ [9,-3,-8,0,-8,-9,-6],
+ [-5,9,-6,-8,2,3,2],
+ [-1,10,-2,4,10,-9,-8],
+ [-5,-10,8,-7,-5,8,-3],
+ [-10,8,-4,8,4,-2,-7],
+ [-8,4,-1,-10,4,3,-4],
+ [9,-1,1,-9,9,4,-1],
+ [-7,-7,9,0,6,-3,-7],
+ [-2,-4,-8,3,9,9,-4],
+ [3,9,9,5,-9,8,7],
+ [-3,-4,-9,1,-2,4,9],
+ [-5,2,3,-6,-10,-9,3],
+ [-3,10,5,7,2,1,0],
+ [3,-10,9,-1,2,-2,-4],
+ [3,3,-6,-1,-6,3,-6],
+ [8,6,-9,7,-3,-4,2],
+ [-5,-1,4,-1,5,-1,-7],
+ [-8,8,0,-7,-8,-10,1],
+ [-1,0,10,2,8,-5,8],
+ [-2,-2,-8,-4,9,-8,-4],
+ [-10,-8,-1,8,2,0,2],
+ [-10,7,-4,3,-4,-6,1],
+ [1,2,-7,-3,-9,-7,2],
+ [-3,10,2,-6,-5,9,-1],
+ [0,4,-9,9,9,7,-6],
+ [5,2,-5,-3,-10,-9,10],
+ [-10,9,3,2,5,-5,3],
+ [-2,5,10,-10,5,3,-6],
+ [-9,1,0,5,8,10,-6],
+ [2,-6,-5,6,2,9,3],
+ [-4,9,-8,-2,0,-10,0],
+ [4,6,-9,-3,-2,-4,2],
+ [-8,-7,2,10,9,-5,4],
+ [-5,4,-5,6,6,2,-1],
+ [3,10,-3,-6,-10,-7,-1],
+ [8,6,-1,-1,9,1,-5],
+ [4,-8,-2,-9,-6,3,-1],
+ [8,-8,-7,-6,-10,3,-6],
+ [-1,6,7,-2,6,-10,1],
+ [5,1,6,-6,-5,-10,-8],
+ [-3,6,9,-3,-4,-9,-3],
+ [3,8,8,-10,-7,-1,3],
+ [10,0,-7,-4,7,-8,9],
+ [2,-9,2,-9,-8,10,8],
+ [-9,8,-3,5,1,-2,9],
+ [-7,4,-7,10,-2,5,-1],
+ [10,-1,-5,3,1,-1,0],
+ [10,3,5,-10,10,7,-1],
+ [2,-5,0,3,4,-7,-3],
+ [-5,-6,-9,7,-1,-9,9],
+ [8,6,-8,-10,3,1,8],
+ [0,1,-1,5,4,0,10],
+ [1,3,-1,-9,9,-6,-3],
+ [-7,0,1,8,7,8,4],
+ [1,8,0,-1,9,-2,4],
+ [-1,-8,-8,-3,-10,3,-10],
+ [7,0,-8,3,0,-6,-4],
+ [5,-3,-10,-10,8,-8,8],
+ [-4,-2,7,-7,0,-4,9],
+ [10,2,5,-2,-4,-7,-2],
+ [1,6,1,1,-1,-6,6],
+ [-1,-7,4,-1,-1,-6,-6],
+ [2,7,-1,1,-2,7,10],
+ [1,-6,2,-3,9,2,-4],
+ [1,10,6,-8,6,7,1],
+ [6,5,-3,-7,-6,8,-2],
+ [6,-1,-10,-4,5,3,3],
+ [10,1,4,-1,7,10,-4],
+ [7,2,-6,-9,6,-8,-6],
+ [-5,8,0,-3,7,5,-10],
+ [-2,0,0,8,-9,-4,0],
+ [9,-8,-6,1,6,10,1],
+ [-6,-7,3,-9,-9,6,-4],
+ [9,-2,-2,6,4,9,10],
+ [8,1,0,-4,-5,-9,10],
+ [-8,-2,6,9,0,6,6],
+ [10,10,4,1,9,-6,4],
+ [1,9,-1,2,10,9,10],
+ [2,-9,-6,-10,0,1,-1],
+ [-2,-1,-3,4,-5,6,-6],
+ [-3,7,9,5,8,9,-2],
+ [-5,5,8,10,8,7,9],
+ [-2,5,-10,3,4,10,-2],
+ [-10,-1,-6,10,6,-1,6],
+ [7,4,-1,-2,9,2,10],
+ [5,-8,-7,4,10,-8,-3],
+ [2,-2,6,4,5,-1,-5],
+ [8,2,-4,6,6,1,2],
+ [-9,2,-3,7,6,2,3],
+ [2,4,1,-2,-8,-10,-2],
+ [-6,-1,5,9,10,5,-1],
+ [-10,9,10,2,-7,-3,-8],
+ [2,-2,2,9,-2,0,-1],
+ [-6,-1,6,7,-8,0,-8],
+ [-4,-1,0,7,1,-4,-2],
+ [10,-6,2,2,10,-5,3],
+ [8,-10,-4,3,2,4,10],
+ [-9,3,1,-5,-10,-6,-5],
+ [9,-3,-2,-2,-8,-1,5],
+ [-5,3,-5,2,3,10,-2],
+ [5,-10,0,-7,1,8,-6],
+ [5,-6,-3,3,-7,0,7],
+ [9,4,1,1,-7,-8,-6],
+ [2,5,-5,-10,1,0,-3],
+ [0,8,9,-5,2,8,10],
+ [9,10,10,2,2,7,-6],
+ [-4,-2,-6,4,10,4,5],
+ [-8,-5,-8,4,8,-10,-6],
+ [-3,-4,7,-1,-10,7,6],
+ [-7,-5,-2,0,-6,-4,7],
+ [-7,5,2,-7,-5,1,-3],
+ [-6,0,8,-4,10,-4,-5],
+ [-10,-5,-8,-1,5,4,6],
+ [-7,-7,2,-5,8,0,2],
+ [9,-7,-7,2,-4,-1,4],
+ [10,-2,-7,-9,7,-5,0],
+ [9,-9,9,7,0,2,5],
+ [2,6,3,10,3,1,8],
+ [7,8,-8,0,-6,-4,1],
+ [-3,7,-10,4,5,4,1],
+ [-9,0,-7,-10,-8,8,9],
+ [3,-8,9,-5,2,-4,1],
+ [-10,-1,5,4,6,-5,4],
+ [-6,-2,-4,8,-1,-8,5],
+ [-7,-7,-7,-4,-7,-4,-8],
+ [-3,-10,0,7,-9,1,-4],
+ [-7,9,-8,3,8,8,-2],
+ [-2,-10,-5,7,-10,-10,-9],
+ [7,-2,9,1,3,4,6],
+ [10,-9,-10,-4,1,6,0],
+ [-4,10,1,9,-7,-4,-8],
+ [9,2,4,-8,-5,-5,1],
+ [-2,-3,-9,1,-6,-4,3],
+ [3,-3,-3,6,2,4,7],
+ [5,-3,-9,2,10,8,-6],
+ [-9,-4,2,0,-10,-8,-9],
+ [3,9,1,3,-8,-5,3],
+ [-2,-6,10,1,9,-4,-8],
+ [-10,3,-10,2,3,8,1],
+ [1,-3,-1,-4,0,-5,2],
+ [-7,5,10,8,2,-7,-3],
+ [-5,-6,9,-3,-10,-8,5],
+ [1,9,3,-4,-5,-2,-6],
+ [1,9,-1,-1,-10,-3,-3],
+ [10,4,3,7,-1,-2,-9],
+ [-4,-7,1,-3,5,6,8],
+ [8,10,-2,4,-6,-5,4],
+ [2,-10,-3,-9,5,6,9],
+ [-3,8,2,-1,-9,0,-2],
+ [1,-6,-8,-1,2,1,-2],
+ [-3,3,5,-4,10,-9,-10],
+ [-8,0,-2,5,-10,0,1],
+ [6,-5,-1,-6,0,0,-10],
+ [-10,0,3,1,0,7,6],
+ [7,-6,1,-6,9,3,2],
+ [4,-8,3,-3,-4,-10,9],
+ [-2,10,-8,-1,1,8,-5],
+ [-3,7,8,-1,-2,3,-3],
+ [-7,-6,-5,4,6,-4,3],
+ [2,9,-1,-1,3,-4,0],
+ [-5,-1,-10,0,-6,-3,-7],
+ [6,1,7,3,1,-9,-3],
+ [9,-7,0,-8,6,10,-3],
+ [-8,0,-3,-2,3,8,3],
+ [-9,-4,8,-2,0,7,-10],
+ [-6,7,-6,0,4,-8,6],
+ [6,-7,0,4,-2,5,-4],
+ [4,-5,5,-9,-7,8,3],
+ [8,-2,-9,-2,9,-5,8],
+ [5,1,10,10,-10,-6,-9],
+ [10,3,9,-4,1,-4,0],
+ [-4,-2,5,2,-7,9,-3],
+ [2,-5,3,8,8,2,9],
+ [-1,-4,10,-4,0,4,9],
+ [5,7,-3,-10,-3,-8,7],
+ [4,10,6,7,1,-3,3],
+ [4,-3,-8,9,-9,5,4],
+ [-8,1,-7,-9,4,-2,-9],
+ [0,5,-9,-8,0,1,-8],
+ [-5,-5,3,6,2,-10,-3],
+ [-4,-1,2,4,7,-4,-7],
+ [8,1,-5,9,1,1,1],
+ [-2,7,9,5,-1,-4,7],
+ [-8,6,8,-9,-7,-5,2],
+ [-10,9,2,-8,-10,5,-5],
+ [1,7,6,-2,9,-8,1],
+ [9,8,8,9,-1,-6,-3],
+ [9,-3,-1,6,7,-3,7],
+ [-8,0,-2,10,-2,-2,-5],
+ [-10,0,9,5,10,-7,-7],
+ [-8,-6,-6,-1,-9,-4,-7],
+ [-5,8,-6,-1,10,-3,2],
+ [3,-8,-6,10,10,-3,0],
+ [-8,6,4,0,-1,-8,-1],
+ [9,2,2,-5,8,-8,-1],
+ [9,-9,-6,8,-10,-10,5],
+ [-6,-3,-7,5,8,-10,6],
+ [-3,-8,-2,-9,1,-6,7],
+ [-3,-7,-1,2,1,7,0],
+ [4,4,-1,8,-9,-7,6],
+ [0,-2,4,0,9,-10,0],
+ [1,-1,5,-2,2,10,-6],
+ [2,5,2,5,-6,5,8],
+ [-6,-1,-3,6,0,2,10],
+ [-1,0,5,9,-4,-5,-3],
+ [1,-8,-4,-3,5,3,-4],
+ [1,9,-7,-9,-10,4,3],
+ [-5,6,-7,-10,7,-5,1],
+ [5,1,-4,-8,-3,10,-1],
+ [-10,8,4,-9,-2,3,-8],
+ [9,6,9,5,8,-6,1],
+ [5,-5,-7,-8,-10,0,-3],
+ [4,-9,-10,2,10,-6,-7],
+ [7,-3,2,-6,0,6,-4],
+ [5,5,1,2,10,7,-7],
+ [1,8,7,8,-3,-1,-9],
+ [-7,9,7,1,-10,7,8],
+ [-3,2,-6,3,3,-1,-7],
+ [0,-2,6,6,3,-9,-10],
+ [4,-5,7,-8,-6,9,3],
+ [9,-10,3,8,6,-2,3],
+ [-7,-6,10,1,-6,-6,1],
+ [10,2,-3,-10,2,-5,10],
+ [10,9,-7,-10,2,8,-6],
+ [-1,6,-5,1,-2,0,-5],
+ [-5,-4,3,3,0,-8,-2],
+ [-4,5,7,6,2,-7,3],
+ [-10,2,4,-4,6,6,0],
+ [2,6,0,3,1,-6,6],
+ [-3,1,-7,2,-9,4,1],
+ [5,7,-3,4,6,-5,-1],
+ [-9,4,9,-10,-7,10,6],
+ [-1,1,-9,-6,-4,-3,-5],
+ [1,-6,2,8,-3,-3,-8],
+ [2,-8,-8,-10,4,-9,-1],
+ [-10,-2,2,-5,5,2,-1],
+ [8,-1,-7,-1,-9,-4,7],
+ [-10,6,-7,-10,6,7,7],
+ [-1,-6,9,-10,-2,-4,-8],
+ [0,5,-9,-10,-7,-5,-7],
+ [9,2,1,-1,0,-4,-5],
+ [3,1,-4,-8,-3,-9,7],
+ [-10,-9,-1,7,-2,-1,-7],
+ [-10,7,2,2,8,2,-7],
+ [2,2,10,-5,-10,6,8],
+ [5,-9,-8,1,-10,-1,6],
+ [0,-5,-3,10,8,-9,-6],
+ [7,-3,2,-10,-8,-7,7],
+ [9,5,6,1,0,0,10],
+ [3,3,-8,7,2,6,4],
+ [2,10,-3,7,4,-4,-8],
+ [-4,6,7,4,0,8,7],
+ [-5,5,8,-7,-9,-2,-1],
+ [-3,-6,-8,-6,4,2,-1],
+ [5,-7,-5,5,8,-3,-10],
+ [6,-1,0,-8,10,7,10],
+ [6,10,3,-7,-8,-6,-2],
+ [-9,-10,4,6,3,-8,6],
+ [-1,8,-1,-2,-10,-6,-1],
+ [-2,2,6,-7,-5,-4,8],
+ [4,-8,5,2,1,-8,5],
+ [8,1,-8,-6,-4,6,7],
+ [8,1,9,3,8,-7,-5],
+ [-8,10,-6,-7,-6,-2,-5],
+ [-5,0,6,-3,5,10,-6],
+ [10,-7,3,-9,-7,-6,-5],
+ [1,-7,-1,4,3,-10,-4],
+ [7,5,2,7,9,-2,-4],
+ [4,-4,4,9,-5,-10,-9],
+ [2,-1,5,1,9,-5,4],
+ [-7,4,-5,-8,-10,6,6],
+ [4,-7,-10,6,-5,-6,-1],
+ [-4,-6,7,4,-8,8,-5],
+ [7,-2,-4,-4,-5,-7,-2],
+ [10,3,-6,3,10,6,-8],
+ [4,1,7,-4,-1,10,-3],
+ [8,0,6,-5,1,7,-9],
+ [-5,-2,8,7,4,1,9],
+ [-5,0,7,9,-3,5,0],
+ [-5,-4,1,-4,-3,-1,7],
+ [-7,5,8,3,-7,-1,-9],
+ [5,-4,4,-2,5,-1,-2],
+ [-7,-3,7,7,1,6,8],
+ [-1,-1,-10,7,-6,7,4],
+ [8,6,8,9,-8,-9,5],
+ [-8,0,8,0,-10,10,7],
+ [-10,2,-9,4,6,-9,-5],
+ [-3,-3,9,8,9,-3,-5],
+ [6,7,10,2,9,6,10],
+ [-1,5,-5,2,2,-4,-7],
+ [5,2,-3,-6,-4,1,7],
+ [5,6,10,-3,8,-5,2],
+ [-7,-6,2,-4,9,10,-3],
+ [-9,3,3,-3,-6,7,-3],
+ [8,8,-7,9,5,-2,-8],
+ [1,-3,-5,-2,-4,0,6],
+ [10,3,6,1,-10,-3,-7],
+ [2,-1,-10,7,-3,9,4],
+ [4,-2,-4,0,-8,5,5],
+ [-3,0,-3,-3,5,-1,3],
+ [5,6,3,7,10,-9,-6],
+ [-5,7,6,7,4,-3,7],
+ [-8,8,-9,4,6,1,2],
+ [-8,-7,-10,2,-5,-2,4],
+ [-2,0,10,4,2,-1,6],
+ [-10,-5,-4,-2,0,10,8],
+ [7,0,7,5,6,-2,-2],
+ [8,-2,5,-10,10,-10,-7],
+ [10,-2,-5,0,-10,3,9],
+ [2,3,10,2,-10,7,-8],
+ [8,-7,5,-9,1,1,-7],
+ [1,1,8,3,9,9,8],
+ [10,10,-8,7,-6,10,10],
+ [-2,-10,-5,-8,-8,6,-6],
+ [-5,-6,3,-8,-2,-5,-5],
+ [-1,-1,3,5,7,2,2],
+ [4,-5,-5,-6,3,5,7],
+ [-1,7,-2,-6,-4,-4,-10],
+ [9,10,2,-6,-3,-8,-7],
+ [3,0,8,9,0,-5,4],
+ [-9,4,-6,10,-7,6,10],
+ [1,-8,-10,-6,5,-1,-6],
+ [-1,-7,0,0,-6,8,-8],
+ [-7,-7,-6,9,-9,-9,6],
+ [8,-2,-2,6,10,6,7],
+ [-1,6,8,-9,8,7,-3],
+ [8,10,-4,9,1,6,9],
+ [3,8,-2,-7,-7,-1,-8],
+ [-6,9,-8,-4,-6,0,2],
+ [-2,-6,-9,7,5,1,9],
+ [-5,1,-7,-1,-10,8,7],
+ [-7,8,-1,1,-10,-2,-9],
+ [6,-5,3,-3,-2,-1,-4],
+ [9,-5,9,10,-1,-10,9],
+ [4,-1,5,3,-7,9,-8],
+ [8,-2,6,-5,3,-2,2],
+ [1,-1,5,4,3,7,-6],
+ [-10,-7,5,4,3,-2,-10],
+ [-4,8,-8,3,-5,-1,4],
+ [-7,-7,0,-2,-4,-7,5],
+ [3,3,7,-7,-10,-9,6],
+ [0,9,7,8,10,-9,6],
+ [-7,1,-9,1,-5,-4,-2],
+ [-1,10,-9,-5,-3,3,9],
+ [-3,2,6,-3,1,-7,-3],
+ [4,10,8,3,2,7,-1],
+ [3,-3,3,-10,-3,8,10],
+ [-3,7,8,-5,3,4,-4],
+ [-5,2,4,-4,-4,9,10],
+ [1,-1,1,-2,4,5,5],
+ [2,9,10,8,-4,-10,2],
+ [3,-9,2,-3,-7,-3,6],
+ [0,-2,4,-6,10,1,6],
+ [2,-3,5,7,0,8,5],
+ [9,7,4,5,-10,8,8],
+ [6,9,8,6,-8,-8,2],
+ [2,-8,3,6,7,10,-9],
+ [9,3,-4,-9,-10,2,2],
+ [-1,5,1,3,6,7,7],
+ [-10,6,-8,10,5,10,-7],
+ [-10,8,2,2,2,-1,2],
+ [1,7,10,-6,-1,-8,-3],
+ [2,-9,8,-4,-9,4,-10],
+ [7,4,-9,7,5,10,9],
+ [4,8,-8,-2,-7,-3,-9],
+ [-5,0,7,-1,-8,-9,7],
+ [-10,-7,-7,1,3,7,-5],
+ [-9,-5,0,-2,-8,2,-2],
+ [9,9,10,5,-6,1,2],
+ [6,10,-2,6,-6,-2,-1],
+ [-4,-5,4,-10,-3,0,9],
+ [-9,-5,-2,-2,7,-7,3],
+ [5,4,-7,7,-5,6,2],
+ [5,-6,-6,-4,-7,-5,-8],
+ [-10,-9,-10,-1,-6,-1,8],
+ [2,-10,2,-6,8,-9,-3],
+ [-8,5,5,-2,2,-7,10],
+ [-3,-3,6,3,-9,-8,3],
+ [2,-6,3,1,8,-8,9],
+ [-3,-7,2,0,-2,10,-8],
+ [9,-9,-2,-8,-4,5,-8],
+ [-5,-4,-10,-5,-1,7,10],
+ [9,-1,8,-8,4,-2,4],
+ [7,0,-2,3,6,2,-8],
+ [-9,7,6,0,3,7,-2],
+ [10,7,-8,-7,-10,1,8],
+ [-9,9,9,1,-7,9,6],
+ [4,8,-8,-9,2,-1,-1],
+ [-8,7,7,-2,-3,0,9],
+ [8,-3,-1,0,-9,5,6],
+ [5,-10,-8,-1,8,-2,-3],
+ [-3,-5,10,5,-5,-8,-3],
+ [10,1,-2,-3,8,-8,-10],
+ [3,9,-9,0,10,0,2],
+ [4,1,-5,-6,-9,-1,6],
+ [-8,4,1,10,1,-3,6],
+ [-7,5,8,-3,10,4,5],
+ [5,-9,10,-3,-3,-2,-8],
+ [4,-10,-2,-6,10,8,9],
+ [-8,5,9,-2,-4,-4,-3],
+ [7,-4,1,-9,-6,-6,-5],
+ [-9,10,-2,4,-4,-8,2],
+ [-4,7,1,-2,-4,-9,10],
+ [-10,-8,10,-7,3,5,3],
+ [7,6,-10,-5,-5,-7,-5],
+ [3,9,6,9,-6,-5,-3],
+ [2,-6,-4,3,-2,-2,6],
+ [-3,-6,4,8,3,0,5],
+ [-5,1,-8,-5,6,5,3],
+ [-7,-3,1,7,2,1,6],
+ [-8,0,8,-8,10,-6,-2],
+ [-4,-6,-10,5,10,-2,10],
+ [8,6,5,9,-2,-7,2],
+ [10,-6,7,-10,1,3,4],
+ [2,-6,-8,-2,2,1,9],
+ [8,-9,-3,-6,-4,6,-3],
+ [7,-2,-10,7,0,-3,-10],
+ [8,1,3,3,-4,2,-6],
+ [-1,-10,9,1,8,10,0],
+ [-6,9,-10,-9,-6,10,-5],
+ [-4,0,-5,3,-9,1,4],
+ [4,-3,-5,-6,0,8,5],
+ [10,5,-6,-9,6,-4,6],
+ [4,7,1,-1,0,-4,7],
+ [-8,-2,-1,-7,7,4,9],
+ [10,-7,-4,2,7,-9,7],
+ [-2,10,-6,-6,-3,1,5],
+ [7,-9,-9,5,-4,5,-5],
+ [-1,2,-5,-5,-8,-1,-3],
+ [-5,-9,4,4,6,-4,1],
+ [-5,-1,-2,5,2,-3,-1],
+ [-5,9,2,-2,7,-6,-1],
+ [-1,-10,-7,6,-9,9,9],
+ [-4,6,-1,-10,10,2,5],
+ [3,7,-4,5,3,-9,2],
+ [4,4,-7,4,0,-9,8],
+ [6,8,-6,-1,10,2,-4],
+ [-3,-6,-7,6,0,-10,8],
+ [-5,6,-10,10,-9,-3,-6],
+ [3,3,-7,2,2,6,0],
+ [2,-3,-7,-10,6,3,-10],
+ [-5,4,-3,-1,8,0,-6],
+ [10,-9,-2,-9,1,-3,2],
+ [7,8,0,-2,-3,0,5],
+ [-6,6,9,-4,1,0,-6],
+ [-3,2,5,-5,3,4,9],
+ [2,-5,-9,-7,0,-7,4],
+ [-9,1,-7,7,0,-1,-3],
+ [9,-8,-10,-8,-5,5,-10],
+ [-9,-5,-4,4,-10,6,-5],
+ [10,-4,-10,-2,-7,7,-2],
+ [3,-9,1,-9,-5,5,0],
+ [-9,3,8,1,9,-2,8],
+ [2,0,-10,2,-9,4,-9],
+ [1,-1,7,7,5,10,-8],
+ [1,0,5,-10,2,2,-1],
+ [-8,4,5,-5,2,-9,5],
+ [9,10,3,-7,8,7,7],
+ [4,-3,-10,2,6,-5,-8],
+ [-2,-4,7,-6,10,-2,-3],
+ [7,0,-1,5,-7,5,3],
+ [-8,3,5,8,5,-7,8],
+ [3,-5,-2,6,2,2,-6],
+ [-10,-7,4,-1,1,-8,-8],
+ [10,-4,9,6,10,-6,6],
+ [-10,-1,2,8,0,-3,-8],
+ [6,2,4,1,7,2,-3],
+ [7,-5,6,6,1,9,0],
+ [-10,6,8,-4,7,0,8],
+ [10,3,7,0,10,8,-10],
+ [2,-9,-6,-6,-6,-3,-7],
+ [-7,9,1,6,8,10,-7],
+ [-10,0,-8,-1,4,8,1],
+ [-2,5,-10,6,-7,3,3],
+ [1,-1,-4,3,7,-7,0],
+ [8,2,3,-2,6,-1,-9],
+ [5,0,0,-7,-1,-9,9],
+ [9,-7,-4,1,-9,-2,-10],
+ [-9,-10,9,-9,-9,4,-7],
+ [8,6,-3,2,-7,1,1],
+ [-10,-4,-6,-10,-5,-4,-8],
+ [-7,5,-7,10,-9,-9,7],
+ [-3,10,-8,-9,-9,4,8],
+ [-10,1,-6,5,7,-6,10],
+ [10,-6,5,-1,-8,3,-5],
+ [-3,2,-3,10,-4,10,-7],
+ [4,4,-9,-7,7,-2,3],
+ [10,2,-6,-1,8,10,1],
+ [0,2,9,-7,1,10,2],
+ [5,0,-8,-7,-10,8,4],
+ [10,7,-4,9,7,6,8],
+ [2,-3,-1,2,9,8,-5],
+ [10,-2,-2,-7,9,-7,0],
+ [0,5,-7,-5,0,-1,-3],
+ [3,1,6,0,-2,-5,0],
+ [5,-2,-10,7,4,5,5],
+ [-2,-10,4,7,10,-7,8],
+ [-7,-2,-1,-8,-3,0,7],
+ [-4,3,-10,1,1,1,10],
+ [1,4,-4,-7,0,1,-8],
+ [3,4,-2,10,-10,7,6],
+ [4,-5,-9,2,10,-1,-2],
+ [-7,6,-5,6,-8,-4,-6],
+ [-3,9,-1,-7,0,5,-2],
+ [-7,0,8,-2,-2,-6,-5],
+ [4,0,1,-7,8,4,0],
+ [-6,8,9,0,-4,-7,-3],
+ [-7,-1,-5,-6,7,1,7],
+ [-4,-9,-1,8,-10,0,3],
+ [7,-1,10,-7,2,-3,-4],
+ [-9,5,-5,-6,9,-1,0],
+ [10,4,5,-9,3,2,-5],
+ [5,-10,9,5,7,0,3],
+ [3,-3,0,-10,0,-6,0],
+ [-8,3,2,-5,-5,2,-1],
+ [-4,-10,-10,10,-4,4,3],
+ [6,1,-9,5,6,8,9],
+ [-8,2,3,5,7,7,5],
+ [-4,0,7,-3,6,-8,-5],
+ [6,9,8,7,9,0,-9],
+ [-9,9,-7,9,-7,-10,0],
+ [-1,7,6,-6,4,-6,6],
+ [-9,-9,-7,8,-6,9,-3],
+ [5,-7,-1,-8,3,6,3],
+ [-1,-2,-10,2,-10,10,3],
+ [-5,-6,3,2,6,5,-7],
+ [-8,1,8,-1,7,-2,3],
+ [1,-6,-4,-8,2,2,6],
+ [3,0,1,6,-1,-1,-2],
+ [-10,0,1,8,2,0,3],
+ [6,-5,-6,7,9,-9,-10],
+ [-10,-10,-10,-5,-5,6,0],
+ [8,-4,-8,1,-10,8,-3],
+ [7,-3,3,-6,1,5,-6],
+ [8,10,10,7,6,-3,0],
+ [4,-6,6,-3,6,0,7],
+ [9,-7,-3,10,-2,-2,-4],
+ [4,-1,-4,-8,-3,7,-7],
+ [-8,2,9,-4,0,-5,-7],
+ [-2,9,7,7,3,7,-5],
+ [8,5,10,-5,-8,-6,-3],
+ [0,6,-10,-5,-6,-7,10],
+ [-9,-2,3,-10,5,-9,-1],
+ [9,-4,-6,-10,10,7,8],
+ [-7,-4,-1,-1,4,7,-7],
+ [-9,-3,-5,-2,-7,-5,7],
+ [-6,8,4,-3,10,-7,9],
+ [-10,-2,4,1,8,-2,-3],
+ [3,-3,-7,-2,9,8,7],
+ [8,9,8,9,9,9,-1],
+ [-3,10,6,-4,3,1,-4],
+ [5,1,1,5,-10,-4,-2],
+ [-7,7,3,2,1,-5,0],
+ [0,-4,-10,-2,9,-3,-3],
+ [2,8,-7,-4,5,7,-9],
+ [4,-4,4,4,0,5,5],
+ [5,-2,2,-3,9,3,8],
+ [4,3,-8,3,-9,4,6]
+];
+  for(i=1,#V,
+    print1(i,":");
+    my(L=hyperellratpoints(Polrev(V[i]),16383));
+    print(apply(affine2proj,L)));
+}