MultiplesOfPi
Introduction
This package exports the type PiExpTimes{N,T}
that satisfies PiExpTimes{N,T}(x) = x*π^N
, and the type PiTimes
that is aliased to PiExpTimes{1}
. It also provides the constant Pi
for convenience, defined as PiTimes(1)
, that behaves like π
except it produces results with higher accuracy in certain trigonometric and algebraic contexts. In most scenarios the numbers Pi
and π
are interchangable.
julia> Pi^2 == π^2
true
Expressing mathematical relations in terms of Pi
instead of the more cumbersome PiExpTimes
is usually cleaner, and is recommended unless it's specifically necessary to do otherwise. One such scenario is exponentiation, more on which is presented below.
Rationale
The number π
is represented as an Irrational
type in julia, and may be computed to an arbitrary degree of precision. In normal course of events it is converted to a float when it encounters another number, for example 2π
is computed by converting both 2
and π
to floats and subsequently carrying out a floating-point multiplication. This is lossy, as both 2
and π
may be represented with arbitrary precision. This package delays the conversion of the π
to a float, treating it as a common factor in algebraic simplifications. This limits floating-point inaccuracies, especially if the terms multiplying π
are exactly representable in binary. As an added advantage, it uses sinpi
and cospi
wherever possible to avoid having to convert π
to a float altogether.
Features
Arithmetic
Delaying the conversion of π
to a float results in satisfying mathematical expressions such as
julia> (1//3)π + (4//3)π == (5//3)π
false
julia> (1//3)Pi + (4//3)Pi == (5//3)Pi
true
We may also simplify algebraic expressions involving powers of Pi
as
julia> (2Pi^2//3) // (4Pi//5)
(5//6)Pi
julia> Pi^-2 / 4Pi^3
0.25*Pi^-5
The powers of Pi
cancel as expected, and Pi^0
is automatically converted to an ordinary real number wherever possible.
julia> Pi^2 / Pi^2
1.0
Expressions involving Pi
are automatically promoted to Complex
as necessary, eg.
julia> (1+im)Pi^3 / 2Pi^2
0.5*Pi + 0.5*Pi*im
julia> (1+im)Pi^3 / 2Pi^2 * 2/Pi
1.0 + 1.0im
Trigonometric functions
The type PiTimes
uses sinpi
and cospi
under the hood when it is used as an argument to sin
and cos
. This results in exact results in several contexts where the inaccuracies arise from floating-point conversions.
julia> cos(3π/2)
-1.8369701987210297e-16
julia> cos(3Pi/2)
0.0
julia> sin(-π)
-1.2246467991473532e-16
julia> sin(-Pi)
-0.0
julia> tan(π/2)
1.633123935319537e16
julia> tan(Pi/2)
Inf
We may compute complex exponential exactly:
julia> exp(im*π/2)
6.123233995736766e-17 + 1.0im
julia> exp(im*Pi/2)
0.0 + 1.0im
# Euler's identity : exp(iπ) + 1 == 0
julia> exp(im*π) + 1
0.0 + 1.2246467991473532e-16im
julia> exp(im*Pi) + 1
0.0 + 0.0im
Hyperbolic functions work as expected:
# cosh(ix) = cos(x)
# Should be exactly zero for x = π/2
julia> cosh(im*π/2)
6.123233995736766e-17 + 0.0im
julia> cosh(im*Pi/2)
0.0 + 0.0im
Algebra
We may convert between types having different exponents without losing accuracy.
julia> convert(PiExpTimes{3},Pi)
Pi^-2*Pi^3
julia> convert(PiExpTimes{3},Pi) == Pi
true
Such an expression may be reduced to a simple form using simplify
.
julia> convert(PiExpTimes{3},Pi) |> MultiplesOfPi.simplify
Pi
Look out
Type-instability
The type PiExpTimes{N}
stores the exponent as a type-parameter, therefore exponentiation is not type-stable in general.
Constructor-abuse to avoid nesting
PiExpTimes{N}(PiExpTimes{M})
is automatically simplified to PiExpTimes{N+M}
. This is an abuse of Julia's constructors as the type of the object changes, however this avoids nested expressions that have performance issues.
julia> PiExpTimes{2}(PiExpTimes{3}(4))
4Pi^5
julia> PiExpTimes{2}(PiExpTimes{3}(4)) |> typeof
PiExpTimes{5,Int64}
Interactions with π
The irrational number π
is usually aggressively converted to PiTimes(1)
, eg:
julia> PiTimes(π)
Pi^2
This ensures that subsequent calculation would not get promoted to a floating-point type. However if this behavior is not desired then one may specify the type explicitly while constructing the object as
julia> PiTimes{Irrational{:π}}(π)
π = 3.1415926535897...*Pi
However, it is not possible to convert a number to the type PiTimes{Irrational{:π}}
.
Floating-point promotion
Addition and subtraction involving mixed exponents of Pi
will involve floating-point conversions, and the resulting expression will have the minimum exponent out of the terms being summed.
julia> Pi + 3Pi
4Pi
julia> Pi + 3Pi^2
10.42477796076938*Pi
This fits with the intuition of the expression being factorized as Pi + 3Pi^2 == Pi*(1 + 3Pi)
.
Note that π
is promoted to Pi
in such operations, so we obtain
julia> Pi + π
2Pi
Conversion vs construction
Constructors for PiExpTimes
act as a wrapper and not as a conversion. Conversion to the type retains the value of a number, whereas construction implies multiplication by an exponent of Pi
.
julia> PiTimes(1)
Pi
julia> convert(PiTimes,1)
Pi^-1*Pi
Promotion of mixed types
promote
and promote_type
work differently with types having different exponents. promote
converts both the types to one that has the minimum exponent, whereas promote_type
leaves the exponent as a free parameter.
julia> promote(Pi,Pi^2) |> typeof
Tuple{PiExpTimes{1,Real},PiExpTimes{1,Real}}
julia> promote_type(typeof(Pi),typeof(Pi^2))
PiExpTimes{N,Int64} where N
This is so that structs of PiTimes
— such as complex numbers — do not lose accuracy in conversion. The behaviour is explained below with some examples:
Arrays of mixed types
Storing different exponents of Pi
in an array will in general lead to conversion to a supertype that can store all the values.
julia> [Pi,Pi^2] # element type is not concrete
2-element Array{PiExpTimes{N,Int64} where N,1}:
Pi
Pi^2
Such an array will not be the most performant, as the element type is not concrete. This may be avoided by specifying a type while creating the array. A concrete type will not be able to store all the numbers losslessly.
julia> PiExpTimes{2,Real}[Pi,Pi^2] # exponent is fixed but Real is not a concrete type
2-element Array{PiExpTimes{2,Real},1}:
Pi^-1*Pi^2
Pi^2
julia> PiExpTimes{2,Float64}[Pi,Pi^2] # eltype is concrete, but result loses accuracy
2-element Array{PiExpTimes{2,Float64},1}:
0.3183098861837907*Pi^2
Pi^2
Complex numbers
Complex numbers rely on promote
to generate the element type
julia> Complex(Pi,Pi^2)
Pi + Pi*Pi*im
julia> Complex(Pi,Pi^2) |> typeof
Complex{PiExpTimes{1,Real}}
In this case converting either the real or imaginary part to a floating-point type would have resulted in a loss of accuracy. Such a type might not be performant, so if a conversion is desired then it might be enforced by specifying the element type while constructing the Complex
struct:
julia> Complex{PiTimes{Float64}}(Pi,Pi^2)
Pi + 3.141592653589793*Pi*im
julia> Complex{PiTimes{Float64}}(Pi,Pi^2) |> typeof
Complex{PiExpTimes{1,Float64}}
Installation
Install the package using
pkg> add MultiplesOfPi