CoCalc -- ternary.pyx

GitHub Repository: sagemath/sagesmc
Path: blob/master/src/sage/quadratic_forms/ternary.pyx
⁸⁸¹⁷ views
1
#*****************************************************************************
2
#       Copyright (C) 2012 Gustavo Rama
3
#
4
#  Distributed under the terms of the GNU General Public License (GPL)
5
#
6
#    This code is distributed in the hope that it will be useful,
7
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
8
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
9
#    General Public License for more details.
10
#
11
#  The full text of the GPL is available at:
12
#
13
#                  http://www.gnu.org/licenses/
14
#*****************************************************************************
15

16

17

18
from sage.rings.integer_ring import ZZ
19
from sage.matrix.constructor import matrix, identity_matrix, diagonal_matrix
20
from sage.modules.free_module_element import vector
21
from sage.rings.arith import inverse_mod
22
from sage.quadratic_forms.extras import extend_to_primitive
23
from sage.rings.finite_rings.integer_mod import mod
24
from sage.misc.prandom import randint
25
from sage.rings.arith import xgcd,gcd
26
from sage.functions.other import ceil, floor
27
from __builtin__ import max
28

29

30

31
def red_mfact(a,b):
32
    """
33
    Auxiliar function for reduction that finds the reduction factor of a, b integers.
34

35
    INPUT:
36
        - a, b integers
37

38
    OUTPUT:
39
        Integer
40

41
    EXAMPLES::
42

43
        sage: from sage.quadratic_forms.ternary import red_mfact
44
        sage: red_mfact(0, 3)
45
        0
46
        sage: red_mfact(-5, 100)
47
        9
48

49
    """
50

51
    if a:
52
      return (-b + abs(a))//(2*a)
53
    else:
54
      return 0
55

56
def _reduced_ternary_form_eisenstein_with_matrix(a1, a2, a3, a23, a13, a12):
57
    """
58
    Find the coefficients of the equivalent unique reduced ternary form according to the conditions
59
    of Dickson's "Studies in the Theory of Numbers", pp164-171, and the tranformation matrix.
60
    See TernaryQF.is_eisenstein_reduced for the conditions.
61

62
    EXAMPLES::
63

64
        sage: from sage.quadratic_forms.ternary import _reduced_ternary_form_eisenstein_with_matrix
65
        sage: Q = TernaryQF([293, 315, 756, 908, 929, 522])
66
        sage: qr, M = _reduced_ternary_form_eisenstein_with_matrix(293, 315, 756, 908, 929, 522)
67
        sage: qr
68
        (1, 2, 2, -1, 0, -1)
69
        sage: M
70
        [ -54  137  -38]
71
        [ -23   58  -16]
72
        [  47 -119   33]
73
        sage: Qr = TernaryQF(qr)
74
        sage: Qr.is_eisenstein_reduced()
75
        True
76
        sage: Q(M) == Qr
77
        True
78

79

80
    """
81

82
    M = identity_matrix(3)
83

84
    loop=1
85

86
    while loop:
87

88

89
        # adjust
90
        v=a1+a2+a23+a13+a12
91
        if (v<0):
92
            M*=matrix(ZZ,3,[1,0,1,0,1,1,0,0,1])
93
            a3+=v
94
            a23+=a12+2*a2
95
            a13+=a12+2*a1
96

97
        # cuadred 12
98
        m=red_mfact(a1,a12)
99
        M*=matrix(ZZ,3,[1,m,0,0,1,0,0,0,1])
100
        t=a1*m
101
        a12+=t
102
        a2+=a12*m
103
        a12+=t
104
        a23+=a13*m
105

106
        # cuadred 23
107
        m=red_mfact(a2,a23)
108
        M*=matrix(ZZ,3,[1,0,0,0,1,m,0,0,1])
109
        t=a2*m
110
        a23+=t
111
        a3+=a23*m
112
        a23+=t
113
        a13+=a12*m
114

115
        # cuadred 13
116
        m=red_mfact(a1,a13)
117
        M*=matrix(ZZ,3,[1,0,m,0,1,0,0,0,1])
118
        t=a1*m
119
        a13+=t
120
        a3+=a13*m
121
        a13+=t
122
        a23+=a12*m
123

124
        # order 12
125
        if a1 > a2 or (a1 == a2 and abs(a23) > abs(a13)):
126
            M*=matrix(ZZ,3,[0,-1,0,-1,0,0,0,0,-1])
127
            [a1,a2]=[a2,a1]
128
            [a13,a23]=[a23,a13]
129

130
        # order 23
131
        if a2 > a3 or (a2 == a3 and abs(a13) > abs(a12)):
132
            M*=matrix(ZZ,3,[-1,0,0,0,0,-1,0,-1,0])
133
            [a2,a3]=[a3,a2]
134
            [a13,a12]=[a12,a13]
135

136
        # order 12
137
        if a1 > a2 or (a1 == a2 and abs(a23) > abs(a13)):
138
            M*=matrix(ZZ,3,[0,-1,0,-1,0,0,0,0,-1])
139
            [a1,a2]=[a2,a1]
140
            [a13,a23]=[a23,a13]
141

142
        # signs
143
        if (a23*a13*a12>0):
144
            # a23, a13, a12 positive
145

146
            if (a23<0):
147
                M*=diagonal_matrix([-1,1,1])
148
                a23=-a23
149
            if (a13<0):
150
                M*=diagonal_matrix([1,-1,1])
151
                a13=-a13
152
            if (a12<0):
153
                M*=diagonal_matrix([1,1,-1])
154
                a12=-a12
155

156
        else:
157
            # a23, a13, a12 nonpositive
158

159
            [s1,s2,s3]=[(a23>0),(a13>0),(a12>0)]
160
            if ((s1+s2+s3)%2):
161
                if (a23==0):
