(unsecure) elliptic numbers is julia

elliptic.jl

#elliptic.jl

#we'll make an "elliptic curve type" that is empty with no parameters.  The advantage to this is that the parameter values
#take no memory in the point representation and are JITted in as assembler immediates.
type EllipticCurve{A, B}; end

#properties of the elliptic curves
discriminant{A,B}(E::Type{EllipticCurve{A,B}}) = -16 * (4*A^3 + 27*B^2)
issmooth(E) = discriminant(E) != 0
haspoint{A,B}(E::Type{EllipticCurve{A,B}}, x, y) = y^2 == x^3 + A*x + B

#some syntactic sugar.
Base.show{A,B}(io::IO, E::Type{EllipticCurve{A,B}}) = print(io, "y^2 = x^3 + $(A)x + $(B)")

#we'll make an "elliptic point type" that is a representation of the elliptic curve with an ideal value.  There
#may actually be a better way to do this using homogeneous coordinates, but we'll forgo that for now.
type EllipticPoint{T, E <: EllipticCurve}
  x::T
  y::T
  ideal::Bool
end

#build a two-parameter constructor that does the obvious thing, but also safely checks to see if the point
#is a valid one on the curve.
function (::Type{EllipticPoint{T,E}}){T,E}(x, y)
  haspoint(E, x, y) || throw(ArgumentError("curve $E doesn't have ($x,$y)"))
  EllipticPoint{T,E}(x, y, false)
end

#create the zero (additive identity) over the abelian group.  And make a check for that, too.
Base.zero{A,B}(T::Type, E::Type{EllipticCurve{A,B}}) = EllipticPoint{T,E}(0, 1, true)
isideal(x::EllipticPoint) = x.ideal

#some syntactic sugar that makes the REPL interactions usable via cut/paste.
function Base.show{T,A,B}(io::IO, e::EllipticPoint{T,EllipticCurve{A,B}})
  write(io,isideal(e) ? "zero($T,EllipticCurve{$A,$B})" :"EllipticPoint{$T,EllipticCurve{$A,$B}}($(e.x),$(e.y))")
end

#key overloaded base functions to make math look right in julia.
Base.:(==)(e1::EllipticPoint, e2::EllipticPoint) = (e1.ideal && e2.ideal) || ((e1.x == e2.x) && (e1.x == e1.y))
function Base.:+{T,A,B}(a::EllipticPoint{T,EllipticCurve{A,B}}, b::EllipticPoint{T,EllipticCurve{A,B}})
  E = EllipticCurve{A,B}
  #make sure zero elements are additive identities
  isideal(a) && return b
  isideal(b) && return a
  #more difficult cases.
  if (a.x, a.y) == (b.x, b.y)
    #use tangent method.  Check for y = 0, which yields the ideal.
    (b.y == zero(T,E)) && return zero(T,E)
    #don't forget to grab the result from E.
    m = (3 * a.x^2 + A) / (2 * a.y)
  elseif (a.x == b.x)
    #this is vertical.
    return zero(T,E)
  else
    #use vieta's formula
    m = (b.y - a.y) / (b.x - a.x)
  end

  x = m^2 - b.x - a.x
  y = m*(x - a.x) + a.y

  return EllipticPoint{T,E}(x, -y)
end

#make subtraction work correctly
Base.:-{T,E}(e::EllipticPoint{T,E}) = isideal(e) ? e : EllipticPoint{T,E}(e.x, -e.y)
Base.:-{T,E}(a::EllipticPoint{T,E}, b::EllipticPoint{T,E}) = a + (-b)

#some syntactic sugar for the integer group action.
function Base.:*{T,E}(n::Integer, x::EllipticPoint{T,E})
  if n < 0
    -sum(x for idx = 1:-n)
  elseif n == 0
    zero(E)
  else
    sum(x for idx = 1:n)
  end
end

#examples.
"""
julia> Q = Rational{Int128}
Rational{Int128}

julia> C = EllipticCurve{-2,4}
y^2 = x^3 + -2x + 4

julia> P1 = EllipticPoint{Q,C}(3,5)
EllipticPoint{Rational{Int128},EllipticCurve{-2,4}}(3//1,5//1)

julia> P2 = EllipticPoint{Q,C}(-2,0)
EllipticPoint{Rational{Int128},EllipticCurve{-2,4}}(-2//1,0//1)

julia> P1 + P1 + P1 + P1 + P1
EllipticPoint{Rational{Int128},EllipticCurve{-2,4}}(2312883//1142761,-3507297955//1221611509)

julia> 5 * P1
EllipticPoint{Rational{Int128},EllipticCurve{-2,4}}(2312883//1142761,-3507297955//1221611509)

julia> Q = Rational{BigInt}
Rational{BigInt}

julia> P1 = EllipticPoint{Q,C}(3,5)
EllipticPoint{Rational{Int128},EllipticCurve{-2,4}}(3//1,5//1)

julia> -20*P1
EllipticPoint{Rational{BigInt},EllipticCurve{-2,4}}(8721716889552403457973789401
45384578112856996417727644408306502486841054959621893457430066791656001//5207831
20481946829397143140761792686044102902921369189488390484560995418035368116532220
330470490000,-274832909312681034314715462652601412804233448172661586199076252096
86954671299076160289194864753864983185162878307166869927581148168092234359162702
751//118846213456054547200920652321763022860552680999545167772762774106916699633
02621761108166472206145876157873100626715793555129780028801183525093000000)
"""
nothing