Is Julia lazy?

New to Julia
61

I’m pretty sure sub-typing Z{P} by Integer is dangerous, as these are not the same algebraic structure. I’m almost sure addition and comparison are both in the interface (or at least deduced in Base from the interface). A function could assume that if n is an instance of Integer, then n < n+1 must hold. This kind of assumption can’t hold for a modular field. The best we can do with Base is I think Number (we would need to check the Interface, and I don’t see any documentation of the contract), but I guess there is almost no useful operations to gather from it. I’m afraid we are forced to re-implement almost every single functions acting on Z{P}, or to find a package which already have an interface more general than Z{P}, and with some basic operations already implemented.

Anyone have an idea why the interfaces of julia are not more documented ?

62

All of Julia’s native integers types are modular, so this is never reliable:

jl> Int8(127) + Int8(1)
-128
1 Like
63

I know, but all the Integer types in Julia almost satisfy the axioms of arithmetic (except for overflow problems). It depends on what users are willing to accept. A user needs to be aware of the risk overflow, and then, if the Int used by the user are small enough, he will be almost sure that all the functions he will called will return coherent result with the axioms of arithmetic. This is the same deal with Float64. Multiplications of Float64 are not necessarily associative, but the user can expect that functions defined for Real are almost associative. If the user is not willing to accept this risk, then he needs to use low level functions specialized for Float64, and manage the rounding errors.
For a modular field, this is a very wild guess to assume Integer functions will work “almost” as we think they will…

64

Then again, without proper documentation of the interface of Integer, every one can project is own idea on how Integer numbers should behave…

65

OK? Can you think of functionality do you foresee could be useful for Z{P}, but could also fail in a way that would not be easy to anticipate? Is this risk very different from the one that exists for other number types?

I have no vested interest in the Integer subtyping. I’m just thinking that Z{P} are integers, so it makes intuitive sense.

66

This is probably a bad example, and I did not looked it in depth, so it might be wrong, but if someone use a range p:q of Z{P} with q < p , it is reasonable to expect it to loop back to 0, then to q.
And if a function defined for Integer use ranges (which is very reasonable), this can give awkward results with Z{P}.

67

I had some issues with ranges, since a range Z13(0):Z13(12) appears to be empty. This had to be fixed manually, but I suspect that it’s an issue with ranges, really, that ought to calculate length with Ints.

But I don’t think p:q with q < p could be anything other than empty.

68

Maybe something like this is what you want?

Minimal definition of a prime field:

struct F{p, T<:Integer} <: Integer
    a::T
    F{p, T}(a::Integer) where {p, T<:Integer} = new{p, T}(mod(T(a), p))
end
F{p}(a::Integer) where p = F{p, typeof(a)}(mod(a, p))
F{p, S}(x::F{p, T}) where {S<:Integer, p, T<:Integer} = F{p, S}(x.a)

Base.promote_rule(::Type{F{p, T}}, ::Type{S}) where {p, T<:Integer, S<:Integer} =
    F{p, promote_type(T, S)}

Base.zero(::Type{F{p, T}}) where {p, T<:Integer} = F{p}(mod(zero(T), p))
Base.one(::Type{F{p, T}}) where {p, T<:Integer} = F{p}(mod(one(T), p))
for op in (:-, :+)
    @eval Base.$op(x::F{p}) where p = F{p}(mod($op(x.a), p))
end
for op in (:-, :+, :*)
    @eval Base.$op(x::F{p}, y::F{p}) where p = F{p}(mod($op(x.a, y.a), p))
end
Base.inv(x::F{p}) where p = F{p}(invmod(x.a, p))
Base.:/(x::F{p}, y::F{p}) where p = x * inv(y)
Base.:\(x::F{p}, y::F{p}) where p = inv(x) * y
Base.:(==)(x::F{p}, y::F{p}) where p = x.a == y.a
Base.:<(x::F{p}, y::F{p}) where p = x.a < y.a

Base.show(io::IO, x::F{p}) where p = print(io, "F", p, '(', x.a, ')')

