Author: Andrew Sutherland
1#################################################################################
2#
3# (c) Copyright 2011 William Stein
4#
5#  This file is part of PSAGE
6#
7#  PSAGE is free software: you can redistribute it and/or modify
9#  the Free Software Foundation, either version 3 of the License, or
10#  (at your option) any later version.
11#
12#  PSAGE is distributed in the hope that it will be useful,
13#  but WITHOUT ANY WARRANTY; without even the implied warranty of
14#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15#  GNU General Public License for more details.
16#
17#  You should have received a copy of the GNU General Public License
18#  along with this program.  If not, see <http://www.gnu.org/licenses/>.
19#
20#################################################################################
21
22"""
23Fast prime ideals of the ring R of integers of Q(sqrt(5)).
24
25This module implements Cython classes for prime ideals of R, and their
26enumeration.  The main entry function is primes_of_bounded_norm::
27
28    sage: from psage.number_fields.sqrt5 import primes_of_bounded_norm
29    sage: v = primes_of_bounded_norm(50); v
30    [2, 5a, 3, 11a, 11b, 19a, 19b, 29a, 29b, 31a, 31b, 41a, 41b, 7]
31
32Note that the output of primes_of_bounded_norm is a list.  Each entry
33is a prime ideal, which prints using a simple label consisting of the
34characteristic of the prime then "a" or "b", where "b" only appears for
35the second split prime.::
36
37    sage: type(v)
38    <type 'psage.number_fields.sqrt5.prime.Prime'>
39    sage: v.sage_ideal()
40    Fractional ideal (a + 5)
41
42AUTHOR:
43    - William Stein
44"""
45
46
47include "stdsage.pxi"
48include "interrupt.pxi"
49
50cdef extern from "pari/pari.h":
51    unsigned long Fl_sqrt(unsigned long, unsigned long)
52    unsigned long Fl_div(unsigned long, unsigned long, unsigned long)
53
54from sage.rings.number_field.number_field_ideal import NumberFieldFractionalIdeal
55
56cdef class Prime:
57    """
58    Nonzero prime ideal of the ring of integers of Q(sqrt(5)).  This
59    is a fast customized Cython class; to get at the corresponding
60    Sage prime ideal use the sage_ideal method.
61    """
62    def __init__(self, p, long r=0, bint first=True):
63        """
64        Create Prime ideal with given residue characteristic, root,
65        and first or not with that characterstic.
66
67        INPUT form 1:
68            - p -- prime
69            - r -- root (or 0)
70            - first -- boolean: True if first prime over p
71
72        INPUT form 2:
73            - p -- prime ideal of integers of Q(sqrt(5)); validity of the
74              input is not checked in any way!
75
76        NOTE: No checking is done to verify that the input is valid.
77
78        EXAMPLES::
79
80            sage: from psage.number_fields.sqrt5.prime import Prime
81            sage: P = Prime(2,0,True); P
82            2
83            sage: type(P)
84            <type 'psage.number_fields.sqrt5.prime.Prime'>
85            sage: P = Prime(5,3,True); P
86            5a
87            sage: P = Prime(11,8,True); P
88            11a
89            sage: P = Prime(11,4,False); P
90            11b
91
92        We can also set P using a prime ideal of the ring of integers::
93
94            sage: from psage.number_fields.sqrt5.prime import Prime; K.<a> = NumberField(x^2-x-1)
95            sage: P1 = Prime(K.primes_above(11)); P2 = Prime(K.primes_above(11)); P1, P2
96            (11b, 11a)
97            sage: P1 > P2
98            True
99            sage: Prime(K.prime_above(2))
100            2
101            sage: P = Prime(K.prime_above(5)); P, P.r
102            (5a, 3)
103            sage: Prime(K.prime_above(3))
104            3
105
106        Test that conversion both ways works for primes up to norm 10^5::
107
108            sage: from psage.number_fields.sqrt5.prime import primes_of_bounded_norm, Prime
109            sage: v = primes_of_bounded_norm(10^5)
110            sage: w = [Prime(z.sage_ideal()) for z in v]
111            sage: v == w
112            True
113        """
114        cdef long t, r1
115        if isinstance(p, NumberFieldFractionalIdeal):
116            # Set self using a prime ideal of Q(sqrt(5)).
117            H = p.pari_hnf()
118            self.p = H[0,0]
119            self.first = True
120            t = self.p % 5
121            if t == 1 or t == 4:
122                self.r = self.p - H[0,1]
123                r1 = self.p + 1 - self.r
124                if self.r > r1:
125                    self.first = False
126            elif t == 0:
127                self.r = 3
128            else:
129                self.r = 0
130        else:
131            self.p = p
132            self.r = r
133            self.first = first
134
135    def __repr__(self):
136        """
137        EXAMPLES::
138
139            sage: from psage.number_fields.sqrt5.prime import Prime
140            sage: Prime(11,4,False).__repr__()
141            '11b'
142            sage: Prime(11,4,True).__repr__()
143            '11a'
144            sage: Prime(7,0,True).__repr__()
145            '7'
146        """
147        if self.r:
148            return '%s%s'%(self.p, 'a' if self.first else 'b')
149        return '%s'%self.p
150
151    def __hash__(self):
152        return self.p*(self.r+1) + int(self.first)
153
154    def _latex_(self):
155        """
156        EXAMPLES::
157
158            sage: from psage.number_fields.sqrt5.prime import Prime
159            sage: Prime(11,8,True)._latex_()
160            '11a'
161            sage: Prime(11,4,False)._latex_()
162            '11b'
163        """
164        return self.__repr__()
165
166    cpdef bint is_split(self):
167        """
168        Return True if this prime is split (and not ramified).
169
170        EXAMPLES::
171
172            sage: from psage.number_fields.sqrt5.prime import Prime
173            sage: Prime(11,8,True).is_split()
174            True
175            sage: Prime(3,0,True).is_split()
176            False
177            sage: Prime(5,3,True).is_split()
178            False
179        """
180        return self.r != 0 and self.p != 5
181
182    cpdef bint is_inert(self):
183        """
184        Return True if this prime is inert.
185
186        EXAMPLES::
187
188            sage: from psage.number_fields.sqrt5.prime import Prime
189            sage: Prime(11,8,True).is_inert()
190            False
191            sage: Prime(3,0,True).is_inert()
192            True
193            sage: Prime(5,3,True).is_inert()
194            False
195        """
196        return self.r == 0
197
198    cpdef bint is_ramified(self):
199        """
200        Return True if this prime is ramified (i.e., the prime over 5).
201
202        EXAMPLES::
203
204            sage: from psage.number_fields.sqrt5.prime import Prime
205            sage: Prime(11,8,True).is_ramified()
206            False
207            sage: Prime(3,0,True).is_ramified()
208            False
209            sage: Prime(5,3,True).is_ramified()
210            True
211        """
212        return self.p == 5
213
214    cpdef long norm(self):
215        """
216        Return the norm of this ideal.
217
218        EXAMPLES::
219
220            sage: from psage.number_fields.sqrt5.prime import Prime
221            sage: Prime(11,4,True).norm()
222            11
223            sage: Prime(7,0,True).norm()
224            49
225        """
226        return self.p if self.r else self.p*self.p
227
228    def __cmp__(self, Prime right):
229        """
230        Compare two prime ideals. First sort by the norm, then in the
231        (only remaining) split case if the norms are the same, compare
232        the residue of (1+sqrt(5))/2 in the interval [0,p).
233
234        WARNING: The ordering is NOT the same as the ordering of
235        fractional ideals in Sage.
236
237        EXAMPLES::
238
239            sage: from psage.number_fields.sqrt5 import primes_of_bounded_norm
240            sage: v = primes_of_bounded_norm(50); v
241            [2, 5a, 3, 11a, 11b, 19a, 19b, 29a, 29b, 31a, 31b, 41a, 41b, 7]
242            sage: v, v
243            (11a, 11b)
244            sage: v < v
245            True
246            sage: v > v
247            True
248
249        We test the ordering a bit by sorting::
250
251            sage: v.sort(); v
252            [2, 5a, 3, 11a, 11b, 19a, 19b, 29a, 29b, 31a, 31b, 41a, 41b, 7]
253            sage: v = list(reversed(v)); v
254            [7, 41b, 41a, 31b, 31a, 29b, 29a, 19b, 19a, 11b, 11a, 3, 5a, 2]
255            sage: v.sort(); v
256            [2, 5a, 3, 11a, 11b, 19a, 19b, 29a, 29b, 31a, 31b, 41a, 41b, 7]
257
258        A bigger test::
259
260            sage: v = primes_of_bounded_norm(10^7)
261            sage: w = list(reversed(v)); w.sort()
262            sage: v == w
263            True
264        """
265        cdef long selfn = self.norm(),  rightn = right.norm()
266        if   selfn > rightn: return 1
267        elif rightn > selfn: return -1
268        elif self.r > right.r: return 1
269        elif right.r > self.r: return -1
270        else: return 0
271
272    def sage_ideal(self):
273        """
274        Return the usual prime fractional ideal associated to this
275        prime.  This is slow, but provides substantial additional
276        functionality.
277
278        EXAMPLES::
279
280            sage: from psage.number_fields.sqrt5 import primes_of_bounded_norm
281            sage: v = primes_of_bounded_norm(20)
282            sage: v.sage_ideal()
283            Fractional ideal (2*a - 1)
284            sage: [P.sage_ideal() for P in v]
285            [Fractional ideal (2), Fractional ideal (2*a - 1), Fractional ideal (3), Fractional ideal (3*a - 1), Fractional ideal (3*a - 2), Fractional ideal (-4*a + 1), Fractional ideal (-4*a + 3)]
286        """
287        from misc import F
288        cdef long p=self.p, r=self.r
289        if r: # split and ramified cases
290            return F.ideal(p, F.gen()-r)
291        else: # inert case
292            return F.ideal(p)
293
294from sage.rings.integer import Integer
295
296def primes_above(long p, bint check=True):
297    """
298    Return ordered list of all primes above p in the ring of integers
299    of Q(sqrt(5)).  See the docstring for primes_of_bounded_norm.
300
301    INPUT:
302        - p -- prime number in integers ZZ (less than 2^{31})
303        - check -- bool (default: True); if True, check that p is prime
304    OUTPUT:
305        - list of 1 or 2 Prime objects
306
307    EXAMPLES::
308
309        sage: from psage.number_fields.sqrt5.prime import primes_above
310        sage: primes_above(2)
311        
312        sage: primes_above(3)
313        
314        sage: primes_above(5)
315        [5a]
316        sage: primes_above(11)
317        [11a, 11b]
318        sage: primes_above(13)
319        
320        sage: primes_above(17)
321        
322        sage: primes_above(4)
323        Traceback (most recent call last):
324        ...
325        ValueError: p must be a prime
326        sage: primes_above(4, check=False)
327        
328    """
329    if check and not Integer(p).is_pseudoprime():
330        raise ValueError, "p must be a prime"
331    cdef long t = p%5
332    if t == 1 or t == 4 or t == 0:
333        return prime_range(p, p+1)
334    else: # inert
335        return [Prime(p, 0, True)]
336
337def primes_of_bounded_norm(bound):
338    """
339    Return ordered list of all prime ideals of the ring of integers of
340    Q(sqrt(5)) of norm less than bound.
341
342    The primes are instances of a special fast Primes class (they are
343    *not* usual Sage prime ideals -- use the sage_ideal() method to
344    get those).  They are sorted first by norm, then in the remaining
345    split case by the integer in the interval [0,p) congruent to
346    (1+sqrt(5))/2.
347
348    INPUT:
349        - bound -- nonnegative integer, less than 2^{31}
350
351    WARNING: The ordering is NOT the same as the ordering of primes by
352    Sage.   Even if you order first by norm, then use Sage's ordering
353    for primes of the same norm, then the orderings do not agree.::
354
355    EXAMPLES::
356
357        sage: from psage.number_fields.sqrt5 import primes_of_bounded_norm
358        sage: primes_of_bounded_norm(0)
359        []
360        sage: primes_of_bounded_norm(10)
361        [2, 5a, 3]
362        sage: primes_of_bounded_norm(50)
363        [2, 5a, 3, 11a, 11b, 19a, 19b, 29a, 29b, 31a, 31b, 41a, 41b, 7]
364        sage: len(primes_of_bounded_norm(10^6))
365        78510
366
367    We grab one of the primes::
368
369        sage: v = primes_of_bounded_norm(100)
370        sage: P = v; type(P)
371        <type 'psage.number_fields.sqrt5.prime.Prime'>
372
373    It prints with a nice label::
374
375        sage: P
376        11a
377
378    You can get the corresponding fractional ideal as a normal Sage ideal::
379
380        sage: P.sage_ideal()
381        Fractional ideal (3*a - 1)
382
383    You can also get the underlying residue characteristic::
384
385        sage: P.p
386        11
387
388    And, the image of (1+sqrt(5))/2 modulo the prime (or 0 in the inert case)::
389
390        sage: P.r
391        4
392        sage: z = P.sage_ideal(); z.residue_field()(z.number_field().gen())
393        4
394
395    Prime enumeration is reasonable fast, even when the input is
396    relatively large (going up to 10^8 takes a few seconds, and up
397    to 10^9 takes a few minutes), and the following should take less
398    than a second::
399
400        sage: len(primes_of_bounded_norm(10^7))  # less than a second
401        664500
402
403    One limitation is that the bound must be less than 2^{31}::
404
405        sage: primes_of_bounded_norm(2^31)
406        Traceback (most recent call last):
407        ...
408        ValueError: bound must be less than 2^31
409    """
410    return prime_range(bound)
411
412def prime_range(long start, stop=None):
413    """
414    Return ordered list of all prime ideals of the ring of integers of
415    Q(sqrt(5)) of norm at least start and less than stop.  If only
416    start is given then return primes with norm less than start.
417
418    The primes are instances of a special fast Primes class (they are
419    *not* usual Sage prime ideals -- use the sage_ideal() method to
420    get those).  They are sorted first by norm, then in the remaining
421    split case by the integer in the interval [0,p) congruent to
422    (1+sqrt(5))/2.  For optimal speed you can use the Prime objects
424    underlying data structure.
425
426    INPUT:
427        - start -- nonnegative integer, less than 2^{31}
428        - stop  -- None or nonnegative integer, less than 2^{31}
429
430    EXAMPLES::
431
432        sage: from psage.number_fields.sqrt5.prime import prime_range
433        sage: prime_range(10, 60)
434        [11a, 11b, 19a, 19b, 29a, 29b, 31a, 31b, 41a, 41b, 7, 59a, 59b]
435        sage: prime_range(2, 11)
436        [2, 5a, 3]
437        sage: prime_range(2, 12)
438        [2, 5a, 3, 11a, 11b]
439        sage: prime_range(3, 12)
440        [2, 5a, 3, 11a, 11b]
441        sage: prime_range(9, 12)
442        [3, 11a, 11b]
443        sage: prime_range(5, 12)
444        [5a, 3, 11a, 11b]
445    """
446    if start >= 2**31 or (stop and stop >= 2**31):
447        raise ValueError, "bound must be less than 2^31"
448
449    cdef long p, p2, sr, r0, r1, t, bound
450    cdef Prime P
451    cdef list v = []
452
453    if stop is None:
454        bound = start
455        start = 2
456    else:
457        bound = stop
458
459    from sage.all import prime_range as prime_range_ZZ
460
461    for p in prime_range_ZZ(bound, py_ints=True):
462        t = p % 5
463        if t == 1 or t == 4:   # split
464            if p >= start:
465                # Compute a square root of 5 modulo p.
466                sr = Fl_sqrt(5, p)
467                # Find the two values of (1+sqrt(5))/2.
468                r0 = Fl_div(1+sr, 2, p)
469                r1 = p+1-r0
470                # Sort
471                if r0 > r1: r0, r1 = r1, r0
472                # Append each prime to the list
473                P = PY_NEW(Prime); P.p = p; P.r = r0; P.first = True; v.append(P)
474                P = PY_NEW(Prime); P.p = p; P.r = r1; P.first = False; v.append(P)
475        elif p == 5:   # ramified
476            if p >= start:
477                v.append(Prime(p, 3, True))
478        else:
479            p2 = p*p
480            if p2 < bound and p2 >= start:
481                v.append(Prime(p, 0, True))
482
483    v.sort()
484    return v
485
486
487
488
489
490