Rationals, complex numbers and maths functions

Simpson’s rule gives the exact integral of a cubic, yet run with 10 steps in Float64, it returns 0.25000000000000006 for the integral of x³ from 0 to 1, whose value is 1/4. Run with Julia’s exact fractions, the same function returns 1//4, and the rounding is gone.

Below, every result is a run recorded on Julia 1.13.1: exact arithmetic with Rational numbers and where it stops, complex numbers for the roots of a quadratic and the impedance of a circuit, functions such as log1p and hypot that keep the digits a direct formula loses, and how the constants π and ℯ turn into numbers. This lesson ends the module, with a map of where you stand in it.

Rationals are exact until their integers overflow

// builds an exact fraction, a rational numberrational number: An exact fraction of two integers, written with // in Julia: 1//3 + 1//6 == 1//2. Julia keeps it in lowest terms and checks its arithmetic for overflow. of two integers. Here are a few, a float turned into one, Simpson’s rule in rationals and in floats, and a sum that gets too long for Int64:

# // builds an exact fraction, kept in lowest terms with the sign on the numerator.
println("1//3 + 1//6 = ", 1//3 + 1//6, ", equal to 1//2: ", 1//3 + 1//6 == 1//2)
println("6//8 = ", 6//8, " (numerator ", numerator(6//8), ", denominator ", denominator(6//8), "), 5//-15 = ", 5//-15)
println("typeof(1//3) = ", typeof(1//3), ", float(1//3) = ", float(1//3))

# Converting floats: rationalize finds a short fraction, Rational gives the exact stored value.
println("rationalize(0.75) = ", rationalize(0.75), ", rationalize(0.1) = ", rationalize(0.1))
println("Rational(0.1) = ", Rational(0.1))
println("0.5 == 1//2: ", 0.5 == 1//2, ", 0.1 == 1//10: ", 0.1 == 1//10)
println("3//5 + 1 = ", 3//5 + 1, ", 3//5 - 0.5 = ", 3//5 - 0.5)

"Composite Simpson's rule with n intervals (n even), in whatever number type a and b have."
function simpson(f, a, b, n)
    h = (b - a) / n
    s = f(a) + f(b)
    for k in 1:n-1
        s += (isodd(k) ? 4 : 2) * f(a + k * h)
    end
    return s * h / 3
end

# Simpson's rule is exact for cubics: with rationals the result shows it.
cube(x) = x^3
println("∫₀¹ x³ dx, Rational: ", simpson(cube, 0//1, 1//1, 10))
println("∫₀¹ x³ dx, Float64:  ", simpson(cube, 0.0, 1.0, 10))
println("sum of the weights:  ", simpson(x -> 1, 0//1, 1//1, 10), " and ", simpson(x -> 1, 0.0, 1.0, 10))

# Rational{Int64} checks its arithmetic, so a long chain stops with an error instead of wrapping.
function harmonic(T, n)
    s = zero(Rational{T})
    for k in 1:n
        s += 1 // T(k)
    end
    return s
end
println("H(46) = ", harmonic(Int, 46))
try
    harmonic(Int, 47)
catch err
    println("H(47) → ", sprint(showerror, err))
end
println("H(47) with BigInt = ", harmonic(BigInt, 47))

println("gcd(12, 18) = ", gcd(12, 18), ", lcm(4, 6) = ", lcm(4, 6))
println("binomial(50, 25) = ", binomial(50, 25))
println("binomial(big(70), 35) = ", binomial(big(70), 35))
println("invmod(3, 7) = ", invmod(3, 7), ", powermod(2, 100, 7) = ", powermod(2, 100, 7))
binomial(70, 35)   # too large for an Int64

from session import files, run

# Exact fractions, Simpson's rule in rationals and in floats, a rational sum that overflows, and integer tools.
files("rationals.jl", "exact_steps.jl", "integer_tools.jl")
run("""
julia rationals.jl
julia exact_steps.jl
!julia integer_tools.jl
""")

$ julia rationals.jl
1//3 + 1//6 = 1//2, equal to 1//2: true
6//8 = 3//4 (numerator 3, denominator 4), 5//-15 = -1//3
typeof(1//3) = Rational{Int64}, float(1//3) = 0.3333333333333333
rationalize(0.75) = 3//4, rationalize(0.1) = 1//10
Rational(0.1) = 3602879701896397//36028797018963968
0.5 == 1//2: true, 0.1 == 1//10: false
3//5 + 1 = 8//5, 3//5 - 0.5 = 0.09999999999999998
$ julia exact_steps.jl
∫₀¹ x³ dx, Rational: 1//4
∫₀¹ x³ dx, Float64:  0.25000000000000006
sum of the weights:  1//1 and 1.0
H(46) = 5943339269060627227//1345655451257488800
H(47) → OverflowError: 5943339269060627227 * 47 overflowed for type Int64
H(47) with BigInt = 280682601097106968469//63245806209101973600
$ julia integer_tools.jl
gcd(12, 18) = 6, lcm(4, 6) = 12
binomial(50, 25) = 126410606437752
binomial(big(70), 35) = 112186277816662845432
invmod(3, 7) = 5, powermod(2, 100, 7) = 2
ERROR: LoadError: OverflowError: binomial(70, 35) overflows
Stacktrace:
 [1] binomial(n::Int64, k::Int64)
   @ Base ./intfuncs.jl:1339
 [2] top-level scope
   @ ~/demo/integer_tools.jl:5
 [3] include(mod::Module, _path::String)
   @ Base ./Base.jl:309
 [4] exec_options(opts::Base.JLOptions)
   @ Base ./client.jl:344
 [5] _start()
   @ Base ./client.jl:577
in expression starting at ~/demo/integer_tools.jl:5
[exit status 1]

1//3 + 1//6 gave exactly 1//2, and 6//8 became 3//4: the manual says rationals are reduced to lowest terms, with the sign on the numerator, so 5//-15 is -1//3. Converting a float goes two ways. rationalize(0.1) finds the short fraction 1//10 that rounds to 0.1, while Rational(0.1) gives the exact value stored in the float, a fraction whose denominator is 2⁵⁵; that is why 0.1 == 1//10 is false while 0.5 == 1//2 is true. A rational with an integer stays rational, 3//5 + 1 is 8//5, and a rational with a float becomes a float.

Simpson’s rule ran unchanged on both types, because its code only adds, multiplies and divides: with 0//1 and 1//1 as limits every step was exact, and the integral of x³ came out as 1//4. Step through one addition:

The documentation states it in one line: rationals are checked for overflow. The harmonic sum 1 + 1/2 + … + 1/46 still fitted, with a denominator of 1,345,655,451,257,488,800, but the 47th term needed a product beyond typemax(Int64), and the run stopped with an OverflowError instead of wrapping. Rational{BigInt} has no such limit.

The same integer toolbox shows up around fractions: gcd and lcm, binomial(50, 25) = 126,410,606,437,752, and modular arithmetic, where invmod(3, 7) is 5 because 3 × 5 = 15 leaves 1 when divided by 7. binomial(70, 35) does not fit an Int64 and throws, so its exact value came from binomial(big(70), 35).

Complex numbers carry a magnitude and a phase

Julia writes the imaginary unit as im, a global constant bound to i, the principal square root of −1; the manual explains that i and j were passed over because loops use them so often as names. A complex numbercomplex number: A number a + b·im with a real part a and an imaginary part b, where im² = −1. Julia's ComplexF64 has two Float64 parts; abs and angle give its size and direction. takes the type of its parts, so 2.0 + 3im is a ComplexF64. This run works through the basics, the roots of x² + 2x + 5 = 0, and then a circuit:

println(typeof(2.0 + 3im), ", ", typeof(2 + 3im), ", im^2 = ", im^2)
z = 3 + 4im
println("abs = ", abs(z), ", angle = ", angle(z), ", conj = ", conj(z), ", real = ", real(z), ", imag = ", imag(z))

# From variables, build a complex number with complex(a, b).
a, b = 1.0, 2.0
println("complex(a, b) = ", complex(a, b))

# The roots of x^2 + 2x + 5 = 0: the discriminant is negative, so take a complex square root.
d = 2^2 - 4 * 1 * 5
println("roots: ", (-2 + sqrt(complex(d))) / 2, " and ", (-2 - sqrt(complex(d))) / 2)
println("exp(im * π) = ", exp(im * π), ", cis(π / 2) = ", cis(π / 2))
println(sqrt(-1))   # a real argument gets a real square root, or an error

using Printf
# A series RC circuit: a 1 kΩ resistor and a 100 nF capacitor.
R = 1_000.0                          # ohms
C = 100e-9                           # farads
impedance(f) = R + 1 / (im * 2π * f * C)
fc = 1 / (2π * R * C)                # where the two parts of Z are equal in size

@printf("cut-off frequency fc = %.1f Hz\n", fc)
@printf("%9s %12s %10s\n", "f (Hz)", "|Z| (ohm)", "phase (°)")
for f in (10.0, 100.0, 1_000.0, fc, 10_000.0, 100_000.0)
    Z = impedance(f)
    @printf("%9.1f %12.1f %10.2f\n", f, abs(Z), rad2deg(angle(Z)))
end
Z = impedance(fc)
println("Z at fc = ", Z, ", |Z| / R = ", abs(Z) / R)

from session import files, run

# Complex arithmetic, the roots of a quadratic, Euler's formula, sqrt of a negative real, then an RC circuit's impedance.
files("complex_basics.jl", "rc_phasor.jl")
run("""
!julia complex_basics.jl
julia rc_phasor.jl
""")

$ julia complex_basics.jl
ComplexF64, Complex{Int64}, im^2 = -1 + 0im
abs = 5.0, angle = 0.9272952180016122, conj = 3 - 4im, real = 3, imag = 4
complex(a, b) = 1.0 + 2.0im
roots: -1.0 + 2.0im and -1.0 - 2.0im
exp(im * π) = -1.0 + 1.2246467991473532e-16im, cis(π / 2) = 6.123233995736766e-17 + 1.0im
ERROR: LoadError: DomainError with -1.0:
sqrt was called with a negative real argument but will only return a complex result if called with a complex argument. Try sqrt(Complex(x)).
Stacktrace:
 [1] throw_complex_domainerror(f::Symbol, x::Float64)
   @ Base.Math ./math.jl:33
 [2] sqrt
   @ ./math.jl:645 [inlined]
 [3] sqrt(x::Int64)
   @ Base.Math ./math.jl:1343
 [4] top-level scope
   @ ~/demo/complex_basics.jl:13
 [5] include(mod::Module, _path::String)
   @ Base ./Base.jl:309
 [6] exec_options(opts::Base.JLOptions)
   @ Base ./client.jl:344
 [7] _start()
   @ Base ./client.jl:577
in expression starting at ~/demo/complex_basics.jl:13
[exit status 1]
$ julia rc_phasor.jl
cut-off frequency fc = 1591.5 Hz
   f (Hz)    |Z| (ohm)  phase (°)
     10.0     159158.1     -89.64
    100.0      15946.9     -86.40
   1000.0       1879.6     -57.86
   1591.5       1414.2     -45.00
  10000.0       1012.6      -9.04
 100000.0       1000.1      -0.91
Z at fc = 1000.0 - 999.9999999999999im, |Z| / R = 1.414213562373095

abs(3 + 4im) is 5.0 and angle gives its direction in radians. complex(a, b) builds a complex number from two variables, the form the manual recommends over a + b*im. The quadratic’s discriminant is −16, and sqrt(complex(d)) gave the roots −1 ± 2im. sqrt(-1) itself stopped the first file with a DomainError whose message says what to do: call sqrt(Complex(x)). A real argument always gets a real result or an error, never a sometimes-complex one, so the type of sqrt(x) stays predictable. exp(im * π) printed −1.0 + 1.22e-16im, not exactly −1, because π became the nearest Float64.

The second program applies phasorphasor: A complex number that stands for a sine wave of one frequency: its absolute value is the amplitude and its angle the phase. Impedances in AC circuits are used the same way. analysis to a circuit. A resistor’s impedanceimpedance: The complex ratio of voltage to current in an AC circuit: R for a resistor and 1/(iωC) for a capacitor of C farads at the angular frequency ω. In series, impedances add. is R and a capacitor’s is 1/(iωC), so impedance(f) is one line, and abs and angle read off the magnitude and the phase. At 10 Hz the capacitor dominates, 159,158.1 Ω at −89.64°; at 100 kHz the resistor does, 1000.1 Ω at −0.91°; and at the cut-off frequency, 1/(2πRC) = 1591.5 Hz, the two parts are equal. Step through that row:

At the cut-off the magnitude is √2 times R and the phase −45°, and 1591.5 Hz is the cut-off frequency that a filter calculator reports for the same two parts. Rounding showed up even here: Z printed as 1000.0 − 999.9999999999999im, and |Z|/R as 1.414213562373095, one unit in the last place from √2.

Special functions keep the digits a direct formula loses

Some formulas lose digits however carefully they are typed, because of where the rounding falls: log(1 + x) rounds 1 + x before the logarithm sees it. Base has a function for each common case. This run sets the direct formula against the special function:

x = 1e-12
println("log(1 + x) = ", log(1 + x), ", log1p(x) = ", log1p(x))
println("exp(x) - 1 = ", exp(x) - 1, ", expm1(x) = ", expm1(x))
println("sqrt(a^2 + b^2) = ", sqrt(1e200^2 + 1e200^2), ", hypot(a, b) = ", hypot(1e200, 1e200))
println("sin(2π) = ", sin(2π), ", sinpi(2) = ", sinpi(2), ", cos(π/2) = ", cos(π / 2), ", cospi(0.5) = ", cospi(0.5))
println("atan(1 / -1) = ", rad2deg(atan(1 / -1)), "°, atan(1, -1) = ", rad2deg(atan(1, -1)), "°")
try
    println("(-27.0)^(1/3) = ", (-27.0)^(1 / 3))
catch err
    println("(-27.0)^(1/3) → ", sprint(showerror, err))
end
println("cbrt(-27) = ", cbrt(-27))

from session import files, run

# Six functions that keep the digits, or the range, that the direct formula loses.
files("small_x.jl")
run("""
julia small_x.jl
""")

$ julia small_x.jl
log(1 + x) = 1.000088900581841e-12, log1p(x) = 9.999999999995e-13
exp(x) - 1 = 1.000088900582341e-12, expm1(x) = 1.0000000000005e-12
sqrt(a^2 + b^2) = Inf, hypot(a, b) = 1.414213562373095e200
sin(2π) = -2.4492935982947064e-16, sinpi(2) = 0.0, cos(π/2) = 6.123233995736766e-17, cospi(0.5) = 0.0
atan(1 / -1) = -45.0°, atan(1, -1) = 135.0°
(-27.0)^(1/3) → DomainError with -27.0:
Exponentiation yielding a complex result requires a complex argument.
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.
cbrt(-27) = -3.0

At x = 1e-12, log(1 + x) printed 1.000088900581841e-12, wrong from the fifth digit, because 1 + 1e-12 had already been rounded to the float grid near 1.0; log1p(x) printed 9.999999999995e-13, correct to the last digit. expm1 does the same for exp(x) − 1. The table says which function replaces which formula, and why:

hypot(1e200, 1e200) returned 1.414213562373095e200 where the direct formula overflowed to Inf. sinpi(2) and cospi(0.5) gave exactly 0.0, while sin(2π) gave −2.4e-16 and cos(π/2) 6.1e-17, because 2π and π/2 had been rounded first. atan(1, -1) gave 135°, the true direction of the point (−1, 1), while atan(1 / -1) gave −45°, the direction of the opposite point. And the cube root written as a power, (-27.0)^(1/3), stopped with a DomainError where cbrt(-27) gave −3.0.

Irrational constants take the precision of the expression they meet

π is not a Float64. It is a value of type Irrational{:π}, an IrrationalIrrational: Julia's type for exact mathematical constants such as π and ℯ. Such a value has no fixed precision: it is rounded to Float32, Float64 or BigFloat when it meets one in an expression. that stands for the exact constant, and the documentation says such a value is rounded to the correct precision in arithmetic with other numbers. Here is π in several types, and the tests where its exactness shows:

println(π, " :: ", typeof(π))
println("Float64(π) = ", Float64(π), ", Float32(π) = ", Float32(π))
println("big(π) = ", big(π))
p128 = setprecision(BigFloat, 128) do
    big(π)
end
println("big(π) at 128 bits = ", p128)
println("precision: ", precision(big(π)), " bits by default, ", precision(p128), " inside setprecision")
println("typeof(2π) = ", typeof(2π), ", typeof(π * 1.0f0) = ", typeof(π * 1.0f0), ", typeof(π * big(1.0)) = ", typeof(π * big(1.0)))
println("π == 3.141592653589793: ", π == 3.141592653589793, ", Float64(π) == 3.141592653589793: ", Float64(π) == 3.141592653589793)
println("sin(π) = ", sin(π), ", sin(Float64(π)) = ", sin(Float64(π)))
println("ℯ = ", ℯ, ", Float64(ℯ) = ", Float64(ℯ), ", exp(1) == Float64(ℯ): ", exp(1) == Float64(ℯ))

from session import files, run

# π and ℯ as exact constants, and the type each takes in an expression.
files("irrationals.jl")
run("""
julia irrationals.jl
""")

$ julia irrationals.jl
π :: Irrational{:π}
Float64(π) = 3.141592653589793, Float32(π) = 3.1415927
big(π) = 3.141592653589793238462643383279502884197169399375105820974944592307816406286198
big(π) at 128 bits = 3.141592653589793238462643383279502884195
precision: 256 bits by default, 128 inside setprecision
typeof(2π) = Float64, typeof(π * 1.0f0) = Float32, typeof(π * big(1.0)) = BigFloat
π == 3.141592653589793: false, Float64(π) == 3.141592653589793: true
sin(π) = 0.0, sin(Float64(π)) = 1.2246467991473532e-16
ℯ = ℯ, Float64(ℯ) = 2.718281828459045, exp(1) == Float64(ℯ): true

2π is a Float64, π * 1.0f0 a Float32, and π * big(1.0) a BigFloat with 256 bits, about 77 digits; inside setprecision(BigFloat, 128) do … end the same big(π) had 128 bits. Writing Float32(π) or big(π) directly gives the best value of that type, so mixing π into single-precision code costs no precision and forces no conversion to Float64.

Because π is exact, π == 3.141592653589793 is false: Julia’s source defines an irrational as never equal to a finite number, while Float64(π) == 3.141592653589793 is true. And sin(π) printed exactly 0.0 because Base defines sin for π itself, while sin(Float64(π)) gave 1.22e-16. That special case covers π alone; sin(2π) was not 0.0, so for multiples of π use sinpi. ℯ, typed as \euler and Tab, works the same way, and exp(1) == Float64(ℯ) printed true.

Practise: exact, complex and accurate numbers

The answers to these questions are in the recorded runs. Practice starts with the idea you are least sure of and explains every choice you make:

Interview questions

Warm-up (fresher to 2–5 years): for which jobs is a Rational the better choice than a Float64 in Julia, and what does it cost?

Use a Rational when the answer must be exact and every value is a fraction of integers: probabilities of discrete events, weights of a quadrature rule, coefficients of a polynomial, ratios in a recipe or a gear train. 1//3 + 1//6 is exactly 1//2, and Simpson’s rule run on rationals gave 1//4 where floats gave 0.25000000000000006. The costs are speed, because every operation finds a gcd to keep the fraction in lowest terms; size, because numerators and denominators can grow quickly in long computations; and range, because Rational{Int64} checks its arithmetic and throws an OverflowError once they no longer fit, as the harmonic sum did at its 47th term. Switch to Rational{BigInt} for long exact chains, and convert with float at the end.

Two harder questions follow, about the impedance of a resonant circuit and about probabilities so small that the direct formulas fail. Note your answer first, then open each first step:

Key takeaways

• // builds exact rationals in lowest terms; rationalize finds a short fraction for a float and Rational its exact stored value, and Rational{Int64} throws an OverflowError rather than wrapping, so long chains need Rational{BigInt}.
• im is the imaginary unit and complex numbers take the type of their parts; abs and angle give a phasor’s magnitude and phase, and sqrt(-1) is an error while sqrt(complex(-1)) is im.
• log1p, expm1, hypot, sinpi, atan(y, x) and cbrt keep the digits, range or quadrant that the direct formulas lose.
• π and ℯ are Irrational values that take the precision of the expression they meet, are never == to a float, and get exact results only where Base defines them, such as sin(π).