Real, complex, and rational numbers

Now we will consider how to handle real and complex numbers in Julia and also introduce an alternate representation of fixed-point reals as a fraction comprising two integers, the Rational datatype.

Further we will discuss the use of the Big() function to handle integers and real numbers which are too large to be represented by the primitive Julia numeric types.

Reals

We have met real numbers a few times already. The generic type is FloatingPoint which is sub-classed from Real:

abstract Real <: Number
abstract FloatingPoint <: Real
bitstype 16 Float16 <: FloatingPoint
bitstype 32 Float32 <: FloatingPoint
bitstype 64 Float64 <: FloatingPoint

A float can be defined as x = 100.0 or x = 1e2 or x = 1f2; all represent the number 100.

The first will be of the type equivalent to WORD_SIZE, the second of type Float64 and the third (using f rather than the e notation) of type Float32.

There is also a p notation which can be used with hexadecimals, that is x = 0x10p2 corresponds to 64.0.

Operators and built-in functions

Julia provides comprehensive operator and function support for real numbers. There is a wealth of mathematical functions built-in. In addition to the 'usual' ones such as exp(), log(), sin(), cos(), and so on, there is support for gamma, bessel, zeta and hankel functions and many others.

The reader should look at the section of the manual on Mathematical Operations (http://docs.julialang.org/en/release-0.3/manual/mathematical-operations/) for a comprehensive list.

One feature to note is that the multiplication operator * can be omitted in places where there is no ambiguity. If x is a variable 2.0x and 2.0*x are both valid. This is useful in cases when dealing with pre-defined constants such as pi, where 2pi is equal to 6.2831.

Special values

In dealing with real numbers Julia defines three special values Inf, -Inf and NaN. Inf and -Inf refer to values greater (or less) than all finite floating-point values and NaN is "not a number" which is a value not equal to any floating-point value (including itself). So 1.0/0.0 is Inf and -1.0/0.0 is -Inf, whereas 0.0/0.0 is NaN, as is 0.0 * Inf.

Observe that typemin(Float64) and typemax(Float64) are defined as -Inf and Inf respectively rather than the minimum/maximum representation.

BigFloats

Earlier, in regard to integers, we met BigInts.

Unsuprisingly, there are also BigFloats which can be used for arbitrary precision arithmetic:

julia>h_atoms_in_universe = 1.0*10.0^82 # => 1.0e82
julia> x = BigFloat(h_atoms_in_universe)
9.999999999999999634067965630886574211027143225273567793680363843427086501542887e+81 with 256 bits of precision

The default precision, nominally 256, and rounding mode of BigFloat can be changed using with_bigfloat_precision() and with_rounding() functions.

Notice that it is also possible to apply the big() function to achieve a similar result and the function returns either a BigInt or BigFloat, corresponding to the type of its argument:

julia> k = big(7^7); typeof(k); # => BigInt
julia> k = big(7.0^7); typeof(k); # => BigFloat

Rationals

Julia has a rational number type to represent 'exact' ratios of integers. A rational is defined by use of the // operator, for example, 5//7. If the numerator and denominator has a common factor then the number is reduced to its simplest form, 21//35 # => 3//5.

Operations on rationals or on mixed rationals and integers return a rational result:

x = 3; y = 5//7;
x*y # => 15//7
y^2 # => 25/49
y/x # => 5//21;

The functions num() and dem() return the numerator and denominator of a rational and float() can be used to convert a rational to a float.

julia> x = 17//100;
num(x) # => 17
den(x) # => 100
float(x) # => 0.17

Complex numbers

There are two ways to define a complex number in Julia. First using the type definition Complex as its associated constructor complex().

c = complex(1, 2); typeof(c) # => Complex{Float64};

Because the complex number consists of an ordered pair of two reals, its size is complex128. Similarly complex64 has Float32 arguments and complex32 has Float16 arguments.

The number Complex(0.0,1.0) corresponds to the imaginary number i, that is.sqrt(-1.0), but Julia uses the symbol im rather the i to avoid confusion with a variable i, frequently used as an index, iterator.

Hence Complex(1, 2) is exactly equivalent to 1 + 2*im, but normally the * operator is omitted and this would be expressed as 1 + 2im.

The complex number supports all the normal arithmetic operations:

julia> c = 1 + 2im; d = 3 + 4im; c*d # => -5 + 10im
julia> c/d # => 0.44 + 0.08im

Also complex functions real(), imag(), conj(), abs(), and angle().

abs and angle can be used to convert the complex arguments to polar form, that is:

julia> c = 1 + 2im; abs(c) # => 2.23606
julia> angle(c) # => 1.107148 (radians)

Complex versions of many mathematical apply:

julia> c = 1 + 2im;
julia> sin(c) = 3.1657 + 1.9596im;
julia>log(c) # => 0.8047 + 1.10715im;
Julia>sqrt(c) # => 1.272 + 0.78615im

Juliasets

The Julia documentation provides the example of generating a Mandelbrot set. Given the name of the language we will instead look at the code to create a related fractal - the Juliaset.

Both the Mandelbrot set and Julia set (for a given constant z0) are the sets of all z (complex numbers) for which the iteration z → z*z + z0 does not perge to infinity. The Mandelbrot set is those z0 for which the Julia set is connected.

We create a file jset.jl and its contents define the function to generate a Julia set:

function juliaset(z, z0, nmax::Int64)
for n = 1:nmax
        if abs(z) > 2 (return n-1) end
        z = z^2 + z0
end
return nmax
end

Here z and z0 are complex values and nmax is the number of trials to make before returning. If the modulus of the complex number z gets above 2 then it can be shown that it will increase without limit.

The function returns the number of iterations until the modulus test succeeds or else nmax.

Also we will write a second file pgmfile.jl to handle displaying the Julia set:

function create_pgmfile(img, outf::String)
  s = open(outf, "w")
   write(s, "P5\n")
   n, m = size(img)
   write(s, "$m $n 255\n")
    for  i=1:n, j=1:m
      p = img[i,j]
    write(s, uint8(p))
   end
   close(s)
end

Although we will not be looking in any depth at graphics until Chapter 7, Graphics it is quite easy to create a simple disk file using the portable bitmap (netpbm) format. This consists of a "magic" number P1 - P6, followed on the next line by the image height, width and a maximum color value which must be greater than 0 and less than 65536; all of these are ASCII values not binary.

Then follows the image values (height x width) which may be ASCII for P1, P2, P3 or binary for P4, P5, P6. There are three different types of portable bitmap; B/W (P1/P4), Grayscale (P2/P5) and Colour (P3/P6).

The function create_pgmfile() creates a binary grayscale file (magic number = P5) from an image matrix where the values are written as Uint8. Notice that the for loop defines the indices i, j in a single statement with correspondingly only one end statement. The image matrix is output in column order which matches the way it is stored in Julia.

So the main program looks like:

include("jset.jl")
include("pgmfile.jl")
h = 400; w = 800; M = Array(Int64, h, w);
c0 = -0.8+0.16im;
pgm_name = "juliaset.pgm";

t0 = time();
for  y=1:h, x=1:w
  c = complex((x-w/2)/(w/2), (y-h/2)/(w/2))
  M[y,x] = juliaset(c, c0, 256)
end
t1 = time();
create_pgmfile(M, pgm_name);
print("Written $pgm_name\nFinished in $(t1-t0) seconds.\n");

Consider how the previous code works:

1. We define an array M of type Int64 to hold the return values from the juliaset function.
2. The constant c0 is arbitrary, different values of c0 will produce different Julia sets.
3. The starting complex number is constructed from the (x,y) coordinates and scaled to the half width.
4. We 'cheat' a little by defining the maximum number of iterations as 256. This is because we are writing byte values (Uint8) and the value which remains bounded will be 256 and since overflow values wrap around will be output as 0 (black).

The script defines a main program in a function jmain():

julia>jmain
Written juliaset.pgm
Finished in 0.458 seconds # => (on my laptop)

The following screenshot shows the resulting image: