Kunerth's Algorithm Example

#!/usr/bin/python3
# Kunerth's algorithm example
# needs python 3.8+ for pow(#,-1,N)
# otherwise make your own modinv

import math
from sympy import legendre_symbol
from sympy.ntheory import sqrt_mod

primes = [1,2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61, 67, 71, 73,
79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163,
 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013, 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069, 1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151, 1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223]

def trysqrt(someN, p):
  if legendre_symbol(someN,p) != 1:
    return(False)
  return sqrt_mod(someN,p)

def tryfactor(A,B,C,N):
  disc = B**2 - 4*A*C
  sqrt_disc = math.sqrt(disc)
  if not sqrt_disc.is_integer():
    return False
  D = int(sqrt_disc)
  f1 = math.gcd(2*A,B+D)
  A1 = 2*A//f1
  B1 = (B+D)//f1
  f2 = math.gcd(2*A,B-D)
  A2 = 2*A//f2
  B2 = (B-D)//f2
  print(" factors: (",A1,"x +",B1,")(",A2,"x +",B2,")")
  return [A1,B2,A2,B2]


def tryKunerth(p,r,a,N):
  aN = a*N
  C = (r**2-aN)//p
  if C < 0:
    return
  sqrt_C = math.sqrt(C)
  if not sqrt_C.is_integer():
    return
  w = int(sqrt_C)
  v = r
  print("a =",a," r =",r," aN =",aN, "sqrt(C) =",w," poly =",p,"z^2 +",2*r,"z +",(r**2-aN)//p)
  for beta in range(-5,5+1):
    alpha = w*(v + w*beta)
    factors = tryfactor(alpha*alpha, 2*alpha*beta + aN, beta*beta - p,N)
    if factors == False:
      continue
    [A1,B1,A2,B2] = factors
    if math.gcd(A1,N) == 1:
      X1 = -pow(A1, -1, N)*B1
      Y1 = (alpha*X1 + beta)%N
      if (Y1*Y1)%N == p:
        print(Y1)
    if math.gcd(A2,N) == 1:
      X2 = -pow(A2, -1, N)*B2
      Y2 = (alpha*X2 + beta)%N
      if (Y2*Y2)%N == p:
        print(Y2)

p = 41
N = 856

for a in primes:
  if a == p:
    next
  r = trysqrt(a*N,p)
  if r != False:
    tryKunerth(p,r,a,N)
  r = trysqrt(-a*N,p)
  if r != False:
    tryKunerth(p,r,-a,N)