scheinerman / Mods.jl

Easy modular arithmetic for Julia
MIT License
33 stars 5 forks source link

Modular exponentiation comes out wrong: Mod{1618034}(1000000)^367730 #19

Closed phma closed 11 months ago

phma commented 11 months ago

I'm working on a repo (not public yet) which is a Julia implementation of wring-twistree. I updated in Pkg, which updated Mods from 1.3.3 to 2.0.0, then ran the tests, and got an early failure: Testing Running tests... Test Failed at /home/phma/src/WringTwistree/test/runtests.jl:3 Expression: findMaxOrder(85) == 54 Evaluated: 1 == 54 You can look at the Rust and Haskell code in https://github.com/phma/wring-twistree . This test is in Rust, but not Haskell.

I then did the next test in the Debugger and got this in the debugging session:

1|debug> so #exiting from the function (==) In findMaxOrder(n) at /home/phma/src/WringTwistree/src/Mix3.jl:57 63 m=start+i÷2+1 64 else 65 m=start-i÷2 66 end

67 if n==1 || isMaxOrder(n,car,fac,m) 68 return m 69 end 70 end 71 one(n) 72 end

About to run: <(WringTwistree.Mix3.isMaxOrder)(1618034, 367730, Primes.Factorization(2 => 1, 5 => 1, 11 => 1, 3343 =...> 1|debug> s In isMaxOrder(modl, car, fac, n) at /home/phma/src/WringTwistree/src/Mix3.jl:43 46 # fac is the factorization of car, 47 # and n is the number being tested. 48 # Returns true if n has maximum order, which implies it's a primitive root 49 # if modulus has any primitive roots.

50 ret=Mod{modl}(n)^car==1 51 for (p,ex) in fac 52 ret=ret&&Mod{modl}(n)^(car÷p)!=1 53 end 54 ret 55 end

About to run: (Core.apply_type)(Mods.Mod, 1618034) 1|debug> fr [1] isMaxOrder(modl, car, fac, n) at /home/phma/src/WringTwistree/src/Mix3.jl:43 | modl::Int64 = 1618034 | car::Int64 = 367730 | fac::Primes.Factorization{Int64} = Primes.Factorization(2 => 1, 5 => 1, 11 => 1, 3343 => 1) | n::Int64 = 1000000 1|julia> Mod{modl}(n)^car Mod{1618034}(376832110955397120) This is invalid; the number in parentheses should be no bigger than the modulus.

Here's the module: module Mix3 using Primes,Mods,Base.Threads export carmichael,findMaxOrder,mix3Parts!

function mix3(a::Integer,b::Integer,c::Integer) mask=(a|b|c)-(a&b&c) (a⊻mask,b⊻mask,c⊻mask) end

function fiboPair(n::Integer) f0=0 f1=1 while f0<=n (f0,f1)=(f1,f1+f0) end (f0,f1) end

function searchDir(n::Integer) (num,den)=fiboPair(n) (q,r)=divrem(nnum,den) if r2<den (q,1) else (q+1,-1) end end

function carmichael(n::Integer) facs=factor(n) ret=one(n) for (p,ex) in facs if p==2 && ex>=3 carfac=p^(ex-2) else carfac=p^(ex-1)*(p-1) end ret=lcm(ret,carfac) end ret end

function isMaxOrder(modl::Integer,car::Integer, fac::Primes.Factorization{<:Integer},n::Integer)

modl is the modulus, car is its Carmichael function,

fac is the factorization of car,

and n is the number being tested.

Returns true if n has maximum order, which implies it's a primitive root

if modulus has any primitive roots.

ret=Mod{modl}(n)^car==1 for (p,ex) in fac ret=ret&&Mod{modl}(n)^(car÷p)!=1 end ret end

function findMaxOrder(n::Integer) car=carmichael(n) fac=factor(car) (start,dir)=searchDir(n) for i in 0:n if (i&1)==1 m=start+i÷2+1 else m=start-i÷2 end if n==1 || isMaxOrder(n,car,fac,m) return m end end one(n) end

function mix3Parts!(buf::Vector{<:Integer},rprime::Integer) len=div(length(buf),3) a=1 b=2len c=2len+1 while len>0 && a<b @inbounds (buf[a],buf[b],buf[c])=mix3(buf[a],buf[b],buf[c]) a+=1 b-=1 c+=rprime if c>3*len c-=len end end end

end

scheinerman commented 11 months ago

Yes. I am so sorry!! I am working on fixing a bunch of problems. I should not have release 2.0.0 yet.

scheinerman commented 11 months ago

@phma I'm not sure where the problem is. Version 2.0.1 gives the correct answer for your question:

julia> m = 1618034
1618034

julia> p = 367730
367730

julia> a = 1000000
1000000

julia> powermod(a,p,m)
809018

julia> Mod{m}(a)^p
Mod{1618034}(809018)

I'm not saying there aren't problems with version 2.0.x, but I think this is OK. Let me know if you can narrow the problem down. Again, I'm sorry for releasing this too soon.

phma commented 11 months ago

I just upgraded to 2.0.2 and the test passes. I also published the repo: https://github.com/phma/WringTwistree.jl .

scheinerman commented 11 months ago

@phma Phew! Glad to hear. I'm working on the re-implementation of the GaussMod type, but I sort of doubt it's of much use to people.