162
                    s1=1
163
                else:
164
                    if (a13==0):
165
                        s2=1
166
                    else:
167
                        if (a12==0):
168
                            s3=1
169
            if s1:
170
                  M*=diagonal_matrix([-1,1,1])
171
                  a23=-a23
172
            if s2:
173
                  M*=diagonal_matrix([1,-1,1])
174
                  a13=-a13
175
            if s3:
176
                  M*=diagonal_matrix([1,1,-1])
177
                  a12=-a12
178

179
        loop = not (abs(a23) <= a2 and abs(a13) <= a1 and abs(a12) <= a1 and a1+a2+a23+a13+a12>=0)
180

181

182
################
183
################
184

185

186

187
    # adj 3
188
    if a1+a2+a23+a13+a12 == 0 and 2*a1+2*a13+a12 > 0:
189
        M*=matrix(ZZ,3,[-1,0,1,0,-1,1,0,0,1])
190
        # a3 += a1+a2+a23+a13+a12
191
        a23=-2*a2-a23-a12
192
        a13=-2*a1-a13-a12
193

194
    # adj 5.12
195
    if a1 == -a12 and a13 != 0:
196
        M*=matrix(ZZ,3,[-1,-1,0,0,-1,0,0,0,1])
197
        #a2 += a1+a12
198
        a23=-a23-a13
199
        a13=-a13
200
        a12=-a12  # = 2*a1+a12
201

202
    # adj 5.13
203
    if a1 == -a13 and a12 != 0:
204
        M*=matrix(ZZ,3,[-1,0,-1,0,1,0,0,0,-1])
205
        # a3 += a1+a13
206
        a23=-a23-a12
207
        a13=-a13  # = 2*a1+a13
208
        a12=-a12
209

210
    # adj 5.23
211
    if a2 == -a23 and a12 != 0:
212
        M*=matrix(ZZ,3,[1,0,0,0,-1,-1,0,0,-1])
213
        # a3 += a2+a23
214
        a23=-a23  # = 2*a2+a23
215
        a13=-a13-a12
216
        a12=-a12
217

218
    # adj 4.12
219
    if a1 == a12 and a13 > 2*a23:
220
        M*=matrix(ZZ,3,[-1,-1,0,0,1,0,0,0,-1])
221
        # a 2 += a1-a12
222
        a23 = -a23 + a13
223
        # a12 = 2*a1 - a12
224

225
    # adj 4.13
226
    if a1 == a13 and a12 > 2*a23:
227
        M*=matrix(ZZ,3,[-1,0,-1,0,-1,0,0,0,1])
228
        # a3 += a1-a13
229
        a23 = -a23 + a12
230
        # a13 = 2*a1 - a13
231

232
    # adj 4.23
233
    if a2 == a23 and a12 > 2*a13:
234
        M*=matrix(ZZ,3,[-1,0,0,0,-1,-1,0,0,1])
235
        # a3 += a2-a23
236
        # a23 = 2*a2 - a23
237
        a13 = -a13 + a12
238

239
    # order 12
240
    if a1 == a2 and abs(a23) > abs(a13):
241
        M*=matrix(ZZ,3,[0,-1,0,-1,0,0,0,0,-1])
242
        [a1,a2]=[a2,a1]
243
        [a13,a23]=[a23,a13]
244

245
    # order 23
246
    if a2 == a3 and abs(a13) > abs(a12):
247
        M*=matrix(ZZ,3,[-1,0,0,0,0,-1,0,-1,0])
248
        [a13,a12]=[a12,a13]
249

250
    # order 12
251
    if a1 == a2 and abs(a23) > abs(a13):
252
        M*=matrix(ZZ,3,[0,-1,0,-1,0,0,0,0,-1])
253
        [a13,a23]=[a23,a13]
254

255
    return((a1,a2,a3,a23,a13,a12),M)
256

257
def _reduced_ternary_form_eisenstein_without_matrix(a1, a2, a3, a23, a13, a12):
258
    """
259
    Find the coefficients of the equivalent unique reduced ternary form according to the conditions
260
    of Dickson's "Studies in the Theory of Numbers", pp164-171.
261
    See TernaryQF.is_eisenstein_reduced for the conditions.
262

263
    EXAMPLES::
264

265
        sage: from sage.quadratic_forms.ternary import _reduced_ternary_form_eisenstein_without_matrix
266
        sage: Q = TernaryQF([293, 315, 756, 908, 929, 522])
267
        sage: qr = _reduced_ternary_form_eisenstein_without_matrix(293, 315, 756, 908, 929, 522)
268
        sage: qr
269
        (1, 2, 2, -1, 0, -1)
270
        sage: Qr = TernaryQF(qr)
271
        sage: Qr.is_eisenstein_reduced()
272
        True
273

274
    """
275

276
    loop=1
277

278
    while loop:
279

280
        # adjust
281
        v=a1+a2+a23+a13+a12
282
        if (v<0):
283
            a3+=v
284
            a23+=a12+2*a2
285
            a13+=a12+2*a1
286

287
        # cuadred 12
288
        m=red_mfact(a1,a12)
289
        t=a1*m
290
        a12+=t
291
        a2+=a12*m
292
        a12+=t
293
        a23+=a13*m
294

295
        # cuadred 23
296
        m=red_mfact(a2,a23)
297
        t=a2*m
298
        a23+=t
299
        a3+=a23*m
300
        a23+=t
301
        a13+=a12*m
302

303
        # cuadred 13
304
        m=red_mfact(a1,a13)
305
        t=a1*m
306
        a13+=t
307
        a3+=a13*m
308
        a13+=t
309
        a23+=a12*m
310

311
        # order 12
312
        if a1 > a2 or (a1 == a2 and abs(a23) > abs(a13)):
313
            [a1,a2]=[a2,a1]