F7 = F{7, BigInt}
x, y = F7(10), F7(-2)
@show(x, y, zero(x), one(x), +x, -x, x + y, x - y, x * y, x / y, x \ y, x^3, x^-5, x == F7(3), x == 3)
println()
A = F7[1 2; 3 4]
@show(A, 4A, A/4)
println()
using LinearAlgebra
L, U = lu(A)
@show(L, U, L * U, det(A), inv(A), inv(A) * A, A * inv(A));

Output:

x = F7(3)
y = F7(5)
zero(x) = F7(0)
one(x) = F7(1)
+x = F7(3)
-x = F7(4)
x + y = F7(1)
x - y = F7(5)
x * y = F7(1)
x / y = F7(2)
x \ y = F7(4)
x ^ 3 = F7(6)
x ^ -5 = F7(3)
x == F7(3) = true
x == 3 = true

A = F{7, BigInt}[F7(1) F7(2); F7(3) F7(4)]
4A = F{7, BigInt}[F7(4) F7(1); F7(5) F7(2)]
A / 4 = F{7, BigInt}[F7(2) F7(4); F7(6) F7(1)]

L = F{7, BigInt}[F7(1) F7(0); F7(5) F7(1)]
U = F{7, BigInt}[F7(3) F7(4); F7(0) F7(3)]
L * U = F{7, BigInt}[F7(3) F7(4); F7(1) F7(2)]
det(A) = F7(5)
inv(A) = F{7, BigInt}[F7(5) F7(1); F7(5) F7(3)]
inv(A) * A = F{7, BigInt}[F7(1) F7(0); F7(0) F7(1)]
A * inv(A) = F{7, BigInt}[F7(1) F7(0); F7(0) F7(1)]

Edit1: x^3 and x^-5 are added.
Edit2: Base.:(==) is defined.
Edit3: NoPivot() is deleted.

69

Great! There they are, some of the missing methods.

I believe these are redundant.

Also, there are quite a few redundant calls to mod, like here:

since mod is guaranteed by the inner constructor.

I’m also not sure about the promote rule:

It promotes the integer type, but it seems a bit pointless to turn F{7, Int8} + Int into F{7, Int}. No reason to promote Int8, I think.

That’s assuming promotion should be to F at all.

BTW. Is there something about the letter F? Why not Z, it seems to be what is used when I look around. Or maybe ℤ?

70

I’m pretty sure sub-typing Z{P} by Integer is dangerous

We have here a real problem. Let’s think a little deeper.
Maths creates objects by extension not by restriction. If one compares the Julia abstract hierarchy we see that top-down specialization and qualification increases.
The higher you are, the less things you can do… With this in mind promotion process is counter-intuitive and completely opposite to usual meaning.
So this organization is based on classical-basic maths number theory starting from naturals and ending with complex numbers.
In this hierarchy, some maths concepts find their place easily without doing anything. It is the case, for example of Gauss integers. Julia has naturally type for them:, and computing with them just consist in making a few tests.
For anyone not convinced try :

z=1+2im
println(typeof(z))

If you want to implement Gaussian euclidean division not necessary to create any type.
Other things are easily created (or need no creation at all).
It’s a well known trap that computations with rationals lead to overflow. Secure computation is provided by rationals where numerator and denominator are big.
Again Julia’s logic works perfectly

R=big(1)//big(2)
println(typeof(R))

p-fields and Hamilton’s quaternions are kind of UFOs, they don’t have natural place in the hierarchy.
So sub-typing is a real problem.
DNF who made a very interesting contributions hesitates between Integer and Number…
The author of the GaloisFields package creates own abstract type sub-typing Number as such

abstract type AbstractGaloisField <: Number end

and this being done

struct PrimeField{I<:Integer, p} <: AbstractGaloisField

