Ticket #13376: smalljac-4.0.23.patch

File smalljac-4.0.23.patch, 84.4 KB (added by pavpanchekha, 9 years ago)

Patch to Sage 5.6

  • module_list.py 

    # HG changeset patch
# User Pavel Panchekha <me@pavpanchekha.com>
# Date 1362951771 14400
# Node ID 3259d69fe5bcd65e01ed73bf8215c091b3f14a74
# Parent  ec1fb07db6e23f9fbd4c34f0d48198d08ec76473
[mq]: current

diff --git a/module_list.py b/module_list.py
    a b  
    19921992                  libraries = ["lrcalc"],
    19931993                  depends = [SAGE_LOCAL + "/include/lrcalc/symfcn.h"]), # should include all .h
    19941994        )
     1995
     1996if is_package_installed('smalljac'):
     1997    ext_modules.append(
     1998        Extension('sage.libs.smalljac.wrapper2',
     1999            sources = ["sage/libs/smalljac/wrapper2.pyx"],
     2000            libraries = ["gmp", "m", "smalljac", "ff_poly"])
     2001        )

  • 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
    - +  
     1"""
     2A wrapper module to import the SmallJac class from the Cython Smalljac interface
     3"""
     4
     5from 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    cdef int SMALLJAC_FILENOTFOUND
     53    cdef int SMALLJAC_BADFILE
     54    cdef int SMALLJAC_NODATA
     55    cdef int SMALLJAC_NOT_OVER_Q
     56
     57    cdef int SMALLJAC_CURVE_STRING_LEN
     58
     59    int smalljac_moments (smalljac_curve_t curve,
     60                          unsigned long start, unsigned long end,
     61                          double moments[], int n, int m,
     62                          char STgroup[16],
     63                          int (*filter_callback)(smalljac_curve_t curve,
     64                                                 unsigned long q,
     65                                                 void *arg),
     66                          void *arg)
     67    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    """
     5    Parse a bounds specification into a start and end of range.  Both
     6    the start and end are inclusive.
     7
     8    EXAMPLES:
     9
     10    Integers represent the range from zero to themselves, inclusive::
     11
     12        >>> parse_limits(7)
     13        (0, 7)
     14        >>> parse_limits(5)
     15        (0, 5)
     16   
     17    Lists represent inclusive ranges::
     18
     19        >>> parse_limits([2, 17])
     20        (2, 17)
     21        >>> parse_limits([10**5, 12**5])
     22        (100000, 248832)
     23
     24    Tuples represent exclusive ranges::
     25
     26        >>> parse_limits((2, 17))
     27        (3, 16)
     28        >>> parse_limits((10**5, 12**5))
     29        (100001, 248831)
     30   
     31    Slices represent standard Python left-inclusive/right-exclusive ranges::
     32
     33        >>> parse_limits(slice(2, 17))
     34        (2, 16)
     35        >>> parse_limits(slice(10**5, 12**5))
     36        (100000, 248831)
     37    """
     38   
     39    start, end = 0, 0
     40    if isinstance(limit, int) or isinstance(limit, Integer):
     41        start = 0
     42        end = int(limit)
     43    elif isinstance(limit, list):
     44        start = limit[0]
     45        end = limit[-1]
     46    elif isinstance(limit, tuple):
     47        start = limit[0] + 1
     48        end = limit[-1] - 1
     49    elif isinstance(limit, slice):
     50        start = limit.start
     51        end = limit.stop - 1
     52    else:
     53        raise ValueError("Cannot interpret limits for trace computation", limit)
     54
     55    if start < 0 or end < 0 or end < start:
     56        raise ValueError("Invalid limits; both must be nonnegative, and end must be not less than start")
     57
     58    return start, end
     59
     60def interpret_smalljac_results(res, correction_fn, construct_fn, raw=False, index="primes", container=list):
     61    """
     62    Interpret results returned by Smalljac into a more usable form.
     63    Results from Smalljac are a list of tuples of prime and value.
     64    Returns a container containing tuples of an index and a values for
     65    that index.
     66
     67    INPUT:
     68
     69    - ``res`` - The results returned by Smalljac
     70    - ``correction_fn`` - A function that takes two inputs: a prime
     71      ``p`` and a list of values at that prime.  Returns a "corrected"
     72      version, where ``None`` values have been replaced.  Used when Sage
     73      already implements fallback routines, such as for trace of Frobenius
     74      calculations at primes that produce bad reductions.
     75    - ``construct_fn`` - A function that takes a computed value and
     76      returns a value to be returned to the user.  This might convert
     77      Python ``int`` to Sage ``Integer`` objects, or compute the
     78      correct format for L-polynomials, or what have you.
     79    - ``raw`` - If true, the original Smalljac results are immediately
     80      returned, with no additional processing.
     81    - ``index`` - A string, either ``"primes"`` or ``"ideals"``.  If
     82      ``"primes"``, the indexes of the returned pairs are primes (as
     83      Sage ``Integer`` objects) and the values are the list of values
     84      for that prime.  If ``"ideals"``, the indexes are a pair of
     85      prime (as a Sage ``Integer``) and a counter value (also as a
     86      Sage ``Integer``), and the values are individual values (in
     87      other words, passing ``"ideals"`` splits up the list associated
     88      with each prime into individual elements of the overall
     89      container.
     90    - ``container`` - The container all these pairs are passed to.
     91      Popular values are ``list`` (the default), in which case you get
     92      back a sorted list of these pairs, and ``dict``, in which case
     93      you get an unsorted hash table of such pairs.  In theory, can be
     94      any function that takes a single list as input.
     95
     96    EXAMPLES:
     97
     98      >>> res = [(2, [5, 1]), (11, None), (13, [10, 7]), (19, [20, 2])]
     99      >>> interpret_smalljac_results(res, lambda p, vals: [2], lambda p, val: str(x))
     100      [(2, ['[5, 1]']), (11, [2]), (13, ['[10, 7]']), (19, ['[20, 2]'])]
     101      >>> interpret_smalljac_results(res, lambda p, vals: [4], lambda p, val: str(x))
     102      [(2, ['[5, 1]']), (11, [4]), (13, ['[10, 7]']), (19, ['[20, 2]'])]
     103      >>> interpret_smalljac_results(res, lambda p, vals: [2], lambda p, val: val[0])
     104      [(2, [5]), (11, [2]), (13, [10]), (19, [20])]
     105      >>> res == interpret_smalljac_results(res, lambda p,v:[2], lambda p,v: v, raw=True)
     106      True
     107      >>> interpret_smalljac_results(res, lambda p,v:[2], lambda p,v: v, index="ideals")
     108      [((2, 0), 5), ((11, 0), 2), ((13, 0), 10), ((19, 0), 20)]
     109      >>> interpret_smalljac_results(res, lambda p,v:[2], lambda p,v: v, container=dict)
     110      {19: [[20, 2]], 2: [[5, 1]], 11: [2], 13: [[10, 7]]}
     111      >>> interpret_smalljac_results(res, lambda p,v:[2], lambda p,v: v, container=lambda x: 1)
     112      1
     113    """
     114   
     115    if raw: return res
     116   
     117    def correct_values(p, values):
     118        """ Correct in case of bad primes """
     119        if not values.count(None): return values
     120        if correction_fn:
     121            correct = correction_fn(p, values)
     122        else:
     123            return values
     124        # Mutating copy
     125        for i in range(len(values)):
     126            values[i] = correct[i]
     127
     128    def make_value(p, value):
     129        if value is None:
     130            return None
     131        else:
     132            return construct_fn(p, value)
     133   
     134    last_p = 0
     135    last_p_lst = []
     136    result = []
     137    for p, val in res:
     138        if p == last_p:
     139            last_p_lst.append(make_value(p, val))
     140        else:
     141            correct_values(last_p, last_p_lst)
     142           
     143            last_p = p
     144            last_p_lst = [make_value(p, val)]
     145            result.append((Integer(p), last_p_lst))
     146   
     147    # Index by ideals, not primes
     148    if index == "ideals":
     149        result2 = []
     150        for p, vals in result:
     151            for i, val in enumerate(vals):
     152                result2.append(((p, Integer(i)), val))
     153        result = result2
     154    elif index == "primes":
     155        pass
     156    else:
     157        raise ValueError("Traces cannot be indexed by `%s`; check your `index` parameter" % index)
     158   
     159    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    """ This callback collects a list of single-valued results in arg.tmp """
     12
     13    if good:
     14        result = (int(p), -a[0])
     15    else:
     16        result = (int(p), None)
     17   
     18    (<SmallJac>arg).tmp.append(result)
     19    return 1
     20
     21cdef int callback_single_single(smalljac_curve_t c, unsigned long p, int good,  long a[],  int n, void *arg):
     22    """ This callback collects a single single-valued result in arg.tmp """
     23
     24    if good:
     25        (<SmallJac>arg).tmp = -a[0]
     26    else:
     27        (<SmallJac>arg).tmp = None
     28    return 1
     29
     30cdef int callback_list_list(smalljac_curve_t c, unsigned long q, int good, long m[], int n, void *arg):
     31    """ This callback collects a list of list-valued results in arg.tmp """
     32
     33    cdef int i
     34
     35    if good:
     36        lst = []
     37        for i in range(n):
     38            lst.append(m[i])
     39        result = (int(q), lst)
     40    else:
     41        result = (int(q), None)
     42
     43    (<SmallJac>arg).tmp.append(result)
     44    return 1
     45
     46cdef int callback_single_list(smalljac_curve_t c, unsigned long q, int good, long m[], int n, void *arg):
     47    """ This callback collects a single list-valued result in arg.tmp """
     48
     49    cdef int i
     50    if good:
     51        lst = []
     52        for i in range(n):
     53            lst.append(m[i])
     54    else:
     55        lst = None
     56
     57    (<SmallJac>arg).tmp = lst
     58    return 1
     59
     60cdef int callback_histogram_single(smalljac_curve_t c, unsigned long q, int good, long a[], int n, void *arg):
     61    """
     62    This callback collects results into a histogram stored in array[2:].
     63
     64    INPUT:
     65
     66    - ``array[0]`` must be the genus of the curve.
     67    - ``array[1]`` must be the number of buckets (the size array, less two).
     68    """
     69
     70    cdef double val
     71    cdef unsigned int bucket
     72    cdef unsigned int buckets
     73    cdef unsigned int *array = <unsigned int *>arg
     74    cdef int genus
     75
     76    genus = array[0]
     77    buckets = array[1]
     78
     79    if good:
     80        val = 0.5 * (0.5 * -a[0] / sqrt(q) / genus + 1)
     81        bucket = (<int>(val * buckets))
     82        if bucket >= buckets:
     83            print "Anomaly at p=%d: a1=%d, so fraction=%f" % (q, a[0], val)
     84        else:
     85            array[bucket + 2] += 1
     86
     87    return 1
     88
     89smalljac_version = SMALLJAC_VERSION_STRING
     90
     91cdef class SmallJac:
     92    cdef smalljac_curve_t c
     93    cdef object tmp
     94    cdef int genus
     95   
     96    def __cinit__(self):
     97        """
     98        Initialize Smalljac curve pointer to the null pointer, for safety.
     99        """
     100       
     101        self.c = <smalljac_curve_t>0
     102       
     103    def __init__(self, v):
     104        """
     105        Create a Smalljac curve for an elliptic curve.
     106
     107        INPUT:
     108
     109        - ``v`` - A string defining the curve.  See Smalljac documentation for format.
     110
     111        EXAMPLE::
     112
     113            >>> SmallJac("[0,1,1,-10,-20]") # optional: smalljac
     114            Smalljac curve defined by [0,1,1,-10,-20]
     115
     116        """
     117       
     118        curve = str(v)
     119        cdef int err
     120
     121        self.c = smalljac_curve_init(curve, &err)
     122        self.genus = smalljac_curve_genus(self.c)
     123        self._check_error(err)
     124
     125    cdef _check_error(self, err):
     126        """
     127        Check whether a value is a error code, and if so raise an
     128        Exception.  If error is one of the standard Smalljac errors,
     129        use a descriptive message.  Some supported errors cannot
     130        actually happen with the functionality currently implemented
     131        (the file-based ones).
     132        """
     133       
     134        if err >= 0:
     135            return
     136        elif err == SMALLJAC_INTERNAL_ERROR:
     137            raise RuntimeError("Internal Smalljac Error")
     138        elif err == SMALLJAC_INVALID_INTERVAL:
     139            raise RuntimeError("Invalid Interval")
     140        elif err == SMALLJAC_PARSE_ERROR:
     141            raise RuntimeError("Cannot parse arguments")
     142        elif err == SMALLJAC_UNSUPPORTED_CURVE:
     143            raise RuntimeError("Unsupported curve; must be degree 3 or 5")
     144        elif err == SMALLJAC_SINGULAR_CURVE:
     145            raise RuntimeError("Curve is singular")
     146        elif err == SMALLJAC_INVALID_PP:
     147            raise RuntimeError("Not prime power, or too large a prime power, passed where prime power is required")
     148        elif err == SMALLJAC_WRONG_GENUS:
     149            raise RuntimeError("Only curves of genus 1 and 2 supported")
     150        elif err == SMALLJAC_INVALID_FLAGS:
     151            raise RuntimeError("Invalid combination of internal bit flags set")
     152        elif err == SMALLJAC_FILENOTFOUND:
     153            # Should never happen
     154            raise RuntimeError("Smalljac couldn't find the specified file")
     155        elif err == SMALLJAC_BADFILE:
     156            # Should never happen
     157            raise RuntimeError("Smalljac couldn't read specified file; is it corrupt?")
     158        elif err == SMALLJAC_NODATA:
     159            # Should never happen
     160            raise RuntimeError("Requested data not in specified file")
     161        elif err == SMALLJAC_NOT_OVER_Q:
     162            raise RuntimeError("Curve is not defined over QQ")
     163        else:
     164            raise RuntimeError("Unknown error %s" % err)
     165
     166    def __repr__(self):
     167        """
     168        Print a descriptive message for a Smalljac object.
     169
     170        EXAMPLE::
     171
     172            >>> SmallJac("[0,1,1,-10,-20]") # optional: smalljac
     173            Smalljac curve defined by [0,1,1,-10,-20]
     174        """
     175       
     176        return "Smalljac curve defined by %s"%smalljac_curve_str(self.c)
     177
     178    def __dealloc__(self):
     179        """ Clean up memory, such as the curve object itself. """
     180       
     181        if self.c:
     182            smalljac_curve_clear(self.c)
     183
     184    def ap(self, unsigned long p, unsigned long b=0):
     185        """
     186        If only p is given, return trace of Frobenius at p.  If both p
     187        and b are given, return a list of (key, value) pairs, with
     188        keys the primes < b, and values the traces of Frobenius at the
     189        good primes, and None at the bad primes, where bad is "bad for
     190        the given model".
     191
     192        EXAMPLE::
     193
     194            >>> s = SmallJac("[0,1,1,-10,-20]") # optional: smalljac
     195            >>> s.ap(7) # optional: smalljac
     196            -1
     197            >>> s.ap(next_prime(2^35)) # optional: smalljac
     198            245262
     199            >>> s.ap(0, 7) # optional: smalljac
     200            [(2, -2), (3, None), (5, -4), (7, -1)]
     201        """
     202
     203        cdef int err
     204
     205        if b == 0:
     206            # Compute for a single p only
     207            err = smalljac_parallel_Lpolys(self.c, p, p, SMALLJAC_A1_ONLY, callback_single_single, <void*>self)
     208        else:
     209            # Compute for a range of primes, starting with p.
     210            self.tmp = []
     211            err = smalljac_parallel_Lpolys(self.c, p, b,  SMALLJAC_A1_ONLY, callback_list_single, <void*>self)
     212
     213        self._check_error(err)
     214        return self.tmp
     215
     216    def sato_tate_buckets(self, unsigned int n, unsigned long a, unsigned long b):
     217        """
     218        Produce a histogram of traces of Frobenius, normalized to be
     219        between zero and one.
     220
     221        INPUT:
     222
     223        - ``n`` - number of buckets
     224        - ``a`` and ``b`` - bounds on the primes for which to compute
     225          traces of Frobenius, inclusive.
     226
     227        EXAMPLE:
     228
     229            >>> s = SmallJac("[0,1,1,-10,-20]") # optional: smalljac
     230            >>> s.sato_tate_buckets(10, 0, 1000) # optional: smalljac
     231            [7L, 12L, 22L, 24L, 18L, 17L, 27L, 14L, 14L, 10L]
     232        """
     233       
     234        cdef unsigned long i
     235        cdef unsigned int *buckets
     236        cdef int err
     237       
     238        buckets = <unsigned int *>(malloc((n + 2) * sizeof(unsigned int)))
     239       
     240        for i in range(n): # Memset would be faster, but it doesn't matter
     241            buckets[i+2] = 0
     242
     243        buckets[0] = self.genus
     244        buckets[1] = n
     245       
     246        err = smalljac_parallel_Lpolys(self.c, a, b,  SMALLJAC_A1_ONLY, callback_histogram_single, <void*>buckets)
     247
     248        self.tmp = []
     249        for i in range(n):
     250            self.tmp.append(buckets[i + 2])
     251       
     252        free(buckets)
     253        self._check_error(err)
     254        return self.tmp
     255
     256    def group(self, unsigned long p, unsigned long b=0):
     257        """
     258        If only p is given, return the group structure of the curve
     259        reduced at p.  If both p and b are given, return a list of
     260        (key, value) pairs, with keys the primes < b, and values are,
     261        at bad primes, None, and at good primes is a list of integers
     262        m1, m2, ..., representing the group Z/m1Z x Z/m2 x ...
     263
     264        EXAMPLE::
     265       
     266            >>> s = SmallJac("[0,1,1,-10,-20]") # optional: smalljac
     267            >>> s.group(7) # optional: smalljac
     268            [9]
     269
     270        This tells us that the curve has group structure Z/9Z when
     271        reduced at the finite field of order 7.
     272       
     273            >>> s.group(next_prime(2^35)) # optional: smalljac
     274            [34359493160]
     275            >>> s.group(0, 7) # optional: smalljac
     276            [(2, [5]), (3, None), (5, [10]), (7, [9])]
     277        """
     278
     279        cdef int err
     280       
     281        if b == 0:
     282            # Compute for a single p only
     283            err = smalljac_parallel_groups(self.c, p, p, 0, callback_single_list, <void*>self)
     284        else:
     285            # Compute for a range of primes, starting with p.
     286            self.tmp = []
     287            err = smalljac_parallel_groups(self.c, p, b, 0, callback_list_list, <void*>self)
     288
     289        self._check_error(err)
     290        return self.tmp
     291
     292    def lpoly(self, unsigned long p, unsigned long b=0):
     293        """
     294        If only p is given, return the L polynomial of the curve
     295        reduced at p.  If both p and b are given, return a list of
     296        (key, value) pairs, with keys the primes < b, and values are,
     297        at bad primes, None, and at good primes is a list of integers
     298        m1, m2, ..., representing the L polynomial
     299        p^n T^2n + m1 p^n-1 T^2n-1 m2 p^n-1 T^2n-2 + ... + m1 T + 1
     300
     301        EXAMPLE::
     302
     303            >>> s = SmallJac("[0,1,1,-10,-20]") # optional: smalljac
     304            >>> s.lpoly(7) # optional: smalljac
     305            [1]
     306
     307        This tells us that the curve has L-polynomial 7 T^2 + T + 1,
     308        where the 1 is the coefficient of T.
     309       
     310            >>> s.lpoly(next_prime(2^35)) # optional: smalljac
     311            [-245262]
     312            >>> s.lpoly(0, 7) # optional: smalljac
     313            [(2, [2]), (3, None), (5, [4]), (7, [1])]
     314        """
     315
     316        cdef int err
     317       
     318        if b == 0:
     319            # Compute for a single p only
     320            err = smalljac_parallel_Lpolys(self.c, p, p, 0, callback_single_list, <void*>self)
     321        # Compute for a range of primes, starting with p.
     322        else:
     323            self.tmp = []
     324            err = smalljac_parallel_Lpolys(self.c, p, b, 0, callback_list_list, <void*>self)
     325
     326        self._check_error(err)
     327        return self.tmp
     328
     329    def moments(self, int moments, unsigned long start, unsigned long end):
     330        """
     331        Computes the first few moments of the Sato-Tate distribution.
     332
     333        INPUT:
     334
     335        - ``moments`` - The number of moments to compute; at most 11.
     336        - ``start`` and ``end`` - the bounds on which reductions to use; inclusive.
     337
     338        EXAMPLE:
     339
     340            >>> s = SmallJac("[0,1,1,-10,-20]") # optional: smalljac
     341            >>> s.moments(5, 0, 1000) # optional: smalljac
     342            ([1.0, 0.002257165809854347, 0.9692801815679776,
     343              -0.06359909058782176, 1.8792555673855094],)
     344
     345        For a higher-genus curve, will return multiple lists of
     346        moments, one for each coefficient of the L-polynomial.
     347
     348            >>> s = SmallJac("x^5+2x-1") # optional: smalljac
     349            >>> s.moments(3, 0, 1000) # optional: smalljac
     350            ([1.0, -0.09346091047494341, 1.052736244717573],
     351             [1.0, 1.0335517476280245, 2.029919838936752])
     352        """
     353
     354        cdef int err
     355
     356        if moments > 11: raise ValueError("Only up to 11 moments supported for now.")
     357
     358        cdef double *moment_arr = <double *>malloc(sizeof(double) * moments * self.genus)
     359
     360        err = smalljac_moments(self.c, start, end, moment_arr, self.genus, moments,
     361                               NULL, NULL, NULL)
     362
     363        if err:
     364            free(moment_arr)
     365            self._check_err(err)
     366            # Should never happen
     367            raise RuntimeError("Error exists but doesn't exist; please report paradox!")
     368
     369        a1moments = []
     370
     371        cdef int i
     372        for i in range(moments):
     373            a1moments.append(moment_arr[i])
     374
     375        if self.genus == 2:
     376            a2moments = []
     377            for i in range(moments):
     378                a2moments.append(moment_arr[i + moments])
     379
     380            free(moment_arr)
     381            return (a1moments, a2moments)
     382        else:
     383            free(moment_arr)
     384            return (a1moments,)
     385
     386    def guess_group(self, unsigned long max):
     387        """
     388        Attempt to guess the Sato-Tate group of a curve, by testing
     389        various characteristics of the distribution of L-polynomial
     390        coefficients for reductions less than a certain size.  May
     391        rely on unproven conjectures.  May yield incorrect answers (it
     392        is, after all, only a guess), but attempts not to.  To
     393        interpret results, see
     394
     395            "Sato-Tate distributions and Galois endomorphism modules
     396            in genus 2", Francesc Fité, Kiran S. Kedlaya, Victor
     397            Rotger, and Andrew V. Sutherland, to appear in Compositio
     398            Mathematica, <http://arxiv.org/abs/1110.6638>.
     399
     400        INPUT:
     401
     402        - ``max`` - Largest norm of prime the reduction at which is tested.
     403
     404        EXAMPLE::
     405
     406            >>> s = SmallJac("[0,1,1,-10,-20]") # optional: smalljac
     407            >>> s.guess_group(1000) # optional: smalljac
     408            >>> s.guess_group(100000) # optional: smalljac
     409            'SU(2)'
     410            >>> s = SmallJac("x^5+2x-1") # optional: smalljac
     411            >>> s.guess_group(1000) # optional: smalljac
     412            >>> s.guess_group(100000) # optional: smalljac
     413            'USp(4)'
     414        """
     415
     416        cdef int i
     417        cdef char STgroup[16]
     418
     419        i = smalljac_identify_STgroup(self.c, STgroup, max)
     420
     421        if i:
     422            return STgroup
     423        else:
     424            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        EXAMPLE::
     1293
     1294            sage: E = EllipticCurve(GF(97), [2,3])
     1295            sage: E._smalljac() # optional: smalljac
     1296            Smalljac curve defined by [0,0,0,2,3]
     1297            sage: E = EllipticCurve('11a').reduction(17)
     1298            sage: E._smalljac() # optional: smalljac
     1299            Smalljac curve defined by [0,16,1,7,14]
     1300        """
     1301
     1302        try:
     1303            return self._smalljac_proxy
     1304        except AttributeError:
     1305            try:
     1306                import sage.libs.smalljac.all as smalljac
     1307            except ImportError:
     1308                raise RuntimeError, "smalljac not installed or inaccessible"
     1309           
     1310            f = "[%s,%s,%s,%s,%s]" % tuple(self.ainvs())
     1311            self._smalljac_proxy = smalljac.SmallJac(f)
     1312            return self._smalljac_proxy
     1313
     1314    def jacobian_structure(self, algorithm="smalljac"):
     1315        r"""
     1316        Returns the group structure of the group of points on this
     1317        elliptic curve.  Unlike ``abelian_group``, it does not compute
     1318        the generators of the group, just the structure, and is thus
     1319        perhaps somewhat less useful.
     1320
     1321        INPUT:
     1322       
     1323        - ``algorithm`` - either "sage", which uses Sage's abelian
     1324          group determination algorithm via ``abelian_group``, or
     1325          ``smalljac``, which uses Drew Sutherland's SmallJac library
     1326          for the computation.  SmallJac is much faster, but can only
     1327          be used on degree-1 fields.  You can also pass
     1328          ``algorithm="all"`` to run both algorithms and verify that
     1329          they produce the same results.
     1330
     1331        OUTPUT:
     1332
     1333        An abelian group isomorphic to the jacobian of this curve.
     1334
     1335        EXAMPLES::
     1336
     1337            sage: E=EllipticCurve(GF(11),[2,5])
     1338            sage: E.jacobian_structure() # optional: smalljac
     1339            Additive abelian group isomorphic to Z/10
     1340            sage: E=EllipticCurve(GF(41),[2,5])
     1341            sage: E.jacobian_structure() # optional: smalljac
     1342            Additive abelian group isomorphic to Z/2 + Z/22
     1343       
     1344        We can use ``algorithm="all"`` to test multiple
     1345        implementations against each other::
     1346       
     1347            sage: E.jacobian_structure(algorithm="all") # optional: smalljac
     1348            Additive abelian group isomorphic to Z/2 + Z/22
     1349        """
     1350
     1351        from sage.groups.additive_abelian.additive_abelian_group import AdditiveAbelianGroup
     1352
     1353        if algorithm == "smalljac":
     1354            K = self.base_field()
     1355            assert K.degree() == 1, "Smalljac can only find group structures in degree-1 fields"
     1356            p = K.cardinality()
     1357            proxy = self._smalljac()
     1358            structure = proxy.group(p)
     1359            return AdditiveAbelianGroup(structure)
     1360        elif algorithm == "sage":
     1361            g = self.abelian_group()
     1362            return AdditiveAbelianGroup(g.generator_orders())
     1363        elif algorithm == "all":
     1364            smalljac_results = self.jacobian_structure(algorithm="smalljac")
     1365            sage_results     = self.jacobian_structure(algorithm="sage")
     1366
     1367            # If the orders were given in different orders, testing
     1368            # two groups for equality can give incorrect results
     1369            smalljac_order = sorted(g.additive_order() for g in smalljac_results.gens())
     1370            sage_order     = sorted(g.additive_order() for g in sage_results.gens())
     1371            assert smalljac_order == sage_order, "An error has occured and smalljac and Sage's results do not match"
     1372            return smalljac_results
     1373        else:
     1374            raise ValueError("Unknown algorithm `%s`" % algorithm)
     1375       
    12871376    def abelian_group(self, debug=False):
    12881377        r"""
    12891378        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  
    21212121            degrees.extend(newdegs)
    21222122            k = k+1
    21232123
    2124         raise NotImplementedError, "Not all isogenies implemented over general number fields."
     2124        raise NotImplementedError, "Not all isogenies implemented over general number fields."
     2125
     2126    def aplists(self, limit, algorithm="smalljac", container=list, index="primes", raw=False):
     2127        """
     2128        Computes the traces of Frobenius at the prime ideals of the
     2129        base number field, sorted by the norm of these ideals.  The
     2130        structure of the return values depends on the ``format``
     2131        argument.  So far, this method requires the Smalljac library.
     2132
     2133        INPUT:
     2134
     2135        - ``limit`` - either an integer, in which case the prime
     2136          ideals will have norm less than this limit, or [a, b], in
     2137          which case the prime ideals with norms between ``a`` and
     2138          ``b`, inclusive, are used.
     2139
     2140        - ``algorithm`` - so far, only "smalljac" is supported, which
     2141          uses Drew Sutherland's smalljac library for elliptic curve
     2142          structure computation.
     2143
     2144        - The remainder (``container``, ``index``, ``raw``)
     2145          define the format of the result.  By default, the return
     2146          value looks like this:
     2147
     2148            [(9, [2]), (11, [4, -4]), (19, [-4, 4]), (31, [4, 8])]
     2149
     2150          This describes the fact that 9 does not split, and its prime
     2151          ideal has trace of Frobeneus 2; that 11 splits, and its two
     2152          prime ideals have traces 4 and -4; and so on.
     2153
     2154          Instead of a list of tuples, one can choose a dictionary,
     2155          mapping norms to traces, by changing ``container`` to
     2156          ``dict`` (in fact, any value for container is just called
     2157          with the list of pairs).
     2158
     2159          Instead of pairs of norms and lists of traces, you can have
     2160          each trace have its own pair, with the first element being a
     2161          pair of norm and index, yielding a result such as:
     2162
     2163            [((9, 0), 2), ((11, 0), 4), ((11, 0), -4), ...]
     2164
     2165          This is chosen with ``index="ideals"``.
     2166
     2167          By default, all integers above are Sage Integers.  If this
     2168          is a performance or memory problem, you likely want
     2169          ``raw=True``.
     2170
     2171          Finally, for performance or memory reasons, you can turn off
     2172          all post-processing by by setting ``raw`` to ``True``.  This
     2173          overrides all other arguments and returns a list formatted
     2174          like so:
     2175
     2176            [(9, 2), (11, None), (11, -4), (19, 4), (19, -4), ...]
     2177
     2178          All integers are python integers.  This last mode is by far
     2179          the fastest, especially for large ranges of prime norms, but
     2180          is also often the hardest to use for other computation.  It
     2181          is often used internally.  Note that, among other things,
     2182          using only raw data means that bad primes will have the
     2183          trace ``None``, instead of the correct trace, since it is
     2184          not SmallJac, but Sage, that does that computation.
     2185
     2186        OUTPUT:
     2187
     2188        A list of traces of Frobenius, with the precise format defined
     2189        by the ``container``, ``index``, and ``raw`` keywords.
     2190
     2191        EXAMPLES:
     2192
     2193        We can find the traces of Frobenius up to a given prime::
     2194
     2195            sage: ec = EllipticCurve('11a')
     2196            sage: ec.aplists(15) # optional: smalljac
     2197            [(2, [-2]), (3, [-1]), (5, [1]), (7, [-2]), (11, [None]), (13, [4])]
     2198
     2199        Or in any arbitrary range::
     2200
     2201            sage: ec.aplists([15, 40]) # optional: smalljac
     2202            [(17, [-2]), (19, [0]), (23, [-1]), (29, [0]), (31, [7]), (37, [3])]
     2203
     2204        Curves over number fields are also supported::
     2205
     2206            sage: K.<a> = NumberField(x^2 - x - 1)
     2207            sage: ec2 = EllipticCurve(K, [a, -a+1, 0, -4, 3*a-5])
     2208            sage: ec2.aplists([15, 50]) # optional: smalljac
     2209            [(19, [-4, 4]), (29, [-2, 6]), (31, [-8, 8]), (41, [2, -6]), (49, [2])]
     2210
     2211        You can choose to separate each ideal, instead of combining them by norm::
     2212
     2213            sage: ec2.aplists([15, 30], index="ideals") # optional: smalljac
     2214            [((19, 0), -4), ((19, 1), 4), ((29, 0), -2), ((29, 1), 6)]
     2215
     2216        You can also choose a dictionary representation, though this
     2217        of course destroys the ordering::
     2218
     2219            sage: ec2.aplists([15, 50], container=dict) # optional: smalljac
     2220            {41: [2, -6], 49: [2], 19: [-4, 4], 29: [-2, 6], 31: [-8, 8]}
     2221            sage: ec2.aplists([15, 30], container=dict, index="ideals") # optional: smalljac
     2222            {(19, 0): -4, (29, 1): 6, (19, 1): 4, (29, 0): -2}
     2223
     2224        Finally, the raw results are useful for large-scale statistics::
     2225
     2226            sage: bsum = __builtins__.sum
     2227            sage: bsum(ap for p, ap in ec.aplists(100000, raw=True) if ap) # optional: smalljac
     2228            4837
     2229            sage: bsum(ap for p, ap in ec2.aplists(100000, raw=True) if ap) # optional: smalljac
     2230            5289
     2231
     2232        Note that the equivalent non-SmallJac code for that last
     2233        example returns a different answer::
     2234
     2235            sage: bsum(ec.aplist(100000))
     2236            4838
     2237
     2238        This is because our curve has a bad reduction at 11, and we
     2239        aren't adding it its trace.
     2240        """
     2241
     2242        import sage.libs.smalljac.util as util
     2243
     2244        start, end = util.parse_limits(limit)
     2245       
     2246        if algorithm == "smalljac":
     2247            proxy = self._smalljac()
     2248            res = proxy.ap(start, end)
     2249
     2250            O = self.base_ring().ring_of_integers()
     2251            def get_aps(p, aps):
     2252                if hasattr(O.ideal(p), "factor"):
     2253                    factors = [ideal for ideal, exponent in O.ideal(p).factor()]
     2254                else:
     2255                    factors = [p]
     2256               
     2257                out = []
     2258                for ideal in factors:
     2259                    try:
     2260                        out.append(self.reduction(ideal).trace_of_frobenius())
     2261                    except (AttributeError, ValueError): # Bad reduction
     2262                        out.append(None)
     2263                return out
     2264            def make_int(p, v):
     2265                return Integer(v)
     2266
     2267            return util.interpret_smalljac_results(
     2268                res, correction_fn=get_aps, construct_fn=make_int,
     2269                raw=raw, index=index, container=container)
     2270        else:
     2271            raise ValueError("Algorithm %s not defined" % algorithm)
     2272
     2273    def aps(self, p):
     2274        """
     2275        Returns a list of the traces of Frobenius of this elliptic
     2276        curve when reduced at prime ideals of norm ``p``.
     2277
     2278        INPUT:
     2279
     2280         - ``p`` - The integer prime at which to reduce.
     2281
     2282        OUTPUT:
     2283
     2284        A list of traces of Frobenius, or None if there was a bad
     2285        reduction.
     2286
     2287        EXAMPLES:
     2288
     2289        We can take a curve and find traces::
     2290
     2291            sage: ec = EllipticCurve('11a')
     2292            sage: ec.aps(next_prime(10^6)) # optional: smalljac
     2293            [284]
     2294            sage: ec.aps(next_prime(10^9)) # optional: smalljac
     2295            [-1962]
     2296
     2297        Curves over quadratic number fields are also supported::
     2298
     2299            sage: K.<a> = NumberField(x^2 - x - 1)
     2300            sage: ec2 = EllipticCurve(K, [a, -a+1, 0, -4, 3*a-5])
     2301            sage: ec2.aps(9949) # optional: smalljac
     2302            [-162, -34]
     2303        """
     2304
     2305        res = self.aplists([p, p])
     2306
     2307        if len(res) == 0 or res[0][0] != p:
     2308            raise ValueError("%d is not a prime (over the chosen field)" % p)
     2309
     2310        return res[0][1] # [0] extracts the prime we care about, [1] extracts the list of traces
     2311
     2312    def grouplists(self, limit, algorithm="smalljac", container=list, index="primes", raw=False):
     2313        r"""
     2314        Computes the group structure of reductions at the prime ideals
     2315        of the base number field, sorted by the norm of these ideals.
     2316        So far, this method requires the Smalljac library.
     2317
     2318        INPUT:
     2319
     2320        - ``limit`` - either an integer, in which case the prime
     2321          ideals will have norm less than this limit, or [a, b], in
     2322          which case the prime ideals with norms between ``a`` and
     2323          ``b`, inclusive, are used.
     2324
     2325        - ``algorithm`` - so far, only "smalljac" is supported, which
     2326          uses Drew Sutherland's smalljac library for elliptic curve
     2327          structure computation.
     2328
     2329        - The remainder (``container``, ``index``, and ``raw``) define
     2330          the format of the result; see ``EllipticCurve.aplists`` for
     2331          a description.
     2332
     2333        OUTPUT:
     2334
     2335        A list of group structures of Jacobians at various primes; the
     2336        format is precisely specified by the ``container``, ``index``,
     2337        and ``raw`` arguments.  If a curve has a bad reduction at a prime,
     2338        ``None`` is instead returned.
     2339
     2340        EXAMPLES::
     2341
     2342            sage: ec = EllipticCurve('11a')
     2343            sage: ec.grouplists(10) # optional: smalljac
     2344            [(2, [Additive abelian group isomorphic to Z/5]),
     2345             (3, [Additive abelian group isomorphic to Z/5]),
     2346             (5, [Additive abelian group isomorphic to Z/5]),
     2347             (7, [Additive abelian group isomorphic to Z/10])]
     2348
     2349            sage: ec.grouplists([50, 60]) # optional: smalljac
     2350            [(53, [Additive abelian group isomorphic to Z/2 + Z/30]), (59, [Additive abelian group isomorphic to Z/55])]
     2351
     2352        The raw results are suitable for larger-scale statistics::
     2353
     2354            sage: def jac_checksum(poly):
     2355            ...       ec = EllipticCurve(poly)
     2356            ...       sum = __builtins__.sum                           
     2357            ...       return sum(sum(gs) for p, gs in ec.grouplists([3, 10000], raw=True) if gs)
     2358            sage: jac_checksum([4, 7]) # optional: smalljac
     2359            5179039
     2360            sage: jac_checksum('11a') # optional: smalljac
     2361            4121837
     2362            sage: jac_checksum('37b') # optional: smalljac
     2363            3455881
     2364        """
     2365
     2366        from sage.groups.additive_abelian.additive_abelian_group import AdditiveAbelianGroup
     2367        import sage.libs.smalljac.util as util
     2368
     2369        start, end = util.parse_limits(limit)
     2370
     2371        if algorithm == "smalljac":
     2372            proxy = self._smalljac()
     2373            res = proxy.group(start, end)
     2374
     2375            O = self.base_ring().ring_of_integers()
     2376            def get_group(p, gps):
     2377                if hasattr(O.ideal(2), "factor"): # Not the rationals
     2378                    factors = [ideal for ideal, exponent in O.ideal(p).factor()]
     2379                else: # Only in the case of QQ
     2380                    factors = [p]
     2381               
     2382                out = []
     2383                for ideal in factors:
     2384                    try:
     2385                        out.append(self.reduction(ideal).jacobian_structure())
     2386                    except AttributeError:
     2387                        out.append(None)
     2388                return out
     2389           
     2390            def make_group(p, v):
     2391                return AdditiveAbelianGroup(v)
     2392
     2393            return util.interpret_smalljac_results(
     2394                res, correction_fn=get_group, construct_fn=make_group,
     2395                raw=raw, index=index, container=container)
     2396        else:
     2397            raise ValueError("Algorithm %s not defined" % algorithm)
     2398
     2399    def lpolylists(self, limit, algorithm="smalljac", container=list, index="primes", normalize=False, raw=False):
     2400        r"""
     2401        Computes the L-polynomials of reductions at the prime ideals
     2402        of the base number field, sorted by the norm of these ideals.
     2403        So far, this method requires the Smalljac library.
     2404
     2405        INPUT:
     2406
     2407        - ``limit`` - either an integer, in which case the prime
     2408          ideals will have norm less than this limit, or [a, b], in
     2409          which case the prime ideals with norms between ``a`` and
     2410          ``b`, inclusive, are used.
     2411
     2412        - ``normalize`` - whether to return normalized L-polynomials.
     2413          By default, non-normalized L-polynomials are returned.
     2414
     2415        - ``algorithm`` - so far, only "smalljac" is supported, which
     2416          uses Drew Sutherland's smalljac library for elliptic curve
     2417          structure computation.
     2418
     2419        - The remainder (``container``, ``index``, and ``raw``) define
     2420          the format of the result; see ``EllipticCurve.aplists`` for
     2421          a description.
     2422
     2423        OUTPUT:
     2424
     2425        A list of L-polynomials of reductions at various primes; the
     2426        format is precisely specified by the ``container``, ``index``,
     2427        and ``raw`` arguments.  If a curve has a bad reduction at a
     2428        prime, ``None`` is instead returned.
     2429
     2430        EXAMPLES:
     2431
     2432        We can construct an elliptic curve and find its L-polynomial::
     2433
     2434            sage: E = EllipticCurve('11a')
     2435            sage: E.lpolylists([89, 89]) # optional: smalljac
     2436            [(89, [89*T^2 - 15*T + 1])]
     2437
     2438        Or we can use a range::
     2439
     2440            sage: E.lpolylists([80, 100]) # optional: smalljac
     2441            [(83, [83*T^2 + 6*T + 1]), (89, [89*T^2 - 15*T + 1]), (97,
     2442            [97*T^2 + 7*T + 1])]
     2443
     2444        We can ask for normalized L-polynomials instead::
     2445
     2446            sage: E.lpolylists([89, 89], normalize=True) # optional: smalljac
     2447            [(89, [T^2 - 1.58999682000954*T + 1.00000000000000])]
     2448
     2449        Number fields are also supported::
     2450
     2451            sage: K.<a> = NumberField(x^2 - x - 1)
     2452            sage: E2 = EllipticCurve(K, [a, -a+1, 0, -4, 3*a-5])
     2453            sage: E2.lpolylists([80, 100]) # optional: smalljac
     2454            [(89, [89*T^2 + 14*T + 1, 89*T^2 - 2*T + 1])]
     2455        """
     2456
     2457        from sage.rings.polynomial.polynomial_ring import PolynomialRing_field
     2458        from sage.rings.real_mpfr import RR
     2459        from sage.rings.rational_field import QQ
     2460        from math import sqrt
     2461        import sage.libs.smalljac.util as util
     2462
     2463        start, end = util.parse_limits(limit)
     2464
     2465        if algorithm == "smalljac":
     2466            proxy = self._smalljac()
     2467            res = proxy.lpoly(start, end)
     2468
     2469            if normalize:
     2470                polynomial_ring = PolynomialRing_field(RR, 'T')
     2471            else:
     2472                polynomial_ring = PolynomialRing_field(QQ, 'T')
     2473            t = polynomial_ring.gens()[0] # I don't want a T
     2474
     2475            def make_lpoly(p, coeffs):
     2476                if not normalize:
     2477                    a1 = coeffs[0]
     2478                    return p*t**2 + a1*t + 1
     2479                else:
     2480                    a1 = coeffs[0]
     2481                    a1 /= sqrt(p)
     2482                    return t**2 + a1*t + 1
     2483
     2484            def get_lpoly(p, coeffs):
     2485                aps = self.aps(p)
     2486
     2487                out = []
     2488                for ap in aps:
     2489                    if not normalize:
     2490                        out.append(p*t**2 - ap*t + 1)
     2491                    else:
     2492                        ap /= sqrt(p)
     2493                        out.append(p*t**2 - ap*t + 1)
     2494                return out
     2495
     2496            return util.interpret_smalljac_results(
     2497                res, correction_fn=get_lpoly, construct_fn=make_lpoly,
     2498                raw=raw, index=index, container=container)
     2499        else:
     2500            raise ValueError("Algorithm %s not defined" % algorithm)
     2501
     2502    def _smalljac(self):
     2503        """
     2504        Computes and caches the smalljac representation of the given
     2505        curve.  For internal use only.
     2506
     2507        EXAMPLES::
     2508
     2509            sage: E = EllipticCurve('11a')
     2510            sage: E._smalljac() # optional: smalljac
     2511            Smalljac curve defined by [0,-1,1,-10,-20]
     2512            sage: K.<a> = NumberField(x^2 - x - 1)
     2513            sage: E = EllipticCurve(K, [a, -a+1, 0, -4, 3*a-5])
     2514            sage: E._smalljac() # optional: smalljac
     2515            Smalljac curve defined by [a, -a + 1, 0, -4, 3*a - 5] / (a^2 - a - 1)
     2516        """
     2517
     2518        try:
     2519            return self._smalljac_proxy
     2520        except AttributeError:
     2521            try:
     2522                import sage.libs.smalljac.all as smalljac
     2523            except ImportError:
     2524                raise RuntimeError, "smalljac not installed or inaccessible"
     2525
     2526            k = self.base_ring()
     2527            gen = k.gen()
     2528            var = str(gen)
     2529            cpoly = gen.charpoly().change_variable_name(var)
     2530            f = "[%s, %s, %s, %s, %s]" % tuple(self.ainvs())
     2531            f = "%s / (%s)" % (f, cpoly)
     2532            self._smalljac_proxy = smalljac.SmallJac(f)
     2533            return self._smalljac_proxy
     2534
     2535    def trace_histogram(self, limit, buckets=None):
     2536        """
     2537        Computes the traces of reductions of this curve at all primes
     2538        in a range, and returns how many, normalized, fell in each
     2539        of ``buckets`` buckets from -2 to 2.
     2540
     2541        INPUT:
     2542
     2543        - ``limits`` - Either an integer, meaning all primes below that
     2544          integer are used, or a list of ``[start, end]``, in which
     2545          case all primes from ``start`` to ``end``, inclusive, are
     2546          used.
     2547
     2548        - ``buckets`` - How many buckets to split the range $[-2, 2]$
     2549          into.  By default, the number of buckets is the square root
     2550          of the range size.
     2551
     2552        OUTPUT:
     2553
     2554        A list of integers, as many as the ``buckets`` parameter,
     2555        denoting the number of trace values fell into that bucket for
     2556        primes in the given range.
     2557
     2558        EXAMPLES:
     2559
     2560        We can take an elliptic curve and compute a histogram for it::
     2561
     2562            sage: E1 = EllipticCurve([4,7])
     2563            sage: E1.trace_histogram(100) # optional: smalljac
     2564            [0, 1, 2, 1, 3, 8, 2, 4, 2, 1]
     2565
     2566        In this case, the default choice of 10 buckets was made.
     2567        Instead, you can specify the number of bucket manually::
     2568
     2569            sage: E1.trace_histogram(100000, buckets=10) # optional: smalljac
     2570            [529, 849, 1059, 1193, 1201, 1153, 1133, 1092, 876, 505]
     2571
     2572        You can choose a range that starts anywhere, not just at 1::
     2573
     2574            sage: E1.trace_histogram([1000, 2000], buckets=10) # optional: smalljac
     2575            [11, 13, 12, 14, 14, 14, 17, 11, 20, 8]
     2576
     2577        This is generally useful for making histograms; a quick way to
     2578        do this is with ``bar_chart``, though custom graphing code
     2579        will look better. ::
     2580
     2581            sage: bar_chart(E1.trace_histogram(100000), width=1) #optional: smalljac
     2582        """
     2583
     2584        import sage.libs.smalljac.util as util
     2585        from math import sqrt
     2586
     2587        start, end = util.parse_limits(limit)
     2588
     2589        if buckets is None:
     2590            range_size = end - start + 1
     2591            buckets = int(sqrt(range_size))
     2592
     2593        return map(Integer,
     2594                   self._smalljac().sato_tate_buckets(buckets, start, end))
     2595
     2596    def moments(self, limits, moments=11):
     2597        """
     2598        Compute the moments of the distributions of the coefficients
     2599        of the L-polynomials of this curve, reduced at primes in a
     2600        range.
     2601
     2602        INPUT:
     2603
     2604        - ``limits`` - Either an integer, meaning all primes below that
     2605          integer are used, or a list of ``[start, end]``, in which
     2606          case all primes from ``start`` to ``end``, inclusive, are
     2607          used.
     2608
     2609        - ``moments`` - How many moments to return; fewer moments
     2610          might be faster to compute.  For now, at most 10 moments can
     2611          be computed.
     2612
     2613        OUTPUT:
     2614
     2615        A tuple whose first element is a list of moments of the
     2616        distribution of the $a_1$ coefficients of L-polynomials of
     2617        this curve, starting from the zeroth moment being 0.  The list
     2618        is as long as ``moments``; that is, if ``moments`` is 5, only
     2619        the zero through fourth moments are returned.
     2620
     2621        .. NOTE::
     2622            A tuple is returned for polymorphism with the
     2623            hyperelliptic curve method of the same name
     2624
     2625        EXAMPLES:
     2626
     2627        You can request the moments of an arbitrary curve::
     2628       
     2629            sage: ec = EllipticCurve('37b')
     2630            sage: ec.moments(100000) # optional: smalljac
     2631            ([1.0, -0.003018470650554107, 0.9999566856391424,
     2632              -0.017815838483148164, 1.993903401871212,
     2633              -0.06113842272925238, 4.963581622954891,
     2634              -0.20681637970442066, 13.83188351590207,
     2635              -0.7065524181435564, 41.306996429058835],)
     2636
     2637        You can request fewer moments if you wish::
     2638
     2639            sage: ec.moments(10000, moments=3) # optional: smalljac
     2640            ([1.0, -0.000812699651702702, 0.9890134939436471],)
     2641
     2642        Or compute in any range::
     2643
     2644            sage: ec.moments([1000, 2000], moments=5) # optional: smalljac
     2645            ([1.0, -0.008380482644508922, 0.9721022301822473,
     2646              -0.16345224489884563, 1.8856843129225396],)
     2647
     2648        Number fields are also supported::
     2649
     2650            sage: K.<a> = NumberField(x^2 - x - 1)
     2651            sage: ec2 = EllipticCurve(K, [a, -a+1, 0, -4, 3*a-5])
     2652            sage: ec2.moments([1000, 2000], moments=5) # optional: smalljac
     2653            ([1.0, -0.047677896726755147, 0.9514627444845026,
     2654              -0.20252616280144933, 1.7720971308862585],)
     2655
     2656        Rounding the moments should produce a recognizable sequence::
     2657
     2658            sage: map(round, ec.moments(100000)[0]) # optional: smalljac
     2659            [1.0, -0.0, 1.0, -0.0, 2.0, -0.0, 5.0, -0.0, 14.0, -1.0, 41.0]
     2660        """
     2661
     2662        import sage.libs.smalljac.util as util
     2663
     2664        start, end = util.parse_limits(limits)
     2665        if moments > 11:
     2666            raise ValueError("Only up to 11 moments supported for now.")
     2667
     2668        return self._smalljac().moments(moments, start, end)
     2669
     2670    def guess_sato_tate_group(self, limit=None):
     2671        """
     2672        Provide a preliminary guess toward the Sato-Tate group
     2673        corresponding to this curve, by considering its reduction at
     2674        primes up to a certain limit.
     2675
     2676        INPUT:
     2677
     2678        - ``limit`` - The maximal prime considered.  The search may
     2679          terminate earlier if an answer is found.  Or, ``None``, to
     2680          use an implementation-defined (but very large) maximum.
     2681
     2682        OUTPUT:
     2683
     2684        A string containing the name of the corresponding Sato-Tate
     2685        group.  If no provisional identification can be made, returns
     2686        ``None``.
     2687
     2688        .. TODO:
     2689           - Should output a textual description of the Sato-Tate group
     2690
     2691        EXAMPLES:
     2692
     2693        Just create a hyperelliptic curve of genus 2 and call
     2694        ``guess_sato_tate_group`` to start it off::
     2695
     2696            sage: ec = EllipticCurve('11a')
     2697            sage: ec.guess_sato_tate_group(100000) # optional: smalljac
     2698            'SU(2)'
     2699            sage: ec = EllipticCurve([0, 1])
     2700            sage: ec.guess_sato_tate_group(100000) # optional: smalljac
     2701            'N(U(1))'
     2702            sage: K.<a> = NumberField(x^2 + 3)
     2703            sage: ec = EllipticCurve(K, [0, 1])
     2704            sage: ec.guess_sato_tate_group(100000) # optional: smalljac
     2705            'U(1)'
     2706
     2707        The above are the only three possibilities (for genus 1).
     2708        """
     2709
     2710        if limit is None:
     2711            limit = 0
     2712
     2713        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  
    750750            return self.__database_curve
    751751           
    752752
    753     def Np(self, p):
     753    def Np(self, p, algorithm="pari"):
    754754        r"""
    755         The number of points on `E` modulo `p`.
     755        The number of points on `E` modulo `p`.  The ``algorithm``
     756        argument is the same as that for ``ap``.
    756757
    757758        INPUT:
    758759
    … …  
    761762
    762763        OUTPUT:
    763764
    764         (int) The number ofpoints on the reduction of `E` modulo `p`
     765        (int) The number of points on the reduction of `E` modulo `p`
    765766        (including the singular point when `p` is a prime of bad
    766767        reduction).
    767768       
    … …  
    784785            1000000000000001426441464441649
    785786        """
    786787        if self.conductor() % p == 0:
    787             return p + 1 - self.ap(p)
    788         return p+1 - self.ap(p)
     788            return p + 1 - self.ap(p, algorithm=algorithm)
     789        return p+1 - self.ap(p, algorithm=algorithm)
    789790
    790791    #def __pari_double_prec(self):
    791792    #    EllipticCurve_number_field._EllipticCurve__pari_double_prec(self)
    … …  
    884885    #  Etc.
    885886    ####################################################################
    886887
    887     def aplist(self, n, python_ints=False):
     888    def aplist(self, n, python_ints=False, algorithm="pari"):
    888889        r"""
    889890        The Fourier coefficients `a_p` of the modular form
    890891        attached to this elliptic curve, for all primes `p\leq n`.
    … …  
    896897       
    897898        -  ``python_ints`` - bool (default: False); if True
    898899           return a list of Python ints instead of Sage integers.
     900
     901        -  ``algorithm`` - std (default: "pari")
     902           
     903           - "pari" - Use PARI's implementation of ellaplist
     904
     905           - "smalljac" - Use the smalljac library, if installed
    899906       
    900907       
    901908        OUTPUT: list of integers
    … …  
    915922            <type 'sage.rings.integer.Integer'>
    916923            sage: type(e.aplist(13, python_ints=True)[0])
    917924            <type 'int'>
    918         """
    919         e = self.pari_mincurve()
    920         v = e.ellaplist(n, python_ints=True)
    921         if python_ints:
     925
     926        We can pass ``algorithm="smalljac"`` to use the smalljac
     927        library; this is generally much faster.
     928
     929            sage: e.aplist(13, algorithm="smalljac") # optional: smalljac
     930            [-2, -3, -2, -1, -5, -2]
     931        """
     932       
     933        if algorithm == "pari":
     934            e = self.pari_mincurve()
     935            v = e.ellaplist(n, python_ints=True)
     936            if python_ints:
     937                return v
     938            else:
     939                return [Integer(a) for a in v]
     940        elif algorithm == "smalljac":
     941            proxy = self._smalljac()
     942            aps = proxy.ap(0, n)
     943           
     944            v = []
     945            for p, ap in aps:
     946                if ap is None:
     947                    # Sage already deals with this case well enough
     948                    ap = self.ap(p, algorithm="pari")
     949                    if not python_ints: ap = int(ap)
     950                    v.append(ap)
     951                else:
     952                    if python_ints: ap = Integer(ap)
     953                    v.append(ap)
     954           
    922955            return v
     956        elif algorithm == "all":
     957            aplist1 = self.aplist(n, python_ints=python_ints, algorithm="pari")
     958            aplist2 = self.aplist(n, python_ints=python_ints, algorithm="smalljac")
     959            if aplist1 != aplist2:
     960                raise ArithmeticError("PARI and smalljac disagree (%s, %s)", aplist1, aplist2)
     961            return aplist1
    923962        else:
    924             return [Integer(a) for a in v]
    925        
    926        
     963            raise RuntimeError, "algorithm '%s' is not known."%algorithm
    927964
    928965    def anlist(self, n, python_ints=False):
    929966        """
    … …  
    13771414        else:
    13781415            raise ValueError, "algorithm %s not defined"%algorithm
    13791416
     1417    def _smalljac(self):
     1418        r"""
     1419        Computes and caches the smalljac representation of the given
     1420        curve.  For internal use only.
     1421
     1422        EXAMPLES::
     1423
     1424            sage: E = EllipticCurve('11a')
     1425            sage: E._smalljac() # optional: smalljac
     1426            Smalljac curve defined by [0,-1,1,-10,-20]
     1427        """
     1428
     1429        try:
     1430            return self._smalljac_proxy
     1431        except AttributeError:
     1432            try:
     1433                import sage.libs.smalljac.all as smalljac
     1434            except ImportError:
     1435                raise RuntimeError, "smalljac not installed or inaccessible"
     1436           
     1437            f = "[%s,%s,%s,%s,%s]" % tuple(self.ainvs())
     1438            self._smalljac_proxy = smalljac.SmallJac(f)
     1439            return self._smalljac_proxy
    13801440
    13811441    def simon_two_descent(self, verbose=0, lim1=5, lim3=50, limtriv=10, maxprob=20, limbigprime=30):
    13821442        r"""
    … …  
    25012561        """
    25022562        return Integer(self.pari_mincurve().ellak(n))
    25032563
    2504     def ap(self, p):
     2564    def ap(self, p, algorithm="pari"):
    25052565        """
    25062566        The p-th Fourier coefficient of the modular form corresponding to
    25072567        this elliptic curve, where p is prime.
    25082568       
     2569        -  ``algorithm`` - std (default: "pari")
     2570           
     2571           - "pari" - Use PARI's implementation of ellaplist
     2572
     2573           - "smalljac" - Use the smalljac library, if installed
     2574       
     2575       
    25092576        EXAMPLES::
    25102577       
    25112578            sage: E=EllipticCurve('37a1')
    25122579            sage: [E.ap(p) for p in prime_range(50)]
    25132580            [-2, -3, -2, -1, -5, -2, 0, 0, 2, 6, -4, -1, -9, 2, -9]
    2514         """
     2581
     2582        If the ``smalljac`` libary is installed, we can use it::
     2583
     2584            sage: [E.ap(p, algorithm="smalljac") for p in prime_range(50)] # optional: smalljac
     2585            [-2, -3, -2, -1, -5, -2, 0, 0, 2, 6, -4, -1, -9, 2, -9]
     2586        """
     2587       
    25152588        if not arith.is_prime(p):
    25162589            raise ArithmeticError, "p must be prime"
    2517         return Integer(self.pari_mincurve().ellap(p))
     2590       
     2591        if algorithm == "pari":
     2592            return Integer(self.pari_mincurve().ellap(p))
     2593        elif algorithm == "smalljac":
     2594            proxy = self._smalljac()
     2595            ap = proxy.ap(p)
     2596            if ap is None:
     2597                # Let PARI handle the icky case
     2598                return self.ap(p, algorithm="pari")
     2599            else:
     2600                return Integer(ap)
     2601        elif algorithm == "all":
     2602            pari_results     = self.ap(p, algorithm="pari")
     2603            smalljac_results = self.ap(p, algorithm="smalljac")
     2604            assert pari_results == smalljac_results, "An error has occured, and SmallJac's and PARI's results are not the same"
     2605            return smalljac_results
    25182606
    25192607    def quadratic_twist(self, D):
    25202608        """
    … …  
    46834771        """
    46844772        return self.conductor().is_squarefree()
    46854773
    4686     def is_ordinary(self, p, ell=None):
     4774    def is_ordinary(self, p, ell=None, algorithm="pari"):
    46874775        """
    46884776        Return True precisely when the mod-p representation attached to
    4689         this elliptic curve is ordinary at ell.
     4777        this elliptic curve is ordinary at ell.  The ``algorithm``
     4778        argument has the same values as that for ``aplist``.
    46904779       
    46914780        INPUT:
    46924781       
    … …  
    47084797        """
    47094798        if ell is None:
    47104799            ell = p
    4711         return self.ap(ell) % p != 0
     4800        return self.ap(ell, algorithm=algorithm) % p != 0
    47124801
    47134802    def is_good(self, p, check=True):
    47144803        """
    … …  
    47684857            ell = p
    47694858        return self.is_good(p) and not self.is_ordinary(p, ell)
    47704859
    4771     def supersingular_primes(self, B):
     4860    def supersingular_primes(self, B, algorithm="pari"):
    47724861        """
    47734862        Return a list of all supersingular primes for this elliptic curve
    4774         up to and possibly including B.
     4863        up to and possibly including B.  The ``algorithm`` argument has
     4864        the same values as that for ``aplist``.
    47754865       
    47764866        EXAMPLES::
    47774867       
    … …  
    47974887            sage: e.supersingular_primes(1)
    47984888            []
    47994889        """
    4800         v = self.aplist(max(B, 3))
     4890        v = self.aplist(max(B, 3), algorithm=algorithm)
    48014891        P = rings.prime_range(max(B,3)+1)
    48024892        N = self.conductor()
    48034893        return [P[i] for i in [0,1] if P[i] <= B and v[i]%P[i]==0 and N%P[i] != 0] + \
    48044894                      [P[i] for i in range(2,len(v)) if v[i] == 0 and N%P[i] != 0]
    48054895
    4806     def ordinary_primes(self, B):
     4896    def ordinary_primes(self, B, algorithm="pari"):
    48074897        """
    48084898        Return a list of all ordinary primes for this elliptic curve up to
    4809         and possibly including B.
     4899        and possibly including B.  The ``algorithm`` argument has the same
     4900        values as that for ``aplist``.
    48104901       
    48114902        EXAMPLES::
    48124903       
    … …  
    48294920            sage: e.ordinary_primes(1)
    48304921            []
    48314922        """
    4832         v = self.aplist(max(B, 3) )
     4923        v = self.aplist(max(B, 3), algorithm=algorithm)
    48334924        P = rings.prime_range(max(B,3) +1)
    48344925        return [P[i] for i in [0,1] if P[i] <= B and v[i]%P[i]!=0] +\
    48354926               [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        EXAMPLES::
     23       
     24            sage: R.<x> = GF(next_prime(10^5))[]
     25            sage: hec = HyperellipticCurve(x^5 + 17*x^4 - 318*x^3 + 2*x - 17)
     26            sage: hec._smalljac() # optional: smalljac
     27            Smalljac curve defined by x^5 + 17*x^4 + 99685*x^3 + 2*x + 99986
     28        """
     29
     30        try:
     31            return self._smalljac_proxy
     32        except AttributeError:
     33            try:
     34                import sage.libs.smalljac.all as smalljac
     35            except ImportError:
     36                raise RuntimeError, "smalljac not installed or inaccessible"
     37           
     38            f, h = self.hyperelliptic_polynomials()
     39            assert h == 0, "SmallJac can only does with hyperelliptic polynomials in Weirstrass form"
     40            s = str(f)
     41            self._smalljac_proxy = smalljac.SmallJac(s)
     42            return self._smalljac_proxy
     43
     44    def jacobian_structure(self, algorithm="smalljac"):
     45        r"""
     46        Returns the group structure of the group of points on this
     47        elliptic curve.  Unlike ``abelian_group``, it does not compute
     48        the generators of the group, just the structure, and is thus
     49        perhaps somewhat less useful.
     50
     51        INPUT:
     52       
     53        - ``algorithm`` - so far, only ``"smalljac"`` is supported,
     54          which uses Drew Sutherland's SmallJac library for the
     55          computation.  SmallJac can only be used on degree-1 fields.
     56
     57        OUTPUT:
     58
     59        An ``AdditiveAbelianGroup`` representing the group structure
     60        of the jacobian of the curve.
     61
     62        EXAMPLES:
     63
     64        We can take an arbitrary curve and ask for the jacobian's
     65        structure::
     66
     67            sage: R.<x> = GF(next_prime(10^5))[]
     68            sage: hec = HyperellipticCurve(x^5 + 17*x^4 - 318*x^3 + 2*x - 17)
     69            sage: hec.jacobian_structure() # optional: smalljac
     70            Additive abelian group isomorphic to Z/2 + Z/4974328512
     71        """
     72       
     73        from sage.groups.additive_abelian.additive_abelian_group import AdditiveAbelianGroup
     74
     75        if algorithm == "smalljac":
     76            K = self.base_ring()
     77            assert K.degree() == 1, "Smalljac can only find group structures in degree-1 fields"
     78            p = K.cardinality()
     79            proxy = self._smalljac()
     80            structure = proxy.group(p)
     81            return AdditiveAbelianGroup(structure)
     82        else:
     83            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        EXAMPLES::
     24       
     25            sage: R.<x> = QQ[]
     26            sage: hec = HyperellipticCurve(x^5 + 2*x^4 - x^3 - 3*x^2 - x)
     27            sage: hec._smalljac() # optional: smalljac
     28            Smalljac curve defined by x^5 + 2*x^4 - x^3 - 3*x^2 - x
     29            sage: hec = HyperellipticCurve(x^5 + 17*x^4 - 318*x^3 + 2*x - 17)
     30            sage: hec._smalljac() # optional: smalljac
     31            Smalljac curve defined by x^5 + 17*x^4 - 318*x^3 + 2*x - 17
     32        """
     33
     34        try:
     35            return self._smalljac_proxy
     36        except AttributeError:
     37            try:
     38                import sage.libs.smalljac.all as smalljac
     39            except ImportError:
     40                raise RuntimeError, "smalljac not installed or inaccessible"
     41           
     42            f, h = self.hyperelliptic_polynomials()
     43            assert h == 0, "SmallJac can only does with hyperelliptic polynomials in Weirstrass form"
     44            s = str(f)
     45            self._smalljac_proxy = smalljac.SmallJac(s)
     46            return self._smalljac_proxy
     47
     48    def aplists(self, limit, algorithm="smalljac", container=list, index="primes", raw=False):
     49        """
     50        Computes the traces of Frobenius at the prime ideals of the
     51        base number field, sorted by the norm of these ideals.  So
     52        far, this method requires the Smalljac library.
     53
     54        INPUT:
     55
     56        - ``limit`` - either an integer, in which case the prime
     57          ideals will have norm less than this limit, or [a, b], in
     58          which case the prime ideals with norms between ``a`` and
     59          ``b`, inclusive, are used.
     60
     61        - ``algorithm`` - so far, only "smalljac" is supported, which
     62          uses Drew Sutherland's smalljac library for elliptic curve
     63          structure computation.
     64
     65        - The remainder (``container``, ``index``, and ``raw``) define
     66          the format of the result; see ``EllipticCurve.aplists`` for
     67          a description.
     68
     69        OUTPUT:
     70
     71        A list of traces of Frobenius at various primes; the format is
     72        precisely specified by the ``container``, ``index``, and
     73        ``raw`` arguments.  If a curve has a bad reduction at a prime,
     74        ``None`` is instead returned.
     75
     76        EXAMPLES:
     77
     78        We can find the trace of Frobenius at a given prime:
     79
     80            sage: R.<x> = QQ[]
     81            sage: hec = HyperellipticCurve(x^5 + 2*x^4 - x^3 - 3*x^2 - x)
     82            sage: hec.aplists([17, 17]) # optional: smalljac
     83            [(17, [-9])]
     84
     85        Or in a range:
     86
     87            sage: hec.aplists(100) # optional: smalljac
     88            [(2, [None]), (3, [-1]), (5, [3]), (7, [None]), (11, [-3]), (13, [0]),
     89            (17, [-9]), (19, [-7]), (23, [15]), (29, [-12]), (31, [-5]), (37, [5]),
     90            (41, [0]), (43, [0]), (47, [-3]), (53, [9]), (59, [-9]), (61, [15]),
     91            (67, [9]), (71, [0]), (73, [3]), (79, [-9]), (83, [24]), (89, [-21]),
     92            (97, [0])]
     93            sage: hec.aplists([50, 100]) # optional: smalljac
     94            [(53, [9]), (59, [-9]), (61, [15]), (67, [9]), (71, [0]), (73, [3]),
     95            (79, [-9]), (83, [24]), (89, [-21]), (97, [0])]
     96
     97        The raw results are suitable for large-scale statistics:
     98
     99            sage: __builtins__.sum(ap for p, ap in hec.aplists(10^4, raw=True) if ap) # optional: smalljac
     100            1900
     101            sage: __builtins__.sum(ap for p, ap in hec.aplists(2 * 10^4, raw=True) if ap) # optional: smalljac
     102            1925
     103            sage: __builtins__.sum(ap for p, ap in hec.aplists(3 * 10^4, raw=True) if ap) # optional: smalljac
     104            5451
     105        """
     106
     107        import sage.libs.smalljac.util as util
     108
     109        start, end = util.parse_limits(limit)
     110
     111        if algorithm == "smalljac":
     112            proxy = self._smalljac()
     113            res = proxy.ap(start, end)
     114
     115            def make_int(p, v):
     116                return Integer(v)
     117
     118            return util.interpret_smalljac_results(
     119                res, correction_fn=None, construct_fn=make_int,
     120                raw=raw, index=index, container=container)
     121        else:
     122            raise ValueError("Algorithm %s not defined" % algorithm)
     123
     124    def aps(self, p):
     125        """
     126        Returns a list of traces of Frobenius when this elliptic curve
     127        is reduced at prime ideals of norm ``p``.
     128
     129        INPUT:
     130
     131         - ``p`` - The integer prime at which to reduce.
     132
     133        OUTPUT:
     134
     135        A list of traces of Frobenius, or None if there was a bad
     136        reduction.
     137
     138        EXAMPLES:
     139
     140        We can take a curve and find arbitrary traces::
     141
     142            sage: R.<x> = QQ[]
     143            sage: hec = HyperellipticCurve(x^5 + 2*x^4 - x^3 - 3*x^2 - x)
     144            sage: hec.aps(17) # optional: smalljac
     145            [-9]
     146            sage: hec.aps(next_prime(10^6)) # optional: smalljac
     147            [909]
     148            sage: hec.aps(next_prime(10^8)) # optional: smalljac
     149            [27651]
     150        """
     151
     152        res = self.aplists([p, p])
     153
     154        return res[0][1] # [0] extracts the prime we care about, [1] extracts the list of traces
     155
     156    def grouplists(self, limit, algorithm="smalljac", container=list, index="primes", raw=False):
     157        r"""
     158        Computes the group structure of reductions at the prime ideals
     159        of the base number field, sorted by the norm of these ideals.
     160        So far, this method requires the Smalljac library.
     161
     162        INPUT:
     163
     164        - ``limit`` - either an integer, in which case the prime
     165          ideals will have norm less than this limit, or [a, b], in
     166          which case the prime ideals with norms between ``a`` and
     167          ``b`, inclusive, are used.
     168
     169        - ``algorithm`` - so far, only "smalljac" is supported, which
     170          uses Drew Sutherland's smalljac library for elliptic curve
     171          structure computation.
     172
     173        - The remainder (``container``, ``index``, and ``raw``) define
     174          the format of the result; see ``EllipticCurve.aplists`` for
     175          a description.
     176
     177        OUTPUT:
     178
     179        A list of group structures of Jacobians at various primes; the
     180        format is precisely specified by the ``container``, ``index``,
     181        and ``raw`` arguments.  If a curve has a bad reduction at a prime,
     182        ``None`` is instead returned.
     183
     184        EXAMPLES:
     185
     186        You can construct a hyperelliptic curve and ask for its
     187        Jacobian's group structure at a certain prime::
     188
     189            sage: R.<x> = QQ[]
     190            sage: hec = HyperellipticCurve(x^5 + 2*x^4 - x^3 - 3*x^2 - x)
     191            sage: hec.grouplists([17, 17]) # optional: smalljac
     192            [(17, [Additive abelian group isomorphic to Z/4 + Z/124])]
     193
     194       
     195        Check for larger group lists; 4 elements
     196
     197            sage: hec2 = HyperellipticCurve(x^5 + x)
     198            sage: hec2.grouplists([17, 17]) # optional: smalljac
     199            [(17, [Additive abelian group isomorphic to Z/2 + Z/2 + Z/12 + Z/12])]
     200
     201        Or in a range::
     202
     203            sage: hec.grouplists([80, 100]) # optional: smalljac
     204            [(83, [Additive abelian group isomorphic to Z/2 + Z/2 + Z/36 + Z/36]),
     205             (89, [Additive abelian group isomorphic to Z/8 + Z/1256]),
     206             (97, [Additive abelian group isomorphic to Z/2 + Z/2 + Z/2 + Z/1158])]
     207           
     208        The raw results are suitable for larger-scale statistics::
     209
     210            sage: def jac_checksum(poly):                                         
     211            ...       hec = HyperellipticCurve(poly)
     212            ...       sum = __builtins__.sum                           
     213            ...       return sum(sum(gs) for p, gs in hec.grouplists([3, 10000], raw=True) if gs)
     214            sage: jac_checksum(x^5-x) # optional: smalljac
     215            2083547730
     216            sage: jac_checksum(x^5+5*x^3+5*x) # optional: smalljac
     217            18393758619
     218            sage: jac_checksum(x^5+1) # optional: smalljac
     219            19340969888
     220            sage: jac_checksum(x^5+x^3+2*x) # optional: smalljac
     221            12945578110
     222            sage: jac_checksum(x^5-10*x^2-9*x+2) # optional: smalljac
     223            12917177585
     224            sage: jac_checksum(x^5+3*x^4+4*x^2-9*x+1) # optional: smalljac
     225            12614952542
     226            sage: jac_checksum(x^5-x+1) # optional: smalljac
     227            28462483128
     228        """
     229
     230        from sage.groups.additive_abelian.additive_abelian_group import AdditiveAbelianGroup
     231        import sage.libs.smalljac.util as util
     232
     233        start, end = util.parse_limits(limit)
     234
     235        if algorithm == "smalljac":
     236            proxy = self._smalljac()
     237            res = proxy.group(start, end)
     238
     239            def make_group(p, v):
     240                return AdditiveAbelianGroup(v)
     241
     242            return util.interpret_smalljac_results(
     243                res, correction_fn=None, construct_fn=make_group,
     244                raw=raw, index=index, container=container)
     245        else:
     246            raise ValueError("Algorithm %s not defined" % algorithm)
     247
     248    def lpolylists(self, limit, normalize=False, algorithm="smalljac", container=list, index="primes", raw=False):
     249        r"""
     250        Computes the L-polynomials of reductions at the prime ideals
     251        of the base number field, sorted by the norm of these ideals.
     252        So far, this method requires the Smalljac library.
     253
     254        INPUT:
     255
     256        - ``limit`` - either an integer, in which case the prime
     257          ideals will have norm less than this limit, or ``[a, b]``,
     258          in which case the prime ideals with norms between ``a`` and
     259          ``b`, inclusive, are used.
     260
     261        - ``normalize`` - whether to return normalized L-polynomials.
     262          By default, non-normalized L-polynomials are returned.
     263
     264        - ``algorithm`` - so far, only "smalljac" is supported, which
     265          uses Drew Sutherland's smalljac library for elliptic curve
     266          structure computation.
     267
     268        - The remainder (``container``, ``index``, and ``raw``) define
     269          the format of the result; see ``EllipticCurve.aplists`` for
     270          a description.
     271
     272        OUTPUT:
     273
     274        A list of L-polynomials of reductions at various primes; the
     275        format is precisely specified by the ``container``, ``index``,
     276        and ``raw`` arguments.  If a curve has a bad reduction at a
     277        prime, ``None`` is instead returned.
     278
     279        EXAMPLES:
     280
     281        You can construct a hyperelliptic curve and ask for its
     282        L-polynomials at a certain prime:
     283
     284            sage: R2.<y> = QQ[]
     285            sage: E = HyperellipticCurve(y^5 + 2*y^4 - y^3 - 3*y^2 - y)
     286            sage: E.lpolylists([89, 89]) # optional: smalljac
     287            [(89, [7921*T^4 + 1869*T^3 + 236*T^2 + 21*T + 1])]
     288
     289        Or in a range:
     290
     291            sage: E.lpolylists([80, 90]) # optional: smalljac
     292            [(83, [6889*T^4 - 1992*T^3 + 310*T^2 - 24*T + 1]),
     293             (89, [7921*T^4 + 1869*T^3 + 236*T^2 + 21*T + 1])]
     294
     295        You can also get normalized L-polynomials:
     296
     297            sage: E.lpolylists([80, 100], normalize=True) # optional: smalljac
     298            [(83, [T^4 - 2.63434223975257*T^3 + 3.00000000000000*T^2 -
     299                   2.63434223975257*T + 1.00000000000000]),
     300             (89, [T^4 + 2.22599554801336*T^3 + 2.00000000000000*T^2 +
     301                   2.22599554801336*T + 1.00000000000000]),
     302             (97, [T^4 - 2.00000000000000*T^2 + 1.00000000000000])]
     303        """
     304
     305        from sage.rings.polynomial.polynomial_ring import PolynomialRing_field
     306        from sage.rings.real_mpfr import RR
     307        from sage.rings.rational_field import QQ
     308        from math import sqrt
     309        import sage.libs.smalljac.util as util
     310
     311        start, end = util.parse_limits(limit)
     312
     313        if algorithm == "smalljac":
     314            proxy = self._smalljac()
     315            res = proxy.lpoly(start, end)
     316
     317            if normalize:
     318                polynomial_ring = PolynomialRing_field(RR, 'T')
     319            else:
     320                polynomial_ring = PolynomialRing_field(QQ, 'T')
     321            t = polynomial_ring.gens()[0] # I don't want a T
     322           
     323            def make_lpoly(p, coeffs):
     324                if not normalize:
     325                    a1, a2 = coeffs
     326                    return p**2*t**4 + a1*p*t**3 + a2*t**2 + a1*t + 1
     327                else:
     328                    a1, a2 = coeffs
     329                    a2 /= p
     330                    a1 /= sqrt(p)
     331                    return t**4 + a1*t**3 + a2*t**2 + a1*t + 1
     332
     333            return util.interpret_smalljac_results(
     334                res, correction_fn=None, construct_fn=make_lpoly,
     335                raw=raw, index=index, container=container)
     336        else:
     337            raise ValueError("Algorithm %s not defined" % algorithm)
     338
     339    def trace_histogram(self, limit, buckets=None):
     340        """
     341        Computes the traces of reductions of this curve at all primes
     342        in a range, and returns how many, normalized, fell in each
     343        of ``buckets`` buckets from -2 to 2.
     344
     345        INPUT:
     346
     347        - ``limits`` - Either an integer, meaning all primes below that
     348          integer are used, or a list of ``[start, end]``, in which
     349          case all primes from ``start`` to ``end``, inclusive, are
     350          used.
     351
     352        - ``buckets`` - How many buckets to split the range $[-2, 2]$
     353          into.  By default, the number of buckets is the square root
     354          of the range size.
     355
     356        OUTPUT:
     357
     358        A list of integers, as many as the ``buckets`` parameter,
     359        denoting the number of trace values fell into that bucket for
     360        primes in the given range.
     361
     362        EXAMPLES:
     363
     364        We can take a hyperelliptic curve and compute a histogram for it::
     365
     366            sage: R.<x> = QQ[]
     367            sage: E3 = HyperellipticCurve(x^5 + 2*x^4 - x^3 - 3*x^2 - x)
     368            sage: E3.trace_histogram(100) # optional: smalljac
     369            [0L, 0L, 4L, 4L, 2L, 6L, 4L, 1L, 2L, 0L]
     370
     371        In this case, the default choice of 10 buckets was made.
     372        Instead, you can specify the number of bucket manually::
     373
     374            sage: E3.trace_histogram(100000, buckets=10) # optional: smalljac
     375            [126L, 401L, 722L, 1317L, 1420L, 3098L, 1255L, 722L, 403L, 126L]
     376
     377        You can choose a range that starts anywhere, not just at 1::
     378
     379            sage: E3.trace_histogram([1000, 2000], buckets=10) # optional: smalljac
     380            [6L, 7L, 5L, 17L, 16L, 45L, 23L, 5L, 8L, 3L]
     381
     382        This is generally useful for making histograms; a quick way to
     383        do this is with ``bar_chart``, though custom graphing code
     384        will look better. ::
     385
     386            sage: bar_chart(E3.trace_histogram(100000), width=1)
     387        """
     388
     389        import sage.libs.smalljac.util as util
     390        from math import sqrt
     391
     392        start, end = util.parse_limits(limit)
     393
     394        if buckets is None:
     395            range_size = end - start + 1
     396            buckets = int(sqrt(range_size))
     397
     398        return map(Integer,
     399                   self._smalljac().sato_tate_buckets(buckets, start, end))
     400
     401    def moments(self, limits, moments=11):
     402        """
     403        Compute the moments of the distributions of the coefficients
     404        of the L-polynomials of this curve, reduced at primes in a
     405        range.
     406
     407        INPUT:
     408
     409        - ``limits`` - Either an integer, meaning all primes below that
     410          integer are used, or a list of ``[start, end]``, in which
     411          case all primes from ``start`` to ``end``, inclusive, are
     412          used.
     413
     414        - ``moments`` - How many moments to return; fewer moments
     415          might be faster to compute.  For now, at most 10 moments can
     416          be computed.
     417
     418        OUTPUT:
     419
     420        A tuple whose first element is a list of moments of the
     421        distribution of the $a_1$ coefficients of L-polynomials of
     422        this curve, starting from the zeroth moment being 0, and whose
     423        second element is a similar list for the $a_2$ coefficients.
     424        The list is as long as ``moments``; that is, if ``moments`` is
     425        5, only the zero through fourth moments are returned.
     426
     427        EXAMPLES:
     428
     429        You can request the moments of an arbitrary curve::
     430       
     431            sage: R.<x> = QQ[]
     432            sage: he = HyperellipticCurve(x^5 + 2*x^4 - x^3 - 3*x^2 - x)
     433            sage: he.moments(10000) # optional: smalljac
     434            ([1.0, -0.012809677791304501, 2.0318601033357764,
     435              -0.022492869661964197, 12.332650329303904,
     436              -0.08293711910285241, 104.79583112542458,
     437              1.4982847405159512, 1046.380403999994, 58.11982235337913,
     438              11453.36487528555],
     439             [0.0, 1.0190697892939953, 4.013283796035834,
     440              11.252490719416159, 45.18973261805002,
     441              179.65488417712342, 796.9073238931065,
     442              3636.7510997563845, 17306.79323001631, 84297.73197858929,
     443              419658.09474476246])
     444
     445        You can request fewer moments if you wish::
     446
     447            sage: he.moments(10000, moments=3) # optional: smalljac
     448            ([1.0, -0.012809677791304501, 2.0318601033357764],
     449             [0.0, 1.0190697892939953, 4.013283796035834])
     450
     451        Or compute in any range::
     452
     453            sage: he.moments([1000, 2000], moments=5) # optional: smalljac
     454            ([1.0, 0.012019115191429789, 2.360758977244879,
     455              0.7305838739637142, 17.033614832007466],
     456             [0.0, 1.0981854467458285, 4.588416007826706,
     457              14.81342939310661, 62.822935207983335])
     458        """
     459
     460        import sage.libs.smalljac.util as util
     461
     462        start, end = util.parse_limits(limits)
     463        if moments > 11:
     464            raise ValueError("Only up to 11 moments supported for now.")
     465
     466        return self._smalljac().moments(moments, start, end)
     467
     468    def guess_sato_tate_group(self, limit=None):
     469        """
     470        Provide a preliminary guess toward the Sato-Tate group
     471        corresponding to this curve, by considering its reduction at
     472        primes up to a certain limit.
     473
     474        INPUT:
     475
     476         - ``limit`` - The maximal prime considered.  The search may
     477           terminate earlier if an answer is found.  Or, ``None``, to
     478           use an implementation-defined (but very large) maximum.
     479
     480        OUTPUT:
     481
     482        A string containing the name of the corresponding Sato-Tate
     483        group.  To interpret these, a copy of "Sato-Tate distributions
     484        and Galois endomorphism modules in genus 2", [fkrs]_, may be
     485        helpful.  If no provisional identification can be made,
     486        returns ``None``
     487
     488        .. TODO:
     489           - Should output a textual description of the Sato-Tate group
     490
     491        EXAMPLES:
     492
     493        Just create a hyperelliptic curve of genus 2 and call
     494        ``guess_sato_tate_group`` to start it off::
     495
     496            sage: R.<x> = QQ[]
     497            sage: he = HyperellipticCurve(x^5 + 2*x^4 - x^3 - 3*x^2 - x)
     498            sage: he.guess_sato_tate_group(100000) # long time - 4 seconds, optional: smalljac
     499            'E_6'
     500            sage: he = HyperellipticCurve(x^5 + x)
     501            sage: he.guess_sato_tate_group(100000) # optional: smalljac
     502            'D_{2,1}'
     503            sage: he = HyperellipticCurve(x^5 + 2*x)
     504            sage: he.guess_sato_tate_group(100000) # optional: smalljac
     505            'D_{4,1}'
     506        """
     507
     508        if limit is None:
     509            limit = 0
     510
     511        return self._smalljac().guess_group(limit)

  • setup.py 

    diff --git a/setup.py b/setup.py
    a b  
    944944                     'sage.libs.cremona',
    945945                     'sage.libs.mpmath',
    946946                     'sage.libs.lcalc',
     947                     'sage.libs.smalljac',
    947948
    948949                     'sage.logic',
    949950