314
            [a13,a23]=[a23,a13]
315

316
        # order 23
317
        if a2 > a3 or (a2 == a3 and abs(a13) > abs(a12)):
318
            [a2,a3]=[a3,a2]
319
            [a13,a12]=[a12,a13]
320

321
        # order 12
322
        if a1 > a2 or (a1 == a2 and abs(a23) > abs(a13)):
323
            [a1,a2]=[a2,a1]
324
            [a13,a23]=[a23,a13]
325

326
        # signs
327
        if a23*a13*a12 > 0:
328
            # a23, a13, a12 positive
329

330
            if (a23<0):
331
                a23=-a23
332
            if (a13<0):
333
                a13=-a13
334
            if (a12<0):
335
                a12=-a12
336

337
        else:
338
            # a23, a13, a12 nonpositive
339

340
            [s1,s2,s3]=[(a23>0),(a13>0),(a12>0)]
341
            if ((s1+s2+s3)%2):
342
                if (a23==0):
343
                    s1=1
344
                else:
345
                    if (a13==0):
346
                        s2=1
347
                    else:
348
                        if (a12==0):
349
                            s3=1
350
            if s1:
351
                  a23=-a23
352
            if s2:
353
                  a13=-a13
354
            if s3:
355
                  a12=-a12
356

357
        loop = not (abs(a23) <= a2 and abs(a13) <= a1 and abs(a12) <= a1 and a1+a2+a23+a13+a12 >= 0)
358

359

360
################
361
################
362

363

364

365
    # adj 3
366
    if a1+a2+a23+a13+a12 == 0 and 2*a1+2*a13+a12 > 0:
367
        # a3 += a1+a2+a23+a13+a12
368
        a23=-2*a2-a23-a12
369
        a13=-2*a1-a13-a12
370

371
    # adj 5.12
372
    if a1 == -a12 and a13 != 0:
373
        #a2 += a1+a12
374
        a23=-a23-a13
375
        a13=-a13
376
        a12=-a12  # = 2*a1+a12
377

378
    # adj 5.13
379
    if a1 == -a13 and a12 != 0:
380
        # a3 += a1+a13
381
        a23=-a23-a12
382
        a13=-a13  # = 2*a1+a13
383
        a12=-a12
384

385
    # adj 5.23
386
    if a2 == -a23 and a12 != 0:
387
        # a3 += a2+a23
388
        a23=-a23  # = 2*a2+a23
389
        a13=-a13-a12
390
        a12=-a12
391

392
    # adj 4.12
393
    if a1 == a12 and a13 > 2*a23:
394
        # a 2 += a1-a12
395
        a23 = -a23 + a13
396
        # a12 = 2*a1 - a12
397

398
    # adj 4.13
399
    if a1 == a13 and a12 > 2*a23:
400
        # a3 += a1-a13
401
        a23 = -a23 + a12
402
        # a13 = 2*a1 - a13
403

404
    # adj 4.23
405
    if a2 == a23 and a12 > 2*a13:
406
        # a3 += a2-a23
407
        # a23 = 2*a2 - a23
408
        a13 = -a13 + a12
409

410
    # order 12
411
    if a1 == a2 and abs(a23) > abs(a13):
412
        [a1,a2]=[a2,a1]
413
        [a13,a23]=[a23,a13]
414

415
    # order 23
416
    if a2 == a3 and abs(a13) > abs(a12):
417
        [a13,a12]=[a12,a13]
418

419
    # order 12
420
    if a1 == a2 and abs(a23) > abs(a13):
421
        [a13,a23]=[a23,a13]
422

423
    return((a1,a2,a3,a23,a13,a12))
424

425

426
def primitivize(long long v0, long long v1, long long v2, p):
427
    """
428
    Given a 3-tuple v not singular mod p, it returns a primitive 3-tuple version of v mod p.
429

430
    EXAMPLES::
431

432
        sage: from sage.quadratic_forms.ternary import primitivize
433
        sage: primitivize(12, 13, 14, 5)
434
        (3, 2, 1)
435
        sage: primitivize(12, 13, 15, 5)
436
        (4, 1, 0)
437

438
    """
439

440
    if v2%p != 0:
441
        v2_inv = inverse_mod(v2, p)
442
        return v2_inv*v0%p, v2_inv*v1%p, 1
443
    elif v1%p != 0:
444
        return inverse_mod(v1, p)*v0%p, 1, 0
445
    else:
446
        return 1, 0, 0
447

448
def evaluate(a, b, c, r, s, t, v):
449
    """
450
    Function to evaluate the ternary quadratic form (a, b, c, r, s, t) in a 3-tuple v.
451

452
    EXAMPLES::
453

454
        sage: from sage.quadratic_forms.ternary import evaluate
455
        sage: Q = TernaryQF([1, 2, 3, -1, 0, 0])
456
        sage: v = (1, -1, 19)
457
        sage: Q(v)
458
        1105
459
        sage: evaluate(1, 2, 3, -1, 0, 0, v)
460
        1105
461

462
    """
463

464
    return a*v[0]**2+b*v[1]**2+c*v[2]**2+r*v[2]*v[1]+s*v[2]*v[0]+t*v[1]*v[0]
465

466
def _find_zeros_mod_p_2(a, b, c, r, s, t):
467
    """
468
    Function to find the zeros mod 2 of a ternary quadratic form.
469

470
    EXAMPLES::
471

472
        sage: Q = TernaryQF([1, 2, 2, -1, 0, 0])
473
        sage: from sage.quadratic_forms.ternary import _find_zeros_mod_p_2
474
        sage: zeros = _find_zeros_mod_p_2(1, 2, 2, -1, 0, 0)
475
        sage: zeros
476
        [(0, 1, 0), (0, 0, 1), (1, 1, 1)]
477
        sage: Q((0, 1, 0))
478
        2
479
        sage: Q((0, 0, 1))
480
        2
481
        sage: Q((1, 1, 1))
482
        4
483

484
    """
