/base/rational.jl

https://github.com/cyrilmhansen/julia · Lisp · 196 lines · 165 code · 31 blank · 0 comment · 38 complexity · 0e44c192e0e1a7c1b4ef4af6f4fd3965 MD5 · raw file

  1. immutable Rational{T<:Integer} <: Real
  2. num::T
  3. den::T
  4. function Rational(num::T, den::T)
  5. if num == 0 && den == 0
  6. error("invalid rational: 0//0")
  7. end
  8. g = gcd(den, num)
  9. new(div(num, g), div(den, g))
  10. end
  11. end
  12. Rational{T<:Integer}(n::T, d::T) = Rational{T}(n,d)
  13. Rational(n::Integer, d::Integer) = Rational(promote(n,d)...)
  14. Rational(n::Integer) = Rational(n,one(n))
  15. //(n::Integer, d::Integer ) = Rational(n,d)
  16. //(x::Rational, y::Integer ) = x.num//(x.den*y)
  17. //(x::Integer, y::Rational) = (x*y.den)//y.num
  18. //(x::Rational, y::Rational) = (x.num*y.den)//(x.den*y.num)
  19. //(x::Complex, y::Real ) = complex(real(x)//y,imag(x)//y)
  20. //(x::Real, y::Complex ) = x*y'//real(y*y')
  21. function //(x::Complex, y::Complex)
  22. xy = x*y'
  23. yy = real(y*y')
  24. complex(real(xy)//yy, imag(xy)//yy)
  25. end
  26. function show(io::IO, x::Rational)
  27. if isinf(x)
  28. print(io, x.num > 0 ? "Inf" : "-Inf")
  29. else
  30. show(io, num(x)); print(io, "//"); show(io, den(x))
  31. end
  32. end
  33. convert{T<:Integer}(::Type{Rational{T}}, x::Rational) = Rational(convert(T,x.num),convert(T,x.den))
  34. convert{T<:Integer}(::Type{Rational{T}}, x::Integer) = Rational(convert(T,x), convert(T,1))
  35. convert(::Type{Rational}, x::Rational) = x
  36. convert(::Type{Rational}, x::Integer) = convert(Rational{typeof(x)},x)
  37. convert(::Type{Bool}, x::Rational) = (x!=0) # to resolve ambiguity
  38. convert{T<:Integer}(::Type{T}, x::Rational) = (isinteger(x) ? convert(T, x.num) : throw(InexactError()))
  39. convert{T<:FloatingPoint}(::Type{T}, x::Rational) = convert(T,x.num)/convert(T,x.den)
  40. function convert{T<:Integer}(::Type{Rational{T}}, x::FloatingPoint)
  41. x == 0 && return zero(T)//one(T)
  42. if !isfinite(x)
  43. x < 0 && return -one(T)//zero(T)
  44. x > 0 && return +one(T)//zero(T)
  45. error(InexactError())
  46. end
  47. # TODO: handle values that can't be converted exactly
  48. p = get_precision(x)-1
  49. n = convert(T, significand(x)*2.0^p)
  50. p -= exponent(x)
  51. z = trailing_zeros(n)
  52. p - z > 0 ? (n>>z)//(one(T)<<(p-z)) :
  53. p > 0 ? (n>>p)//one(T) :
  54. (n<<-p)//one(T)
  55. end
  56. convert(::Type{Rational}, x::Float64) = convert(Rational{Int64}, x)
  57. convert(::Type{Rational}, x::Float32) = convert(Rational{Int}, x)
  58. promote_rule{T<:Integer}(::Type{Rational{T}}, ::Type{T}) = Rational{T}
  59. promote_rule{T<:Integer,S<:Integer}(::Type{Rational{T}}, ::Type{S}) = Rational{promote_type(T,S)}
  60. promote_rule{T<:Integer,S<:Integer}(::Type{Rational{T}}, ::Type{Rational{S}}) = Rational{promote_type(T,S)}
  61. promote_rule{T<:Integer,S<:FloatingPoint}(::Type{Rational{T}}, ::Type{S}) = promote_type(T,S)
  62. function rationalize{T<:Integer}(::Type{T}, x::FloatingPoint; tol::Real=eps(x))
  63. if isnan(x); return zero(T)//zero(T); end
  64. if x < typemin(T); return -one(T)//zero(T); end
  65. if typemax(T) < x; return one(T)//zero(T); end
  66. tm = x < 0 ? typemin(T) : typemax(T)
  67. z = x*tm
  68. if z <= 0.5 return zero(T)//one(T) end
  69. if z <= 1.0 return one(T)//tm end
  70. y = x
  71. a = d = 1
  72. b = c = 0
  73. while true
  74. f = itrunc(y); y -= f
  75. p, q = f*a+c, f*b+d
  76. typemin(T) <= p <= typemax(T) &&
  77. typemin(T) <= q <= typemax(T) || break
  78. a, b, c, d = p, q, a, b
  79. if y == 0 || abs(a/b-x) <= tol
  80. break
  81. end
  82. y = inv(y)
  83. end
  84. return convert(T,a)//convert(T,b)
  85. end
  86. rationalize(x::Union(Float64,Float32); tol::Real=eps(x)) = rationalize(Int, x, tol=tol)
  87. num(x::Integer) = x
  88. den(x::Integer) = one(x)
  89. num(x::Rational) = x.num
  90. den(x::Rational) = x.den
  91. sign(x::Rational) = sign(x.num)
  92. signbit(x::Rational) = signbit(x.num)
  93. copysign(x::Rational, y::Real) = copysign(x.num,y) // x.den
  94. copysign(x::Rational, y::Rational) = copysign(x.num,y.num) // x.den
  95. isnan(x::Rational) = false
  96. isinf(x::Rational) = x.den == 0
  97. isfinite(x::Rational) = x.den != 0
  98. typemin{T<:Integer}(::Type{Rational{T}}) = -one(T)//zero(T)
  99. typemax{T<:Integer}(::Type{Rational{T}}) = one(T)//zero(T)
  100. isinteger(x::Rational) = x.den == 1
  101. isfloat64(x::Rational) = ispow2(x.den) & (abs(x.num) <= x.den*maxintfloat(Float64))
  102. hash(x::Rational) = isinteger(x) ? hash(x.num) :
  103. isfloat64(x) ? hash(float64(x)) :
  104. bitmix(hash(x.num), hash(x.den))
  105. -(x::Rational) = (-x.num) // x.den
  106. for op in (:+, :-, :rem, :mod)
  107. @eval begin
  108. function ($op)(x::Rational, y::Rational)
  109. g = gcd(x.den, y.den)
  110. Rational(($op)(x.num * div(y.den, g), y.num * div(x.den, g)), x.den * div(y.den, g))
  111. end
  112. end
  113. end
  114. *(x::Rational, y::Rational) = (x.num*y.num) // (x.den*y.den)
  115. /(x::Rational, y::Rational) = (x.num*y.den) // (x.den*y.num)
  116. /(x::Rational, z::Complex ) = inv(z/x)
  117. ==(x::Rational, y::Rational) = (x.den == y.den) & (x.num == y.num)
  118. ==(x::Rational, y::Integer ) = (x.den == 1) & (x.num == y)
  119. ==(x::Integer , y::Rational) = y == x
  120. # needed to avoid ambiguity between ==(x::Real, z::Complex) and ==(x::Rational, y::Number)
  121. ==(z::Complex , x::Rational) = isreal(z) & (real(z) == x)
  122. ==(x::Rational, z::Complex ) = isreal(z) & (real(z) == x)
  123. ==(x::FloatingPoint, q::Rational) = ispow2(q.den) & (x == q.num/q.den) & (x*q.den == q.num)
  124. ==(q::Rational, x::FloatingPoint) = ispow2(q.den) & (x == q.num/q.den) & (x*q.den == q.num)
  125. # TODO: fix inequalities to be in line with equality check
  126. < (x::Rational, y::Rational) = x.den == y.den ? x.num < y.num :
  127. widemul(x.num,y.den) < widemul(x.den,y.num)
  128. < (x::Rational, y::Integer ) = x.num < widemul(x.den,y)
  129. < (x::Rational, y::Real ) = x.num < x.den*y
  130. < (x::Integer , y::Rational) = widemul(x,y.den) < y.num
  131. < (x::Real , y::Rational) = x*y.den < y.num
  132. <=(x::Rational, y::Rational) = x.den == y.den ? x.num <= y.num :
  133. widemul(x.num,y.den) <= widemul(x.den,y.num)
  134. <=(x::Rational, y::Integer ) = x.num <= widemul(x.den,y)
  135. <=(x::Rational, y::Real ) = x.num <= x.den*y
  136. <=(x::Integer , y::Rational) = widemul(x,y.den) <= y.num
  137. <=(x::Real , y::Rational) = x*y.den <= y.num
  138. div(x::Rational, y::Rational) = div(x.num*y.den, x.den*y.num)
  139. div(x::Rational, y::Real ) = div(x.num, x.den*y)
  140. div(x::Real , y::Rational) = div(x*y.den, y.num)
  141. fld(x::Rational, y::Rational) = fld(x.num*y.den, x.den*y.num)
  142. fld(x::Rational, y::Real ) = fld(x.num, x.den*y)
  143. fld(x::Real , y::Rational) = fld(x*y.den, y.num)
  144. itrunc(x::Rational) = div(x.num,x.den)
  145. ifloor(x::Rational) = fld(x.num,x.den)
  146. iceil (x::Rational) = -fld(-x.num,x.den)
  147. iround(x::Rational) = div(x.num*2 + copysign(x.den,x.num), x.den*2)
  148. trunc(x::Rational) = Rational(itrunc(x))
  149. floor(x::Rational) = Rational(ifloor(x))
  150. ceil (x::Rational) = Rational(iceil(x))
  151. round(x::Rational) = Rational(iround(x))
  152. ## rational to int coercion ##
  153. for f in (:int8, :int16, :int32, :int64, :int128,
  154. :uint8, :uint16, :uint32, :uint64, :uint128,
  155. :signed, :integer, :unsigned, :int, :uint)
  156. @eval ($f)(x::Rational) = ($f)(iround(x))
  157. end
  158. ^(x::Rational, y::Integer) = y < 0 ?
  159. Rational(x.den^-y, x.num^-y) : Rational(x.num^y, x.den^y)
  160. ^(x::Number, y::Rational) = x^(y.num/y.den)
  161. ^{T<:FloatingPoint}(x::T, y::Rational) = x^(convert(T,y.num)/y.den)
  162. ^{T<:FloatingPoint}(x::Complex{T}, y::Rational) = x^(convert(T,y.num)/y.den)
  163. ^{T<:Rational}(z::Complex{T}, n::Bool) = n ? z : one(z) # to resolve ambiguity
  164. ^{T<:Rational}(z::Complex{T}, n::Integer) =
  165. n>=0 ? power_by_squaring(z,n) : power_by_squaring(inv(z),-n)