Kongruenssiyhtalo.jl

# Multiplicative order Julia
# Juhani Kaukoranta 21.5.2020
using Primes

function factors(n)
    f = [one(n)]
    for (p,e) in factor(n)
        f = reduce(vcat, [f*p^j for j in 1:e], init=f)
    end
    return length(f) == 1 ? [one(n), n] : sort!(f)
end

function multorder(a, m)
    gcd(a,m) == 1 || error("$a and $m are not coprime")
    res = one(m)
    for (p,e) in factor(m)
        m = p^e
        t = div(m, p) * (p-1)
        for f in factors(t)
            if powermod(a, f, m) == 1
                res = lcm(res, f)
                break
            end
        end
    end
    return res
end
# multiplicative order of 2 modulo primes
b= primes(1,10000000);
max = length(b);

function sopivat_alkuluvut(max)
    properprimes = Int64[]  # kongruenssiyhtälölle kelvolliset kanditaatit
    s = 0
    for i in 2 : max
        if multorder(2,b[i]) % 2 != 0 && (b[i] % 8 == 1 || b[i] % 8 == 7)
            #println(b[i])
            s += 1
            push!(properprimes,b[i])
        end
    end
    print("sopivia alkulukuja oli ",s)
    return properprimes
end
# testaan kelvollisilla alkulukukanditaateilla p kongruenssiyhtälö
# testataan n = ((2k+1)*ord(2,p))+2)*p yhtälölle 2^(n-2) = 1 (mod n)
properprimes = sopivat_alkuluvut(max);
nnumbers = Int64[]
nmaxproper =length(properprimes)
for i in 1 : nmaxproper
    mult = multorder(2,properprimes[i])
    for k in 1 : 100000
        n = ((2*k+1)*mult+2)*properprimes[i]
        if powermod(2,n-2,n)==1
            #println(n," factors: ",factor(n))
            push!(nnumbers,n)
        end
    end
end
nsorted = sort(nnumbers);
nlukuja = length(nsorted)
for i in 1 : nlukuja
    println(nsorted[i],"  ",factor(nsorted[i]))
end