| 504 | cdef bint Zp_soluble_siksek_large_p(mpz_t a, mpz_t b, mpz_t c, mpz_t d, mpz_t e, mpz_t pp, |
| 505 | fmpz_poly_t f1, fmpz_poly_t linear): |
| 506 | """ |
| 507 | Uses the approach of Algorithm 5.3.1 of Siksek's thesis to test for |
| 508 | solubility of y^2 == ax^4 + bx^3 + cx^2 + dx + e over Zp. |
| 509 | """ |
| 510 | cdef unsigned long v_min, v |
| 511 | cdef mpz_t roots[4] |
| 512 | cdef int i, j, has_roots, has_single_roots |
| 513 | cdef bint result |
| 514 | |
| 515 | cdef mpz_t aa, bb, cc, dd, ee |
| 516 | cdef mpz_t aaa, bbb, ccc, ddd, eee |
| 517 | cdef mpz_t qq, rr, ss, tt |
| 518 | cdef Integer A,B,C,D,E,P |
| 519 | |
| 520 | # Step 0: divide out all common p from the quartic |
| 521 | v_min = valuation(a, pp) |
| 522 | if mpz_cmp_ui(b, ui0) != 0: |
| 523 | v = valuation(b, pp) |
| 524 | if v < v_min: v_min = v |
| 525 | if mpz_cmp_ui(c, ui0) != 0: |
| 526 | v = valuation(c, pp) |
| 527 | if v < v_min: v_min = v |
| 528 | if mpz_cmp_ui(d, ui0) != 0: |
| 529 | v = valuation(d, pp) |
| 530 | if v < v_min: v_min = v |
| 531 | if mpz_cmp_ui(e, ui0) != 0: |
| 532 | v = valuation(e, pp) |
| 533 | if v < v_min: v_min = v |
| 534 | for 0 <= v < v_min: |
| 535 | mpz_divexact(a, a, pp) |
| 536 | mpz_divexact(b, b, pp) |
| 537 | mpz_divexact(c, c, pp) |
| 538 | mpz_divexact(d, d, pp) |
| 539 | mpz_divexact(e, e, pp) |
| 540 | |
| 541 | if not v_min%2: |
| 542 | # Step I in Alg. 5.3.1 of Siksek's thesis |
| 543 | A = Integer(0); B = Integer(0); C = Integer(0); D = Integer(0); E = Integer(0); P = Integer(0) |
| 544 | mpz_set(A.value, a); mpz_set(B.value, b); mpz_set(C.value, c); mpz_set(D.value, d); mpz_set(E.value, e); mpz_set(P.value, pp) |
| 545 | f = ntl.ZZ_pX([E,D,C,B,A], P) |
| 546 | f /= ntl.ZZ_pX([A], P) # now f is monic, and we are done with A,B,C,D,E |
| 547 | mpz_set(qq, A.value) # qq is the leading coefficient of the polynomial |
| 548 | f_factzn = f.factor() |
| 549 | result = 0 |
| 550 | for factor, exponent in f_factzn: |
| 551 | if exponent&1: |
| 552 | result = 1 |
| 553 | break |
| 554 | if result == 0 and mpz_legendre(qq, pp) == 1: |
| 555 | result = 1 |
| 556 | if result: |
| 557 | return 1 |
| 558 | |
| 559 | f = ntl.ZZ_pX([1], P) |
| 560 | for factor, exponent in f_factzn: |
| 561 | for j from 0 <= j < (exponent/2): |
| 562 | f *= factor |
| 563 | |
| 564 | f /= f.leading_coefficient() |
| 565 | f_factzn = f.factor() |
| 566 | |
| 567 | has_roots = 0 |
| 568 | j = 0 |
| 569 | for factor, exponent in f_factzn: |
| 570 | if factor.degree() == 1: |
| 571 | has_roots = 1 |
| 572 | A = P - Integer(factor[0]) |
| 573 | mpz_set(roots[j], A.value) |
| 574 | j += 1 |
| 575 | if not has_roots: |
| 576 | return 0 |
| 577 | |
| 578 | i = f.degree() |
| 579 | mpz_init(aaa) |
| 580 | mpz_init(bbb) |
| 581 | mpz_init(ccc) |
| 582 | mpz_init(ddd) |
| 583 | mpz_init(eee) |
| 584 | |
| 585 | if i == 0: # g == 1 |
| 586 | mpz_set(aaa, a) |
| 587 | mpz_set(bbb, b) |
| 588 | mpz_set(ccc, c) |
| 589 | mpz_set(ddd, d) |
| 590 | mpz_sub(eee, e, qq) |
| 591 | elif i == 1: # g == x + rr |
| 592 | mpz_set(aaa, a) |
| 593 | mpz_set(bbb, b) |
| 594 | mpz_sub(ccc, c, qq) |
| 595 | A = Integer(f[0]) |
| 596 | mpz_set(rr, A.value) |
| 597 | mpz_mul(ss, rr, qq) |
| 598 | mpz_set(ddd,d) |
| 599 | mpz_sub(ddd, ddd, ss) |
| 600 | mpz_sub(ddd, ddd, ss) |
| 601 | mpz_set(eee,e) |
| 602 | mpz_mul(ss, ss, rr) |
| 603 | mpz_sub(eee, eee, ss) |
| 604 | mpz_divexact(ss, ss, rr) |
| 605 | elif i == 2: # g == x^2 + rr*x + ss |
| 606 | mpz_sub(aaa, a, qq) |
| 607 | A = Integer(f[1]) |
| 608 | mpz_set(rr, A.value) |
| 609 | mpz_init(tt) |
| 610 | mpz_mul(tt, rr, qq) |
| 611 | mpz_set(bbb,b) |
| 612 | mpz_submul_ui(bbb, tt, ui2) |
| 613 | mpz_set(ccc,c) |
| 614 | mpz_submul(ccc, tt, rr) |
| 615 | A = Integer(f[0]) |
| 616 | mpz_set(ss, A.value) |
| 617 | mpz_mul(tt, ss, qq) |
| 618 | mpz_set(eee,e) |
| 619 | mpz_submul(eee, tt, ss) |
| 620 | mpz_mul_ui(tt, tt, ui2) |
| 621 | mpz_sub(ccc, ccc, tt) |
| 622 | mpz_set(ddd,d) |
| 623 | mpz_submul(ddd, tt, rr) |
| 624 | mpz_clear(tt) |
| 625 | mpz_divexact(aaa, aaa, pp) |
| 626 | mpz_divexact(bbb, bbb, pp) |
| 627 | mpz_divexact(ccc, ccc, pp) |
| 628 | mpz_divexact(ddd, ddd, pp) |
| 629 | mpz_divexact(eee, eee, pp) |
| 630 | # now aaa,bbb,ccc,ddd,eee represents h(x) |
| 631 | |
| 632 | result = 0 |
| 633 | mpz_init(tt) |
| 634 | for i from 0 <= i < j: |
| 635 | mpz_mul(tt, aaa, roots[i]) |
| 636 | mpz_add(tt, tt, bbb) |
| 637 | mpz_mul(tt, tt, roots[i]) |
| 638 | mpz_add(tt, tt, ccc) |
| 639 | mpz_mul(tt, tt, roots[i]) |
| 640 | mpz_add(tt, tt, ddd) |
| 641 | mpz_mul(tt, tt, roots[i]) |
| 642 | mpz_add(tt, tt, eee) |
| 643 | # tt == h(r) mod p |
| 644 | mpz_mod(tt, tt, pp) |
| 645 | if mpz_sgn(tt) == 0: |
| 646 | fmpz_poly_zero(f1) |
| 647 | fmpz_poly_zero(linear) |
| 648 | fmpz_poly_set_coeff_mpz(f1, 0, e) |
| 649 | fmpz_poly_set_coeff_mpz(f1, 1, d) |
| 650 | fmpz_poly_set_coeff_mpz(f1, 2, c) |
| 651 | fmpz_poly_set_coeff_mpz(f1, 3, b) |
| 652 | fmpz_poly_set_coeff_mpz(f1, 4, a) |
| 653 | fmpz_poly_set_coeff_mpz(linear, 0, roots[i]) |
| 654 | fmpz_poly_set_coeff_mpz(linear, 1, pp) |
| 655 | fmpz_poly_compose(f1, f1, linear) |
| 656 | fmpz_poly_scalar_div_mpz(f1, f1, pp) |
| 657 | fmpz_poly_scalar_div_mpz(f1, f1, pp) |
| 658 | |
| 659 | mpz_init(aa) |
| 660 | mpz_init(bb) |
| 661 | mpz_init(cc) |
| 662 | mpz_init(dd) |
| 663 | mpz_init(ee) |
| 664 | fmpz_poly_get_coeff_mpz(aa, f1, 4) |
| 665 | fmpz_poly_get_coeff_mpz(bb, f1, 3) |
| 666 | fmpz_poly_get_coeff_mpz(cc, f1, 2) |
| 667 | fmpz_poly_get_coeff_mpz(dd, f1, 1) |
| 668 | fmpz_poly_get_coeff_mpz(ee, f1, 0) |
| 669 | result = Zp_soluble_siksek_large_p(aa, bb, cc, dd, ee, pp, f1, linear) |
| 670 | mpz_clear(aa) |
| 671 | mpz_clear(bb) |
| 672 | mpz_clear(cc) |
| 673 | mpz_clear(dd) |
| 674 | mpz_clear(ee) |
| 675 | if result == 1: |
| 676 | break |
| 677 | mpz_clear(aaa) |
| 678 | mpz_clear(bbb) |
| 679 | mpz_clear(ccc) |
| 680 | mpz_clear(ddd) |
| 681 | mpz_clear(eee) |
| 682 | mpz_clear(tt) |
| 683 | return result |
| 684 | else: |
| 685 | # Step II in Alg. 5.3.1 of Siksek's thesis |
| 686 | A = Integer(0); B = Integer(0); C = Integer(0); D = Integer(0); E = Integer(0); P = Integer(0) |
| 687 | mpz_set(A.value, a); mpz_set(B.value, b); mpz_set(C.value, c); mpz_set(D.value, d); mpz_set(E.value, e); mpz_set(P.value, pp) |
| 688 | f = ntl.ZZ_pX([E,D,C,B,A], P) |
| 689 | f /= ntl.ZZ_pX([A], P) # now f is monic |
| 690 | f_factzn = f.factor() |
| 691 | |
| 692 | has_roots = 0 |
| 693 | has_single_roots = 0 |
| 694 | j = 0 |
| 695 | for factor, exponent in f_factzn: |
| 696 | if factor.degree() == 1: |
| 697 | has_roots = 1 |
| 698 | if exponent == 1: |
| 699 | has_single_roots = 1 |
| 700 | break |
| 701 | A = P - Integer(factor[0]) |
| 702 | mpz_set(roots[j], A.value) |
| 703 | j += 1 |
| 704 | |
| 705 | if not has_roots: return 0 |
| 706 | if has_single_roots: return 1 |
| 707 | |
| 708 | result = 0 |
| 709 | if j > 0: |
| 710 | mpz_init(aa) |
| 711 | mpz_init(bb) |
| 712 | mpz_init(cc) |
| 713 | mpz_init(dd) |
| 714 | mpz_init(ee) |
| 715 | for i from 0 <= i < j: |
| 716 | fmpz_poly_zero(f1) |
| 717 | fmpz_poly_zero(linear) |
| 718 | fmpz_poly_set_coeff_mpz(f1, 0, e) |
| 719 | fmpz_poly_set_coeff_mpz(f1, 1, d) |
| 720 | fmpz_poly_set_coeff_mpz(f1, 2, c) |
| 721 | fmpz_poly_set_coeff_mpz(f1, 3, b) |
| 722 | fmpz_poly_set_coeff_mpz(f1, 4, a) |
| 723 | fmpz_poly_set_coeff_mpz(linear, 0, roots[i]) |
| 724 | fmpz_poly_set_coeff_mpz(linear, 1, pp) |
| 725 | fmpz_poly_compose(f1, f1, linear) |
| 726 | fmpz_poly_scalar_div_mpz(f1, f1, pp) |
| 727 | fmpz_poly_get_coeff_mpz(aa, f1, 4) |
| 728 | fmpz_poly_get_coeff_mpz(bb, f1, 3) |
| 729 | fmpz_poly_get_coeff_mpz(cc, f1, 2) |
| 730 | fmpz_poly_get_coeff_mpz(dd, f1, 1) |
| 731 | fmpz_poly_get_coeff_mpz(ee, f1, 0) |
| 732 | result = Zp_soluble_siksek_large_p(aa, bb, cc, dd, ee, pp, f1, linear) |
| 733 | if result == 1: |
| 734 | break |
| 735 | if j > 0: |
| 736 | mpz_clear(aa) |
| 737 | mpz_clear(bb) |
| 738 | mpz_clear(cc) |
| 739 | mpz_clear(dd) |
| 740 | mpz_clear(ee) |
| 741 | return result |
| 742 | |