Author Archives: Mosè Giordano

Returned value of assignment in Julia

By: Mosè Giordano

Re-posted from: https://giordano.github.io/blog/2017-12-02-julia-assignment/

Today I realized by chance that in
the Julia programming language, like in C/C++,
assignment returns the assigned value.

To be more precise, assignment returns the right-hand side of the assignment
expression. This doesn’t appear to be documented in the current (as of December
2017) stable version of the manual, but
a
FAQ has
been added to its development version. Be aware of this subtlety, because you
can get unexpected results if you rely on the type of the returned value of an
assignment:

julia> let a
         a = 3.0
       end
3.0

julia> let a
         a::Int = 3.0
       end
3.0

In both cases the let blocks return 3.0, even though in the former case a
is really 3.0 and in the latter case it’s been marked to be an Int.

Assignment in the REPL

The fact that assignment returns the assigned value is a very handy feature.
First of all, this is probably the reason why in the Julia’s REPL the assigned
value is printed after an assignment, so that you don’t have to type again the
name of the variable to view its value:

julia> x = sin(2.58) ^ 2 + cos(2.58) ^ 2
1.0

Instead in Python, for example, assignment returns None. Thus, in the
interactive mode of the interpreter you have to type again the name of the
variable to check its value:

>>> from math import *
>>> x = sin(2.58) ** 2 + cos(2.58) ** 2
>>> x
1.0

Using assignment as condition

Another useful consequence of this feature is that you can use assignment as
condition in, e.g., if and while, making them more concise. For example,
consider this brute-force method to compute
the Euler’s number
stopping when the desired precision is reached:

euler = 0.0
i = 0
while true
    tmp = inv(factorial(i))
    euler += tmp
    tmp < eps(euler) && break
    i += 1
end

By using this feature the while can be simplified to:

euler = 0.0
i = 0
while (tmp = inv(factorial(i))) >= eps(euler)
    euler += tmp
    i += 1
end

Can this feature be a problem in Julia?

A common objection against having this feature in Python is that it is
error-prone. Quoting from
a Python tutorial:

Note that in Python, unlike C, assignment cannot occur inside expressions. C
programmers may grumble about this, but it avoids a common class of problems
encountered in C programs: typing = in an expression when == was intended.

This is indeed an issue because in C and Python conditions can be pretty much
anything that can be converted to plain-bit 0 and 1, including, e.g.,
characters. In addition, in Python a variable can freely change the type, so
that

Consider the following C program:

#include <stdio.h>

int main()
{
  char a = '\0';

  if (a = 0)
    printf("True\n");
  else
    printf("False\n");

  if (a == 0)
    printf("True\n");
  else
    printf("False\n");

  return 0;
}

which prints

False
True

Mistyping the condition (= instead of == or vice-versa) in this language can
lead to unexpected results.

Instead, this is mostly a non-issue in Julia because conditions can only be
instances of Bool type. One way to make an error in Julia is the following:

julia> a = false
false

julia> if (a = false)
           println("True")
       else
           println("False")
       end
False

julia> if a == false
           println("True")
       else
           println("False")
       end
True

but it is kind of useless to test equality with a Bool as a condition, since
you can use the Bool itself. The only other possible error you can make in
Julia I came up with is the following:

julia> a = false
false

julia> b = false
false

julia> if (a = b)
           println("True")
       else
           println("False")
       end
False

julia> if a == b
           println("True")
       else
           println("False")
       end
True

For this case I don’t have an easy replacement that could prevent you from an
error, but this is also an unlikely situation. The Julia style guide suggests
to
not parenthesize conditions.
Following this recommendation is a good way to catch the error in this case
because the assignment requires parens:

julia> if a = b
           println("True")
       else
           println("False")
       end

ERROR: syntax: unexpected "="



Computing the hexadecimal value of pi

By: Mosè Giordano

Re-posted from: https://giordano.github.io/blog/2017-11-21-hexadecimal-pi/

The
Bailey–Borwein–Plouffe formula is
one of the
several
algorithms to compute .
Here it is:

What makes this formula stand out among other approximations of is that
it allows one to directly extract the -th fractional digit of the
hexadecimal value of without computing the preceding ones.

Digits of pi

Image credit: Cormullion, Julia
code
here.

The Wikipedia article about the Bailey–Borwein–Plouffe formula explains that the
-th fractional digit (well, actually it’s the -th) is
given by

where

Only the fractional part of expression in square brackets on the right side of
is relevant, thus, in order to avoid rounding errors, when we compute
each term of the finite sum above we can take only the fractional part. This
allows us to always use ordinary double precision floating-point arithmetic,
without resorting to arbitrary-precision numbers. In addition note that the
terms of the infinite sum get quickly very small, so we can stop the summation
when they become negligible.

Implementation in Julia language

