Julia's built-in numeric types include a wide range of
π
π = 3.1415926535897...
typeof(π)
Irrational{:π}
Float64(π)
3.141592653589793
pi = BigFloat(π)
3.141592653589793238462643383279502884197169399375105820974944592307816406286198
exp(pi*1im)
-1.000000000000000000000000000000000000000000000000000000000000000000000000000000 + 1.096917440979352076742130626395698021050758236508687951179005716992142688513354e-77im
exp(π*1im)
-1.0 + 1.2246467991473532e-16im
typeof(ans)
Complex{BigFloat}
Rational(Float32(pi))
13176795//4194304
@show Float32(13176795/4194304);
@show Float32(pi);
Float32(13176795 / 4194304) = 3.1415927f0 Float32(pi) = 3.1415927f0
5.0 + 3im
5.0 + 3.0im
typeof(ans)
Complex{Float64}
5//3 + 7//9im
5//3 - 7//9*im
typeof(ans)
Complex{Rational{Int64}}
subtypes(Number)
2-element Array{Any,1}: Complex{T<:Real} Real
subtypes(Real)
4-element Array{Any,1}: AbstractFloat Integer Irrational{sym} Rational{T<:Integer}
subtypes(AbstractFloat)
4-element Array{Any,1}: BigFloat Float16 Float32 Float64
Int32 <: Number # are Int64's numbers?
true
Int64 <: Real # are Int64's a subset of the reals?
true
Int64 <: Float64 # are Int64's a subset of 64-bit floating-point numbers?
false
17 ∈ 0:3:21 # is 17 an element of the set 0,3,6,9,12,15,18,21? (a set operation)
false
typeof(ans)
Bool
illustrates linear algebra over arbitrary-precision rationals
A = convert(Matrix{Rational{BigInt}}, rand(1:10,5,5))/10
5x5 Array{Rational{BigInt},2}: 2//5 9//10 4//5 3//5 4//5 3//10 1//1 4//5 1//10 3//5 3//10 1//2 3//10 7//10 3//10 2//5 1//1 3//10 7//10 9//10 1//2 1//10 4//5 3//5 3//10
x = [3//4; 17//3; -1//13; -7//11 ; 3//19]
b = A*x
5-element Array{Rational{BigInt},1}: 69052//13585 382199//65208 859823//326040 229868//40755 177913//326040
x̂ = A\b
5-element Array{Rational{BigInt},1}: 3//4 17//3 -1//13 -7//11 3//19
x - x̂
5-element Array{Rational{BigInt},1}: 0//1 0//1 0//1 0//1 0//1
Note that Julia's backslash operator (and LU decomp) works over all its numeric types
Float32 and Float64 LU, QR, SVD, etc. are calls to LAPACK
example from Andreas Noack, CSAIL MIT http://andreasnoack.github.io/talks/2015AprilStanford_AndreasNoack.ipynb
# Scalar finite fields. P is the modulus, T is the integer type (Int16, Int32, ...)
immutable GF{P,T<:Integer} <: Number
data::T
function GF(x::Integer)
return new(mod(x, P))
end
end
immutable: In Julia variables change values, but numbers do not.
# basic methods for scalar finite field
import Base: convert, inv, one, promote_rule, show, zero
function call{P}(::Type{GF{P}}, x::Integer)
if !isprime(P)
throw(ArgumentError("P must be a prime"))
end
return GF{P,typeof(x)}(mod(x, P))
end
convert{P,T}(::Type{GF{P,T}}, x::Integer) = GF{P}(x)
convert{P}(::Type{GF{P}}, x::Integer) = GF{P}(x)
convert{P,T}(::Type{GF{P,T}}, x::GF{P}) = GF{P,T}(x.data)
promote_rule{P,T1,T2<:Integer}(::Type{GF{P,T1}}, ::Type{T2}) = GF{P,promote_type(T1,T2
)}
show(io::IO, x::GF) = show(io, x.data)
show (generic function with 107 methods)
import Base: +, -, *, /
for op in (:+, :-, :*)
@eval begin
($op){P,T}(x::GF{P,T}, y::GF{P,T}) = GF{P,T}($(op)(x.data, y.data))
end
end
x, y = GF{5}(9), GF{5}(8)
@show x
@show y
@show x + y
@show x - y
@show x * y
;
x = 4 y = 3 x + y = 2 x - y = 1 x * y = 2
# Division requires slightly more care
function inv{P,T}(x::GF{P,T})
if x == zero(x)
throw(DivideError())
end
r, u, v = gcdx(x.data, P)
GF{P,T}(u)
end
(/){P}(x::GF{P}, y::GF{P}) = x*inv(y)
;
@show x / y
@show x \ y # backslash on any Number is defined in terms of /, so we get it autmomatically
;
x / y = 3 x \ y = 2
# create 4x4 matrix of random GF(5) elems
srand(1234)
A = [GF{5}(rand(0:4)) for i = 1:4, j = 1:4]
4x4 Array{GF{5,Int64},2}: 3 1 1 0 1 3 4 3 4 3 2 2 0 0 1 1
b = [GF{5}(rand(0:4)) for i = 1:4]
4-element Array{GF{5,Int64},1}: 4 1 2 0
x = A\b
4-element Array{GF{5,Int64},1}: 2 3 0 0
A*x - b
4-element Array{GF{5,Int64},1}: 0 0 0 0
typeof(x)
Array{GF{5,Int64},1}
methods(+)
@which x + y
@which GF{7}(4) + GF{7}(1)
@which 4+8
matrix-matrix mult
A1, A2 = rand(1:100, 100, 100), rand(1:100, 100, 100)
A1*A2 # warm up to be sure function is compiled
print("int mat mult ")
@time A1*A2
AF1, AF2 = map(GF{5}, A1), map(GF{5}, A2)
AF1*AF2
print("GF(p) mat mult ")
@time AF1*AF2
;
int mat mult 0.001053 seconds (157 allocations: 88.651 KB) GF(p) mat mult 0.009569 seconds (13 allocations: 78.656 KB)
LU factorization: Float64 via LAPACK, GF(p) via generic LU algorithm
lufact(A1)
print("Float64 mat lufact ")
@time lufact(A1) # Promoted to Float64 and calls LAPACK
F = lufact(AF1,Val{false})
while F.info != 0
AF1[F.info, F.info] += 1
F = lufact(AF1, Val{false})
end
lufact(AF1)
print("GF(p) mat lufact ")
@time lufact(AF1) # Non-blocked generic LU implemented in Julia
;
Float64 mat lufact 0.000415 seconds (8 allocations: 79.250 KB) GF(p) mat lufact 0.004178 seconds (8 allocations: 79.250 KB)