Primality and Factoring

import numpy
from EuclideanGCD import *
from math import ceil, floor, exp, log, sqrt
from decimal import *
from random import *

getcontext().prec = 800

# Simbolo di jacobi


def Jacobi(a, m):
    a = a % m
    t = 1
    while a != 0:
        while a % 2 == 0:
            a = a // 2
            if m % 8 in [3, 5]:
                t = -t
        (a, m) = (m, a)
        if (a % 4 == m % 4) and (m % 4 == 3):
            t = -t
        a = a % m
    if m == 1:
        return t
    return 0


# TEST DI PRIMALITA'


def Fermat(n):
    for i in range(1, n):
        if pow(i, n - 1, n) != 1:
            print("Witness=", i, "\n")
            return False
    return True


def SolvayStrassen(n):
    if n == 2:
        return True
    if n > 2 and (n % 2 == 0):
        return False
    for i in range(2, int((n + 3) // 2)):
        if (
            pow(i, (n - 1) // 2, n)
            != (n + Jacobi(i, n)) % n
        ):
            print("Witness=", i)
            return False
    return True


def SPprime(n, a):
    # Strong Pseudoprime mod a
    t = n - 1
    s = 0
    while t % 2 == 0:
        s += 1
        t = t // 2
    # print "n=1+",t,"*(2**",s,")"
    b = pow(int(a), int(t), int(n))
    if (b == 1) or (b == n - 1):
        return True
    for j in range(1, s):
        b = pow(b, 2, n)
        if b == n - 1:
            return True
    return False


def MillerTest(n):
    # Questo test di primalita' dipende da ERH
    # altrimenti e' probabilistico
    W = min(round(2 * log(n * 1.0) ** 2), n - 1)
    for a in range(2, int(W + 1)):
        if not (SPprime(n, a)):
            #       print "a=",a
            return False
    return True


PrimeP = MillerTest

# Costruzione primi


def NextPrime(N):
    i = N
    while i < 2 * N + 2:
        i += 1
        if PrimeP(i):
            return i


def PrimesUpTo(B):
    if B < 2:
        return ()
    i = 2
    r = []
    while i < B + 1:
        r.append(i)
        i = NextPrime(i)
    return r


# Fattorizzazione con Pollard Rho


def PollardRho(n: int) -> int:
    maxiter = 100000
    for t in range(0, maxiter):
        print("iter=", t, "\n")
        a = int(randrange(1, n - 3 + 1))
        s = int(randrange(0, n - 1 + 1))
        (U, V) = (s, s)
        print("a=", a, " s=", s, "\n")

        def F(x):
            return (pow(int(x), 2) + a) % n

        g = 1
        while g == 1:
            U = F(U)
            V = F(V)
            V = F(V)
            T = abs(int(U - V))
            print(
                "U=",
                U,
                "\n",
                "V=",
                V,
                "\n",
                "U-V=",
                T,
                "\n",
            )
            h = XEuclidean(T, n)
            print(h)
            g = h[2]
        if g != n:
            return g
    return False