485

486
    zeros=[]
487
    v=(1,0,0)
488
    if evaluate(a,b,c,r,s,t,v)%2==0:
489
        zeros.append(v)
490
    for i in range(2):
491
        v=(i,1,0)
492
        if evaluate(a,b,c,r,s,t,v)%2==0:
493
            zeros.append(v)
494
    for i in range(2):
495
        for j in range(2):
496
            v=(i,j,1)
497
            if evaluate(a,b,c,r,s,t,v)%2==0:
498
                zeros.append(v)
499
    return zeros
500

501
def pseudorandom_primitive_zero_mod_p(a, b, c, r, s, t, p):
502
    """
503
    Find a zero of the form (a, b, 1) of the ternary quadratic form given by the coefficients (a, b, c, r, s, t)
504
    mod p, where p is a odd prime that doesn't divide the discriminant.
505

506
    EXAMPLES::
507

508
        sage: Q = TernaryQF([1, 2, 2, -1, 0, 0])
509
        sage: p = 1009
510
        sage: from sage.quadratic_forms.ternary import pseudorandom_primitive_zero_mod_p
511
        sage: v = pseudorandom_primitive_zero_mod_p(1, 2, 2, -1, 0, 0, p)
512
        sage: v[2]
513
        1
514
        sage: Q(v)%p
515
        0
516

517
    """
518

519
    #[a,b,c,r,s,t] = Q.coefficients()
520
    while True:
521

522
        r1 = randint(0,p-1)
523
        r2 = randint(0,p-1)
524
        alpha = (b*r1**2+t*r1+a)%p
525
        if alpha != 0:
526

527
            beta = (2*b*r1*r2+t*r2+r*r1+s)%p
528
            gamma = (b*r2**2+r*r2+c)%p
529
            disc = beta**2-4*alpha*gamma
530
            if mod(disc,p).is_square():
531

532
                z = (-beta+mod(disc,p).sqrt().lift())*(2*alpha).inverse_mod(p)
533
                #return vector((z,r1*z+r2,1))%p
534
                return z%p, (r1*z+r2)%p, 1
535

536
def _find_zeros_mod_p_odd(long long a, long long b, long long c, long long r, long long s, long long t, long long p, v):
537
    """
538
    Find the zeros mod p, where p is an odd prime, of a ternary quadratic form given by its coefficients and a given zero of the form v.
539
    The prime p doesn't divides the discriminant of the form.
540

541
    EXAMPLES::
542

543
        sage: from sage.quadratic_forms.ternary import _find_zeros_mod_p_odd
544
        sage: Q = TernaryQF([1, 2, 2, -1, 0, 0])
545
        sage: p = 1009
546
        sage: v = (817, 974, 1)
547
        sage: Q(v)%1009
548
        0
549
        sage: zeros_1009 = _find_zeros_mod_p_odd(1, 2, 2, -1, 0, 0, 1009, v)
550
        sage: len(zeros_1009)
551
        1010
552
        sage: zeros_1009.sort()
553
        sage: zeros_1009[0]
554
        (0, 32, 1)
555
        sage: Q((0, 32, 1))
556
        2018
557

558
    """
559

560
    cdef long long a_i
561
    cdef long long c_i
562
    cdef long long a_inf
563
    cdef long long c_inf
564
    cdef long long i
565
    cdef long long l
566
    cdef long long x0
567
    cdef long long y0
568

569
    zeros=[v]
570
    x0=v[0]
571
    y0=v[1]
572
    more=False
573
    for i in xrange(p):
574
        a_i=(a*x0**2+b*i**2-2*b*i*y0+b*y0**2-t*i*x0+t*x0*y0-2*a*x0+t*i-t*y0+s*x0-r*i+r*y0+a+c-s)%p
575
        #b_i=(-2*b*i**2+2*b*i*y0+t*i*x0+2*a*x0-2*t*i+t*y0+r*i-2*a+s)%p
576
        if a_i==0:
577
            w=((x0-1)%p,(y0-i)%p,1)
578
            if w==v:
579
                more=True
580
            else:
581
                zeros.append(w)
582
        else:
583
            c_i=(b*i**2+t*i+a)%p
584
            l=c_i*ZZ(a_i).inverse_mod(p)
585
            w = primitivize(l*(x0-1)+1, l*(y0-i)+i, l , p)
586
            if (w[0]==v[0] and w[1]==v[1] and w[2]==v[2]):
587
                more=True
588
            else:
589
                zeros.append(w)
590
    if more:
591
        a_inf=(a*x0**2+b*y0**2+t*x0*y0-2*b*y0-t*x0+s*x0+r*y0+b+c-r)%p
592
        #b_inf=(2*b*y0+t*x0-2*b+r)%p
593
        if a_inf==0:
594
            w=(x0%p,(y0-1)%p,1)
595
            zeros.append(w)
596
        else:
597
            c_inf=b%p
598
            l=c_inf*ZZ(a_inf).inverse_mod(p)
599
            w = primitivize(l*x0, l*(y0-1)+1, l, p)
600
            zeros.append(w)
601

602
    return zeros
603

604

605

606
def _find_zeros_mod_p(a, b, c, r, s, t, p):
607
    """
608
    Finds the zeros mod p of the ternary quadratic form given by the coefficients (a, b, c, r, s, t), where p is
609
    a prime that doesn't divides the discriminant of the form.
610

611
    EXAMPLES::
612

613
        sage: from sage.quadratic_forms.ternary import _find_zeros_mod_p
614
        sage: Q = TernaryQF([1, 2, 2, -1, 0, 0])
615
        sage: p = 1009
616
        sage: zeros_1009 = _find_zeros_mod_p(1, 2, 2, -1, 0, 0, p)
617
        sage: len(zeros_1009)
618
        1010
619
        sage: zeros_1009.sort()
620
        sage: zeros_1009[0]
621
        (0, 32, 1)
622
        sage: Q((0, 32, 1))
623
        2018
624
        sage: zeros_2 = _find_zeros_mod_p(1, 2, 2, -1, 0, 0, 2)
625
        sage: zeros_2
626
        [(0, 1, 0), (0, 0, 1), (1, 1, 1)]
627

628
    """