Here is a Julia implementation of the algorithm to
extract the -th fractional digit of :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
# Return the fractional part of x, modulo 1, always positive
fpart(x) = mod(x, one(x))

function Σ(n, j)
    # Compute the finite sum
    s = 0.0
    denom = j
    for k in 0:n
        s = fpart(s + powermod(16, n - k, denom) / denom)
        denom += 8
    end
    # Compute the infinite sum
    num = 1 / 16
    frac = num / denom
    while frac > eps(s)
        s     += frac
        num   /= 16
        denom += 8
        frac   = num / denom
    end
    return fpart(s)
end

pi_digit(n) =
    floor(Int, 16 * fpart(4Σ(n-1, 1) - 2Σ(n-1, 4) - Σ(n-1, 5) - Σ(n-1, 6)))

pi_string(n) = "0x3." * join(hex.(pi_digit.(1:n))) * "p0"

The pi_digit function gives the -th hexadecimal fractional digit of
as a base-10 integer, and the pi_string function returns the first
hexadecimal digits of as a valid hexadecimal floating-point
literal:

julia> pi_digit(1)
2

julia> pi_digit(6)
10

julia> pi_string(1000)
"0x3.243f6a8885a308d313198a2e03707344a4093822299f31d0082efa98ec4e6c89452821e638d01377be5466cf34e90c6cc0ac29b7c97c50dd3f84d5b5b54709179216d5d98979fb1bd1310ba698dfb5ac2ffd72dbd01adfb7b8e1afed6a267e96ba7c9045f12c7f9924a19947b3916cf70801f2e2858efc16636920d871574e69a458fea3f4933d7e0d95748f728eb658718bcd5882154aee7b54a41dc25a59b59c30d5392af26013c5d1b023286085f0ca417918b8db38ef8e79dcb0603a180e6c9e0e8bb01e8a3ed71577c1bd314b2778af2fda55605c60e65525f3aa55ab945748986263e8144055ca396a2aab10b6b4cc5c341141e8cea15486af7c72e993b3ee1411636fbc2a2ba9c55d741831f6ce5c3e169b87931eafd6ba336c24cf5c7a325381289586773b8f48986b4bb9afc4bfe81b6628219361d809ccfb21a991487cac605dec8032ef845d5de98575b1dc262302eb651b8823893e81d396acc50f6d6ff383f442392e0b4482a484200469c8f04a9e1f9b5e21c66842f6e96c9a670c9c61abd388f06a51a0d2d8542f68960fa728ab5133a36eef0b6c137a3be4ba3bf0507efb2a98a1f1651d39af017666ca593e82430e888cee8619456f9fb47d84a5c33b8b5ebee06f75d885c12073401a449f56c16aa64ed3aa62363f77061bfedf72429b023d37d0d724d00a1248db0fead3p0"

While I was preparing this post I found an unregistered
package PiBBP.jl that implements the
Bailey–Borwein–Plouffe formula. This is faster than my code above, mostly
thanks to a function
for
modular exponentiation
more efficient than that available in Julia standard library.

Test results

Let’s check if the function is working correctly. We can use
the
parse function
to convert the string to a decimal floating point
number. IEEE 754 double precision
floating-point numbers have a 53-bit mantissa, amounting to hexadecimal digits:

julia> pi_string(13)
"0x3.243f6a8885a30p0"

julia> parse(Float64, pi_string(13))
3.141592653589793

julia> Float64(π) == parse(Float64, pi_string(13))
true

Generator expressions allow
us to obtain the decimal value of the number in a very simple way, without using
the hexadecimal string:

julia> 3 + sum(pi_digit(n)/16^n for n in 1:13)
3.141592653589793

We can use
the
arbitrary-precision BigFloat
to check the correctness of the result for even more digits. By default,
BigFloat numbers in Julia have a 256-bit mantissa:

julia> precision(BigFloat)
256

The result is correct for the first hexadecimal
digits:

julia> pi_string(64)
"0x3.243f6a8885a308d313198a2e03707344a4093822299f31d0082efa98ec4e6c89p0"

julia> BigFloat(π) == parse(BigFloat, pi_string(64))
true

julia> 3 + sum(pi_digit(n)/big(16)^n for n in 1:64)
3.141592653589793238462643383279502884197169399375105820974944592307816406286198

It’s possible to increase the precision of BigFloat numbers, to further test
the accuracy of the Bailey–Borwein–Plouffe formula:

julia> setprecision(BigFloat, 4000) do
           BigFloat(π) == parse(BigFloat, pi_string(1000))
       end
true

Multi-threading

Since the Bailey–Borwein–Plouffe formula extracts the -th digit of
without computing the other ones, we can write a multi-threaded version of
pi_string, taking advantage of native support
for
multi-threading
in Julia:

1
2
3
4
5
6
7
function pi_string_threaded(N)
    digits = Vector{Int}(N)
    Threads.@threads for n in eachindex(digits)
        digits[n] = pi_digit(n)
    end
    return "0x3." * join(hex.(digits)) * "p0"
end

For example, running Julia with 4 threads gives a 2× speed-up:

julia> Threads.nthreads()
4

julia> using BenchmarkTools

julia> pi_string(1000) == pi_string_threaded(1000)
true

julia> @benchmark pi_string(1000)
BenchmarkTools.Trial:
  memory estimate:  105.28 KiB
  allocs estimate:  2016
  --------------
  minimum time:     556.228 ms (0.00% GC)
  median time:      559.198 ms (0.00% GC)
  mean time:        559.579 ms (0.00% GC)
  maximum time:     564.502 ms (0.00% GC)
  --------------
  samples:          9
  evals/sample:     1

julia> @benchmark pi_string_threaded(1000)
BenchmarkTools.Trial:
  memory estimate:  113.25 KiB
  allocs estimate:  2018
  --------------
  minimum time:     270.577 ms (0.00% GC)
  median time:      271.075 ms (0.00% GC)
  mean time:        271.598 ms (0.00% GC)
  maximum time:     278.350 ms (0.00% GC)
  --------------
  samples:          19
  evals/sample:     1



Drawing the analemma with Julia

By: Mosè Giordano

Re-posted from: https://giordano.github.io/blog/2017-11-12-analemma/

You may know that if you check the position of the Sun every day in the same
place at the same time (accounting for daylight saving time if necessary),
you’ll find that it slightly moves. This is a combination of the tilt of the
Earth’s axis and the Earth’s orbital eccentricity. The path traced out by the
position in the sky of the Sun during its wandering is
called analemma.

Analemma Murray Hill

Afternoon analemma taken in 1998–99 by Jack Fishburn in Murray Hill, New
Jersey, USA. Image
credit:
Jfishburn, Wikimedia Commons, GFDL 1.2+ and CC-BY-SA 3.0.

We can use Julia to plot the analemma. In particular,
we’ll employ AstroLib.jl to do the
needed calculations. Throughout this post I’ll assume you have installed
the latest stable version of Julia and the
necessary packages with
the
built-in package manager.

What we want to do is to determine
the position of the Sun for
a specific time every day in a year, say at noon for the whole 2018. This is
the recipe:

  1. compute the Julian dates of all
    the wanted times
  2. calculate
    the
    equatorial coordinates for
    the given Julian dates
  3. convert the equatorial coordinates
    to
    horizontal coordinates in
    the desired place. For example, we
    choose Heidelberg, in Germany,
    which has coordinates 49°25′N 08°43′E and elevation of 114 m.

The trickiest part is to get the right Julian dates.
The
jdcnv function
in AstroLib.jl assumes that times are given
in UTC standard, but
Heidelberg is one hour ahead of Greenwich. In order to work around this issue
we can use the TimeZones.zdt2julian provided by
the TimeZones.jl package which
takes care of the time zones. In addition, Germany adopts daylight saving time
from March to October, thus noon on May 15th is not actually the
same time of day as
noon on November 7th. However, noon on January 1st is the same time of day as
noon on December 31st, so we can create a range between these two times with
step one (Julian) day.

using AstroLib, TimeZones

function analemma(start_date, end_date,
                  latitude, longitude, elevation)
    julian_dates = TimeZones.zdt2julian(start_date):TimeZones.zdt2julian(end_date)
    right_ascension, declination = sunpos(julian_dates)
    altaz = eq2hor.(right_ascension, declination,
	                julian_dates, latitude, longitude, elevation)
    altitude = getindex.(altaz, 1)
    azimuth  = getindex.(altaz, 2)
    return azimuth, altitude
end

We have
used
sunpos to
get the position of the Sun in equatorial coordinates and converted them
with
eq2hor to
horizontal coordinates, specifying the coordinates of Heidelberg.
The broadcast version of this
function returns an array of 2-tuples, being the first element the altitude of
the Sun and the second element its azimuth. We’ve used getindex.(altaz, i) to
obtain the arrays with the i-th elements of the tuples. Now we can draw the
analemma. I recommend using
the Plots.jl package, which provides
a single interface to several different back-ends (GR, PyPlot, PGFPlots,
etc…).

using Plots, Base.Dates

azimuth, altitude =
    analemma(ZonedDateTime(2018,  1,  1, 12, tz"Europe/Berlin"),
             ZonedDateTime(2018, 12, 31, 12, tz"Europe/Berlin"),
             ten(49, 25), ten(8, 43), 114)

scatter(azimuth, altitude, aspect_ratio = :equal,
        xlabel = "Azimuth (°)", ylabel = "Altitude (°)")

Analemma Heidelberg

You can check with the JPL HORIZONS System
that this is accurate within a
few arcminutes.