## MultiplesOfPi.jl

Numbers that produce accurate results when used as arguments to trigonometric functions
Author jishnub
Popularity
4 Stars
Updated Last
2 Years Ago
Started In
March 2020

# 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`