CoCalc Shared Filescode / prime.pyxOpen in CoCalc with one click!
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
8
# it under the terms of the GNU General Public License as published by
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
"""
23
Fast prime ideals of the ring R of integers of Q(sqrt(5)).
24
25
This module implements Cython classes for prime ideals of R, and their
26
enumeration. 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
32
Note that the output of primes_of_bounded_norm is a list. Each entry
33
is a prime ideal, which prints using a simple label consisting of the
34
characteristic of the prime then "a" or "b", where "b" only appears for
35
the second split prime.::
36
37
sage: type(v[8])
38
<type 'psage.number_fields.sqrt5.prime.Prime'>
39
sage: v[8].sage_ideal()
40
Fractional ideal (a + 5)
41
42
AUTHOR:
43
- William Stein
44
"""
45
46
47
include "stdsage.pxi"
48
include "interrupt.pxi"
49
50
cdef 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
54
from sage.rings.number_field.number_field_ideal import NumberFieldFractionalIdeal
55
56
cdef 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)[0]); P2 = Prime(K.primes_above(11)[1]); 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[3], v[4]
243
(11a, 11b)
244
sage: v[3] < v[4]
245
True
246
sage: v[4] > v[3]
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[1].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
294
from sage.rings.integer import Integer
295
296
def 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
[2]
312
sage: primes_above(3)
313
[3]
314
sage: primes_above(5)
315
[5a]
316
sage: primes_above(11)
317
[11a, 11b]
318
sage: primes_above(13)
319
[13]
320
sage: primes_above(17)
321
[17]
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
[2]
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
337
def 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[3]; 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
412
def 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
423
directly from Cython, which provides direct C-level access to the
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