I didn’t make any intervention about the order question in Fp-fields. For me the problem is quite non existing on a mathematical point of view. I’ve never seen anybody using any order on any Zp for any purpose.
I can see some interest on a programmatic point of view. If for example you want to test a perfectly symmetric property on all couples you can save time using a filter for … if k<h with the explicit convention that k<h means k.x<h.x but it’s just programming ‘cuisine’, nothing more than a cooking-book recipe.
To finish with this post I disagree with people who call C++ ‘insane’. C++ is like an old lady and deserves some respect. It is from another age and tried with some success to bring new programming concepts as OOP in a widely used language. C++ created some standards to inspire later creations. C++ is not outdated to solve problems where direct access to memory is needed or automatic garbage collecting dangerous.

71

100% agree. All due respect. But C++ template metaprogramming truly is insane and was discovered by accident. Of course, high performance numerical libraries in C++ wouldn’t be possible without it, and the newer libraries make SFINAE/etc. more palatable, but I think the inventors would probably agree it is nuts as it wasn’t designed into the language intentionally. Speaking of C++: I think the best book for this sort of thing is Stepanov and Rose “From Mathematics to Generic Programming” (https://www.fm2gp.com/) . I think virtually all of the lessons would apply to Julia, you just use different dispatching techniques (i.e. overloading-style dispatching more than template dispatching). If you don’t have it, my guess is that you would really enjoy it given your perspective.

3 Likes
72

Not to reinvent the wheel. Checking the library :

using GaloisFields
const F = @GaloisField 5  
squareroots(k)=[x for x in F if x*x==k] 
#just to check, that for example 3 is not a perfect square
for k in F
    println(squareroots(k))
end
#so we will construct F[√3]
const G, β =GaloisField(5, :β => [-3, 0, 1]) #minimal polynomial X^2-3

#checking that we have a field with 25 elements by finding the inverse of every non nul élément
for g in G
    if !iszero(g)
        print(g); print("  ×  "); print(g^-1);print("  =  ");println(g*g^-1)
    end
end
73

Don’t you need det, inv, etc. of a matrix for math educational purpose?

using GaloisFields, LinearAlgebra
GF7 = @GaloisField 7
B = GF7[1 2; 3 4]
det(B)

MethodError: no method matching abs(::𝔽₇)
...

inv(B)

MethodError: no method matching abs(::𝔽₇)
...

In the sense of my minimal definition example, minimality requires that inv and lu in LinearAlgebra.jl be available. I have tried to make the code small so that you can easily read it for using Julia in math education.

P.S. Iterator example (continuation of the minimal definition example above)

"""See https://docs.julialang.org/en/v1/manual/interfaces/"""
Base.iterate(Fp::Type{F{p, T}}) where {p, T<:Integer} = (zero(Fp), zero(T))
function Base.iterate(Fp::Type{F{p, T}}, state) where {p, T<:Integer}
    nextstate = state + 1
    nextstate < p ? (Fp(nextstate), nextstate) : nothing
end
Base.IteratorSize(Fp::Type{F{p, T}}) where {p, T<:Integer} = Base.HasLength()
Base.length(Fp::Type{F{p, T}}) where {p, T<:Integer} = p
Base.eltype(Fp::Type{F{p, T}}) where {p, T<:Integer} = Fp

squares(Fp) = Fp[x^2 for x in Fp]
squareroots(k, Fp) = Fp[x for x in Fp if x^2 == k]
@show(collect(F7), squares(F7), squareroots.(0:6, Ref(F7)));

Output:

collect(F7) = F{7, BigInt}[F7(0), F7(1), F7(2), F7(3), F7(4), F7(5), F7(6)]
squares(F7) = F{7, BigInt}[F7(0), F7(1), F7(4), F7(2), F7(2), F7(4), F7(1)]
squareroots.(0:6, Ref(F7)) = Vector{F{7, BigInt}}[[F7(0)], [F7(1), F7(6)], [F7(3), F7(4)], [], [F7(2), F7(5)], [], []]
74

You can use AbsteactAlgebra.jl to deal with elementary abstract algebra:

using AbstractAlgebra

@show GF7 = GF(7)
x, y = GF7(10), GF7(-2)
@show(x, y, zero(x), one(x), -x, x + y, x - y, x * y, x^3, x^-5, x == GF7(3), x == 3)
println()

@show C = GF7[1 2; 3 4]
@show P, x = PolynomialRing(GF7, "x")
@show(det(C), inv(C), lu(C), charpoly(P, C))
println()

squares(Fp) = [x^2 for x in Fp]
squareroots(k, Fp) = [x for x in Fp if x^2 == k]
@show(collect(GF7), squares(GF7), squareroots.(0:6, Ref(GF7)));

Output:

GF7 = GF(7) = Finite field F_7
x = 3
y = 5
zero(x) = 0
one(x) = 1
-x = 4
x + y = 1
x - y = 5
x * y = 1
x ^ 3 = 6
x ^ -5 = 3
x == GF7(3) = true
x == 3 = true

C = GF7[1 2; 3 4] = [1 2; 3 4]
(P, x) = PolynomialRing(GF7, "x") = (Univariate Polynomial Ring in x over Finite field F_7, x)
det(C) = 5
inv(C) = [5 1; 5 3]
lu(C) = (2, (), [1 0; 3 1], [1 2; 0 5])
charpoly(P, C) = x^2 + 2*x + 5

collect(GF7) = AbstractAlgebra.GFElem{Int64}[0, 1, 2, 3, 4, 5, 6]
squares(GF7) = AbstractAlgebra.GFElem{Int64}[0, 1, 4, 2, 2, 4, 1]
squareroots.(0:6, Ref(GF7)) = Vector{AbstractAlgebra.GFElem{Int64}}[[0], [1, 6], [3, 4], [], [2, 5], [], []]
2 Likes
75

I already use Julia in math education with 230 programs in my site . This job is not finished yet concerning only for now Bases (Sets, relations, applications, operations, binary logic) , numbers and linear algebra.
I give examples of vector spaces with base field R, Q, C and Zp.
Interactive Javascripts applets are proposed allowing to fix base field and order as well
Of course I need det and inv for linear systems.
As I said earlier the best models for teaching are Zp and Q. But since all students are not familiar with modular arithmetic, I choose preferably Q.
For det, for example I show programs using the definition with permutations, development following one line or one column, Gauss Pivot etc. (this page)
For sample programs till now I use mainly Floats and Rationals but I plan to give exercises with Zp.
My purpose is to write functions ad-hoc to interpret algorithms and understand process and to use libraries as well not to reinvent the wheel.
For ‘usual’ matrices i already use Julia for LU, QR, Cholesky, Jacobi, Gauss-Seidel algorithms both with hand-made functions and LinearAlgebra.
GaloisFields is really powerful, I’m impressed. See above how it built decomposition field of some irreducible second degree polynomial over Z5.
In any case thank you for your help. I will take your suggestions into account in my future contributions.
Have a good week-end and thank you again for your help.

76

None of them, you should subtype Number.

The hierarchy for Integer goes Integer <: Real <: Number <: Any so any Integer is also Real and therefore Integer represents the rational integers, ℤ.

So while your type is for a ring of integers in the broader number theoretic sense, they are not Real.

As an example of what goes wrong, I think it’s reasonable to assume that < and abs are part of the interface for Real, but the integers mod P are neither an ordered ring nor an Archimedean valued ring. So code that tried to use this interface on your type would break.

77

What about Unsigned?

Edit: You imply ‘not Unsigned, because it’s an Integer’, I guess.

78

Unsigned <: Integer so the same comments apply.

79

I thought so. I just missed it in the hierarchy you drew up.

80

I read carefully your code. It was very useful and instructive.
Official doc about “value types” is less than minimal :slight_smile:
Fortunately some articles exist on the subject which were helpful even not so easy to read.
Finally I dunno if Julia is lazy, but one thing is sure, I am !!!
So I decided once for all to use GaloisFields library it’s very well done.
I began to update a few of my pages to include Julia programs computing linear combinations, products of Fp matrices etc…
example here, and there
I will continue next week with inv, det and system resolutions.
Later one we’ll begin with affine geometry with Fp fields. It is strange, segments without middle points for example, taxi-cab distance, and so on… I must combine with graphics.
Thank you so much again for your help