modint example is broken

[StefanKarpinski]

Some recent change to linear algebra dispatch broke the classic ModInt example. We should probably go ahead and make that a test just so that we don't break it. It seems that now multiplication of ModInt matrices ends up using generic_matmatmul which expects zero to be defined for the element type.

[JeffBezanson]

The fallback zero expects to be able to convert 0.0, and ModInt only implements conversion from integers. Maybe we should add another fallback zero for integer types that converts 0.

[StefanKarpinski]

When did that change? Why not convert the integer 0 instead of the float 0.0?

[JeffBezanson]

@andreasnoackjensen changed it since calls like zero(Any) caused problems. We figured Float64 is a better default type than Int.

[StefanKarpinski]

Limiting zero to numbers seems reasonable but changing it to Float64 seems unnecessary.

[StefanKarpinski]

Logically, if you can convert a Float64 to some type, then you can convert an Int to that type too, but the converse does not hold at all.

[JeffBezanson]

Yeah, I think I'm ok with changing it back to an Int.

[andreasnoack]

JuliaLang/julia#6183 specifically changed the default to Float64. As is clear from the discussion there, you can argue for both, but I think it is a little fast to revert a pull request within an hour after an issue has been opened. (It could be that it was five in the morning for someone interested in the issue.)

Float64 arrays accept integer values, but Int arrays don't accept Float64 values and therefore I argued for using Float64 as the default for one(Real). In this case it wouldn't work to use float(one(Real)) because then I would get floating point results from integer arithmetic. The whole business of zero/one(AbstractType) is a bit problematic and I think it is reasonable to require that you define zero and one when defining a new number type. That would make the ModInt multiplication work.

[ivarne]

See also JuliaLang/julia#4808

[StefanKarpinski]

The previous behavior that we've had forever worked very well for almost all scalar situations. This was introduced and started breaking code within an hour. Let's make an actual discussion with examples of various problematic situations for defaulting to Int versus Float64 – that other issue just has a bunch of links to other even longer and less clear discussions and no concrete examples anywhere in sight.

[andreasnoack]

The motivating example was something like

x = Real[1.1, 1.2]

and I want to find the return type of x[1]*x[2] + x[1]*x[2] because I need to store the value in an array. Now

julia> y=x[1]*x[2]+x[1]*x[2];
julia> typeof(y)
Float64

but

julia> T=typeof(one(eltype(x))*one(eltype(x))+one(eltype(x))*one(eltype(x)))
Int64

and therefore

julia> T[y]
ERROR: InexactError()
 in getindex at array.jl:113

In the pull request, I advertised that the change could break code with type inference problems, but there I mainly thought about the Any problem.

An alternative to reverting could have been

zero{T<:Integer}(::Type{T}) = oftype(T, 0)

[StefanKarpinski]

That may be ok, but I'm not entirely convinced that Float64 is really the right default for non-integer reals – that's problematic for Rationals, for example.

[andreasnoack]

I have thought a little more about this. Using Float64 as default for zero(Real) and one(Real) is not a general solution for handling type inference for non-contrete arrays. Furthermore it required specialized method of zero and one for rationals and, as discussed here, integers are also problematic. The methods for rationals are now redundant, so I'll remove them.

[StefanKarpinski]

I agree that the current situation is problematic but I don't think defaulting to Float64 instead of Int is a real solution (we're in agreement). Since that doesn't really fix the matter, I'd rather not break code.