Prime factorization Elms

Some time ago I played with the idea of "plant-like" prime factorization visualization, whereas prime decomposition of a number is rendered as a plant, with branching that corresponds to its prime factors.

There was a little C++/OpenGL demo producing amusing pictures and all, however one of the obvious problems was that one had to download and compile the code to try it.

So I figured I'd give it another shot, from a different angle.

The Elm programming language looked like an interesting candidate for a few reasons:

• It's a web-frontend language, so there we go - the demo accessibility problem is sorted out
• It's high level, succinct, focused on visuals and interaction, so presumably fewer lines of code to do the same thing
• It's built upon the Functional Reactive Programming paradigm, and I heard that learning new paradigms is a good thing to do
• It's kind of like Haskell, and Haskell is hard, but these are more like nice pictures and Marios, so I could do those instead and look a little bit like the cool kids who can do Haskell!

Anyway, here's the online demo of the result (SPACE toggles pause mode, arrows allow to swap numbers back and forth):

The code itself is on GitHub, it's pastable verbatim into "try elm", so feel free to... well, try :)

(update: also on share-elm)

Now, step by step, about how the code works.

Drawing a trunk and a fruit

Let's start from drawing a "trunk" and a "fruit" on top (the numbers are hardcoded to make thing simpler):

import Color
import Graphics.Collage exposing (..)
import Graphics.Element exposing (..)
 
-- draws a trunk
trunk: Form 
trunk =
  polygon [(-5, 0), (-2, 80), (2, 80), (5, 0)]
      |> filled (Color.hsl 0.3 0.3 0.5)
 
-- draws a fruit
fruit: Form 
fruit =
  let grad =
    Color.radial
        (0, 0) 5
        (4, -4) 16
        [(0, Color.hsl 1.5 0.4 0.9),
        (1, Color.hsl 1.5 0.4 0.7)]
  in
    circle 15
    |> gradient grad
 
main: Element 
main = collage 400 300 [trunk, fruit |> moveY 80]

There are two functions "trunk" and "fruit" that create corresponding graphical elements ("Forms").

Trunk is just a polygon (trapezoid), filled with a brown color. Fruit is a circle, filled with a green-ish gradient (specified by an inner center/radius, outer center/radius and two gradient stops).

These forms are combined into an "Element" (like a DOM element) using "collage".

The fruit is moved 80 pixels upwards. Note that the center of coordinates is in the middle of the collage, and Y axis is looking upwards.

This collage is what the main program produces and what we see on the screen.

Elm borrowed the pipe function ("|>") from F#, and it's used quite a few times in the code to "thread" an argument through several modifier functions.

Drawing a tree

Now let's repeat this trunk/fruit drawing recursively, traversing a list of integers.

Each integer is interpreted as a branching factor on every level:

--  finds a rotation angle for i-th branch from n,
-- so that the branches are evenly fanned out in a cone
fanOutAngle: Float -> Int -> Int -> Float 
fanOutAngle cone i n =
  if n <= 1 then 0
  else degrees ((toFloat i)/((toFloat n) - 1) - 0.5)*cone
 
-- draws a subtree recursively,
-- according to the list of integer factors
subTree: Float -> List Int -> Form 
subTree cone s =
  case s of
    n::xs ->
      let fork = \k ->
            [ trunk
            , subTree (cone*0.95) xs
                 |> scale 0.8
                 |> moveY 80
            ]
            |> group
            |> rotate (fanOutAngle cone k n)
      in
      [0..(n - 1)]
      |> List.map fork
      |> group
    _ -> fruit
 
main: Element 
main = collage 400 300 [subTree 100 [1, 2, 3, 3, 5, 7] |> moveY -150]

That makes the recursive magic happen:

The gist of the recursion step for a "branch" is:

• draw first whatever is supposed to be on the tip of the branch (recursively)
• scale the result down a bit
• move it upwards to give some space for the parent trunk
• draw the parent trunk itself
• rotate all together with the trunk by a given angle (fanOutAngle)
• repeat for all the children (the current number n from the list)

That's it. In order to understand recursion, one needs to first understand recursion.

The fanOutAngle is a little helper function which finds how much to rotate a branch with index i, so that the total of n branches are evenly distributed inside a cone. For example, that's how it looks for fanOutAngle(cone,1,6):

Drawing a bud

We've got to tackle the situation when a number becomes too big, and drawing it as a fan of branches produces too cluttered image.

The solution is to draw it as a "bud" using a Fermat spiral layout. This both gives a nice evenly looking distribution and actually does happen in nature (the so-called phyllotaxis, or an arrangement of leaves in e.g. sunflowers), so it's a good plant metaphor.

Code-wise, the change is minimal. We introduce an alternative layout function that instead of branch-fan-out layout does the Fermat Spiral one, and switch between two based on a number threshold:

