Ticket #13376: smalljac-4.0.14.patch

File smalljac-4.0.14.patch, 70.9 KB (added by pavpanchekha, 10 years ago)

Patch to Sage 5.2 to support interfacing with SmallJac?

  • module_list.py

    diff --git a/module_list.py b/module_list.py
    a b  
    19191919                  depends = [SAGE_LOCAL + "/include/lrcalc/symfcn.h"]), # should include all .h
    19201920        )
    19211921
     1922if is_package_installed('smalljac'):
     1923    ext_modules.append(
     1924            Extension('sage.libs.smalljac.wrapper2',
     1925                sources = ["sage/libs/smalljac/wrapper2.pyx"],
     1926                libraries = ["gmp", "m", "smalljac", "ff_poly"])
     1927            )
  • new file sage/libs/smalljac/__init__.py

    diff --git a/sage/libs/smalljac/__init__.py b/sage/libs/smalljac/__init__.py
    new file mode 100644
    - +  
     1# smalljac -- Drew Sutherland's code for point counting on Jacobians of curves over finite fields
     2# This file is only here to make this directory a Python package
  • new file sage/libs/smalljac/all.py

    diff --git a/sage/libs/smalljac/all.py b/sage/libs/smalljac/all.py
    new file mode 100644
    - +  
     1from wrapper2 import SmallJac
  • new file sage/libs/smalljac/defs.pxd

    diff --git a/sage/libs/smalljac/defs.pxd b/sage/libs/smalljac/defs.pxd
    new file mode 100644
    - +  
     1cdef extern from "smalljac.h":
     2    ctypedef void *smalljac_curve_t
     3   
     4    smalljac_curve_t smalljac_curve_init (char *str, int *err)
     5    void  smalljac_curve_clear (smalljac_curve_t c)
     6    char *smalljac_curve_str   (smalljac_curve_t c)   
     7    int   smalljac_curve_genus (smalljac_curve_t c)
     8    char *smalljac_curve_nf    (smalljac_curve_t c)
     9    char *smalljac_curve_cm    (smalljac_curve_t c)
     10   
     11    long smalljac_Lpolys (smalljac_curve_t c,
     12                          unsigned long start,
     13                          unsigned long end,                           
     14                          unsigned long flags,
     15                          int (*callback)(smalljac_curve_t c, unsigned long q, int good, long a[], int n, void *arg),
     16                          void *arg)
     17    long smalljac_parallel_Lpolys (smalljac_curve_t c,
     18                                   unsigned long start,
     19                                   unsigned long end,                           
     20                                   unsigned long flags,
     21                                   int (*callback)(smalljac_curve_t c, unsigned long q, int good, long a[], int n, void *arg),
     22                                   void *arg)
     23    int smalljac_Lpoly (long a[], char *curve_str, unsigned long q, unsigned long flags)
     24   
     25    long smalljac_groups (smalljac_curve_t curve,
     26                          unsigned long start,
     27                          unsigned long end,                           
     28                          unsigned long flags,
     29                          int (*callback)(smalljac_curve_t curve, unsigned long q, int good, long m[], int n, void *arg),
     30                          void *arg)
     31    long smalljac_parallel_groups (smalljac_curve_t curve,
     32                                   unsigned long start,
     33                                   unsigned long end,                           
     34                                   unsigned long flags,
     35                                   int (*callback)(smalljac_curve_t curve, unsigned long q, int good, long m[], int n, void *arg),
     36                                   void *arg)
     37    int smalljac_group (long m[], char *curve_str, unsigned long q, unsigned long flags)
     38
     39    cdef int SMALLJAC_A1_ONLY
     40    cdef int SMALLJAC_DEGREE1_ONLY
     41    cdef char *SMALLJAC_VERSION_STRING
     42
     43    cdef int SMALLJAC_INTERNAL_ERROR
     44    cdef int SMALLJAC_INVALID_INTERVAL
     45    cdef int SMALLJAC_PARSE_ERROR
     46    cdef int SMALLJAC_UNSUPPORTED_CURVE
     47    cdef int SMALLJAC_SINGULAR_CURVE
     48    cdef int SMALLJAC_INVALID_PP
     49    cdef int SMALLJAC_WRONG_GENUS
     50    cdef int SMALLJAC_INVALID_FLAGS
     51    cdef int SMALLJAC_NOT_OVER_Q
     52
     53    cdef int SMALLJAC_CURVE_STRING_LEN
     54
     55    int smalljac_moments (smalljac_curve_t curve,
     56                          unsigned long start, unsigned long end,
     57                          double moments[], int n, int m,
     58                          char STgroup[16],
     59                          int (*filter_callback)(smalljac_curve_t curve,
     60                                                 unsigned long q,
     61                                                 void *arg),
     62                          void *arg)
     63    int smalljac_identify_STgroup (smalljac_curve_t c, char STgroup[16], long maxp)
  • new file sage/libs/smalljac/util.py

    diff --git a/sage/libs/smalljac/util.py b/sage/libs/smalljac/util.py
    new file mode 100644
    - +  
     1from sage.rings.integer import Integer
     2
     3def parse_limits(limit):
     4    start, end = 0, 0
     5    if isinstance(limit, int) or isinstance(limit, Integer):
     6        start = 0
     7        end = int(limit)
     8    elif isinstance(limit, list):
     9        start = limit[0]
     10        end = limit[-1]
     11    elif isinstance(limit, tuple):
     12        start = limit[0] + 1
     13        end = limit[-1] - 1
     14    else:
     15        raise ValueError("Cannot interpret limits for trace computation", limit)
     16
     17    if start < 0 or end < 0 or end < start:
     18        raise ValueError("Invalid limits; both must be nonnegative, and end must be not less than start")
     19
     20    return start, end
     21
     22def interpret_smalljac_results(res, correction_fn, construct_fn, raw=False, index="primes", container=list):
     23    if raw: return res
     24   
     25    def correct_values(p, values):
     26        """ Correct in case of bad primes """
     27        if not values.count(None): return values
     28        if correction_fn:
     29            correct = correction_fn(p, values)
     30        else:
     31            return values
     32        # Mutating copy
     33        for i in range(len(values)):
     34            values[i] = correct[i]
     35
     36    def make_value(p, value):
     37        if value is None:
     38            return None
     39        else:
     40            return construct_fn(p, value)
     41   
     42    last_p = 0
     43    last_p_lst = []
     44    result = []
     45    for p, val in res:
     46        if p == last_p:
     47            last_p_lst.append(make_value(p, val))
     48        else:
     49            correct_values(last_p, last_p_lst)
     50           
     51            last_p = p
     52            last_p_lst = [make_value(p, val)]
     53            result.append((Integer(p), last_p_lst))
     54   
     55    # Index by ideals, not primes
     56    if index == "ideals":
     57        result2 = []
     58        for p, vals in result:
     59            for i, val in enumerate(vals):
     60                result2.append(((p, Integer(i)), val))
     61        result = result2
     62    elif index == "primes":
     63        pass
     64    else:
     65        raise ValueError("Traces cannot be indexed by `%s`; check your `index` parameter" % index)
     66   
     67    return container(result)
  • new file sage/libs/smalljac/wrapper2.pyx

    diff --git a/sage/libs/smalljac/wrapper2.pyx b/sage/libs/smalljac/wrapper2.pyx
    new file mode 100644
    - +  
     1# -*- mode: python -*-
     2from defs cimport *
     3
     4cdef extern from "math.h":
     5    double sqrt(double x)
     6cdef extern from "stdlib.h":
     7    void *malloc(int size)
     8    void free(void *)
     9
     10cdef int callback_list_single(smalljac_curve_t c, unsigned long p, int good,  long a[],  int n, void *arg):
     11
     12    if good:
     13        result = (int(p), -a[0])
     14    else:
     15        result = (int(p), None)
     16   
     17    (<SmallJac>arg).tmp.append(result)
     18    return 1
     19
     20cdef int callback_single_single(smalljac_curve_t c, unsigned long p, int good,  long a[],  int n, void *arg):
     21    if good:
     22        (<SmallJac>arg).tmp = -a[0]
     23    else:
     24        (<SmallJac>arg).tmp = None
     25    return 1
     26
     27cdef int callback_list_list(smalljac_curve_t c, unsigned long q, int good, long m[], int n, void *arg):
     28    cdef int i
     29
     30    if good:
     31        lst = []
     32        for i in range(n):
     33            lst.append(m[i])
     34        result = (int(q), lst)
     35    else:
     36        result = (int(q), None)
     37
     38    (<SmallJac>arg).tmp.append(result)
     39    return 1
     40
     41cdef int callback_single_list(smalljac_curve_t c, unsigned long q, int good, long m[], int n, void *arg):
     42    cdef int i
     43    if good:
     44        lst = []
     45        for i in range(n):
     46            lst.append(m[i])
     47    else:
     48        lst = None
     49
     50    (<SmallJac>arg).tmp = lst
     51    return 1
     52
     53cdef int callback_histogram_single(smalljac_curve_t c, unsigned long q, int good, long a[], int n, void *arg):
     54    cdef double val
     55    cdef unsigned int bucket
     56    cdef unsigned int buckets
     57    cdef unsigned int *array = <unsigned int *>arg
     58    cdef int genus
     59
     60    genus = array[0]
     61    buckets = array[1]
     62
     63    if good:
     64        val = 0.5 * (0.5 * -a[0] / sqrt(q) / genus + 1)
     65        bucket = (<int>(val * buckets))
     66        if bucket >= buckets:
     67            print "Anomaly at p=%d: a1=%d, so fraction=%f" % (q, a[0], val)
     68        else:
     69            array[bucket + 2] += 1
     70
     71    return 1
     72
     73smalljac_version = SMALLJAC_VERSION_STRING
     74
     75cdef class SmallJac:
     76    cdef smalljac_curve_t c
     77    cdef object tmp
     78    cdef int genus
     79   
     80    def __cinit__(self):
     81        self.c = <smalljac_curve_t>0
     82       
     83    def __init__(self, v):
     84        curve = str(v)
     85        cdef int err
     86
     87        self.c = smalljac_curve_init(curve, &err)
     88        self.genus = smalljac_curve_genus(self.c)
     89
     90        if not err:
     91            return
     92        elif err == SMALLJAC_INTERNAL_ERROR:
     93            raise RuntimeError("Internal Smalljac Error")
     94        elif err == SMALLJAC_INVALID_INTERVAL:
     95            raise RuntimeError("Invalid Interval")
     96        elif err == SMALLJAC_PARSE_ERROR:
     97            raise RuntimeError("Cannot parse arguments")
     98        elif err == SMALLJAC_UNSUPPORTED_CURVE:
     99            raise RuntimeError("Unsupported curve; must be degree 3 or 5")
     100        elif err == SMALLJAC_SINGULAR_CURVE:
     101            raise RuntimeError("Curve is singular")
     102        elif err == SMALLJAC_INVALID_PP:
     103            raise RuntimeError("Not prime power, or too large a prime power, passed where prime power is required")
     104        elif err == SMALLJAC_WRONG_GENUS:
     105            raise RuntimeError("Only curves of genus 1 and 2 supported")
     106        elif err == SMALLJAC_INVALID_FLAGS:
     107            raise RuntimeError("Invalid combination of bit flags set")
     108        elif err == SMALLJAC_NOT_OVER_Q:
     109            raise RuntimeError("Curve is not defined over QQ")
     110        else:
     111            raise RuntimeError("Unknown error %s" % err)
     112
     113    def __repr__(self):
     114        return "Smalljac curve defined by %s"%smalljac_curve_str(self.c)
     115
     116    def __dealloc__(self):
     117        if self.c:
     118            smalljac_curve_clear(self.c)
     119
     120    def ap(self, unsigned long p, unsigned long b=0, parallel=True):
     121        """
     122        If only p is given, return trace of Frobenius at p.  If both p
     123        and b are given, return a list of (key, value) pairs, with
     124        keys the primes < b, and values the traces of Frobenius at the
     125        good primes, and None at the bad primes, where bad is "bad for
     126        the given model".
     127        """
     128
     129        if b == 0:
     130            # Compute for a single p only
     131            smalljac_parallel_Lpolys(self.c, p, p, SMALLJAC_A1_ONLY, callback_single_single, <void*>self)
     132            return self.tmp
     133        # Compute for a range of primes, starting with p.
     134        self.tmp = []
     135        smalljac_parallel_Lpolys(self.c, p, b,  SMALLJAC_A1_ONLY, callback_list_single, <void*>self)
     136        return self.tmp
     137
     138    def sato_tate_buckets(self, unsigned int n, unsigned long a, unsigned long b):
     139        cdef unsigned long i
     140        cdef unsigned int *buckets
     141        buckets = <unsigned int *>(malloc((n + 2) * sizeof(unsigned int)))
     142       
     143        for i in range(n): # Memset would be faster, but it doesn't matter
     144            buckets[i+2] = 0
     145
     146        buckets[0] = self.genus
     147        buckets[1] = n
     148       
     149        smalljac_parallel_Lpolys(self.c, a, b,  SMALLJAC_A1_ONLY, callback_histogram_single, <void*>buckets)
     150
     151        self.tmp = []
     152        for i in range(n):
     153            self.tmp.append(buckets[i + 2])
     154       
     155        free(buckets)
     156        return self.tmp
     157
     158    def group(self, unsigned long p, unsigned long b=0):
     159        """
     160        If only p is given, return the group structure of the curve
     161        reduced at p.  If both p and b are given, return a list of
     162        (key, value) pairs, with keys the primes < b, and values are,
     163        at bad primes, None, and at good primes is a list of integers
     164        m1, m2, ..., representing the group Z/m1Z x Z/m2 x ...
     165        """
     166        if b == 0:
     167            # Compute for a single p only
     168            smalljac_parallel_groups(self.c, p, p, 0, callback_single_list, <void*>self)
     169            return self.tmp
     170        # Compute for a range of primes, starting with p.
     171        self.tmp = []
     172        smalljac_parallel_groups(self.c, p, b, 0, callback_list_list, <void*>self)
     173        return self.tmp
     174
     175    def lpoly(self, unsigned long p, unsigned long b=0):
     176        """
     177        If only p is given, return the L polynomial of the curve
     178        reduced at p.  If both p and b are given, return a list of
     179        (key, value) pairs, with keys the primes < b, and values are,
     180        at bad primes, None, and at good primes is a list of integers
     181        m1, m2, ..., representing the L polynomial
     182        p^n T^2n + m1 p^n-1 T^2n-1 m2 p^n-1 T^2n-2 + ... + m1 T + 1
     183        """
     184        if b == 0:
     185            # Compute for a single p only
     186            smalljac_parallel_Lpolys(self.c, p, p, 0, callback_single_list, <void*>self)
     187            return self.tmp
     188        # Compute for a range of primes, starting with p.
     189        self.tmp = []
     190        smalljac_parallel_Lpolys(self.c, p, b, 0, callback_list_list, <void*>self)
     191        return self.tmp
     192
     193    def moments(self, int moments, unsigned long start, unsigned long end):
     194        """
     195        Stuff
     196        """
     197
     198        if moments > 11: raise ValueError("Only up to 11 moments supported for now.")
     199
     200        cdef double *moment_arr = <double *>malloc(sizeof(double) * moments * self.genus)
     201
     202        smalljac_moments(self.c, start, end, moment_arr, self.genus, moments,
     203                         NULL, NULL, NULL)
     204
     205        a1moments = []
     206
     207        cdef int i
     208        for i in range(moments):
     209            a1moments.append(moment_arr[i])
     210
     211        if self.genus == 2:
     212            a2moments = []
     213            for i in range(moments):
     214                a2moments.append(moment_arr[i + moments])
     215
     216            free(moment_arr)
     217            return (a1moments, a2moments)
     218        else:
     219            free(moment_arr)
     220            return (a1moments,)
     221
     222    def guess_group(self, unsigned long max):
     223        """
     224        More stuff
     225        """
     226
     227        cdef int i
     228        cdef char STgroup[16]
     229
     230        i = smalljac_identify_STgroup(self.c, STgroup, max)
     231
     232        if i:
     233            return STgroup
     234        else:
     235            return None
  • sage/schemes/elliptic_curves/ell_finite_field.py

    diff --git a/sage/schemes/elliptic_curves/ell_finite_field.py b/sage/schemes/elliptic_curves/ell_finite_field.py
    a b  
    12841284        """
    12851285        return self.points()[n]
    12861286
     1287    def _smalljac(self):
     1288        r"""
     1289        Computes and caches the smalljac representation of the given
     1290        curve.  For internal use only.
     1291        """
     1292
     1293        try:
     1294            return self._smalljac_proxy
     1295        except AttributeError:
     1296            try:
     1297                import sage.libs.smalljac.all as smalljac
     1298            except ImportError:
     1299                raise RuntimeError, "smalljac not installed or inaccessible"
     1300           
     1301            f = "[%s,%s,%s,%s,%s]" % tuple(self.ainvs())
     1302            self._smalljac_proxy = smalljac.SmallJac(f)
     1303            return self._smalljac_proxy
     1304
     1305    def jacobian_structure(self, algorithm="smalljac"):
     1306        r"""
     1307        Returns the group structure of the group of points on this
     1308        elliptic curve.  Unlike ``abelian_group``, it does not compute
     1309        the generators of the group, just the structure, and is thus
     1310        perhaps somewhat less useful.
     1311
     1312        INPUT:
     1313       
     1314        - ``algorithm`` - either "sage", which uses Sage's abelian
     1315          group determination algorithm via ``abelian_group``, or
     1316          ``smalljac``, which uses Drew Sutherland's SmallJac library
     1317          for the computation.  SmallJac is much faster, but can only
     1318          be used on degree-1 fields.  You can also pass
     1319          ``algorithm="all"`` to run both algorithms and verify that
     1320          they produce the same results.
     1321
     1322        OUTPUT:
     1323
     1324        An abelian group isomorphic to the jacobian of this curve.
     1325
     1326        EXAMPLES::
     1327
     1328            sage: E=EllipticCurve(GF(11),[2,5])
     1329            sage: E.jacobian_structure() # optional: smalljac
     1330            Additive abelian group isomorphic to Z/10
     1331            sage: E=EllipticCurve(GF(41),[2,5])
     1332            sage: E.jacobian_structure() # optional: smalljac
     1333            Additive abelian group isomorphic to Z/2 + Z/22
     1334       
     1335        We can use ``algorithm="all"`` to test multiple
     1336        implementations against each other::
     1337       
     1338            sage: E.jacobian_structure(algorithm="all") # optional: smalljac
     1339            Additive abelian group isomorphic to Z/2 + Z/22
     1340        """
     1341
     1342        from sage.groups.additive_abelian.additive_abelian_group import AdditiveAbelianGroup
     1343
     1344        if algorithm == "smalljac":
     1345            K = self.base_field()
     1346            assert K.degree() == 1, "Smalljac can only find group structures in degree-1 fields"
     1347            p = K.cardinality()
     1348            proxy = self._smalljac()
     1349            structure = proxy.group(p)
     1350            return AdditiveAbelianGroup(structure)
     1351        elif algorithm == "sage":
     1352            g = self.abelian_group()
     1353            return AdditiveAbelianGroup(g.generator_orders())
     1354        elif algorithm == "all":
     1355            smalljac_results = self.jacobian_structure(algorithm="smalljac")
     1356            sage_results     = self.jacobian_structure(algorithm="sage")
     1357
     1358            # If the orders were given in different orders, testing
     1359            # two groups for equality can give incorrect results
     1360            smalljac_order = sorted(g.additive_order() for g in smalljac_results.gens())
     1361            sage_order     = sorted(g.additive_order() for g in sage_results.gens())
     1362            assert smalljac_order == sage_order, "An error has occured and smalljac and Sage's results do not match"
     1363            return smalljac_results
     1364        else:
     1365            raise ValueError("Unknown algorithm `%s`" % algorithm)
     1366       
    12871367    def abelian_group(self, debug=False):
    12881368        r"""
    12891369        Returns the abelian group structure of the group of points on this
  • sage/schemes/elliptic_curves/ell_number_field.py

    diff --git a/sage/schemes/elliptic_curves/ell_number_field.py b/sage/schemes/elliptic_curves/ell_number_field.py
    a b  
    20982098            degrees.extend(newdegs)
    20992099            k = k+1
    21002100
    2101         raise NotImplementedError, "Not all isogenies implemented over general number fields."
     2101        raise NotImplementedError, "Not all isogenies implemented over general number fields."
     2102
     2103    def aplists(self, limit, algorithm="smalljac", container=list, index="primes", raw=False):
     2104        """
     2105        Computes the traces of Frobenius at the prime ideals of the
     2106        base number field, sorted by the norm of these ideals.  The
     2107        structure of the return values depends on the ``format``
     2108        argument.  So far, this method requires the Smalljac library.
     2109
     2110        INPUT:
     2111
     2112        - ``limit`` - either an integer, in which case the prime
     2113          ideals will have norm less than this limit, or [a, b], in
     2114          which case the prime ideals with norms between ``a`` and
     2115          ``b`, inclusive, are used.
     2116
     2117        - ``algorithm`` - so far, only "smalljac" is supported, which
     2118          uses Drew Sutherland's smalljac library for elliptic curve
     2119          structure computation.
     2120
     2121        - The remainder (``container``, ``index``, ``raw``)
     2122          define the format of the result.  By default, the return
     2123          value looks like this:
     2124
     2125            [(9, [2]), (11, [4, -4]), (19, [-4, 4]), (31, [4, 8])]
     2126
     2127          This describes the fact that 9 does not split, and its prime
     2128          ideal has trace of Frobeneus 2; that 11 splits, and its two
     2129          prime ideals have traces 4 and -4; and so on.
     2130
     2131          Instead of a list of tuples, one can choose a dictionary,
     2132          mapping norms to traces, by changing ``container`` to
     2133          ``dict`` (in fact, any value for container is just called
     2134          with the list of pairs).
     2135
     2136          Instead of pairs of norms and lists of traces, you can have
     2137          each trace have its own pair, with the first element being a
     2138          pair of norm and index, yielding a result such as:
     2139
     2140            [((9, 0), 2), ((11, 0), 4), ((11, 0), -4), ...]
     2141
     2142          This is chosen with ``index="ideals"``.
     2143
     2144          By default, all integers above are Sage Integers.  If this
     2145          is a performance or memory problem, you likely want
     2146          ``raw=True``.
     2147
     2148          Finally, for performance or memory reasons, you can turn off
     2149          all post-processing by by setting ``raw`` to ``True``.  This
     2150          overrides all other arguments and returns a list formatted
     2151          like so:
     2152
     2153            [(9, 2), (11, None), (11, -4), (19, 4), (19, -4), ...]
     2154
     2155          All integers are python integers.  This last mode is by far
     2156          the fastest, especially for large ranges of prime norms, but
     2157          is also often the hardest to use for other computation.  It
     2158          is often used internally.  Note that, among other things,
     2159          using only raw data means that bad primes will have the
     2160          trace ``None``, instead of the correct trace, since it is
     2161          not SmallJac, but Sage, that does that computation.
     2162
     2163        OUTPUT:
     2164
     2165        A list of traces of Frobenius, with the precise format defined
     2166        by the ``container``, ``index``, and ``raw`` keywords.
     2167
     2168        EXAMPLES:
     2169
     2170        We can find the traces of Frobenius up to a given prime::
     2171
     2172            sage: ec = EllipticCurve('11a')
     2173            sage: ec.aplists(15)
     2174            [(2, [-2]), (3, [-1]), (5, [1]), (7, [-2]), (11, [None]), (13, [4])]
     2175
     2176        Or in any arbitrary range::
     2177
     2178            sage: ec.aplists([15, 40])
     2179            [(17, [-2]), (19, [0]), (23, [-1]), (29, [0]), (31, [7]), (37, [3])]
     2180
     2181        Curves over number fields are also supported::
     2182
     2183            sage: K.<a> = NumberField(x^2 - x - 1)
     2184            sage: ec2 = EllipticCurve(K, [a, -a+1, 0, -4, 3*a-5])
     2185            sage: ec2.aplists([15, 50])
     2186            [(19, [-4, 4]), (29, [-2, 6]), (31, [-8, 8]), (41, [2, -6]), (49, [2])]
     2187
     2188        You can choose to separate each ideal, instead of combining them by norm::
     2189
     2190            sage: ec2.aplists([15, 30], index="ideals")
     2191            [((19, 0), -4), ((19, 1), 4), ((29, 0), -2), ((29, 1), 6)]
     2192
     2193        You can also choose a dictionary representation, though this
     2194        of course destroys the ordering::
     2195
     2196            sage: ec2.aplists([15, 50], container=dict)
     2197            {41: [2, -6], 49: [2], 19: [-4, 4], 29: [-2, 6], 31: [-8, 8]}
     2198            sage: ec2.aplists([15, 30], container=dict, index="ideals")
     2199            {(19, 0): -4, (29, 1): 6, (19, 1): 4, (29, 0): -2}
     2200
     2201        Finally, the raw results are useful for large-scale statistics::
     2202
     2203            sage: bsum = __builtins__.sum
     2204            sage: bsum(ap for p, ap in ec.aplists(100000, raw=True) if ap)
     2205            4837
     2206            sage: bsum(ap for p, ap in ec2.aplists(100000, raw=True) if ap)
     2207            5289
     2208
     2209        Note that the equivalent non-SmallJac code for that last
     2210        example returns a different answer::
     2211
     2212            sage: bsum(ec.aplist(100000))
     2213            4838
     2214
     2215        This is because our curve has a bad reduction at 11, and we
     2216        aren't adding it its trace.
     2217        """
     2218
     2219        import sage.libs.smalljac.util as util
     2220
     2221        start, end = util.parse_limits(limit)
     2222       
     2223        if algorithm == "smalljac":
     2224            proxy = self._smalljac()
     2225            res = proxy.ap(start, end)
     2226
     2227            O = self.base_ring().ring_of_integers()
     2228            def get_aps(p, aps):
     2229                if hasattr(O.ideal(p), "factor"):
     2230                    factors = [ideal for ideal, exponent in O.ideal(p).factor()]
     2231                else:
     2232                    factors = [p]
     2233               
     2234                out = []
     2235                for ideal in factors:
     2236                    try:
     2237                        out.append(self.reduction(ideal).trace_of_frobenius())
     2238                    except AttributeError: # Bad reduction
     2239                        out.append(None)
     2240                return out
     2241            def make_int(p, v):
     2242                return Integer(v)
     2243
     2244            return util.interpret_smalljac_results(
     2245                res, correction_fn=get_aps, construct_fn=make_int,
     2246                raw=raw, index=index, container=container)
     2247        else:
     2248            raise ValueError("Algorithm %s not defined" % algorithm)
     2249
     2250    def aps(self, p):
     2251        """
     2252        Returns a list of the traces of Frobenius of this elliptic
     2253        curve when reduced at prime ideals of norm ``p``.
     2254
     2255        INPUT:
     2256
     2257         - ``p`` - The integer prime at which to reduce.
     2258
     2259        OUTPUT:
     2260
     2261        A list of traces of Frobenius, or None if there was a bad
     2262        reduction.
     2263
     2264        EXAMPLES:
     2265
     2266        We can take a curve and find traces::
     2267
     2268            sage: ec = EllipticCurve('11a')
     2269            sage: ec.aps(next_prime(10^6))
     2270            [284]
     2271            sage: ec.aps(next_prime(10^9))
     2272            [-1962]
     2273
     2274        Curves over quadratic number fields are also supported::
     2275
     2276            sage: K.<a> = NumberField(x^2 - x - 1)
     2277            sage: ec2 = EllipticCurve(K, [a, -a+1, 0, -4, 3*a-5])
     2278            sage: ec2.aps(9949)
     2279            [-162, -34]
     2280        """
     2281
     2282        res = self.aplists([p, p])
     2283
     2284        if len(res) == 0 or res[0][0] != p:
     2285            raise ValueError("%d is not a prime (over the chosen field)" % p)
     2286
     2287        return res[0][1] # [0] extracts the prime we care about, [1] extracts the list of traces
     2288
     2289    def grouplists(self, limit, algorithm="smalljac", container=list, index="primes", raw=False):
     2290        r"""
     2291        Computes the group structure of reductions at the prime ideals
     2292        of the base number field, sorted by the norm of these ideals.
     2293        So far, this method requires the Smalljac library.
     2294
     2295        INPUT:
     2296
     2297        - ``limit`` - either an integer, in which case the prime
     2298          ideals will have norm less than this limit, or [a, b], in
     2299          which case the prime ideals with norms between ``a`` and
     2300          ``b`, inclusive, are used.
     2301
     2302        - ``algorithm`` - so far, only "smalljac" is supported, which
     2303          uses Drew Sutherland's smalljac library for elliptic curve
     2304          structure computation.
     2305
     2306        - The remainder (``container``, ``index``, and ``raw``) define
     2307          the format of the result; see ``EllipticCurve.aplists`` for
     2308          a description.
     2309
     2310        OUTPUT:
     2311
     2312        A list of group structures of Jacobians at various primes; the
     2313        format is precisely specified by the ``container``, ``index``,
     2314        and ``raw`` arguments.  If a curve has a bad reduction at a prime,
     2315        ``None`` is instead returned.
     2316
     2317        EXAMPLES::
     2318
     2319            sage: ec = EllipticCurve('11a')
     2320            sage: ec.grouplists(10)
     2321            [(2, [Additive abelian group isomorphic to Z/5]),
     2322             (3, [Additive abelian group isomorphic to Z/5]),
     2323             (5, [Additive abelian group isomorphic to Z/5]),
     2324             (7, [Additive abelian group isomorphic to Z/10])]
     2325
     2326            sage: ec.grouplists([50, 60])
     2327            [(53, [Additive abelian group isomorphic to Z/2 + Z/30]), (59, [Additive abelian group isomorphic to Z/55])]
     2328
     2329        The raw results are suitable for larger-scale statistics::
     2330
     2331            sage: def jac_checksum(poly):                                         
     2332            ...       ec = EllipticCurve(poly)
     2333            ...       sum = __builtins__.sum                           
     2334            ...       return sum(sum(gs) for p, gs in ec.grouplists([3, 10000], raw=True) if gs)
     2335            sage: jac_checksum([4, 7])
     2336            5179039
     2337            sage: jac_checksum('11a')
     2338            4121837
     2339            sage: jac_checksum('37b')
     2340            3455881
     2341        """
     2342
     2343        from sage.groups.additive_abelian.additive_abelian_group import AdditiveAbelianGroup
     2344        import sage.libs.smalljac.util as util
     2345
     2346        start, end = util.parse_limits(limit)
     2347
     2348        if algorithm == "smalljac":
     2349            proxy = self._smalljac()
     2350            res = proxy.group(start, end)
     2351
     2352            O = self.base_ring().ring_of_integers()
     2353            def get_group(p, gps):
     2354                if hasattr(O.ideal(2), "factor"): # Not the rationals
     2355                    factors = [ideal for ideal, exponent in O.ideal(p).factor()]
     2356                else: # Only in the case of QQ
     2357                    factors = [p]
     2358               
     2359                out = []
     2360                for ideal in factors:
     2361                    try:
     2362                        out.append(self.reduction(ideal).jacobian_structure())
     2363                    except AttributeError:
     2364                        out.append(None)
     2365                return out
     2366           
     2367            def make_group(p, v):
     2368                return AdditiveAbelianGroup(v)
     2369
     2370            return util.interpret_smalljac_results(
     2371                res, correction_fn=get_group, construct_fn=make_group,
     2372                raw=raw, index=index, container=container)
     2373        else:
     2374            raise ValueError("Algorithm %s not defined" % algorithm)
     2375
     2376    def lpolylists(self, limit, algorithm="smalljac", container=list, index="primes", normalize=False, raw=False):
     2377        r"""
     2378        Computes the L-polynomials of reductions at the prime ideals
     2379        of the base number field, sorted by the norm of these ideals.
     2380        So far, this method requires the Smalljac library.
     2381
     2382        INPUT:
     2383
     2384        - ``limit`` - either an integer, in which case the prime
     2385          ideals will have norm less than this limit, or [a, b], in
     2386          which case the prime ideals with norms between ``a`` and
     2387          ``b`, inclusive, are used.
     2388
     2389        - ``normalize`` - whether to return normalized L-polynomials.
     2390          By default, non-normalized L-polynomials are returned.
     2391
     2392        - ``algorithm`` - so far, only "smalljac" is supported, which
     2393          uses Drew Sutherland's smalljac library for elliptic curve
     2394          structure computation.
     2395
     2396        - The remainder (``container``, ``index``, and ``raw``) define
     2397          the format of the result; see ``EllipticCurve.aplists`` for
     2398          a description.
     2399
     2400        OUTPUT:
     2401
     2402        A list of L-polynomials of reductions at various primes; the
     2403        format is precisely specified by the ``container``, ``index``,
     2404        and ``raw`` arguments.  If a curve has a bad reduction at a
     2405        prime, ``None`` is instead returned.
     2406
     2407        EXAMPLES:
     2408
     2409        We can construct an elliptic curve and find its L-polynomial::
     2410
     2411            sage: E = EllipticCurve('11a')
     2412            sage: E.lpolylists([89, 89])
     2413            [(89, [89*T^2 - 15*T + 1])]
     2414
     2415        Or we can use a range::
     2416
     2417            sage: E.lpolylists([80, 100])
     2418            [(83, [83*T^2 + 6*T + 1]), (89, [89*T^2 - 15*T + 1]), (97,
     2419            [97*T^2 + 7*T + 1])]
     2420
     2421        We can ask for normalized L-polynomials instead::
     2422
     2423            sage: E.lpolylists([89, 89], normalize=True)
     2424            [(89, [T^2 - 1.58999682000954*T + 1.00000000000000])]
     2425
     2426        Number fields are also supported::
     2427
     2428            sage: K.<a> = NumberField(x^2 - x - 1)
     2429            sage: E2 = EllipticCurve(K, [a, -a+1, 0, -4, 3*a-5])
     2430            sage: E2.lpolylists([80, 100])
     2431            [(89, [89*T^2 + 14*T + 1, 89*T^2 - 2*T + 1])]
     2432        """
     2433
     2434        from sage.rings.polynomial.polynomial_ring import PolynomialRing_field
     2435        from sage.rings.real_mpfr import RR
     2436        from sage.rings.rational_field import QQ
     2437        from math import sqrt
     2438        import sage.libs.smalljac.util as util
     2439
     2440        start, end = util.parse_limits(limit)
     2441
     2442        if algorithm == "smalljac":
     2443            proxy = self._smalljac()
     2444            res = proxy.lpoly(start, end)
     2445
     2446            if normalize:
     2447                polynomial_ring = PolynomialRing_field(RR, 'T')
     2448            else:
     2449                polynomial_ring = PolynomialRing_field(QQ, 'T')
     2450            t = polynomial_ring.gens()[0] # I don't want a T
     2451
     2452            def make_lpoly(p, coeffs):
     2453                if not normalize:
     2454                    a1 = coeffs[0]
     2455                    return p*t**2 + a1*t + 1
     2456                else:
     2457                    a1 = coeffs[0]
     2458                    a1 /= sqrt(p)
     2459                    return t**2 + a1*t + 1
     2460
     2461            def get_lpoly(p, coeffs):
     2462                aps = self.aps(p)
     2463
     2464                out = []
     2465                for ap in aps:
     2466                    if not normalize:
     2467                        out.append(p*t**2 - ap*t + 1)
     2468                    else:
     2469                        ap /= sqrt(p)
     2470                        out.append(p*t**2 - ap*t + 1)
     2471                return out
     2472
     2473            return util.interpret_smalljac_results(
     2474                res, correction_fn=get_lpoly, construct_fn=make_lpoly,
     2475                raw=raw, index=index, container=container)
     2476        else:
     2477            raise ValueError("Algorithm %s not defined" % algorithm)
     2478
     2479    def _smalljac(self):
     2480        """
     2481        Computes and caches the smalljac representation of the given
     2482        curve.  For internal use only.
     2483        """
     2484
     2485        try:
     2486            return self._smalljac_proxy
     2487        except AttributeError:
     2488            try:
     2489                import sage.libs.smalljac.all as smalljac
     2490            except ImportError:
     2491                raise RuntimeError, "smalljac not installed or inaccessible"
     2492
     2493            k = self.base_ring()
     2494            gen = k.gen()
     2495            var = str(gen)
     2496            cpoly = gen.charpoly().change_variable_name(var)
     2497            f = "[%s, %s, %s, %s, %s]" % tuple(self.ainvs())
     2498            f = "%s / (%s)" % (f, cpoly)
     2499            self._smalljac_proxy = smalljac.SmallJac(f)
     2500            return self._smalljac_proxy
     2501
     2502    def trace_histogram(self, limit, buckets=None):
     2503        """
     2504        Computes the traces of reductions of this curve at all primes
     2505        in a range, and returns how many, normalized, fell in each
     2506        of ``buckets`` buckets from -2 to 2.
     2507
     2508        INPUT:
     2509
     2510        - ``limits`` - Either an integer, meaning all primes below that
     2511          integer are used, or a list of ``[start, end]``, in which
     2512          case all primes from ``start`` to ``end``, inclusive, are
     2513          used.
     2514
     2515        - ``buckets`` - How many buckets to split the range $[-2, 2]$
     2516          into.  By default, the number of buckets is the square root
     2517          of the range size.
     2518
     2519        OUTPUT:
     2520
     2521        A list of integers, as many as the ``buckets`` parameter,
     2522        denoting the number of trace values fell into that bucket for
     2523        primes in the given range.
     2524
     2525        EXAMPLES:
     2526
     2527        We can take an elliptic curve and compute a histogram for it::
     2528
     2529            sage: E1 = EllipticCurve([4,7])
     2530            sage: E1.trace_histogram(100) # optional: smalljac
     2531            [0, 1, 2, 1, 3, 8, 2, 4, 2, 1]
     2532
     2533        In this case, the default choice of 10 buckets was made.
     2534        Instead, you can specify the number of bucket manually::
     2535
     2536            sage: E1.trace_histogram(100000, buckets=10) # optional: smalljac
     2537            [529, 849, 1059, 1193, 1201, 1153, 1133, 1092, 876, 505]
     2538
     2539        You can choose a range that starts anywhere, not just at 1::
     2540
     2541            sage: E1.trace_histogram([1000, 2000], buckets=10) # optional: smalljac
     2542            [11, 13, 12, 14, 14, 14, 17, 11, 20, 8]
     2543
     2544        This is generally useful for making histograms; a quick way to
     2545        do this is with ``bar_chart``, though custom graphing code
     2546        will look better. ::
     2547
     2548            sage: bar_chart(E1.trace_histogram(100000), width=1) #optional: smalljac
     2549        """
     2550
     2551        import sage.libs.smalljac.util as util
     2552        from math import sqrt
     2553
     2554        start, end = util.parse_limits(limit)
     2555
     2556        if buckets is None:
     2557            range_size = end - start + 1
     2558            buckets = int(sqrt(range_size))
     2559
     2560        return map(Integer,
     2561                   self._smalljac().sato_tate_buckets(buckets, start, end))
     2562
     2563    def moments(self, limits, moments=11):
     2564        """
     2565        Compute the moments of the distributions of the coefficients
     2566        of the L-polynomials of this curve, reduced at primes in a
     2567        range.
     2568
     2569        INPUT:
     2570
     2571        - ``limits`` - Either an integer, meaning all primes below that
     2572          integer are used, or a list of ``[start, end]``, in which
     2573          case all primes from ``start`` to ``end``, inclusive, are
     2574          used.
     2575
     2576        - ``moments`` - How many moments to return; fewer moments
     2577          might be faster to compute.  For now, at most 10 moments can
     2578          be computed.
     2579
     2580        OUTPUT:
     2581
     2582        A tuple whose first element is a list of moments of the
     2583        distribution of the $a_1$ coefficients of L-polynomials of
     2584        this curve, starting from the zeroth moment being 0.  The list
     2585        is as long as ``moments``; that is, if ``moments`` is 5, only
     2586        the zero through fourth moments are returned.
     2587
     2588        .. NOTE::
     2589            A tuple is returned for polymorphism with the
     2590            hyperelliptic curve method of the same name
     2591
     2592        EXAMPLES:
     2593
     2594        You can request the moments of an arbitrary curve::
     2595       
     2596            sage: ec = EllipticCurve('37b')
     2597            sage: ec.moments(100000)
     2598            ([0.0, -0.0030184706505541115, 0.9999566856391422,
     2599              -0.017815838483148088, 1.993903401871212,
     2600              -0.06113842272925202, 4.963581622954875,
     2601              -0.2068163797044211, 13.831883515902087,
     2602              -0.7065524181435547, 41.30699642905881],)
     2603
     2604        You can request fewer moments if you wish::
     2605
     2606            sage: ec.moments(10000, moments=3)
     2607            ([0.0, -0.0008126996517027156, 0.9890134939436451],)
     2608
     2609        Or compute in any range::
     2610
     2611            sage: ec.moments([1000, 2000], moments=5)
     2612            ([0.0, -0.008380482644508892, 0.9721022301822476,
     2613              -0.16345224489884563, 1.8856843129225394],)
     2614
     2615        Number fields are also supported::
     2616
     2617            sage: K.<a> = NumberField(x^2 - x - 1)
     2618            sage: ec2 = EllipticCurve(K, [a, -a+1, 0, -4, 3*a-5])
     2619            sage: ec2.moments([1000, 2000], moments=5)
     2620            ([0.0, -0.04767789672675515, 0.9514627444845026,
     2621              -0.20252616280144925, 1.7720971308862588],)
     2622
     2623        Rounding the moments should produce a recognizable sequence::
     2624
     2625            sage: map(round, ec.moments(100000)[0])
     2626            [0.0, -0.0, 1.0, -0.0, 2.0, -0.0, 5.0, -0.0, 14.0, -1.0, 41.0]
     2627        """
     2628
     2629        import sage.libs.smalljac.util as util
     2630
     2631        start, end = util.parse_limits(limits)
     2632        if moments > 11:
     2633            raise ValueError("Only up to 11 moments supported for now.")
     2634
     2635        return self._smalljac().moments(moments, start, end)
     2636
     2637    def guess_sato_tate_group(self, limit=None):
     2638        """
     2639        Provide a preliminary guess toward the Sato-Tate group
     2640        corresponding to this curve, by considering its reduction at
     2641        primes up to a certain limit.
     2642
     2643        INPUT:
     2644
     2645        - ``limit`` - The maximal prime considered.  The search may
     2646          terminate earlier if an answer is found.  Or, ``None``, to
     2647          use an implementation-defined (but very large) maximum.
     2648
     2649        OUTPUT:
     2650
     2651        A string containing the name of the corresponding Sato-Tate
     2652        group.  If no provisional identification can be made, returns
     2653        ``None``.
     2654
     2655        .. TODO:
     2656           - Should output a textual description of the Sato-Tate group
     2657
     2658        EXAMPLES:
     2659
     2660        Just create a hyperelliptic curve of genus 2 and call
     2661        ``guess_sato_tate_group`` to start it off::
     2662
     2663            sage: ec = EllipticCurve('11a')
     2664            sage: ec.guess_sato_tate_group(100000)
     2665            'SU(2)'
     2666            sage: ec = EllipticCurve([0, 1])
     2667            sage: ec.guess_sato_tate_group(100000)
     2668            'N(U(1))'
     2669            sage: K.<a> = NumberField(x^2 + 3)
     2670            sage: ec = EllipticCurve(K, [0, 1])
     2671            sage: ec.guess_sato_tate_group(100000)
     2672            'U(1)'
     2673
     2674        The above are the only three possibilities (for genus 1).
     2675        """
     2676
     2677        if limit is None:
     2678            limit = 0
     2679
     2680        return self._smalljac().guess_group(limit)
  • sage/schemes/elliptic_curves/ell_rational_field.py

    diff --git a/sage/schemes/elliptic_curves/ell_rational_field.py b/sage/schemes/elliptic_curves/ell_rational_field.py
    a b  
    746746            return self.__database_curve
    747747           
    748748
    749     def Np(self, p):
     749    def Np(self, p, algorithm="pari"):
    750750        r"""
    751         The number of points on `E` modulo `p`.
     751        The number of points on `E` modulo `p`.  The ``algorithm``
     752        argument is the same as that for ``ap``.
    752753
    753754        INPUT:
    754755
     
    757758
    758759        OUTPUT:
    759760
    760         (int) The number ofpoints on the reduction of `E` modulo `p`
     761        (int) The number of points on the reduction of `E` modulo `p`
    761762        (including the singular point when `p` is a prime of bad
    762763        reduction).
    763764       
     
    780781            1000000000000001426441464441649
    781782        """
    782783        if self.conductor() % p == 0:
    783             return p + 1 - self.ap(p)
    784         return p+1 - self.ap(p)
     784            return p + 1 - self.ap(p, algorithm=algorithm)
     785        return p+1 - self.ap(p, algorithm=algorithm)
    785786
    786787    #def __pari_double_prec(self):
    787788    #    EllipticCurve_number_field._EllipticCurve__pari_double_prec(self)
     
    880881    #  Etc.
    881882    ####################################################################
    882883
    883     def aplist(self, n, python_ints=False):
     884    def aplist(self, n, python_ints=False, algorithm="pari"):
    884885        r"""
    885886        The Fourier coefficients `a_p` of the modular form
    886887        attached to this elliptic curve, for all primes `p\leq n`.
     
    892893       
    893894        -  ``python_ints`` - bool (default: False); if True
    894895           return a list of Python ints instead of Sage integers.
     896
     897        -  ``algorithm`` - std (default: "pari")
     898           
     899           - "pari" - Use PARI's implementation of ellaplist
     900
     901           - "smalljac" - Use the smalljac library, if installed
    895902       
    896903       
    897904        OUTPUT: list of integers
     
    911918            <type 'sage.rings.integer.Integer'>
    912919            sage: type(e.aplist(13, python_ints=True)[0])
    913920            <type 'int'>
    914         """
    915         e = self.pari_mincurve()
    916         v = e.ellaplist(n, python_ints=True)
    917         if python_ints:
     921
     922        We can pass ``algorithm="smalljac"`` to use the smalljac
     923        library; this is generally much faster.
     924
     925            sage: e.aplist(13, algorithm="smalljac") # optional: smalljac
     926            [-2, -3, -2, -1, -5, -2]
     927        """
     928       
     929        if algorithm == "pari":
     930            e = self.pari_mincurve()
     931            v = e.ellaplist(n, python_ints=True)
     932            if python_ints:
     933                return v
     934            else:
     935                return [Integer(a) for a in v]
     936        elif algorithm == "smalljac":
     937            proxy = self._smalljac()
     938            aps = proxy.ap(0, n)
     939           
     940            v = []
     941            for p, ap in aps:
     942                if ap is None:
     943                    # Sage already deals with this case well enough
     944                    ap = self.ap(p, algorithm="pari")
     945                    if not python_ints: ap = int(ap)
     946                    v.append(ap)
     947                else:
     948                    if python_ints: ap = Integer(ap)
     949                    v.append(ap)
     950           
    918951            return v
     952        elif algorithm == "all":
     953            aplist1 = self.aplist(n, python_ints=python_ints, algorithm="pari")
     954            aplist2 = self.aplist(n, python_ints=python_ints, algorithm="smalljac")
     955            if aplist1 != aplist2:
     956                raise ArithmeticError("PARI and smalljac disagree (%s, %s)", aplist1, aplist2)
     957            return aplist1
    919958        else:
    920             return [Integer(a) for a in v]
    921        
    922        
     959            raise RuntimeError, "algorithm '%s' is not known."%algorithm
    923960
    924961    def anlist(self, n, python_ints=False):
    925962        """
     
    13731410        else:
    13741411            raise ValueError, "algorithm %s not defined"%algorithm
    13751412
     1413    def _smalljac(self):
     1414        r"""
     1415        Computes and caches the smalljac representation of the given
     1416        curve.  For internal use only.
     1417        """
     1418
     1419        try:
     1420            return self._smalljac_proxy
     1421        except AttributeError:
     1422            try:
     1423                import sage.libs.smalljac.all as smalljac
     1424            except ImportError:
     1425                raise RuntimeError, "smalljac not installed or inaccessible"
     1426           
     1427            f = "[%s,%s,%s,%s,%s]" % tuple(self.ainvs())
     1428            self._smalljac_proxy = smalljac.SmallJac(f)
     1429            return self._smalljac_proxy
    13761430
    13771431    def simon_two_descent(self, verbose=0, lim1=5, lim3=50, limtriv=10, maxprob=20, limbigprime=30):
    13781432        r"""
     
    24962550        """
    24972551        return Integer(self.pari_mincurve().ellak(n))
    24982552
    2499     def ap(self, p):
     2553    def ap(self, p, algorithm="pari"):
    25002554        """
    25012555        The p-th Fourier coefficient of the modular form corresponding to
    25022556        this elliptic curve, where p is prime.
    25032557       
     2558        -  ``algorithm`` - std (default: "pari")
     2559           
     2560           - "pari" - Use PARI's implementation of ellaplist
     2561
     2562           - "smalljac" - Use the smalljac library, if installed
     2563       
     2564       
    25042565        EXAMPLES::
    25052566       
    25062567            sage: E=EllipticCurve('37a1')
    25072568            sage: [E.ap(p) for p in prime_range(50)]
    25082569            [-2, -3, -2, -1, -5, -2, 0, 0, 2, 6, -4, -1, -9, 2, -9]
    2509         """
     2570
     2571        If the ``smalljac`` libary is installed, we can use it::
     2572
     2573            sage: [E.ap(p, algorithm="smalljac") for p in prime_range(50)]
     2574            [-2, -3, -2, -1, -5, -2, 0, 0, 2, 6, -4, -1, -9, 2, -9]
     2575        """
     2576       
    25102577        if not arith.is_prime(p):
    25112578            raise ArithmeticError, "p must be prime"
    2512         return Integer(self.pari_mincurve().ellap(p))
     2579       
     2580        if algorithm == "pari":
     2581            return Integer(self.pari_mincurve().ellap(p))
     2582        elif algorithm == "smalljac":
     2583            proxy = self._smalljac()
     2584            ap = proxy.ap(p)
     2585            if ap is None:
     2586                # Let PARI handle the icky case
     2587                return self.ap(p, algorithm="pari")
     2588            else:
     2589                return Integer(ap)
     2590        elif algorithm == "all":
     2591            pari_results     = self.ap(p, algorithm="pari")
     2592            smalljac_results = self.ap(p, algorithm="smalljac")
     2593            assert pari_results == smalljac_results, "An error has occured, and SmallJac's and PARI's results are not the same"
     2594            return smalljac_results
    25132595
    25142596    def quadratic_twist(self, D):
    25152597        """
     
    46834765        """
    46844766        return self.conductor().is_squarefree()
    46854767
    4686     def is_ordinary(self, p, ell=None):
     4768    def is_ordinary(self, p, ell=None, algorithm="pari"):
    46874769        """
    46884770        Return True precisely when the mod-p representation attached to
    4689         this elliptic curve is ordinary at ell.
     4771        this elliptic curve is ordinary at ell.  The ``algorithm``
     4772        argument has the same values as that for ``aplist``.
    46904773       
    46914774        INPUT:
    46924775       
     
    47084791        """
    47094792        if ell is None:
    47104793            ell = p
    4711         return self.ap(ell) % p != 0
     4794        return self.ap(ell, algorithm=algorithm) % p != 0
    47124795
    47134796    def is_good(self, p, check=True):
    47144797        """
     
    47684851            ell = p
    47694852        return self.is_good(p) and not self.is_ordinary(p, ell)
    47704853
    4771     def supersingular_primes(self, B):
     4854    def supersingular_primes(self, B, algorithm="pari"):
    47724855        """
    47734856        Return a list of all supersingular primes for this elliptic curve
    4774         up to and possibly including B.
     4857        up to and possibly including B.  The ``algorithm`` argument has
     4858        the same values as that for ``aplist``.
    47754859       
    47764860        EXAMPLES::
    47774861       
     
    47974881            sage: e.supersingular_primes(1)
    47984882            []
    47994883        """
    4800         v = self.aplist(max(B, 3))
     4884        v = self.aplist(max(B, 3), algorithm=algorithm)
    48014885        P = rings.prime_range(max(B,3)+1)
    48024886        N = self.conductor()
    48034887        return [P[i] for i in [0,1] if P[i] <= B and v[i]%P[i]==0 and N%P[i] != 0] + \
    48044888                      [P[i] for i in range(2,len(v)) if v[i] == 0 and N%P[i] != 0]
    48054889
    4806     def ordinary_primes(self, B):
     4890    def ordinary_primes(self, B, algorithm="pari"):
    48074891        """
    48084892        Return a list of all ordinary primes for this elliptic curve up to
    4809         and possibly including B.
     4893        and possibly including B.  The ``algorithm`` argument has the same
     4894        values as that for ``aplist``.
    48104895       
    48114896        EXAMPLES::
    48124897       
     
    48294914            sage: e.ordinary_primes(1)
    48304915            []
    48314916        """
    4832         v = self.aplist(max(B, 3) )
     4917        v = self.aplist(max(B, 3), algorithm=algorithm)
    48334918        P = rings.prime_range(max(B,3) +1)
    48344919        return [P[i] for i in [0,1] if P[i] <= B and v[i]%P[i]!=0] +\
    48354920               [P[i] for i in range(2,len(v)) if v[i] != 0]
  • sage/schemes/hyperelliptic_curves/hyperelliptic_g2_finite_field.py

    diff --git a/sage/schemes/hyperelliptic_curves/hyperelliptic_g2_finite_field.py b/sage/schemes/hyperelliptic_curves/hyperelliptic_g2_finite_field.py
    a b  
    1313class HyperellipticCurve_g2_finite_field(
    1414    hyperelliptic_g2_generic.HyperellipticCurve_g2_generic,
    1515    hyperelliptic_finite_field.HyperellipticCurve_finite_field):
    16     pass
    1716
     17    def _smalljac(self):
     18        r"""
     19        Computes and caches the smalljac representation of the given
     20        curve.  For internal use only.
     21        """
     22
     23        try:
     24            return self._smalljac_proxy
     25        except AttributeError:
     26            try:
     27                import sage.libs.smalljac.all as smalljac
     28            except ImportError:
     29                raise RuntimeError, "smalljac not installed or inaccessible"
     30           
     31            f, h = self.hyperelliptic_polynomials()
     32            assert h == 0, "SmallJac can only does with hyperelliptic polynomials in Weirstrass form"
     33            s = str(f)
     34            self._smalljac_proxy = smalljac.SmallJac(s)
     35            return self._smalljac_proxy
     36
     37    def jacobian_structure(self, algorithm="smalljac"):
     38        r"""
     39        Returns the group structure of the group of points on this
     40        elliptic curve.  Unlike ``abelian_group``, it does not compute
     41        the generators of the group, just the structure, and is thus
     42        perhaps somewhat less useful.
     43
     44        INPUT:
     45       
     46        - ``algorithm`` - so far, only ``"smalljac"`` is supported,
     47          which uses Drew Sutherland's SmallJac library for the
     48          computation.  SmallJac can only be used on degree-1 fields.
     49
     50        OUTPUT:
     51
     52        An ``AdditiveAbelianGroup`` representing the group structure
     53        of the jacobian of the curve.
     54
     55        EXAMPLES:
     56
     57        We can take an arbitrary curve and ask for the jacobian's
     58        structure::
     59
     60            sage: R.<x> = GF(next_prime(10^5))[]
     61            sage: hec = HyperellipticCurve(x^5 + 17*x^4 - 318*x^3 + 2*x - 17)
     62            sage: hec.jacobian_structure() # optional: smalljac
     63            Additive abelian group isomorphic to Z/2 + Z/4974328512
     64        """
     65       
     66        from sage.groups.additive_abelian.additive_abelian_group import AdditiveAbelianGroup
     67
     68        if algorithm == "smalljac":
     69            K = self.base_ring()
     70            assert K.degree() == 1, "Smalljac can only find group structures in degree-1 fields"
     71            p = K.cardinality()
     72            proxy = self._smalljac()
     73            structure = proxy.group(p)
     74            return AdditiveAbelianGroup(structure)
     75        else:
     76            raise ValueError("Unknown algorithm `%s`" % algorithm)
  • sage/schemes/hyperelliptic_curves/hyperelliptic_g2_rational_field.py

    diff --git a/sage/schemes/hyperelliptic_curves/hyperelliptic_g2_rational_field.py b/sage/schemes/hyperelliptic_curves/hyperelliptic_g2_rational_field.py
    a b  
    99#*****************************************************************************
    1010
    1111import hyperelliptic_g2_generic, hyperelliptic_rational_field
     12from sage.rings.integer import Integer
    1213
    1314class HyperellipticCurve_g2_rational_field(
    1415    hyperelliptic_g2_generic.HyperellipticCurve_g2_generic,
    1516    hyperelliptic_rational_field.HyperellipticCurve_rational_field):
    16     pass
     17
     18    def _smalljac(self):
     19        r"""
     20        Computes and caches the smalljac representation of the given
     21        curve.  For internal use only.
     22        """
     23
     24        try:
     25            return self._smalljac_proxy
     26        except AttributeError:
     27            try:
     28                import sage.libs.smalljac.all as smalljac
     29            except ImportError:
     30                raise RuntimeError, "smalljac not installed or inaccessible"
     31           
     32            f, h = self.hyperelliptic_polynomials()
     33            assert h == 0, "SmallJac can only does with hyperelliptic polynomials in Weirstrass form"
     34            s = str(f)
     35            self._smalljac_proxy = smalljac.SmallJac(s)
     36            return self._smalljac_proxy
     37
     38    def aplists(self, limit, algorithm="smalljac", container=list, index="primes", raw=False):
     39        """
     40        Computes the traces of Frobenius at the prime ideals of the
     41        base number field, sorted by the norm of these ideals.  So
     42        far, this method requires the Smalljac library.
     43
     44        INPUT:
     45
     46        - ``limit`` - either an integer, in which case the prime
     47          ideals will have norm less than this limit, or [a, b], in
     48          which case the prime ideals with norms between ``a`` and
     49          ``b`, inclusive, are used.
     50
     51        - ``algorithm`` - so far, only "smalljac" is supported, which
     52          uses Drew Sutherland's smalljac library for elliptic curve
     53          structure computation.
     54
     55        - The remainder (``container``, ``index``, and ``raw``) define
     56          the format of the result; see ``EllipticCurve.aplists`` for
     57          a description.
     58
     59        OUTPUT:
     60
     61        A list of traces of Frobenius at various primes; the format is
     62        precisely specified by the ``container``, ``index``, and
     63        ``raw`` arguments.  If a curve has a bad reduction at a prime,
     64        ``None`` is instead returned.
     65
     66        EXAMPLES:
     67
     68        We can find the trace of Frobenius at a given prime:
     69
     70            sage: R.<x> = QQ[]
     71            sage: hec = HyperellipticCurve(x^5 + 2*x^4 - x^3 - 3*x^2 - x)
     72            sage: hec.aplists([17, 17]) # optional: smalljac
     73            [(17, [-9])]
     74
     75        Or in a range:
     76
     77            sage: hec.aplists(100) # optional: smalljac
     78            [(2, [None]), (3, [-1]), (5, [3]), (7, [None]), (11, [-3]), (13, [0]),
     79            (17, [-9]), (19, [-7]), (23, [15]), (29, [-12]), (31, [-5]), (37, [5]),
     80            (41, [0]), (43, [0]), (47, [-3]), (53, [9]), (59, [-9]), (61, [15]),
     81            (67, [9]), (71, [0]), (73, [3]), (79, [-9]), (83, [24]), (89, [-21]),
     82            (97, [0])]
     83            sage: hec.aplists([50, 100]) # optional: smalljac
     84            [(53, [9]), (59, [-9]), (61, [15]), (67, [9]), (71, [0]), (73, [3]),
     85            (79, [-9]), (83, [24]), (89, [-21]), (97, [0])]
     86
     87        The raw results are suitable for large-scale statistics:
     88
     89            sage: __builtins__.sum(ap for p, ap in hec.aplists(10^4, raw=True) if ap) # optional: smalljac
     90            1900
     91            sage: __builtins__.sum(ap for p, ap in hec.aplists(2 * 10^4, raw=True) if ap) # optional: smalljac
     92            1925
     93            sage: __builtins__.sum(ap for p, ap in hec.aplists(3 * 10^4, raw=True) if ap) # optional: smalljac
     94            5451
     95        """
     96
     97        import sage.libs.smalljac.util as util
     98
     99        start, end = util.parse_limits(limit)
     100
     101        if algorithm == "smalljac":
     102            proxy = self._smalljac()
     103            res = proxy.ap(start, end)
     104
     105            def make_int(p, v):
     106                return Integer(v)
     107
     108            return util.interpret_smalljac_results(
     109                res, correction_fn=None, construct_fn=make_int,
     110                raw=raw, index=index, container=container)
     111        else:
     112            raise ValueError("Algorithm %s not defined" % algorithm)
     113
     114    def aps(self, p):
     115        """
     116        Returns a list of traces of Frobenius when this elliptic curve
     117        is reduced at prime ideals of norm ``p``.
     118
     119        INPUT:
     120
     121         - ``p`` - The integer prime at which to reduce.
     122
     123        OUTPUT:
     124
     125        A list of traces of Frobenius, or None if there was a bad
     126        reduction.
     127
     128        EXAMPLES:
     129
     130        We can take a curve and find arbitrary traces::
     131
     132            sage: R.<x> = QQ[]
     133            sage: hec = HyperellipticCurve(x^5 + 2*x^4 - x^3 - 3*x^2 - x)
     134            sage: hec.aps(17) # optional: smalljac
     135            [-9]
     136            sage: hec.aps(next_prime(10^6)) # optional: smalljac
     137            [909]
     138            sage: hec.aps(next_prime(10^8)) # optional: smalljac
     139            [27651]
     140        """
     141
     142        res = self.aplists([p, p])
     143
     144        return res[0][1] # [0] extracts the prime we care about, [1] extracts the list of traces
     145
     146    def grouplists(self, limit, algorithm="smalljac", container=list, index="primes", raw=False):
     147        r"""
     148        Computes the group structure of reductions at the prime ideals
     149        of the base number field, sorted by the norm of these ideals.
     150        So far, this method requires the Smalljac library.
     151
     152        INPUT:
     153
     154        - ``limit`` - either an integer, in which case the prime
     155          ideals will have norm less than this limit, or [a, b], in
     156          which case the prime ideals with norms between ``a`` and
     157          ``b`, inclusive, are used.
     158
     159        - ``algorithm`` - so far, only "smalljac" is supported, which
     160          uses Drew Sutherland's smalljac library for elliptic curve
     161          structure computation.
     162
     163        - The remainder (``container``, ``index``, and ``raw``) define
     164          the format of the result; see ``EllipticCurve.aplists`` for
     165          a description.
     166
     167        OUTPUT:
     168
     169        A list of group structures of Jacobians at various primes; the
     170        format is precisely specified by the ``container``, ``index``,
     171        and ``raw`` arguments.  If a curve has a bad reduction at a prime,
     172        ``None`` is instead returned.
     173
     174        EXAMPLES:
     175
     176        You can construct a hyperelliptic curve and ask for its
     177        Jacobian's group structure at a certain prime::
     178
     179            sage: R.<x> = QQ[]
     180            sage: hec = HyperellipticCurve(x^5 + 2*x^4 - x^3 - 3*x^2 - x)
     181            sage: hec.grouplists([17, 17]) # optional: smalljac
     182            [(17, [Additive abelian group isomorphic to Z/4 + Z/124])]
     183
     184        Or in a range::
     185
     186            sage: hec.grouplists([80, 100]) # optional: smalljac
     187            [(83, [Additive abelian group isomorphic to Z/2 + Z/2 + Z/36 + Z/36]),
     188             (89, [Additive abelian group isomorphic to Z/8 + Z/1256]),
     189             (97, [Additive abelian group isomorphic to Z/2 + Z/2 + Z/2 + Z/1158])]
     190           
     191        The raw results are suitable for larger-scale statistics::
     192
     193            sage: def jac_checksum(poly):                                         
     194            ...       hec = HyperellipticCurve(poly)
     195            ...       sum = __builtins__.sum                           
     196            ...       return sum(sum(gs) for p, gs in hec.grouplists([3, 10000], raw=True) if gs)
     197            sage: jac_checksum(x^5-x) # optional: smalljac
     198            2083547730
     199            sage: jac_checksum(x^5+5*x^3+5*x) # optional: smalljac
     200            18393758619
     201            sage: jac_checksum(x^5+1) # optional: smalljac
     202            19340969888
     203            sage: jac_checksum(x^5+x^3+2*x) # optional: smalljac
     204            12945578110
     205            sage: jac_checksum(x^5-10*x^2-9*x+2) # optional: smalljac
     206            12917177585
     207            sage: jac_checksum(x^5+3*x^4+4*x^2-9*x+1) # optional: smalljac
     208            12614952542
     209            sage: jac_checksum(x^5-x+1) # optional: smalljac
     210            28462483128
     211        """
     212
     213        from sage.groups.additive_abelian.additive_abelian_group import AdditiveAbelianGroup
     214        import sage.libs.smalljac.util as util
     215
     216        start, end = util.parse_limits(limit)
     217
     218        if algorithm == "smalljac":
     219            proxy = self._smalljac()
     220            res = proxy.group(start, end)
     221
     222            def make_group(p, v):
     223                return AdditiveAbelianGroup(v)
     224
     225            return util.interpret_smalljac_results(
     226                res, correction_fn=None, construct_fn=make_group,
     227                raw=raw, index=index, container=container)
     228        else:
     229            raise ValueError("Algorithm %s not defined" % algorithm)
     230
     231    def lpolylists(self, limit, normalize=False, algorithm="smalljac", container=list, index="primes", raw=False):
     232        r"""
     233        Computes the L-polynomials of reductions at the prime ideals
     234        of the base number field, sorted by the norm of these ideals.
     235        So far, this method requires the Smalljac library.
     236
     237        INPUT:
     238
     239        - ``limit`` - either an integer, in which case the prime
     240          ideals will have norm less than this limit, or ``[a, b]``,
     241          in which case the prime ideals with norms between ``a`` and
     242          ``b`, inclusive, are used.
     243
     244        - ``normalize`` - whether to return normalized L-polynomials.
     245          By default, non-normalized L-polynomials are returned.
     246
     247        - ``algorithm`` - so far, only "smalljac" is supported, which
     248          uses Drew Sutherland's smalljac library for elliptic curve
     249          structure computation.
     250
     251        - The remainder (``container``, ``index``, and ``raw``) define
     252          the format of the result; see ``EllipticCurve.aplists`` for
     253          a description.
     254
     255        OUTPUT:
     256
     257        A list of L-polynomials of reductions at various primes; the
     258        format is precisely specified by the ``container``, ``index``,
     259        and ``raw`` arguments.  If a curve has a bad reduction at a
     260        prime, ``None`` is instead returned.
     261
     262        EXAMPLES:
     263
     264        You can construct a hyperelliptic curve and ask for its
     265        L-polynomials at a certain prime:
     266
     267            sage: R2.<y> = QQ[]
     268            sage: E = HyperellipticCurve(y^5 + 2*y^4 - y^3 - 3*y^2 - y)
     269            sage: E.lpolylists([89, 89])
     270            [(89, [7921*T^4 + 1869*T^3 + 236*T^2 + 21*T + 1])]
     271
     272        Or in a range:
     273
     274            sage: E.lpolylists([80, 90])
     275            [(83, [6889*T^4 - 1992*T^3 + 310*T^2 - 24*T + 1]),
     276             (89, [7921*T^4 + 1869*T^3 + 236*T^2 + 21*T + 1])]
     277
     278        You can also get normalized L-polynomials:
     279
     280            sage: E.lpolylists([80, 100], normalized=True)
     281            [(83, [T^4 - 2.63434223975257*T^3 + 3.00000000000000*T^2 -
     282                   2.63434223975257*T + 1.00000000000000]),
     283             (89, [T^4 + 2.22599554801336*T^3 + 2.00000000000000*T^2 +
     284                   2.22599554801336*T + 1.00000000000000]),
     285             (97, [T^4 - 2.00000000000000*T^2 + 1.00000000000000])]
     286        """
     287
     288        from sage.rings.polynomial.polynomial_ring import PolynomialRing_field
     289        from sage.rings.real_mpfr import RR
     290        from sage.rings.rational_field import QQ
     291        from math import sqrt
     292        import sage.libs.smalljac.util as util
     293
     294        start, end = util.parse_limits(limit)
     295
     296        if algorithm == "smalljac":
     297            proxy = self._smalljac()
     298            res = proxy.lpoly(start, end)
     299
     300            if normalize:
     301                polynomial_ring = PolynomialRing_field(RR, 'T')
     302            else:
     303                polynomial_ring = PolynomialRing_field(QQ, 'T')
     304            t = polynomial_ring.gens()[0] # I don't want a T
     305           
     306            def make_lpoly(p, coeffs):
     307                if not normalize:
     308                    a1, a2 = coeffs
     309                    return p**2*t**4 + a1*p*t**3 + a2*t**2 + a1*t + 1
     310                else:
     311                    a1, a2 = coeffs
     312                    a2 /= p
     313                    a1 /= sqrt(p)
     314                    return t**4 + a1*t**3 + a2*t**2 + a1*t + 1
     315
     316            return util.interpret_smalljac_results(
     317                res, correction_fn=None, construct_fn=make_lpoly,
     318                raw=raw, index=index, container=container)
     319        else:
     320            raise ValueError("Algorithm %s not defined" % algorithm)
     321
     322    def trace_histogram(self, limit, buckets=None):
     323        """
     324        Computes the traces of reductions of this curve at all primes
     325        in a range, and returns how many, normalized, fell in each
     326        of ``buckets`` buckets from -2 to 2.
     327
     328        INPUT:
     329
     330        - ``limits`` - Either an integer, meaning all primes below that
     331          integer are used, or a list of ``[start, end]``, in which
     332          case all primes from ``start`` to ``end``, inclusive, are
     333          used.
     334
     335        - ``buckets`` - How many buckets to split the range $[-2, 2]$
     336          into.  By default, the number of buckets is the square root
     337          of the range size.
     338
     339        OUTPUT:
     340
     341        A list of integers, as many as the ``buckets`` parameter,
     342        denoting the number of trace values fell into that bucket for
     343        primes in the given range.
     344
     345        EXAMPLES:
     346
     347        We can take a hyperelliptic curve and compute a histogram for it::
     348
     349            sage: R.<x> = QQ[]
     350            sage: E3 = HyperellipticCurve(x^5 + 2*x^4 - x^3 - 3*x^2 - x)
     351            sage: E3.trace_histogram(100) # optional: smalljac
     352            [0L, 0L, 4L, 4L, 2L, 6L, 4L, 1L, 2L, 0L]
     353
     354        In this case, the default choice of 10 buckets was made.
     355        Instead, you can specify the number of bucket manually::
     356
     357            sage: E3.trace_histogram(100000, buckets=10) # optional: smalljac
     358            [126L, 401L, 722L, 1317L, 1420L, 3098L, 1255L, 722L, 403L, 126L]
     359
     360        You can choose a range that starts anywhere, not just at 1::
     361
     362            sage: E3.trace_histogram([1000, 2000], buckets=10) # optional: smalljac
     363            [6L, 7L, 5L, 17L, 16L, 45L, 23L, 5L, 8L, 3L]
     364
     365        This is generally useful for making histograms; a quick way to
     366        do this is with ``bar_chart``, though custom graphing code
     367        will look better. ::
     368
     369            sage: bar_chart(E3.trace_histogram(100000), width=1)
     370        """
     371
     372        import sage.libs.smalljac.util as util
     373        from math import sqrt
     374
     375        start, end = util.parse_limits(limit)
     376
     377        if buckets is None:
     378            range_size = end - start + 1
     379            buckets = int(sqrt(range_size))
     380
     381        return map(Integer,
     382                   self._smalljac().sato_tate_buckets(buckets, start, end))
     383
     384    def moments(self, limits, moments=11):
     385        """
     386        Compute the moments of the distributions of the coefficients
     387        of the L-polynomials of this curve, reduced at primes in a
     388        range.
     389
     390        INPUT:
     391
     392        - ``limits`` - Either an integer, meaning all primes below that
     393          integer are used, or a list of ``[start, end]``, in which
     394          case all primes from ``start`` to ``end``, inclusive, are
     395          used.
     396
     397        - ``moments`` - How many moments to return; fewer moments
     398          might be faster to compute.  For now, at most 10 moments can
     399          be computed.
     400
     401        OUTPUT:
     402
     403        A tuple whose first element is a list of moments of the
     404        distribution of the $a_1$ coefficients of L-polynomials of
     405        this curve, starting from the zeroth moment being 0, and whose
     406        second element is a similar list for the $a_2$ coefficients.
     407        The list is as long as ``moments``; that is, if ``moments`` is
     408        5, only the zero through fourth moments are returned.
     409
     410        EXAMPLES:
     411
     412        You can request the moments of an arbitrary curve::
     413       
     414            sage: R.<x> = QQ[]
     415            sage: he = HyperellipticCurve(x^5 + 2*x^4 - x^3 - 3*x^2 - x)
     416            sage: he.moments(10000)
     417            ([0.0, -0.01280967779130451, 2.031860103335779,
     418              -0.022492869661964253, 12.3326503293039,
     419              -0.08293711910285151, 104.79583112542454,
     420              1.4982847405160258, 1046.380403999995, 58.119822353379035,
     421              11453.36487528556],
     422             [0.0, 1.0190697892939964, 4.01328379603583,
     423             11.252490719416176, 45.189732618049945,
     424             179.65488417712353, 796.9073238931069,
     425             3636.7510997563822, 17306.793230016287, 84297.7319785892,
     426             419658.09474476264])
     427
     428        You can request fewer moments if you wish::
     429
     430            sage: he.moments(10000, moments=3)
     431            ([0.0, -0.01280967779130451, 2.031860103335779],
     432             [0.0, 1.0190697892939964, 4.01328379603583])
     433
     434        Or compute in any range::
     435
     436            sage: he.moments([1000, 2000], moments=5)
     437            ([0.0, 0.012019115191429827, 2.360758977244878,
     438              0.7305838739637147, 17.03361483200747],
     439             [0.0, 1.0981854467458285, 4.588416007826705,
     440              14.813429393106606, 62.82293520798331])
     441        """
     442
     443        import sage.libs.smalljac.util as util
     444
     445        start, end = util.parse_limits(limits)
     446        if moments > 11:
     447            raise ValueError("Only up to 11 moments supported for now.")
     448
     449        return self._smalljac().moments(moments, start, end)
     450
     451    def guess_sato_tate_group(self, limit=None):
     452        """
     453        Provide a preliminary guess toward the Sato-Tate group
     454        corresponding to this curve, by considering its reduction at
     455        primes up to a certain limit.
     456
     457        INPUT:
     458
     459         - ``limit`` - The maximal prime considered.  The search may
     460           terminate earlier if an answer is found.  Or, ``None``, to
     461           use an implementation-defined (but very large) maximum.
     462
     463        OUTPUT:
     464
     465        A string containing the name of the corresponding Sato-Tate
     466        group.  To interpret these, a copy of "Sato-Tate distributions
     467        and Galois endomorphism modules in genus 2", [fkrs]_, may be
     468        helpful.  If no provisional identification can be made,
     469        returns ``None``
     470
     471        .. TODO:
     472           - Should output a textual description of the Sato-Tate group
     473
     474        EXAMPLES:
     475
     476        Just create a hyperelliptic curve of genus 2 and call
     477        ``guess_sato_tate_group`` to start it off::
     478
     479            sage: R.<x> = QQ[]
     480            sage: he = HyperellipticCurve(x^5 + 2*x^4 - x^3 - 3*x^2 - x)
     481            sage: he.guess_sato_tate_group(100000) # long time - 4 seconds
     482            'E_6'
     483            sage: he = HyperellipticCurve(x^5 + x)
     484            sage: he.guess_sato_tate_group(100000)
     485            'D_{2,1}'
     486            sage: he = HyperellipticCurve(x^5 + 2*x)
     487            sage: he.guess_sato_tate_group(100000)
     488            'D_{4,1}'
     489        """
     490
     491        if limit is None:
     492            limit = 0
     493
     494        return self._smalljac().guess_group(limit)
  • setup.py

    diff --git a/setup.py b/setup.py
    a b  
    935935                     'sage.libs.cremona',
    936936                     'sage.libs.mpmath',
    937937                     'sage.libs.lcalc',
     938                     'sage.libs.smalljac',
    938939
    939940                     'sage.logic',
    940941