Playing around with IntMod

pitao · December 15, 2016, 9:27pm

I was having a play with https://github.com/JuliaLang/julia/blob/master/examples/modint.jl

I tried adding a few goodies (rand, ^ and /):

import Base: +, -, *, ^, /

immutable ModInt{n} <: Integer
    k::Int
    ModInt(k) = new( k ≥ 0 ? k%n : (k%n + n) ) # fixed so always ≥ 0
end

-{n}(a::ModInt{n}) = ModInt{n}(-a.k)
+{n}(a::ModInt{n}, b::ModInt{n}) = ModInt{n}(a.k+b.k)
-{n}(a::ModInt{n}, b::ModInt{n}) = ModInt{n}(a.k-b.k)
*{n}(a::ModInt{n}, b::ModInt{n}) = ModInt{n}(a.k*b.k)
^{n}(a::ModInt{n}, b::ModInt{n}) = ModInt{n}( b.k == 0 ? 1 : a.k^(b.k-1) )


Base.convert{n}(::Type{ModInt{n}}, i::Int) = ModInt{n}(i)
Base.promote_rule{n}(::Type{ModInt{n}}, ::Type{Int}) = ModInt{n}

Base.show{n}(io::IO, k::ModInt{n}) = print(io, "$(k.k) mod $n")
Base.showcompact(io::IO, k::ModInt) = print(io, k.k)

Base.inv{n}(a::ModInt{n}) = ModInt{n}(invmod(a.k, n))

/{n}(a::ModInt{n}, b::ModInt{n}) = ModInt{n}( a*inv(b) )

Base.rand{N}(rng::AbstractRNG, T::Type{ModInt{N}}) = T(rand(rng, 0:N-1))

However,

k = ModInt{7}(3)
1*k  # works
1/k  # ... fails with:

LoadError: MethodError: Cannot `convert` an object of type ModInt{7}
    to an object of type AbstractFloat
This may have arisen from a call to the constructor AbstractFloat(...),
since type constructors fall back to convert methods.
while loading In[4], in expression starting on line 1

 in /(::Int64, ::ModInt{7}) at ./int.jl:35

I can’t see why it successfully handles 1 * k but fails with 1 / k. There is a promotion_rule in place! How can it possibly succeed in one case and fail in the other?

PS Other tests:

A = ModInt{7}[1 2; 3 4]
A^2

B = rand(ModInt{7}, 3, 3)
B * inv(B)

kristoffer.carlsson · December 16, 2016, 8:16am

Because of:

@edit 1/k

takes you to https://github.com/JuliaLang/julia/blob/2d8f5bfc05b1894b160b9fa3597f0e1fe54d1b93/base/int.jl#L35

pitao · December 16, 2016, 9:19am

I can’t understand why it is attempting to use this overload.

I have defined a convert and a promote_rule from Integer → ModInt{k} but there is nothing for going the other way.

So in the case of 1/k how come it is not promoting (Integer,ModInt{k}) to (ModInt{k},ModInt{k}) and finding the / overload I have provided?

Instead it is finding:

/(x::Integer, y::Integer) = float(x)/float(y)  # in int.jl

… which means it must be promoting (demoting?!) k from ModInt{k} to int.

kristoffer.carlsson · December 16, 2016, 9:27am

There is a method defined for / where both types are a subtype of Integer. Both 1 and k are of types that subtype Integer and no more specific method exists, hence the method I linked is the one that will be executed.

Julia does not do implicit conversions. You need to add your own explicit method,/(a::Int, b::ModInt) = /(promote(a,b)...) or something like that.

pitao · December 16, 2016, 9:48am

I see, thanks!

Looking in conversion.jl line 190,

+(x::Number, y::Number) = +(promote(x,y)...)
*(x::Number, y::Number) = *(promote(x,y)...)  # A
-(x::Number, y::Number) = -(promote(x,y)...)
/(x::Number, y::Number) = /(promote(x,y)...)  # B
^(x::Number, y::Number) = ^(promote(x,y)...)

So 1*k finds (A) whereas 1/k finds:

# int.jl line 35
/(x::Integer, y::Integer) = float(x)/float(y)

i.e. a specialisation of (B) that exists because Integer / Integer typically produces a nonInteger.

PS @edit is really useful!

kristoffer.carlsson · December 16, 2016, 10:18am

Your analysis is spot on.

pitao · December 17, 2016, 4:13pm

I have extended this file to handle ^ and / and would like to offer it back as a PR. However, I’m getting a vexing intermittent fault: see (*).

It is happening about 50% the time. Can anyone see what’s going wrong?

# This file is a part of Julia. License is MIT: http://julialang.org/license

module ModInts
export ModInt, test

import Base: +, -, *, ^, /

immutable ModInt{n} <: Integer
    k::Int
    ModInt(k) = new( k ≥ 0 ? k%n : (k%n + n) ) # fixed so always ≥ 0
    ModInt(k::ModInt{n}) = new( k.k )
end