goldenPhi: Float 
goldenPhi = 137.5/180.0*pi
 
-- draws a subtree recursively,
-- according to the list of integer factors
subTree: Float -> List Int -> Form 
subTree cone s =
  case s of
    n::xs ->
      let fork = \k ->
            [ trunk
            , subTree (cone*0.95) xs
                 |> scale 0.8
                 |> moveY 80
            ]
            |> group
            |> rotate (fanOutAngle cone k n)
          bud = \k ->
            [subTree cone xs
                |> moveY (sqrt(goldenPhi*(toFloat k))*3)
                |> scale 0.4
            ]
            |> group
            |> rotate (goldenPhi*(toFloat k))
            |> moveY 30
      in
      [0..(n - 1)]
      |> List.map (if n < 17 then fork else bud)
      |> group
    _ -> fruit
 
main: Element 
main = collage 400 300 [subTree 100 [1, 2, 3, 79] |> moveY -150]

Prime number decomposition

Now let's step aside from drawing for a moment and look at the "data model" of what we draw.

The used algorithm for prime factorization is simple and rather inefficient, but it does the job just fine:

• divide by 2 until it divides (gathering the corresponding 2s as factors)
• continue from trying to divide by 3 and upwards, with step 2
• repeat until the current divisor is bigger than square root of the original number

There is also a bit of code to format the description string for the prime decomposition:

import Text
import String
import Color
import Graphics.Element exposing (..)
 