629

630
    if p==2:
631
        return _find_zeros_mod_p_2(a, b, c, r, s, t)
632
    else:
633
        v = pseudorandom_primitive_zero_mod_p(a, b, c, r, s, t, p)
634
        return _find_zeros_mod_p_odd(a, b, c, r, s, t, p, v)
635

636

637
def _find_all_ternary_qf_by_level_disc(long long N, long long d):
638
    """
639
    Find the coefficients of all the reduced ternary quadratic forms given its discriminant d and level N.
640
    If N|4d and d|N^2, then it may be some forms with that discriminant and level.
641

642
    EXAMPLES::
643

644
        sage: from sage.quadratic_forms.ternary import _find_all_ternary_qf_by_level_disc
645
        sage: _find_all_ternary_qf_by_level_disc(44, 11)
646
        [(1, 1, 3, 0, -1, 0), (1, 1, 4, 1, 1, 1)]
647
        sage: _find_all_ternary_qf_by_level_disc(44, 11^2 * 16)
648
        [(3, 15, 15, -14, -2, -2), (4, 11, 12, 0, -4, 0)]
649
        sage: Q = TernaryQF([1, 1, 3, 0, -1, 0])
650
        sage: Q.is_eisenstein_reduced()
651
        True
652
        sage: Q.reciprocal_reduced()
653
        Ternary quadratic form with integer coefficients:
654
        [4 11 12]
655
        [0 -4 0]
656
        sage: _find_all_ternary_qf_by_level_disc(44, 22)
657
        []
658
        sage: _find_all_ternary_qf_by_level_disc(44, 33)
659
        Traceback (most recent call last):
660
        ...
661
        ValueError: There are no ternary forms of this level and discriminant
662

663

664
    """
665

666
    cdef long long a
667
    cdef long long b
668
    cdef long long c
669
    cdef long long r
670
    cdef long long s
671
    cdef long long t
672
    cdef long long m
673
    cdef long long mu
674
    cdef long long g
675
    cdef long long u
676
    cdef long long v
677
    cdef long long g1
678
    cdef long long m_q
679
    cdef double a_max
680
    cdef double r_max
681
    cdef double b_max
682
    cdef long long stu2
683
    cdef long long mg
684

685

686

687
    l=[]
688

689
    if (4*d)%N!=0:
690
         raise ValueError, "There are no ternary forms of this level and discriminant"
691
    else:
692
        m=4*d//N
693

694
    if (N**2)%d!=0:
695
        raise ValueError, "There are no ternary forms of this level and discriminant"
696
    else:
697
        mu=N*N//d
698

699
    m_2=m%2
700
    mu_2=mu%2
701

702
    a=1
703
    a_max=(d/2.)**(1/3.)
704
    while a<=a_max:
705

706
        [g,u,v]=xgcd(4*a,m)
707
        g1=(ZZ(g).squarefree_part()*g).sqrtrem()[0]
708
        t=0
709
        while t<=a:
710

711
            alpha=max(a,m/4/a)