-{n}(a::ModInt{n}) = ModInt{n}(-a.k)
+{n}(a::ModInt{n}, b::ModInt{n}) = ModInt{n}(a.k+b.k)
-{n}(a::ModInt{n}, b::ModInt{n}) = ModInt{n}(a.k-b.k)
*{n}(a::ModInt{n}, b::ModInt{n}) = ModInt{n}(a.k*b.k)
^{n}(a::ModInt{n}, b::ModInt{n}) = ModInt{n}( b.k == 0 ? 1 : a.k^(b.k-1) )


Base.convert{     n}(::Type{ModInt{n}}, i::     Int ) = ModInt{n}(i)
Base.promote_rule{n}(::Type{ModInt{n}},  ::Type{Int}) = ModInt{n}

Base.show{n}(    io::IO, k::ModInt{n}) = print(io, "$(k.k) mod $n")
Base.showcompact(io::IO, k::ModInt   ) = print(io, k.k)


# https://discourse.julialang.org/t/playing-around-with-intmod/980
Base.inv{n}(a::ModInt{n}) = ModInt{n}(invmod(a.k, n))

/{n}(a::ModInt{n}, b::ModInt{n}) = ModInt{n}( a*inv(b) )
/{n}(a::ModInt{n}, b::Integer  ) = /(promote(a,b)...)
/{n}(a::Integer  , b::ModInt{n}) = /(promote(a,b)...)


Base.rand{N}(rng::AbstractRNG, T::Type{ModInt{N}}) = T(rand(rng, 0:N-1))

macro disp(x)
    quote
        $(string(x)) |> display
        $(   esc(x)) |> display  # <-- line 41
    end
end


function test()
    @disp k = ModInt{7}(3)
    @disp 1/k
    @disp k * (1/k)
    
    @disp A = ModInt{7}[1 2; 3 4]
    @disp A^2

    @disp B = rand(ModInt{7}, 3, 3)
    @disp B * inv(B)                     # (*) fails often!
    @disp B * B^(-1)
end

ModInts.test()

end # module

Here’s the error:

Line 41 is marked above.

pfitzseb · December 17, 2016, 4:37pm

Well, inv calls the generic fallback for lu factorization, which does pivoting by default, which in turn requires comparison of elements via < (actually, via > in lu.jl:36, but that falls back on <). The code should therefore work if you define <{n}(a::ModInt{n}, b::ModInt{n}).

pitao · December 17, 2016, 5:21pm

Thanks! I can see why that might give rise to an inconsistent failure – it is a random matrix, so it is probably choosing a different pivot point each time.

Now for the tricky part: what on earth is a < operator going to do using clock arithmetic!

Is p-1 < 0? They are right next to one another, and p-1 is to the left of 0 just as 0 is to the left of 1.

Maybe I should do it by shortest distance…

pitao · December 18, 2016, 11:38pm

ok I’ve revised the code (now works):

import Base: /

immutable ModInt{n} <: Integer
    k::Int
    ModInt(k::Int) = new( k ≥ 0 ? k%n : (k%n + n) ) # fixed so always ≥ 0
    ModInt(k::ModInt{n}) = new( k.k ) # (B)
end

Base.convert{     n}(::Type{ModInt{n}}, i::     Int ) = ModInt{n}(i)
Base.promote_rule{n}(::Type{ModInt{n}},  ::Type{Int}) = ModInt{n}

Base.inv{n}(a::ModInt{n}) = ModInt{n}(invmod(a.k, n))

/{n}(a::ModInt{n}, b::ModInt{n}) = ModInt{n}( a*inv(b) )
/{n}(a::ModInt{n}, b::Integer  ) = /(promote(a,b)...)
/{n}(a::Integer  , b::ModInt{n}) = /(promote(a,b)...)

# test
k = ModInt{7}(3)
1 / k

^ strangely constructor (B) is required.
Which means that promote acting on (a::Integer , b::ModInt{n}) is having to ‘promote’ ModInt{n} -> ModInt{n} which looks inefficient.
Why not just use the original value?
(it seems that Julia is performing a needless object copy)

Also it doesn’t quite seem aesthetically right that I need that Base.convert line. Shouldn’t Julia be able to infer from the fact that ModInt{n} has a constructor that takes an Int?

StefanKarpinski · December 19, 2016, 11:35am

I wonder if we should add that fallback to Base. It would probably produce less confusing behavior in just such cases – i.e. a definition is provided for dividing a type by itself.

StefanKarpinski · December 25, 2016, 6:29am

Issue opened: https://github.com/JuliaLang/julia/issues/19714. @pitao, if you feel like taking a crack at it, this is a pretty straightforward change.

pitao · December 25, 2016, 9:05am

Thanks – I’m on it! This is just the kind of stuff I need right now.

StefanKarpinski · January 3, 2017, 4:18pm

Fixed. This ended up being trickier than expected because of needing to avoid performance regressions for built-in types but in the end various division operations involving BigInts actually ended up getting faster. Sorry to have snagged the issue, @pitao – I wanted to get it fixed in time for 0.6.