Kongruenssiyhtalo3.jl

# Kongruenssiyhtälön 2^(n-2) = 1 (mod n) ratkaisut
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

#b= primes(1,1000000);
b = primes(1,100000000);
maksimi = length(b);
# multiplicative order of 2 modulo primes


function sopivat_alkuluvut(maksimi)
    properprimes = Int64[]  # kongruenssiyhtälölle kelvolliset kanditaatit
    multord = Int64[]   # kelvollisten multiplicative order of 2
    # multorder(2,p) pariton, p jaettuna 8:lla antaa jakojäännöksen 1 tai 7
    s = 0
    for i in 2 : maksimi
        #if multorder(2,b[i]) % 2 != 0 && (b[i] % 8 == 1 || b[i] % 8 == 7)
        temp1 = b[i] % 8
        if temp1 == 1 || temp1 == 7
            temp2 = multorder(2,b[i])
            if temp2 % 2 != 0
                s += 1
                push!(properprimes,b[i])
                push!(multord,temp2)
            end
        end
    end
    println("sopivia alkulukuja oli ",s)
    return properprimes,multord
end
# testaan kelvollisilla alkulukukanditaateilla p kongruenssiyhtälö
# testataan n = q*p =(k*ord(2,p))+2)*p yhtälölle 2^(n-2) = 1 (mod n)
properprimes,multord = sopivat_alkuluvut(maksimi);
nnumbers = Int64[]
nmaxproper =length(properprimes)
half1 = div(nmaxproper,2)+1  # alkulukujen puoliväli
half2 = half1-1   # toinen puolisko
# lasketaan ensin suurimpien alkulukutekiojöiden avulla
# näin saadaan sisäsilmukan k:t pieniksi
for i in nmaxproper : -1 : half1      # aloitetaan suurimmista alkuluvuista
    mult = multord[i]
    for k in 1 : 2: 10000
        n = (k*mult+2)*properprimes[i]
        if powermod(2,n-2,n)==1
            #println(n," factors: ",factor(n))
            push!(nnumbers,n)
        end
    end
end
for i in half2 : -1 : 1  # jatketaan puolivalista alkuun
    mult = multord[i]
    for k in 1 : 2: 10000
        n = (k*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)

# poistetaan n-tuloslistan toistot
# puhdistettu lista nlukujavalmis
nlukujavalmis = Int64[nsorted[1]]
for i in 1 : nlukuja-1
    if nsorted[i] != nsorted[i+1]
        push!(nlukujavalmis,nsorted[i+1])
    end
end
nlukuja_valmis = length(nlukujavalmis);
# tulostetaan valmiiden lista
for i in 1 : length(nlukujavalmis)
    println(i,"  ",nlukujavalmis[i],"  ",factor(nlukujavalmis[i]))
end