primeFactors: Int -> Int -> List Int 
primeFactors n s =
  if | s*s > n -> [n]
     | n%s == 0 -> s::(primeFactors (n//s) s)
     | s == 2 -> primeFactors n (s + 1)
     | otherwise -> primeFactors n (s + 2)
 
primes: Int -> List Int 
primes n = 1::(primeFactors n 2)
 
primeDescription: Int -> List Int -> Element 
primeDescription n factors =
    let fstr = factors
        |> List.drop 1
        |> List.map toString
        |> String.join "x"
    in
      toString n ++
         (if List.length factors == 2
          then " - PRIME!"
          else " = " ++ fstr)
        |> Text.fromString
        |> centered
 
main =
  flow down (List.map (\k -> primeDescription k (primes k)) [2..99])

which displays

2 - PRIME!
3 - PRIME!
4 = 2x2
5 - PRIME!
6 = 2x3
7 - PRIME!
8 = 2x2x2
9 = 3x3
10 = 2x5
11 - PRIME!
12 = 2x2x3
...

(by the way, one is not a prime number, but we want that one in the list anyway)

Now, we can get all of the earlier code together and display 42, decomposed into prime factors:

primeView: Int -> Element 
primeView n =
  let factors = primes n
  in
    flow down
    [ primeDescription n factors |> width  400
    , collage 400 300 [subTree 100 factors |> moveY -150]
    ]
 
main: Element 
main = primeView 42

Signals plumbing

So far we've been producing a static image as the program output. Now, getting to the actual "Functional Reactive Programming" part.

We'd like to do at least a few things:

• animate the current number, starting from 2 upwards, with a certain delay between the consecutive numbers
• allow to pause the animation with space
• allow to step the current number forward/backward with the arrow buttons
• have the graphics adjust to the current window size

This is done via "signals" in Elm - streams of data (such as time and user input) that get transformed, filtered, merged etc. And in the end, the result of the whole program is still in a sense a static image, or rather a sequence of them, built as a function of these signals.

The signals we employ:

• Time.every - fires periodically, with specified period
• Keyboard.arrows - fires whenever arrow buttons are pressed
• Keyboard.presses - general "button pressed" signal (to know when SPACE is pressed)
• Window.dimensions - when window size is changed

Essentially, our plant image gets updated correspondingly whenever any of the signals gets fired:

The code looks like:

primeView: (Int, Bool) -> (Int, Int) -> Element 
primeView (n, paused) (w, h) =
  -- ...
  -- mostly as before, only dependent on "paused" and window dimensions now
 
-- the stream of updates from timer and keyboard
type Update = Tick Float | Arrows {x: Int, y:Int} | Press Keyboard.KeyCode
updates: Signal Update 
updates =
  mergeMany
  [ map Tick (Time.every (Time.second*3))
  , map Arrows Keyboard.arrows
  , map Press Keyboard.presses
  ]
 
-- evaluate the current number and paused status based on the updates stream
foldUpdates: Update -> (Int, Bool) -> (Int, Bool) 
foldUpdates update (i, paused) =
  case update of
    Arrows {x, y} -> (max 2 ((i + x + y)%1000), paused)
    Press 32 -> (i, not paused)
    Tick _ -> ((i + (if paused then 0 else 1))%1000, paused)
    _ -> (i, paused)
 
 
main: Signal Element 
main =
  primeView
  <~ foldp foldUpdates (2, False) updates
   ~ Window.dimensions

The first three signals get merged into a single signal updates via Signal.mergeMany function. A discriminated union Update is used to identify which event in the merged signal belongs to which original signal.

Then, this merged signal is threaded through Signal.foldp, which in a sense allows to maintain a current "state", represented by tuple (i, paused), where i is the current number to display, and paused tells whether the animation is currently paused.

So after the foldp we get a signal which is a sequence of states (i, paused) to display, together with Window.dimensions signal it gets mapped into the primeView function (it has a corresponding signature):

Other tidbits

A couple of other things to mention.

One is that the we don't really want to have all these magic numbers in the code, so it's better to have them as named "constants" (which is a kind of oxymoron in our case, since in Elm everything is, you know, "constant").

In the demo code I grouped all these "tweakable constants" into a record, having in mind trying to later expose them to JavaScript via ports, so we could have an additional "tweakables" signal coming from an external JavaScript part, using e.g. dat.GUI.

Another thing is the fruit coloring: rotating hue by some constant factor depending on the branch index gives some nice visual variety:

So there we go: the beauty of mathematics combined with the beauty of functional (reactive) programming. All in your web-browser.

10/10 will try Elm again.

A test: which programming language would you be?

Growing prime factorization plants

Inspired by a wonderful  Factorization Diagrams and Dancing Factors, I've come up with a little experiment in visualization of prime factorization for integers. 

A la "math is beautiful" and all that.

The idea is that prime factors of a number are drawn as "plants". Every prime factor contributes to a branching point in the "plant", the branches have "fruits" at the tips, and the total amount of fruits is thus equal to the number itself.

For example, here's how factorization of 20 (2x2x5) and 54 (2x3x3x3) would look like:

I've hacked together a little program in C++,  it's available on github. Tested on Ubuntu Linux and Windows 7 (there is a project file for Visual Studio 2010, but you'll need to download and compile freeglut yourself in order for it to build).

The code there is as inefficient as it can get (so please don't take it as an example of how to write proper OpenGL programs), but for the sake of experiment it worked just fine.

Going a bit more into depth, at the heart of the whole thing is, of course, the integer factorization itself.

To keep things simple, it was done here with a very basic Trial Division method:

template <typename T>
void prime_factors(T val, std::vector& factors, std::vector* primes_cache = NULL)
{
    std::vector _primes;
    std::vector& primes = primes_cache ? *primes_cache : _primes;

    factors.clear();

    T max_to_check = (T)sqrt((double)val) + 1;
    find_primes(max_to_check, primes);
    const int np = primes.size();

    //  do trial division on the primes one by one
    T cval = val;
    int pidx = 0;
    while (pidx < np)
    {
        const T cp = primes[pidx];
        if (cp > max_to_check) break;
        while (cval % cp == 0) 
        {
            factors.push_back(cp);
            cval /= cp;
        }
        pidx++;
    }

    if (cval > 1) factors.push_back(cval);
}

The list of primes (find_primes) is obtained with an Eratosthenes Sieve. With the latter I went a little bit overboard (alright... it just was a kind of fun), and added a few optimizations, even though it was not needed for the particular problem at hand:

• the primes less than 1000 are precomputed
• only odd numbers get sieved out
• a bit array is used for the "sieved out" status
• sieving is done in "segments", in order to use less memory and improve cache coherence
• for the next prime cp all the numbers less than square of cp will not get checked, since they have been already sieved out by smaller primes before

So here's the function itself:

template <typename T>
void find_primes(T up_bound, std::vector& out_primes)
{
    if (out_primes.size() > 0 && out_primes.back() >= up_bound) 
    {
        //  already have the needed prime range
        return;
    }

    const T BAKED_PRIMES[] = {
        2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 
        101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 
        211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 
        307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 
        401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 
        503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 
        601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701, 
        709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 
        809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 
        907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997 };
    const int NUM_BAKED_PRIMES = sizeof(BAKED_PRIMES)/sizeof(BAKED_PRIMES[0]);

    //  search for the upper bound inside the baked range
    const T* ptop = std::lower_bound(BAKED_PRIMES, BAKED_PRIMES + NUM_BAKED_PRIMES, up_bound + 1);
        
    if (ptop < BAKED_PRIMES + NUM_BAKED_PRIMES)
    {
        //  got all the needed primes from the baked range
        out_primes = std::vector(BAKED_PRIMES, ptop);
        return;
    }

    if (out_primes.size() < NUM_BAKED_PRIMES)
    {
        //  take all the backed primes range before started sieving
        out_primes = std::vector(BAKED_PRIMES, BAKED_PRIMES + NUM_BAKED_PRIMES);
    }

    //  do the Sieve of Eratosthenes on the rest of the range, segment by segment
    const T SEGMENT_SIZE = 512;
    const T SIEVED_ARR_SIZE = (SEGMENT_SIZE + 7)/8;
    char sieved[SIEVED_ARR_SIZE]; // bit array of the "sieved" flags

    T last_sieved = out_primes.back();
    while (last_sieved < up_bound - 1)
    {
        const T num_sieved   = std::min(SEGMENT_SIZE, (up_bound - last_sieved)/2);
        const T first_sieved = last_sieved + 2;
        last_sieved = (num_sieved - 1)*2 + first_sieved;
    
        //  reset the "sieved" flags
        const T sieved_bytes = (num_sieved + 7)/8;
        std::fill(sieved, sieved + sieved_bytes, 0);

        //  sieve out the composites from the given range using the previously found primes 
        //  (only the odd numbers are considered)
        const T max_to_check = (T)sqrt((double)last_sieved) + 1;
        for (int i = 1, e = out_primes.size(); i < e; i++)
        {
            const T cp = out_primes[i];
            if (cp > max_to_check) break;
            const T fs   = std::max(first_sieved, cp*cp); // number to start sieving from
            const T fd   = ((fs + cp - 1)/cp)*cp;         // first dividend in the range
            const T fdo  = fd + cp*(1 - (fd&1));          // first odd dividend in the range
            const T di   = (fdo - first_sieved)/2;        // index of the first odd dividend
            for (T j = di; j < num_sieved; j += cp) 
            {
                //  set the corresponding "sieved" bit
                sieved[j/8] |= (1 << (j%8));
            }
        }

        //  collect all the primes (the numbers with the "sieved" bit not set)
        for (T i = 0; i < sieved_bytes; i++) 
        {
            if (sieved[i] == 0xFF) continue; // no primes in the whole byte
            for (int j = 0; j < 8; j++)
            {
                const T n = (i*8 + j);
                if (!(sieved[i]&(1 << j)) && n < num_sieved)
                {
                    out_primes.push_back(n*2 + first_sieved);
                }
            }
        }
    }
}

      

It's a good example of how a simple piece of logic can eventually turn into an incomprehensible abomination. With all the evolution, history, bizarre optimizations and personality factors.

Anyway, getting to the visualization part, here's how the numbers from 1 to 16 look like:

 

From the prime_factors function we get the list of numbers, which is already sorted from the smallest primes to the biggest ones (which naturally helps in our visualization, by avoiding too much branching early on).

Then, we are going recursively from "leaves" (that is, from the last number in the primes list towards the first one), starting with a simple branch and a "fruit" at its tip:

With every prime number P from the list we repeat the same operation: 

• all the branches we've got so far are cloned P times
• a new common root branch is added
• the parameters for the P cloned branches' are computed, such that their bases are neatly arranged along the semi-circle at the top of the newly added root branch
• each one from the P cloned branches is transformed by a composition of:
• a non-uniform scale: the scale about X axis is such that the base can fit together with other bases on the root's top semi-circle; scale about Y axis is fixed
• a shift upwards, so that the base lies on the root's top semi-circle
• rotation with the corresponding angle, so that all the P cloned branches are laid out evenly
• another shift upwards, to give space for the new root branch

It's worth noticing, that going leaves-to-root gives the tree building algorithm complexity of O(N*logN), while going in opposite direction could probably made with O(N), since we'd only need to transform every branch once. But that would only matter if we needed to optimize the algorithm, which at the moment we don't.

Because the scale that we apply on every step is non-uniform, the branches get "skewed", and this effect gets more profound towards the leaves. This may give the whole thing a rather bizarre look (which, let's say, is fine from the point of view of "aesthetics"):

In case when the prime number is too big, the "branches" are represented just with "fruits", and they are arranged into a "bud" (or "inflorescence"), which is using Fermat's spiral layout. Here's how 347 (a prime number) looks like:

The buds will be also drawn at the higher levels, albeit the non-uniform scale makes them look somewhat weird (and it's not obvious how exactly to scale the fruit radii, so I am just doing some kind of an average scale). That's what 1908 (2x2x3x3x53) turns into:

The fruits are colored randomly, using a "bipolar" green-purple color map... just because it looked funky this way (and also because I think that Kenneth Moreland's work is cool).

In the test program, it can be toggled on/off via 'c' button:

So that's about it, we've got our "factorization plants".

Turns out, prime factorization is a fascinating problem, which is important in the computing world, being, for example, at the core of widely used cryptography algorithms. It's a complex problem as well.

Getting a little bit of a visual insight into it feels like a little tiny step into understanding the problem. Even though in reality it's just some amusing pictures.

But hey, it looks nice, so who cares.