By Arch D. Robison (Intel), Added

Julia is a new language for technical computing. It combines the convenience of dynamic languages with good performance -- about twice as slow as C. The Julia home page shows some benchmark times compared to C. The key is a dynamic type system designed to let a compiler infer most types accurately. I won't go into the type system more, but instead will show Julia from a beginner's perspective. I've provided hyperlinks to background mathematical material, though even if you don't follow the math, the example should make Julia mechanics clear.

## Horsengoggle

For the sake of learning Julia, I analyzed a round in the game of Horsengoggle, which can be technically summarized as:

*n*players are arranged in a circle, with players numbered consecutively from 0 to*n*-1.- Each player secretly picks an integer between 0 and 5 inclusive.
- All the players announce their chosen integers (by a show of fingers) simultaneously.
- The sum
*s*of all those integers is computed.. - The player numbered
*s*mod*n*is declared the winner.

My son had a class where they played multiple rounds in elimination form, to determine the order that they had to do presentations. The "winner" of the kth round became responsible for the kth presentation. After each round, the winner was removed from subsequent rounds, and the next player after the winner became the new 0 position in the circle. My son intuited the final result of the analysis below and managed to "lose" every round, thus getting to present last.

My question was whether a round is fair, or favors certain players. It's fairly easy to see that for *n*=1, *n*=2, *n*=3, or *n*=6 the game is fair since any single player's choice distributes the possible outcomes equally. But what about 10 players or 100 players?

## Analysis with Julia

Let's assume each player chooses their integer randomly. (Of course if they've done this analysis or intuited the result, they might do otherwise.) The probability distribution for one player's choices will be 1/6 for each choice 0 through 5. This distribution can be represented as a one-dimensional array containing 6 elements, each equal to 1/6. Here's the code:

ones(6)/6

That's all. The predefined function ones(*n*) creates an array of *n* ones of type **Float64** (double precision). Other precisions have to be specified explicitly, for example **ones(Int32,6)** creates an array of 6 single-precision ones. Division by a scalar implicitly maps over an array.

The probability distribution for the sum of 2 players is the convolution of the probability distributions for each player. The predefined function **conv**(*a*,*b*) computes the convolution of two sequences *a* and *b*. The probability distribution *n* players requires *n*-1 convolutions. (See end of this article for a Fourier shortcut.) Here's the Julia code for a recursive function **convPower** that does that.

# Compute a convolved n-1 times with itself function convPower(a,n) @assert n>=1 n<=1 ? a : conv(a,convPower(a,n-1)) end

The **#** introduces a comment. The **@assert** denotes an assertion. The ternary operator **?:** should be familiar to most programmers. There is no explicit return in the code because the last line of a function is implicitly the return value.

The lack of type declarations contributes to brevity and polymorphism. The Julia JIT compiler applies whole-program analysis to eliminate most of the overhead associated with dynamic typing. Julia nonetheless does allow declaring types, even parametric types, which can be useful for overloading function definitions.

To account for counting off players in a circular arrangement, the probability distribution of the sum needs to be wrapped cyclically (modulo *n*) to generate the probability distribution for who wins. Here is the wrapping function:

# Compute sum of array a wrapped modulo n function wrap(a,n) b = zeros(n) for i=1:length(a) b[mod1(i,n)] += a[i]; end b end

Local variable b is initialized with a one-dimensional array of n zeros, and then values of a are summed into it. The for loop operates over the closed range 1 through **length(a)**. The function finishes by returning b.

Brackets denote array subscripts. In Julia, as in Fortran and MATLAB, array subscripts start at 1. (The Julia creators know their audience!) Hence the function **mod1**(*i*,*n*) comes in handy. It returns a origin-one variant of the modulo function. Specifically, it returns an integer *j* congruent to *i* modulo *n*, where *j* is between 1 and *n* inclusive. The remainder operation is also available, so I could have written the subscript as (i+n-1)%n+1 instead.

Now to combine the two functions described so far. The code below uses a Julia shorthand for definining a one-line function **game(n)**:

# Compute probability distribution for one round with n players game(n) = wrap(convPower(ones(6)/6,n),n)

Though I could have typed the functions into an interactive session directly, I typed them into a file **game.jl** and loaded them into an interactive session. Here is a transcript of a session:

julia> include ("game.jl") # methods for generic function game game(n) at C:UsersDadDocumentsJuliagame-analysisgame.jl:17 julia> game(10) 10-element Array{Float64,1}: 0.0997478 0.099796 0.0999221 0.100078 0.100204 0.100252 0.100204 0.100078 0.0999221 0.099796

So yes, the game is slightly unfair, otherwise the values would all be 0.1. With more players it gets worse. With 100 players, the ratio of the maximum value to the minimum value is >36.

The increasing unfairness with increasing N is implied by the Central Limit Theorem. For large *n*, the probability distribution will be a normal distribution centered at 2.5*n* with a standard deviation that grows proportionately as sqrt(N). Hence the players opposite the starting player.become more favored compared to players near the starting player as N grows.

## For Readers Familiar with Fourier Transforms

Julia comes with many predefined functions for technical computing, including Fourier transforms (using the venerable FFTW). Here is an asymptotically faster version of my game function that employs Fourier transforms:

# Fast version of game(n) function gamef(n) @assert n>=6 real(ifft(fft([i<=6 for i=1:n]/6).^n)) end

The code transforms the problem into the Fourier domain, where circular convolution become elementwise multiplication, and repeated circular convolutions become exponentiation. The expression **[i<=6 for i=1:n]** denotes a sequence with n elements where the first 6 are 1 and the rest are 0. Predefined function **fft** takes the Fourier transform of a sequence; **ifft** is its inverse. The **.^n** raises each element of a sequence to the nth power. Most basic operations seem to map across sequences automatically, but exponentiation of complex with **^** did not, so I used the Julia dot notation to force it. Predefined function **real** takes the real parts. A somewhat subtle point is that using a transform of exactly length **n** (not necessarily a power of 2) implicitly handles the wrap-around arising from the circular arrangement of Horsengoggle players.

## Summary

Julia is an intriguing new language designed for technical computing. It has the convenience of dynamic typing with respectable speed. Its standard library comes with the tools commonly needed for technical computing. Visit the Julia web site to learn more.

[Updated 2010-Oct-10 to remove superfluous keyword **local **from declaration of **b** in function **wrap, **after I reread scope rules for Julia**.]**

## Add a Comment

Top(For technical discussions visit our developer forums. For site or software product issues contact support.)

Please sign in to add a comment. Not a member? Join today