712
            b=alpha+((((u*t*t//g)-alpha))%(m//g))
713
            b_max=(d/2.0/a)**(1/2.)
714
            while b<=b_max:
715

716
                s=0
717
                beta=g//gcd(g,2*t)
718
                while s<=a:
719

720
                    r = -b+((((2*s*t*u//g)+b))%(m//g))
721
                    if s*t==0:
722
                        r_max=0
723
                    else:
724
                        r_max=b
725
                    while r<=r_max:
726

727
                        if (d-r*s*t+a*r*r+b*s*s)%(4*a*b-t**2)==0:
728

729
                            c=(d-r*s*t+a*r*r+b*s*s)//(4*a*b-t**2)
730
                            if r <= 0:
731
                                is_reduced = True
732
                                if r < -b:
733
                                    is_reduced = False
734
                                elif not (b <= c and 0 <= a+b+r-s-t):
735
                                    is_reduced = False
736
                                elif a == b and abs(r) > abs(s):
737
                                    is_reduced = False
738
                                elif b == c and abs(s) > abs(t):
739
                                    is_reduced = False
740
                                elif a+b+r-s-t == 0 and 2*a-2*s-t > 0:
741
                                    is_reduced = False
742
                                elif a == t and s != 0:
743
                                    is_reduced = False
744
                                elif a == s and t != 0:
745
                                    is_reduced = False
746
                                elif b == -r and t != 0:
747
                                    is_reduced = False
748
                                if is_reduced:
749
                                    m_q=gcd((4*b*c-r**2, 4*a*c-s**2, 4*a*b-t**2, 2*s*t-4*a*r, -2*r*t+4*b*s, -2*r*s+4*c*t))
750
                                    if m_q==m:
751
                                        l.append((ZZ(a), ZZ(b), ZZ(c), ZZ(r), ZZ(-s), ZZ(-t)))
752
                            else:
753
                                is_reduced=True
754
                                if not (b <= c and 0 <= a+b+r+s+t):
755
                                    is_reduced = False
756
                                elif a == b and abs(r) > abs(s):
757
                                    is_reduced = False
758
                                elif b == c and abs(s) > abs(t):
759
                                    is_reduced = False
760
                                elif a+b+r+s+t == 0 and 2*a+2*s+t > 0:
761
                                    is_reduced = False
762
                                elif a == t and s > 2*r:
763
                                    is_reduced = False
764
                                elif a == s and t > 2*r:
765
                                    is_reduced = False
766
                                elif b == r and t > 2*s:
767
                                    is_reduced = False
768
                                if is_reduced:
769
                                    m_q=gcd((4*b*c-r**2, 4*a*c-s**2, 4*a*b-t**2, 2*s*t-4*a*r, 2*r*t-4*b*s, 2*r*s-4*c*t))
770
                                    if m_q==m:
771
                                        l.append((ZZ(a), ZZ(b), ZZ(c), ZZ(r), ZZ(s), ZZ(t)))
772
                        r+=(m//g)
773
                    s+=beta
774
                b+=m//g
775
            t+=g1
776

777
        a+=1
778
    return l
779

780

781

782
def _find_a_ternary_qf_by_level_disc(long long N, long long d):
783
    """
784
    Find the coefficients of a reduced ternary quadratic form given its discriminant d and level N.
785
    If N|4d and d|N^2, then it may be a form with that discriminant and level.
786

787
    EXAMPLES::
788

789
        sage: from sage.quadratic_forms.ternary import _find_a_ternary_qf_by_level_disc
790
        sage: _find_a_ternary_qf_by_level_disc(44, 11)
791
        (1, 1, 3, 0, -1, 0)
792
        sage: _find_a_ternary_qf_by_level_disc(44, 11^2 * 16)
793
        (3, 15, 15, -14, -2, -2)
794
        sage: Q = TernaryQF([1, 1, 3, 0, -1, 0])
795
        sage: Q.is_eisenstein_reduced()
796
        True
797
        sage: Q.level()
798
        44
799
        sage: Q.disc()
800
        11
801
        sage: _find_a_ternary_qf_by_level_disc(44, 22)
802
        sage: _find_a_ternary_qf_by_level_disc(44, 33)
803
        Traceback (most recent call last):
804
        ...
805
        ValueError: There are no ternary forms of this level and discriminant
806

807
    """
808

809
    cdef long long a
810
    cdef long long b
811
    cdef long long c
812
    cdef long long r
813
    cdef long long s
814
    cdef long long t
815
    cdef long long m
816
    cdef long long mu
817
    cdef long long g
818
    cdef long long u
819
    cdef long long v
820
    cdef long long g1
821
    cdef long long m_q
822
    cdef double a_max
823
    cdef double r_max
824
    cdef double b_max
825
    cdef long long stu2
826
    cdef long long mg
827

828

829

830
    if (4*d)%N!=0:
831
         raise ValueError, "There are no ternary forms of this level and discriminant"
832
    else:
833
        m=4*d//N
834

835
    if (N**2)%d!=0:
836
        raise ValueError, "There are no ternary forms of this level and discriminant"
837
    else:
838
        mu=N*N//d
839

840
    m_2=m%2
841
    mu_2=mu%2
842

843
    a=1
844
    a_max=(d/2.)**(1/3.)
845
    while a<=a_max:
846

847
        [g,u,v]=xgcd(4*a,m)
848
        g1=(ZZ(g).squarefree_part()*g).sqrtrem()[0]
849
        t=0
850
        while t<=a:
851

852
            alpha=max(a,m/4/a)
853
            b=alpha+((((u*t*t//g)-alpha))%(m//g))
854
            b_max=(d/2.0/a)**(1/2.)
855
            while b<=b_max:
856

857
                s=0
858
                beta=g//gcd(g,2*t)
859
                while s<=a:
860

861
                    r = -b+((((2*s*t*u//g)+b))%(m//g))
862
                    if s*t==0:
863
                        r_max=0
864
                    else:
865
                        r_max=b
866
                    while r<=r_max:
867

868
                        if (d-r*s*t+a*r*r+b*s*s)%(4*a*b-t**2)==0:
869

870
                            c=(d-r*s*t+a*r*r+b*s*s)//(4*a*b-t**2)
871
                            if r <= 0:
872
                                is_reduced = True
873
                                if r < -b:
874
                                    is_reduced = False
875
                                elif not (b <= c and 0 <= a+b+r-s-t):
876
                                    is_reduced = False
877
                                elif a == b and abs(r) > abs(s):
878
                                    is_reduced = False
879
                                elif b == c and abs(s) > abs(t):
880
                                    is_reduced = False
881
                                elif a+b+r-s-t == 0 and 2*a-2*s-t > 0:
882
                                    is_reduced = False
883
                                elif a == t and s != 0:
884
                                    is_reduced = False
885
                                elif a == s and t != 0:
886
                                    is_reduced = False
887
                                elif b == -r and t != 0:
888
                                    is_reduced = False
889
                                if is_reduced:
890
                                    m_q=gcd((4*b*c-r**2, 4*a*c-s**2, 4*a*b-t**2, 2*s*t-4*a*r, -2*r*t+4*b*s, -2*r*s+4*c*t))
891
                                    if m_q==m:
892
                                        return ZZ(a), ZZ(b), ZZ(c), ZZ(r), ZZ(-s), ZZ(-t)
893
                            else:
894
                                is_reduced = True
895
                                if not (b <= c and 0 <= a+b+r+s+t):
896
                                    is_reduced = False
897
                                elif a == b and abs(r) > abs(s):
898
                                    is_reduced = False
899
                                elif b == c and abs(s) > abs(t):
900
                                    is_reduced = False
901
                                elif a+b+r+s+t == 0 and 2*a+2*s+t > 0:
902
                                    is_reduced = False
903
                                elif a==t and s > 2*r:
904
                                    is_reduced = False
905
                                elif a == s and t>2*r:
906
                                    is_reduced = False
907
                                elif b == r and t > 2*s:
908
                                    is_reduced = False
909
                                if is_reduced:
910
                                    m_q=gcd((4*b*c-r**2, 4*a*c-s**2, 4*a*b-t**2, 2*s*t-4*a*r, 2*r*t-4*b*s, 2*r*s-4*c*t))
911
                                    if m_q==m:
912
                                        return ZZ(a), ZZ(b), ZZ(c), ZZ(r), ZZ(s), ZZ(t)
913
                        r+=(m//g)
914
                    s+=beta
915
                b+=m//g
916
            t+=g1
917

918
        a+=1
919

920

921

922
def extend(v):
923
    """
924
    Return the coefficients of a matrix M such that M has determinant gcd(v) and the first column is v.
925

926
    EXAMPLES::
927

928
        sage: from sage.quadratic_forms.ternary import extend
929
        sage: v = (6, 4, 12)
930
        sage: m = extend(v)
931
        sage: M = matrix(3, m)
932
        sage: M
933
        [ 6  1  0]
934
        [ 4  1  0]
935
        [12  0  1]
936
        sage: M.det()
937
        2
938
        sage: v = (-12, 20, 30)
939
        sage: m = extend(v)
940
        sage: M = matrix(3, m)
941
        sage: M
942
        [-12   1   0]
943
        [ 20  -2   1]
944
        [ 30   0  -7]
945
        sage: M.det()
946
        2
947

948
    """
949

950
    b1 = xgcd(v[0], v[1])
951
    b2 = xgcd(b1[1], b1[2])
952
    b3 = xgcd(b1[0], v[2])
953

954
    return v[0], -b1[2], -b2[1]*b3[2], v[1], b1[1], -b2[2]*b3[2], v[2], 0, b3[1]
955

956

957
def _find_p_neighbor_from_vec(a, b, c, r, s, t, p, v, mat = False):
958
    """
959
    Finds the coefficients of the reduced equivalent of the p-neighbor
960
    of the ternary quadratic given by Q = (a, b, c, r, s, t)  form associated
961
    to a given vector v satisfying:
962

963
    1. Q(v) = 0  mod p
964

965
    2. v is a non-singular point of the conic Q(v) = 0 mod p.
966

967
    Reference:  Gonzalo Tornaria's Thesis, Thrm 3.5, p34.
968

969
    EXAMPLES:
970

971
        sage: from sage.quadratic_forms.ternary import _find_p_neighbor_from_vec
972
        sage: Q = TernaryQF([1, 3, 3, -2, 0, -1])
973
        sage: Q
974
        Ternary quadratic form with integer coefficients:
975
        [1 3 3]
976
        [-2 0 -1]
977
        sage: Q.disc()
978
        29
979
        sage: v = (9, 7, 1)
980
        sage: v in Q.find_zeros_mod_p(11)
981
        True
982
        sage: q11, M = _find_p_neighbor_from_vec(1, 3, 3, -2, 0, -1, 11, v, mat = True)
983
        sage: Q11 = TernaryQF(q11)
984
        sage: Q11
985
        Ternary quadratic form with integer coefficients:
986
        [1 2 4]
987
        [-1 -1 0]
988
        sage: M = matrix(3, M)
989
        sage: M
990
        [    -1  -5/11   7/11]
991
        [     0 -10/11   3/11]
992
        [     0  -3/11  13/11]
993
        sage: Q(M) == Q11
994
        True
995

996
    """
997

998
    v0, w0, u0, v1, w1, u1, v2, w2, u2 = extend(v)
999

1000
    m00 = 2*a*v0**2 + 2*b*v1**2 + 2*c*v2**2 + 2*r*v1*v2 + 2*s*v0*v2 + 2*t*v0*v1
1001
    m11 = 2*a*w0**2 + 2*b*w1**2 + 2*t*w0*w1
1002
    m22 = 2*a*u0**2 + 2*b*u1**2 + 2*c*u2**2 + 2*r*u1*u2 + 2*s*u0*u2 + 2*t*u0*u1
1003
    m01 = 2*a*v0*w0 + 2*b*v1*w1 + r*v2*w1 + s*v2*w0 + t*v0*w1 + t*v1*w0
1004
    m02 = 2*a*u0*v0 + 2*b*u1*v1 + 2*c*u2*v2 + r*u1*v2 + r*u2*v1 + s*u0*v2 + s*u2*v0 + t*u0*v1 + t*u1*v0
1005
    m12 = 2*a*u0*w0 + 2*b*u1*w1 + r*u2*w1 + s*u2*w0 + t*u0*w1 + t*u1*w0
1006

1007

1008
    if m02%p!=0:
1009

1010
        m0 = (-m00/m02/2) % p**2
1011
        m1 = (-m01/m02) % p
1012

1013
        b00 = m0**2*m22/p**2 + 2*m0*m02/p**2 + m00/p**2
1014
        b11 = m1**2*m22 + 2*m1*m12 + m11
1015
        b22 = m22*p**2
1016
        b01 = m0*m1*m22/p + m0*m12/p + m02*m1/p + m01/p
1017
        b02 = m0*m22 + m02
1018
        b12 = m1*m22*p + m12*p
1019

1020
        if mat:
1021
            q, Mr = _reduced_ternary_form_eisenstein_with_matrix(ZZ(b00/2), ZZ(b11/2), ZZ(b22/2), ZZ(b12), ZZ(b02), ZZ(b01))
1022
            r00, r01, r02, r10, r11, r12, r20, r21, r22 = Mr.list()
1023
            t00 = p*r20*u0 + (m0*u0/p + v0/p)*r00 + (m1*u0 + w0)*r10
1024
            t01 = p*r21*u0 + (m0*u0/p + v0/p)*r01 + (m1*u0 + w0)*r11
1025
            t02 = p*r22*u0 + (m0*u0/p + v0/p)*r02 + (m1*u0 + w0)*r12
1026
            t10 = p*r20*u1 + (m0*u1/p + v1/p)*r00 + (m1*u1 + w1)*r10
1027
            t11 = p*r21*u1 + (m0*u1/p + v1/p)*r01 + (m1*u1 + w1)*r11
1028
            t12 = p*r22*u1 + (m0*u1/p + v1/p)*r02 + (m1*u1 + w1)*r12
1029
            t20 = p*r20*u2 + (m0*u2/p + v2/p)*r00 + (m1*u2 + w2)*r10
1030
            t21 = p*r21*u2 + (m0*u2/p + v2/p)*r01 + (m1*u2 + w2)*r11
1031
            t22 = p*r22*u2 + (m0*u2/p + v2/p)*r02 + (m1*u2 + w2)*r12
1032
            return q, (t00, t01, t02, t10, t11, t12, t20, t21, t22)
1033
        else:
1034
            return _reduced_ternary_form_eisenstein_without_matrix(ZZ(b00/2), ZZ(b11/2), ZZ(b22/2), ZZ(b12), ZZ(b02), ZZ(b01))
1035

1036
    if m01%p!=0:
1037

1038
        m0 = (-m00/m01/2) % p**2
1039
        m1 = (-m02/m01) % p
1040

1041
        b00 = m0**2*m11/p**2 + 2*m0*m01/p**2 + m00/p**2
1042
        b11 = m1**2*m11 + 2*m1*m12 + m22
1043
        b22 = m11*p**2
1044
        b01 = m0*m1*m11/p + m0*m12/p + m01*m1/p + m02/p
1045
        b02 = m0*m11 + m01
1046
        b12 = m1*m11*p + m12*p
1047

1048
        if mat:
1049
            q, Mr = _reduced_ternary_form_eisenstein_with_matrix(ZZ(b00/2), ZZ(b11/2), ZZ(b22/2), ZZ(b12), ZZ(b02), ZZ(b01))
1050
            r00, r01, r02, r10, r11, r12, r20, r21, r22 = Mr.list()
1051
            t00 = p*r20*w0 + (m0*w0/p + v0/p)*r00 + (m1*w0 + u0)*r10
1052
            t01 = p*r21*w0 + (m0*w0/p + v0/p)*r01 + (m1*w0 + u0)*r11
1053
            t02 = p*r22*w0 + (m0*w0/p + v0/p)*r02 + (m1*w0 + u0)*r12
1054
            t10 = p*r20*w1 + (m0*w1/p + v1/p)*r00 + (m1*w1 + u1)*r10
1055
            t11 = p*r21*w1 + (m0*w1/p + v1/p)*r01 + (m1*w1 + u1)*r11
1056
            t12 = p*r22*w1 + (m0*w1/p + v1/p)*r02 + (m1*w1 + u1)*r12
1057
            t20 = p*r20*w2 + (m0*w2/p + v2/p)*r00 + (m1*w2 + u2)*r10
1058
            t21 = p*r21*w2 + (m0*w2/p + v2/p)*r01 + (m1*w2 + u2)*r11
1059
            t22 = p*r22*w2 + (m0*w2/p + v2/p)*r02 + (m1*w2 + u2)*r12
1060
            return q, (t00, t01, t02, t10, t11, t12, t20, t21, t22)
1061
        else:
1062
            return _reduced_ternary_form_eisenstein_without_matrix(ZZ(b00/2), ZZ(b11/2), ZZ(b22/2), ZZ(b12), ZZ(b02), ZZ(b01))
1063

1064

1065
def _basic_lemma_vec(a, b, c, r, s, t, n):
1066
    """
1067
    Find a vector v such that the ternary quadratic form given by (a, b, c, r, s, t) evaluated at v is
1068
    coprime with n a prime or 1.
1069

1070
    EXAMPLES::
1071

1072
        sage: from sage.quadratic_forms.ternary import _basic_lemma_vec
1073
        sage: Q = TernaryQF([5, 2, 3, -1, 0, 0])
1074
        sage: v = _basic_lemma_vec(5, 2, 3, -1, 0, 0, 5)
1075
        sage: v
1076
        (0, 1, 0)
1077
        sage: Q(v)
1078
        2
1079

1080
    """
1081

1082
    if n == 1:
1083
        return 0, 0, 0
1084

1085
    if a%n != 0:
1086
        return 1, 0, 0
1087
    elif b%n != 0:
1088
        return 0, 1, 0
1089
    elif c%n != 0:
1090
        return 0, 0, 1
1091

1092
    if r%n != 0:
1093
        return 0, 1, 1
1094
    elif s%n != 0:
1095
        return 1, 0, 1
1096
    elif t%n != 0:
1097
        return 1, 1, 0
1098

1099
    raise ValueError, "not primitive form"
1100

1101
def _basic_lemma(a, b, c, r, s, t, n):
1102
    """
1103
    Finds a number represented by the ternary quadratic form given by the coefficients (a, b, c, r, s, t)
1104
    and coprime to the prime n.
1105

1106
    EXAMPLES::
1107

1108
        sage: from sage.quadratic_forms.ternary import _basic_lemma
1109
        sage: Q = TernaryQF([5, 2, 3, -1, 0, 0])
1110
        sage: _basic_lemma(5, 2, 3, -1, 0, 0, 5)
1111
        2
1112

1113
    """
1114

1115
    if n == 1:
1116
        return 0
1117

1118
    if a%n != 0:
1119
        return a
1120
    elif b%n != 0:
1121
        return b
1122
    elif c%n != 0:
1123
        return c
1124

1125
    if r%n != 0:
1126
        return b + c + r
1127
    elif s%n != 0:
1128
        return a + c + s
1129
    elif t%n != 0:
1130
        return a + b + t
1131

1132
    raise ValueError, "not primitive form"
1133

1134

1135
Product

Resources

Company
