Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
sagemath
GitHub Repository: sagemath/sagelib
Path: blob/master/sage/libs/mwrank/mwrank.pyx
4097 views
1
"""
2
Cython interface to Cremona's ``eclib`` library (also known as ``mwrank``)
3
4
EXAMPLES::
5
6
sage: from sage.libs.mwrank.mwrank import _Curvedata, _mw
7
sage: c = _Curvedata(1,2,3,4,5)
8
9
sage: print c
10
[1,2,3,4,5]
11
b2 = 9 b4 = 11 b6 = 29 b8 = 35
12
c4 = -183 c6 = -3429
13
disc = -10351 (# real components = 1)
14
#torsion not yet computed
15
16
sage: t= _mw(c)
17
sage: t.search(10)
18
sage: t
19
[[1:2:1]]
20
"""
21
22
import os
23
import sys
24
25
include '../../ext/interrupt.pxi'
26
include '../../ext/stdsage.pxi'
27
28
# Need to permit tabs in order to doctest verbose output.
29
"""
30
SAGE_DOCTEST_ALLOW_TABS
31
"""
32
33
cdef extern from "wrap.h":
34
### misc functions ###
35
long mwrank_get_precision()
36
void mwrank_set_precision(long n)
37
void mwrank_initprimes(char* pfilename, int verb)
38
39
40
### bigint ###
41
struct bigint
42
bigint* new_bigint()
43
void del_bigint(bigint* x)
44
bigint* str_to_bigint(char* s)
45
char* bigint_to_str(bigint* x)
46
47
### Curvedata ###
48
struct Curvedata
49
Curvedata* Curvedata_new(bigint* a1, bigint* a2,
50
bigint* a3, bigint* a4,
51
bigint* a6, int min_on_init)
52
void Curvedata_del(Curvedata* curve)
53
char* Curvedata_repr(Curvedata* curve)
54
double Curvedata_silverman_bound(Curvedata* curve)
55
double Curvedata_cps_bound(Curvedata* curve)
56
double Curvedata_height_constant(Curvedata* curve)
57
char* Curvedata_getdiscr(Curvedata* curve)
58
char* Curvedata_conductor(Curvedata* m)
59
char* Curvedata_isogeny_class(Curvedata* E, int verbose)
60
61
## mw ##
62
struct mw
63
mw* mw_new(Curvedata* curve, int verb, int pp, int maxr)
64
void mw_del(mw* m)
65
int mw_process(Curvedata* curve, mw* m,
66
bigint* x, bigint* y,
67
bigint* z, int sat)
68
char* mw_getbasis(mw* m)
69
char* mw_regulator(mw* m)
70
int mw_rank(mw* m)
71
int mw_saturate(mw* m, bigint* index, char** unsat,
72
long sat_bd, int odd_primes_only)
73
void mw_search(mw* m, char* h_lim, int moduli_option, int verb)
74
75
### two_descent ###
76
struct two_descent
77
two_descent* two_descent_new(Curvedata* curve,
78
int verb, int sel,
79
long firstlim, long secondlim,
80
long n_aux, int second_descent)
81
82
void two_descent_del(two_descent* t)
83
int two_descent_ok(two_descent* t)
84
long two_descent_get_certain(two_descent* t)
85
char* two_descent_get_basis(two_descent* t)
86
char* two_descent_regulator(two_descent* t)
87
long two_descent_get_rank(two_descent* t)
88
long two_descent_get_rank_bound(two_descent* t)
89
long two_descent_get_selmer_rank(two_descent* t)
90
void two_descent_saturate(two_descent* t, long sat_bd)
91
92
93
94
cdef object string_sigoff(char* s):
95
sig_off()
96
# Makes a python string and deletes what is pointed to by s.
97
t = str(s)
98
sage_free(s)
99
return t
100
101
class __init:
102
pass
103
104
_INIT = __init()
105
# set the default
106
mwrank_set_precision(50)
107
108
def get_precision():
109
"""
110
Returns the working floating point precision of mwrank.
111
112
OUTPUT:
113
114
(int) The current precision in decimal digits.
115
116
EXAMPLE::
117
118
sage: from sage.libs.mwrank.mwrank import get_precision
119
sage: get_precision()
120
50
121
122
"""
123
return mwrank_get_precision()
124
125
def set_precision(n):
126
"""
127
Sets the working floating point precision of mwrank.
128
129
INPUT:
130
131
- ``n`` (int) -- a positive integer: the number of decimal digits.
132
133
OUTPUT:
134
135
None.
136
137
EXAMPLE::
138
139
sage: from sage.libs.mwrank.mwrank import set_precision
140
sage: set_precision(50)
141
142
"""
143
mwrank_set_precision(n)
144
145
def initprimes(filename, verb=False):
146
"""
147
Initialises mwrank/eclib's internal prime list.
148
149
INPUT:
150
151
- ``filename`` (string) -- the name of a file of primes.
152
153
- ``verb`` (bool: default ``False``) -- verbose or not?
154
155
EXAMPLES::
156
157
sage: file= SAGE_TMP + '/PRIMES'
158
sage: open(file,'w').write(' '.join([str(p) for p in prime_range(10^7,10^7+20)]))
159
sage: mwrank_initprimes(file, verb=True)
160
161
The previous command does produce the following output when run
162
from the command line, but not during a doctest::
163
164
Computed 78519 primes, largest is 1000253
165
reading primes from file ...
166
read extra prime 10000019
167
finished reading primes from file ...
168
Extra primes in list: 10000019
169
170
sage: mwrank_initprimes("x" + file, True)
171
Traceback (most recent call last):
172
...
173
IOError: No such file or directory: ...
174
"""
175
if not os.path.exists(filename):
176
raise IOError, 'No such file or directory: %s'%filename
177
mwrank_initprimes(filename, verb)
178
if verb:
179
sys.stdout.flush()
180
sys.stderr.flush()
181
182
############# bigint ###########################################
183
#
184
# In mwrank (and eclib) bigint is synonymous with NTL's ZZ class.
185
#
186
################################################################
187
188
cdef class _bigint:
189
"""
190
Cython class wrapping eclib's bigint class.
191
"""
192
cdef bigint* x
193
194
def __init__(self, x="0"):
195
"""
196
Constructor for bigint class.
197
198
INPUT:
199
200
- ``x`` -- string or int: a string representing a decimal
201
integer, or a Sage integer
202
203
EXAMPLES::
204
205
sage: from sage.libs.mwrank.mwrank import _bigint
206
sage: _bigint(123)
207
123
208
sage: _bigint('123')
209
123
210
sage: type(_bigint(123))
211
<type 'sage.libs.mwrank.mwrank._bigint'>
212
"""
213
if not (x is _INIT):
214
s = str(x)
215
if s.isdigit() or s[0] == "-" and s[1:].isdigit():
216
self.x = str_to_bigint(s)
217
else:
218
raise ValueError, "invalid _bigint: %s"%x
219
220
def __dealloc__(self):
221
"""
222
Destructor for bigint class (releases memory).
223
"""
224
del_bigint(self.x)
225
226
def __repr__(self):
227
"""
228
String representation of bigint.
229
230
OUTPUT:
231
232
(string) the bigint as a string (decimal)
233
234
EXAMPLES::
235
236
sage: from sage.libs.mwrank.mwrank import _bigint
237
sage: a = _bigint('123')
238
sage: a.__repr__()
239
'123'
240
sage: a = _bigint('-456')
241
sage: a.__repr__()
242
'-456'
243
"""
244
sig_on()
245
return string_sigoff(bigint_to_str(self.x))
246
247
248
cdef make_bigint(bigint* x):
249
cdef _bigint y
250
sig_off()
251
y = _bigint(_INIT)
252
y.x = x
253
return y
254
255
############# Curvedata #################
256
257
cdef class _Curvedata: # cython class wrapping eclib's Curvedata class
258
cdef Curvedata* x
259
260
def __init__(self, a1, a2, a3,
261
a4, a6, min_on_init=0):
262
"""
263
Constructor for Curvedata class.
264
265
INPUT:
266
267
- ``a1``,``a2``,``a3``,``a4``,``a6`` (int) -- integer
268
coefficients of a Weierstrass equation (must be
269
nonsingular).
270
271
- ``min_on_init`` (int, default 0) -- flag controlling whether
272
the constructed curve is replaced by a global minimal model.
273
If nonzero then this minimisation does take place.
274
275
EXAMPLES::
276
277
sage: from sage.libs.mwrank.mwrank import _Curvedata
278
sage: _Curvedata(1,2,3,4,5)
279
[1,2,3,4,5]
280
b2 = 9 b4 = 11 b6 = 29 b8 = 35
281
c4 = -183 c6 = -3429
282
disc = -10351 (# real components = 1)
283
#torsion not yet computed
284
285
A non-minimal example::
286
287
sage: _Curvedata(0,0,0,0,64)
288
[0,0,0,0,64]
289
b2 = 0 b4 = 0 b6 = 256 b8 = 0
290
c4 = 0 c6 = -55296
291
disc = -1769472 (# real components = 1)
292
#torsion not yet computed
293
294
sage: _Curvedata(0,0,0,0,64,min_on_init=1)
295
[0,0,0,0,1] (reduced minimal model)
296
b2 = 0 b4 = 0 b6 = 4 b8 = 0
297
c4 = 0 c6 = -864
298
disc = -432 (# real components = 1)
299
#torsion not yet computed
300
"""
301
cdef _bigint _a1, _a2, _a3, _a4, _a6
302
_a1 = _bigint(a1)
303
_a2 = _bigint(a2)
304
_a3 = _bigint(a3)
305
_a4 = _bigint(a4)
306
_a6 = _bigint(a6)
307
self.x = Curvedata_new(_a1.x, _a2.x, _a3.x, _a4.x, _a6.x, min_on_init)
308
if self.discriminant() == 0:
309
raise ArithmeticError, "Invariants (= %s) do not describe an elliptic curve."%([a1,a2,a3,a4,a6])
310
311
def __dealloc__(self):
312
"""
313
Destructor for Curvedata class.
314
"""
315
Curvedata_del(self.x)
316
317
def __repr__(self):
318
"""
319
String representation of Curvedata
320
321
OUTPUT:
322
323
(string) the Curvedata as a string
324
325
EXAMPLES::
326
327
sage: from sage.libs.mwrank.mwrank import _Curvedata
328
sage: E = _Curvedata(1,2,3,4,5)
329
sage: E.__repr__()
330
'[1,2,3,4,5]\nb2 = 9\t b4 = 11\t b6 = 29\t b8 = 35\nc4 = -183\t\tc6 = -3429\ndisc = -10351\t(# real components = 1)\n#torsion not yet computed'
331
sage: E
332
[1,2,3,4,5]
333
b2 = 9 b4 = 11 b6 = 29 b8 = 35
334
c4 = -183 c6 = -3429
335
disc = -10351 (# real components = 1)
336
#torsion not yet computed
337
"""
338
sig_on()
339
return string_sigoff(Curvedata_repr(self.x))[:-1]
340
341
def silverman_bound(self):
342
"""
343
The Silverman height bound for this elliptic curve.
344
345
OUTPUT:
346
347
(float) A non-negative real number `B` such that for every
348
rational point on this elliptic curve `E`, `h(P)\le\hat{h}(P)
349
+ B`, where `h(P)` is the naive height and `\hat{h}(P)` the
350
canonical height.
351
352
TODO:
353
354
Since eclib can compute this to arbitrary precision it would
355
make sense to return a Sage real.
356
357
EXAMPLES::
358
359
sage: from sage.libs.mwrank.mwrank import _Curvedata
360
sage: E = _Curvedata(1,2,3,4,5)
361
sage: E.silverman_bound()
362
6.52226179519101...
363
sage: type(E.silverman_bound())
364
<type 'float'>
365
"""
366
return Curvedata_silverman_bound(self.x)
367
368
def cps_bound(self):
369
"""
370
The Cremona-Prickett-Siksek height bound for this elliptic curve.
371
372
OUTPUT:
373
374
(float) A non-negative real number `B` such that for every
375
rational point on this elliptic curve `E`, `h(P)\le\hat{h}(P)
376
+ B`, where `h(P)` is the naive height and `\hat{h}(P)` the
377
canonical height.
378
379
TODO:
380
381
Since eclib can compute this to arbitrary precision it would
382
make sense to return a Sage real.
383
384
EXAMPLES::
385
386
sage: from sage.libs.mwrank.mwrank import _Curvedata
387
sage: E = _Curvedata(1,2,3,4,5)
388
sage: E.cps_bound()
389
0.11912451909250982
390
391
Note that this is a better bound than Silverman's in this case::
392
393
sage: E.silverman_bound()
394
6.52226179519101...
395
"""
396
cdef double x
397
# We declare x so there are *no* Python library
398
# calls within the sig_on()/sig_off().
399
sig_on()
400
x = Curvedata_cps_bound(self.x)
401
sig_off()
402
return x
403
404
def height_constant(self):
405
"""
406
A height bound for this elliptic curve.
407
408
OUTPUT:
409
410
(float) A non-negative real number `B` such that for every
411
rational point on this elliptic curve `E`, `h(P)\le\hat{h}(P)
412
+ B`, where `h(P)` is the naive height and `\hat{h}(P)` the
413
canonical height. This is the minimum of the Silverman and
414
Cremona_Prickett-Siksek height bounds.
415
416
TODO:
417
418
Since eclib can compute this to arbitrary precision it would
419
make sense to return a Sage real.
420
421
EXAMPLES::
422
423
sage: from sage.libs.mwrank.mwrank import _Curvedata
424
sage: E = _Curvedata(1,2,3,4,5)
425
sage: E.height_constant()
426
0.11912451909250982
427
"""
428
return Curvedata_height_constant(self.x)
429
430
def discriminant(self):
431
"""
432
The discriminant of this elliptic curve.
433
434
OUTPUT:
435
436
(Integer) The discriminant.
437
438
EXAMPLES::
439
440
sage: from sage.libs.mwrank.mwrank import _Curvedata
441
sage: E = _Curvedata(1,2,3,4,5)
442
sage: E.discriminant()
443
-10351
444
sage: E = _Curvedata(100,200,300,400,500)
445
sage: E.discriminant()
446
-1269581104000000
447
sage: ZZ(E.discriminant())
448
-1269581104000000
449
"""
450
sig_on()
451
from sage.rings.all import Integer
452
return Integer(string_sigoff(Curvedata_getdiscr(self.x)))
453
454
def conductor(self):
455
"""
456
The conductor of this elliptic curve.
457
458
OUTPUT:
459
460
(Integer) The conductor.
461
462
463
EXAMPLES::
464
465
sage: from sage.libs.mwrank.mwrank import _Curvedata
466
sage: E = _Curvedata(1,2,3,4,5)
467
sage: E.discriminant()
468
-10351
469
sage: E = _Curvedata(100,200,300,400,500)
470
sage: E.conductor()
471
126958110400
472
"""
473
sig_on()
474
from sage.rings.all import Integer
475
return Integer(string_sigoff(Curvedata_conductor(self.x)))
476
477
def isogeny_class(self, verbose=False):
478
"""
479
The isogeny class of this elliptic curve.
480
481
OUTPUT:
482
483
(tuple) A tuple consisting of (1) a list of the curves in the
484
isogeny class, each as a list of its Weierstrass coefficients;
485
(2) a matrix of the degrees of the isogenies between the
486
curves in the class (prime degrees only).
487
488
.. warning::
489
490
The list may not be complete, if the precision is too low.
491
Use ``mwrank_set_precision()`` to increase the precision.
492
493
EXAMPLES::
494
495
sage: from sage.libs.mwrank.mwrank import _Curvedata
496
sage: E = _Curvedata(1,0,1,4,-6)
497
sage: E.conductor()
498
14
499
sage: E.isogeny_class()
500
([[1, 0, 1, 4, -6], [1, 0, 1, -36, -70], [1, 0, 1, -1, 0], [1, 0, 1, -171, -874], [1, 0, 1, -11, 12], [1, 0, 1, -2731, -55146]], [[0, 2, 3, 3, 0, 0], [2, 0, 0, 0, 3, 3], [3, 0, 0, 0, 2, 0], [3, 0, 0, 0, 0, 2], [0, 3, 2, 0, 0, 0], [0, 3, 0, 2, 0, 0]])
501
502
"""
503
sig_on()
504
s = string_sigoff(Curvedata_isogeny_class(self.x, verbose))
505
if verbose:
506
sys.stdout.flush()
507
sys.stderr.flush()
508
from sage.misc.all import preparse
509
from sage.rings.all import Integer
510
return eval(s)
511
512
513
############# _mw #################
514
515
cdef class _mw:
516
"""
517
Cython class wrapping eclib's mw class.
518
"""
519
cdef mw* x
520
cdef Curvedata* curve
521
cdef int verb
522
523
def __init__(self, _Curvedata curve, verb=False, pp=1, maxr=999):
524
"""
525
Constructor for mw class.
526
527
INPUT:
528
529
- ``curve`` (_Curvedata) -- an elliptic curve
530
531
- ``verb`` (bool, default False) -- verbosity flag (controls
532
amount of output produced in point searches)
533
534
- ``pp`` (int, default 1) -- process points flag (if nonzero,
535
the points found are processed, so that at all times only a
536
`\ZZ`-basis for the subgroup generated by the points found
537
so far is stored; if zero, no processing is done and all
538
points found are stored).
539
540
- ``maxr`` (int, default 999) -- maximum rank (quit point
541
searching once the points found generate a subgroup of this
542
rank; useful if an upper bound for the rank is already
543
known).
544
545
EXAMPLE::
546
547
sage: from sage.libs.mwrank.mwrank import _mw
548
sage: from sage.libs.mwrank.mwrank import _Curvedata
549
sage: E = _Curvedata(1,0,1,4,-6)
550
sage: EQ = _mw(E)
551
sage: EQ
552
[]
553
sage: type(EQ)
554
<type 'sage.libs.mwrank.mwrank._mw'>
555
556
sage: E = _Curvedata(0,0,1,-7,6)
557
sage: EQ = _mw(E)
558
sage: EQ.search(2)
559
sage: EQ
560
[[1:-1:1], [-2:3:1], [-14:25:8]]
561
562
Example to illustrate the verbose parameter::
563
564
sage: from sage.libs.mwrank.mwrank import _mw
565
sage: from sage.libs.mwrank.mwrank import _Curvedata
566
sage: E = _Curvedata(0,0,1,-7,6)
567
sage: EQ = _mw(E, verb=False)
568
sage: EQ.search(1)
569
sage: EQ = _mw(E, verb=True)
570
sage: EQ.search(1)
571
572
The previous command produces the following output::
573
574
P1 = [0:1:0] is torsion point, order 1
575
P1 = [-3:0:1] is generator number 1
576
saturating up to 20...Checking 2-saturation
577
Points have successfully been 2-saturated (max q used = 7)
578
Checking 3-saturation
579
Points have successfully been 3-saturated (max q used = 7)
580
Checking 5-saturation
581
Points have successfully been 5-saturated (max q used = 23)
582
Checking 7-saturation
583
Points have successfully been 7-saturated (max q used = 41)
584
Checking 11-saturation
585
Points have successfully been 11-saturated (max q used = 17)
586
Checking 13-saturation
587
Points have successfully been 13-saturated (max q used = 43)
588
Checking 17-saturation
589
Points have successfully been 17-saturated (max q used = 31)
590
Checking 19-saturation
591
Points have successfully been 19-saturated (max q used = 37)
592
done
593
P2 = [-2:3:1] is generator number 2
594
saturating up to 20...Checking 2-saturation
595
possible kernel vector = [1,1]
596
This point may be in 2E(Q): [14:-52:1]
597
...and it is!
598
Replacing old generator #1 with new generator [1:-1:1]
599
Points have successfully been 2-saturated (max q used = 7)
600
Index gain = 2^1
601
Checking 3-saturation
602
Points have successfully been 3-saturated (max q used = 13)
603
Checking 5-saturation
604
Points have successfully been 5-saturated (max q used = 67)
605
Checking 7-saturation
606
Points have successfully been 7-saturated (max q used = 53)
607
Checking 11-saturation
608
Points have successfully been 11-saturated (max q used = 73)
609
Checking 13-saturation
610
Points have successfully been 13-saturated (max q used = 103)
611
Checking 17-saturation
612
Points have successfully been 17-saturated (max q used = 113)
613
Checking 19-saturation
614
Points have successfully been 19-saturated (max q used = 47)
615
done (index = 2).
616
Gained index 2, new generators = [ [1:-1:1] [-2:3:1] ]
617
P3 = [-14:25:8] is generator number 3
618
saturating up to 20...Checking 2-saturation
619
Points have successfully been 2-saturated (max q used = 11)
620
Checking 3-saturation
621
Points have successfully been 3-saturated (max q used = 13)
622
Checking 5-saturation
623
Points have successfully been 5-saturated (max q used = 71)
624
Checking 7-saturation
625
Points have successfully been 7-saturated (max q used = 101)
626
Checking 11-saturation
627
Points have successfully been 11-saturated (max q used = 127)
628
Checking 13-saturation
629
Points have successfully been 13-saturated (max q used = 151)
630
Checking 17-saturation
631
Points have successfully been 17-saturated (max q used = 139)
632
Checking 19-saturation
633
Points have successfully been 19-saturated (max q used = 179)
634
done (index = 1).
635
P4 = [-1:3:1] = -1*P1 + -1*P2 + -1*P3 (mod torsion)
636
P4 = [0:2:1] = 2*P1 + 0*P2 + 1*P3 (mod torsion)
637
P4 = [2:13:8] = -3*P1 + 1*P2 + -1*P3 (mod torsion)
638
P4 = [1:0:1] = -1*P1 + 0*P2 + 0*P3 (mod torsion)
639
P4 = [2:0:1] = -1*P1 + 1*P2 + 0*P3 (mod torsion)
640
P4 = [18:7:8] = -2*P1 + -1*P2 + -1*P3 (mod torsion)
641
P4 = [3:3:1] = 1*P1 + 0*P2 + 1*P3 (mod torsion)
642
P4 = [4:6:1] = 0*P1 + -1*P2 + -1*P3 (mod torsion)
643
P4 = [36:69:64] = 1*P1 + -2*P2 + 0*P3 (mod torsion)
644
P4 = [68:-25:64] = -2*P1 + -1*P2 + -2*P3 (mod torsion)
645
P4 = [12:35:27] = 1*P1 + -1*P2 + -1*P3 (mod torsion)
646
sage: EQ
647
[[1:-1:1], [-2:3:1], [-14:25:8]]
648
649
Example to illustrate the process points ``pp`` parameter::
650
651
sage: from sage.libs.mwrank.mwrank import _mw
652
sage: from sage.libs.mwrank.mwrank import _Curvedata
653
sage: E = _Curvedata(0,0,1,-7,6)
654
sage: EQ = _mw(E, pp=1)
655
sage: EQ.search(1); EQ
656
[[1:-1:1], [-2:3:1], [-14:25:8]]
657
sage: EQ = _mw(E, pp=0)
658
sage: EQ.search(1); EQ
659
[[-3:0:1], [-2:3:1], [-14:25:8], [-1:3:1], [0:2:1], [2:13:8], [1:0:1], [2:0:1], [18:7:8], [3:3:1], [4:6:1], [36:69:64], [68:-25:64], [12:35:27]]
660
661
"""
662
self.curve = curve.x
663
self.x = mw_new(curve.x, verb, pp, maxr)
664
self.verb = verb
665
666
def __dealloc__(self):
667
"""
668
Destructor for mw class.
669
"""
670
mw_del(self.x)
671
672
def __repr__(self):
673
"""
674
String representation of the current basis of this mw group.
675
676
OUTPUT:
677
678
(string) the current basis of this mw as a string
679
680
EXAMPLES::
681
682
sage: from sage.libs.mwrank.mwrank import _Curvedata
683
sage: from sage.libs.mwrank.mwrank import _mw
684
sage: E = _Curvedata(0,0,1,-7,6)
685
sage: EQ = _mw(E)
686
sage: EQ # indirect doctest
687
[]
688
sage: EQ.search(2)
689
sage: EQ
690
[[1:-1:1], [-2:3:1], [-14:25:8]]
691
"""
692
sig_on()
693
return string_sigoff(mw_getbasis(self.x))
694
695
696
def process(self, point, sat=0):
697
"""
698
Processes the given point, adding it to the mw group.
699
700
INPUT:
701
702
- ``point`` (tuple or list) -- tuple or list of 3 integers.
703
An ``ArithmeticError`` is raised if the point is not on the
704
curve.
705
706
- ``sat`` (int, default 0) --saturate at primes up to ``sat``.
707
No saturation is done if ``sat=0``. (Note that it is more
708
efficient to add several points at once and then saturate
709
just once at the end).
710
711
.. note::
712
713
The eclib function which implements this only carries out
714
any saturation if the rank of the points increases upon
715
adding the new point. This is because it is assumed that
716
one saturates as ones goes along.
717
718
EXAMPLES::
719
720
sage: from sage.libs.mwrank.mwrank import _Curvedata
721
sage: from sage.libs.mwrank.mwrank import _mw
722
sage: E = _Curvedata(0,1,1,-2,0)
723
sage: EQ = _mw(E)
724
725
Initially the list of gens is empty::
726
727
sage: EQ
728
[]
729
730
We process a point of infinite order::
731
732
sage: EQ.process([-1,1,1])
733
sage: EQ
734
[[-1:1:1]]
735
736
And another independent one::
737
738
sage: EQ.process([0,-1,1])
739
sage: EQ
740
[[-1:1:1], [0:-1:1]]
741
742
Processing a point dependent on the current basis will not
743
change the basis::
744
745
sage: EQ.process([4,8,1])
746
sage: EQ
747
[[-1:1:1], [0:-1:1]]
748
749
"""
750
if not isinstance(point, (tuple, list)) and len(point) == 3:
751
raise TypeError, "point must be a list or tuple of length 3."
752
cdef _bigint x,y,z
753
sig_on()
754
x,y,z = _bigint(point[0]), _bigint(point[1]), _bigint(point[2])
755
r = mw_process(self.curve, self.x, x.x, y.x, z.x, sat)
756
sig_off()
757
if r != 0:
758
raise ArithmeticError, "point (=%s) not on curve."%point
759
760
def getbasis(self):
761
"""
762
Returns the current basis of the mw structure.
763
764
OUTPUT:
765
766
(string) String representation of the points in the basis of
767
the mw group.
768
769
EXAMPLES::
770
771
sage: from sage.libs.mwrank.mwrank import _Curvedata
772
sage: from sage.libs.mwrank.mwrank import _mw
773
sage: E = _Curvedata(0,1,1,-2,0)
774
sage: EQ = _mw(E)
775
sage: EQ.search(3)
776
sage: EQ.getbasis()
777
'[[0:-1:1], [-1:1:1]]'
778
sage: EQ.rank()
779
2
780
"""
781
sig_on()
782
s = string_sigoff(mw_getbasis(self.x))
783
return s
784
785
def regulator(self):
786
"""
787
Returns the regulator of the current basis of the mw group.
788
789
OUTPUT:
790
791
(float) The current regulator.
792
793
TODO:
794
795
``eclib`` computes the regulator to arbitrary precision, and
796
the full precision value should be returned.
797
798
EXAMPLES::
799
800
sage: from sage.libs.mwrank.mwrank import _Curvedata
801
sage: from sage.libs.mwrank.mwrank import _mw
802
sage: E = _Curvedata(0,1,1,-2,0)
803
sage: EQ = _mw(E)
804
sage: EQ.search(3)
805
sage: EQ.getbasis()
806
'[[0:-1:1], [-1:1:1]]'
807
sage: EQ.rank()
808
2
809
sage: EQ.regulator()
810
0.15246017277240753
811
"""
812
cdef float f
813
sig_on()
814
f = float(string_sigoff(mw_regulator(self.x)))
815
return f
816
817
def rank(self):
818
"""
819
Returns the rank of the current basis of the mw group.
820
821
OUTPUT:
822
823
(Integer) The current rank.
824
825
EXAMPLES::
826
827
sage: from sage.libs.mwrank.mwrank import _Curvedata
828
sage: from sage.libs.mwrank.mwrank import _mw
829
sage: E = _Curvedata(0,1,1,-2,0)
830
sage: EQ = _mw(E)
831
sage: EQ.search(3)
832
sage: EQ.getbasis()
833
'[[0:-1:1], [-1:1:1]]'
834
sage: EQ.rank()
835
2
836
"""
837
sig_on()
838
r = mw_rank(self.x)
839
sig_off()
840
from sage.rings.all import Integer
841
return Integer(r)
842
843
def saturate(self, int sat_bd=-1, int odd_primes_only=0):
844
"""
845
Saturates the current subgroup of the mw group.
846
847
INPUT:
848
849
- ``sat_bnd`` (int, default -1) -- bound on primes at which to
850
saturate. If -1 (default), compute a bound for the primes
851
which may not be saturated, and use that.
852
853
- ``odd_primes_only`` (bool, default False) -- only do
854
saturation at odd primes. (If the points have been found
855
via 2-descent they should already be 2-saturated.)
856
857
OUTPUT:
858
859
(tuple) (success flag, index, list) The success flag will be 1
860
unless something failed (usually an indication that the points
861
were not saturated but the precision is not high enough to
862
divide out successfully). The index is the index of the mw
863
group before saturation in the mw group after. The list is a
864
string representation of the primes at which saturation was
865
not proved or achieved.
866
867
EXAMPLES::
868
869
sage: from sage.libs.mwrank.mwrank import _Curvedata
870
sage: from sage.libs.mwrank.mwrank import _mw
871
sage: E = _Curvedata(0,1,1,-2,0)
872
sage: EQ = _mw(E)
873
sage: EQ.process([494, -5720, 6859]) # 3 times another point
874
sage: EQ
875
[[494:-5720:6859]]
876
sage: EQ.saturate()
877
(1, 3, '[ ]')
878
sage: EQ
879
[[-1:1:1]]
880
881
If we set the saturation bound at 2, then saturation will fail::
882
883
sage: EQ = _mw(E)
884
sage: EQ.process([494, -5720, 6859]) # 3 times another point
885
sage: EQ.saturate(sat_bd=2)
886
(0, 1, '[ ]')
887
sage: EQ
888
[[494:-5720:6859]]
889
890
The following output is also seen in the preceding example::
891
892
Saturation index bound = 10
893
WARNING: saturation at primes p > 2 will not be done;
894
points may be unsaturated at primes between 2 and index bound
895
Failed to saturate MW basis at primes [ ]
896
897
898
"""
899
cdef _bigint index
900
cdef char* s
901
cdef int ok
902
sig_on()
903
index = _bigint()
904
ok = mw_saturate(self.x, index.x, &s, sat_bd, odd_primes_only)
905
unsat = string_sigoff(s)
906
return ok, index, unsat
907
908
def search(self, h_lim, int moduli_option=0, int verb=0):
909
"""
910
Search for points in the mw group.
911
912
INPUT:
913
914
- ``h_lim`` (int) -- bound on logarithmic naive height of points
915
916
- ``moduli_option`` (int, default 0) -- option for sieving
917
strategy. The default (0) uses an adapted version of
918
Stoll's ratpoints code and is recommended.
919
920
- ``verb`` (int, default 0) -- level of verbosity. If 0, no
921
output. If positive, the points are output as found and
922
some details of the processing, finding linear relations,
923
and partial saturation are output.
924
925
.. note::
926
927
The effect of the search is also governed by the class
928
options, notably whether the points found are processed:
929
meaning that linear relations are found and saturation is
930
carried out, with the result that the list of generators
931
will always contain a `\ZZ`-span of the saturation of the
932
points found, modulo torsion.
933
934
OUTPUT:
935
936
None. The effect of the search is to update the list of
937
generators.
938
939
EXAMPLE::
940
941
sage: from sage.libs.mwrank.mwrank import _Curvedata
942
sage: from sage.libs.mwrank.mwrank import _mw
943
sage: E = _Curvedata(0,0,1,-19569,-4064513) # 873c1
944
sage: EQ = _mw(E)
945
sage: EQ = _mw(E)
946
sage: for i in [1..11]: print i, EQ.search(i), EQ
947
1 None []
948
2 None []
949
3 None []
950
4 None []
951
5 None []
952
6 None []
953
7 None []
954
8 None []
955
9 None []
956
10 None []
957
11 None [[3639568:106817593:4096]]
958
"""
959
cdef char* _h_lim
960
961
h_lim = str(h_lim)
962
_h_lim = h_lim
963
964
sig_on()
965
mw_search(self.x, _h_lim, moduli_option, verb)
966
if verb:
967
sys.stdout.flush()
968
sys.stderr.flush()
969
sig_off()
970
971
972
############# two_descent #################
973
cdef class _two_descent:
974
"""
975
Cython class wrapping eclib's two_descent class.
976
"""
977
cdef two_descent* x
978
979
def __init__(self):
980
"""
981
Constructor for two_descent class.
982
983
EXAMPLES::
984
985
sage: from sage.libs.mwrank.mwrank import _two_descent
986
sage: D2 = _two_descent()
987
"""
988
self.x = <two_descent*> 0
989
990
def __dealloc__(self):
991
"""
992
Destructor for two_descent class.
993
"""
994
if self.x:
995
two_descent_del(self.x)
996
997
def do_descent(self, _Curvedata curve,
998
int verb = 1,
999
int sel = 0,
1000
int firstlim = 20,
1001
int secondlim = 8,
1002
int n_aux = -1,
1003
int second_descent = 1):
1004
"""
1005
Carry out a 2-descent.
1006
1007
INPUT:
1008
1009
- ``curvedata`` (_Curvedata) -- the curve on which to do descent.
1010
1011
- ``verb`` (int, default 1) -- verbosity level.
1012
1013
- ``sel`` (int, default 0) -- Selmer-only flag. If 1, only
1014
the 2-Selmer group will be computed, with no rational
1015
points. Useful as a faster way of getting an upper bound on
1016
the rank.
1017
1018
- ``firstlim`` (int, default 20) -- naive height bound on
1019
first point search on quartic homogeneous spaces (before
1020
testing local solubility; very simple search with no
1021
overheads).
1022
1023
- ``secondlim`` (int, default 8) -- naive height bound on
1024
second point search on quartic homogeneous spaces (after
1025
testing local solubility; sieve-assisted search)
1026
1027
- ``n_aux`` (int, default -1) -- If positive, the number of
1028
auxiliary primes used in sieve-assisted search for quartics.
1029
If -1 (the default) use a default value (set in the eclib
1030
code in ``src/qrank/mrank1.cc`` in DEFAULT_NAUX: currently 8).
1031
Only relevant for curves with no 2-torsion, where full
1032
2-descent is carried out. Worth increasing for curves
1033
expected to be of of rank>6 to one or two more than the
1034
expected rank.
1035
1036
- ``second_descent`` (int, default 1) -- flag specifying
1037
whether or not a second descent will be carried out (yes if
1038
1, the default; no if 0). Only relevant for curves with
1039
2-torsion. Recommended left as the default except for
1040
experts interested in details of Selmer groups.
1041
1042
OUTPUT:
1043
1044
None
1045
1046
EXAMPLES::
1047
1048
sage: from sage.libs.mwrank.mwrank import _Curvedata
1049
sage: CD = _Curvedata(0,0,1,-7,6)
1050
sage: from sage.libs.mwrank.mwrank import _two_descent
1051
sage: D2 = _two_descent()
1052
sage: D2.do_descent(CD)
1053
sage: D2.getrank()
1054
3
1055
sage: D2.getcertain()
1056
1
1057
sage: D2.ok()
1058
1
1059
"""
1060
sig_on()
1061
self.x = two_descent_new(curve.x, verb, sel, firstlim, secondlim, n_aux, second_descent)
1062
if verb:
1063
sys.stdout.flush()
1064
sys.stderr.flush()
1065
sig_off()
1066
1067
def getrank(self):
1068
"""
1069
Returns the rank (after doing a 2-descent).
1070
1071
OUTPUT:
1072
1073
(Integer) the rank (or an upper bound).
1074
1075
EXAMPLES::
1076
1077
sage: from sage.libs.mwrank.mwrank import _Curvedata
1078
sage: CD = _Curvedata(0,0,1,-7,6)
1079
sage: from sage.libs.mwrank.mwrank import _two_descent
1080
sage: D2 = _two_descent()
1081
sage: D2.do_descent(CD)
1082
sage: D2.getrank()
1083
3
1084
"""
1085
cdef int r
1086
sig_on()
1087
r = two_descent_get_rank(self.x)
1088
sig_off()
1089
from sage.rings.all import Integer
1090
return Integer(r)
1091
1092
def getrankbound(self):
1093
"""
1094
Returns the rank upper bound (after doing a 2-descent).
1095
1096
OUTPUT:
1097
1098
(Integer) an upper bound on the rank.
1099
1100
EXAMPLES::
1101
1102
sage: from sage.libs.mwrank.mwrank import _Curvedata
1103
sage: CD = _Curvedata(0,0,1,-7,6)
1104
sage: from sage.libs.mwrank.mwrank import _two_descent
1105
sage: D2 = _two_descent()
1106
sage: D2.do_descent(CD)
1107
sage: D2.getrankbound()
1108
3
1109
"""
1110
cdef int r
1111
sig_on()
1112
r = two_descent_get_rank_bound(self.x)
1113
sig_off()
1114
from sage.rings.all import Integer
1115
return Integer(r)
1116
1117
def getselmer(self):
1118
"""
1119
Returns the 2-Selmer rank (after doing a 2-descent).
1120
1121
OUTPUT:
1122
1123
(Integer) The 2-Selmer rank.
1124
1125
EXAMPLES::
1126
1127
sage: from sage.libs.mwrank.mwrank import _Curvedata
1128
sage: CD = _Curvedata(0,0,1,-7,6)
1129
sage: from sage.libs.mwrank.mwrank import _two_descent
1130
sage: D2 = _two_descent()
1131
sage: D2.do_descent(CD)
1132
sage: D2.getselmer()
1133
3
1134
"""
1135
sig_on()
1136
r = two_descent_get_selmer_rank(self.x)
1137
sig_off()
1138
from sage.rings.all import Integer
1139
return Integer(r)
1140
1141
def ok(self):
1142
"""
1143
Returns the success flag (after doing a 2-descent).
1144
1145
OUTPUT:
1146
1147
(bool) Flag indicating whether or not 2-descent was successful.
1148
1149
EXAMPLES::
1150
1151
sage: from sage.libs.mwrank.mwrank import _Curvedata
1152
sage: CD = _Curvedata(0,0,1,-7,6)
1153
sage: from sage.libs.mwrank.mwrank import _two_descent
1154
sage: D2 = _two_descent()
1155
sage: D2.do_descent(CD)
1156
sage: D2.ok()
1157
1
1158
"""
1159
return two_descent_ok(self.x)
1160
1161
def getcertain(self):
1162
"""
1163
Returns the certainty flag (after doing a 2-descent).
1164
1165
OUTPUT:
1166
1167
(bool) True if the rank upper and lower bounds are equal.
1168
1169
EXAMPLES::
1170
1171
sage: from sage.libs.mwrank.mwrank import _Curvedata
1172
sage: CD = _Curvedata(0,0,1,-7,6)
1173
sage: from sage.libs.mwrank.mwrank import _two_descent
1174
sage: D2 = _two_descent()
1175
sage: D2.do_descent(CD)
1176
sage: D2.getcertain()
1177
1
1178
"""
1179
return two_descent_get_certain(self.x)
1180
1181
def saturate(self, saturation_bound=0):
1182
"""
1183
Carries out saturation of the points found by a 2-descent.
1184
1185
OUTPUT:
1186
1187
None.
1188
1189
EXAMPLES::
1190
1191
sage: from sage.libs.mwrank.mwrank import _Curvedata
1192
sage: CD = _Curvedata(0,0,1,-7,6)
1193
sage: from sage.libs.mwrank.mwrank import _two_descent
1194
sage: D2 = _two_descent()
1195
sage: D2.do_descent(CD)
1196
sage: D2.saturate()
1197
sage: D2.getbasis()
1198
'[[1:-1:1], [-2:3:1], [-14:25:8]]'
1199
"""
1200
sig_on()
1201
two_descent_saturate(self.x, saturation_bound)
1202
sig_off()
1203
1204
def getbasis(self):
1205
"""
1206
Returns the basis of points found by doing a 2-descent.
1207
1208
If the success and certain flags are 1, this will be a
1209
`\ZZ/2\ZZ`-basis for `E(\QQ)/2E(\QQ)` (modulo torsion),
1210
otherwise possibly only for a proper subgroup.
1211
1212
.. note::
1213
1214
You must call ``saturate()`` first, or a RunTimeError will be raised.
1215
1216
OUTPUT:
1217
1218
(string) String representation of the list of points after
1219
saturation.
1220
1221
EXAMPLES::
1222
1223
sage: from sage.libs.mwrank.mwrank import _Curvedata
1224
sage: CD = _Curvedata(0,0,1,-7,6)
1225
sage: from sage.libs.mwrank.mwrank import _two_descent
1226
sage: D2 = _two_descent()
1227
sage: D2.do_descent(CD)
1228
sage: D2.saturate()
1229
sage: D2.getbasis()
1230
'[[1:-1:1], [-2:3:1], [-14:25:8]]'
1231
"""
1232
sig_on()
1233
return string_sigoff(two_descent_get_basis(self.x))
1234
1235
def regulator(self):
1236
"""
1237
Returns the regulator of the points found by doing a 2-descent.
1238
1239
OUTPUT:
1240
1241
(float) The regulator (of the subgroup found by 2-descent).
1242
1243
EXAMPLES::
1244
1245
sage: from sage.libs.mwrank.mwrank import _Curvedata
1246
sage: CD = _Curvedata(0,0,1,-7,6)
1247
sage: from sage.libs.mwrank.mwrank import _two_descent
1248
sage: D2 = _two_descent()
1249
sage: D2.do_descent(CD)
1250
1251
If called before calling ``saturate()``, a bogus value of 1.0
1252
is returned::
1253
1254
sage: D2.regulator()
1255
1.0
1256
1257
After saturation, both ``getbasis()`` and ``regulator()``
1258
return the basis and regulator of the subgroup found by
1259
2-descent::
1260
1261
sage: D2.saturate()
1262
sage: D2.getbasis()
1263
'[[1:-1:1], [-2:3:1], [-14:25:8]]'
1264
sage: D2.regulator()
1265
0.417143558758384
1266
"""
1267
sig_on()
1268
return float(string_sigoff(two_descent_regulator(self.x)))
1269
1270