Archive for the ‘Uncategorised’ Category

Bird Flocking on GPU - CUDA

Bird Flocking Simulation

Bird flocking is an extremely interesting natural phenomenon that have been widely studied as witnessed by the number of papers in literature). I will present here a work on aggregate motion of large number of boids in a virtual environment with the presence of predators using CUDA as computational framework.

Beautiful Flocking motion

Collective motion or flocking appears at different fields and scales in nature and several mathematical tools have been developed for analyzing such motions

  1. organism are threaten as particles in Brownian motion combined with attraction/repulsion forces
  2. Differential equation models
  3. Agent based models.

The model presented here is based on a work of Reynolds (1987) which is based on three key behavioral rules:

  • Cohesion: to attempt to stay close to nearby flock-mates;
  • Collision avoidance: to evade objects that are too close;
  • Velocity/Heading Matching: to head in the same direction of nearby flock-mates

and extends it adding predator avoidance and multiple species and bird groups interaction

Bird Flocking Model

The environment is parameterized using the following a set of parameters that describe the size of the virtual environment and the duration of a timestep. What is really interesting is the set of bird's parameter that describe how a bird behave and react to event in its surrounding. Some notable parameters  include the Fyeld Of View (FOV), peak velocity v_p, thrust a and others (see figure).

Bird Parameters

Bird Parameters


The environment is partially observable and the portion of space that is visible to each bird is defined by its FOV. FOV is defined as follows:

Let p'_n= \{p^x_n-v^x_o,p^y_n-v^y_o,p^z_n-v^z_o\} the position vector of the object n in the o's frame of reference, then n is o's neighbor if and only if the followings hold:

    \[\delta_s=|| p_o - p_n ||,\; \delta_s \leq d_s\]

    \[-\; \frac{s_h}{2} \leq \theta \leq \frac{s_h}{2}, \ \ \ \ - \; \frac{s_v}{2} \leq \phi \leq \frac{s_v}{2}\]

where s_h is the maximum horizontal range of view, s_v is the maximum vertical range of view and

    \[\phi = \arccos \left(\frac{p'^z_n}{\sqrt{(p'^x_n)^2 + (p'^y_n)^2 + (p'^z_n)^2}}\right)\]

    \[\theta = atan2 \left(\frac{p'^y_n}{p'^x_n} \right) \]


In formal terms {C}_b^i, the bird b's centroid at time i, is given by:

    \[C_b^i = \frac{1}{|\mathcal{N}_b|}\sum_{n=1}^{|\mathcal{N}_b|}{\vec{p}_n \frac{d_{i,j}}{d_s}} \]



which is basically a weighted average of neighbors position.



A bird try to keep certain distance between itself and
its neighbors. Bird b's separation vector {S}_b^i at time i is given

    \[{S}_b^i =\begin{cases}\left[\sum_{j \in \mathcal{N}_b}{\frac{\vec{p}_b -\vec{p}_j}{||\vec{p}_b - \vec{p}_j||}} \;f_s\right] + a ,&\mbox{if } 0 < |\vec{S}_i|\leq v_p\\v_p, &\mbox{ otherwise }\end{cases} \]

where fs determines how strong is the repulsion against to the neighbor j.



Bird's alignment is computed as follows

    \[\vec{A}_i = \left[\sum_{j \in \mathcal{N}'_b}{\vec{v_j}}\;f_a\right] + a, \;\;0 < |\vec{A}_i| \leq v_p\]

It is a weighted average of the neighbors's heading direction.

Other species and predator avoidance
Other species avoidance is a behavior pretty much similar to the separation. The only difference is that only birds that belong to other species contribute to the result.

Predator avoidance is also a "flee or separation" behavior, but what happens here is that we do not take into account the current predator position, but instead, birds try to "separate" from predators's next position (the prediction is made from current position and velocity and acceleration of the predator).

Predator Avoidance Flee/Separation Behavior

The predator avoidance vector \vec{\Gamma}_b^i is defined as follows:

    \[\vec{\Gamma}_b^i = \left[\sum_{j \in \mathcal{P}_b }{\frac{\vec{p}_i -(\vec{p}_j+\vec{v}_j)}{||\vec{p}_i - (\vec{p}_j+\vec{v}_j)||}} \;f_{p}\right] +a, \;\;0 < |\vec{\Gamma}_i| \leq v_p \notag \]


  •  \mathcal{P}_b is the b's set of predators
  •  f_{p} = \begin{cases} 0 &\mbox{ if } d_{i,j} > r_p\\1 -\frac{d_{i,j}}{r_p} & \mbox{ otherwise}\end{cases}
    is the predator avoidance coefficient, where r_p is the minimum distance bird avoid predator.

The model has been implemented in CUDA to speedup the simulation. The following is a short video which I used during my presentation at PDP16 conference.  The model and the implementation described in this article are much more greatly described in the following slides (download here).

Largest Prime Number

Largest Prime Number

Recently a team at the University of Central Missouri, headed by Curtis Cooper has announced, via press release from the Mersenne organization

to have discovered a new largest prime. The number have more than 22M digits, to be precise. It's so large that writing 4 digits per centimeter one would be able to cover the entire distance between Edinburgh and Glasgow! Here the full number.

The following is an Haskell micro-script for computing and writing the number to a file (using Haskell because it's lazy input/output allows to write big files to disk without any concern about memory).

import Data.Time
import System.Environment (getArgs)
main :: IO ()
main = do
let n = 2^74207281 -1
startTime <- getCurrentTime
writeFile path (show n)
stopTime <-getCurrentTime
putStrLn ("ElapsedTime: " ++ show (diffUTCTime stopTime startTime))

Compile using using : ghc --make

and execute passing the output filename as the only parameter.

The full number is ~22MB, not surprisingly as one char occupies one byte in memory.

Here the first digits of the number:

268344381703937005859988258738844104703265786972872467031538046586054465054455 ....


Programming Question - Compute GCD

Programming Question: Compute the Greatest Common Divisor of two integers (easy)

This question is divided in two parts: you will first asked to write a function that computes the  GCD withouth any space/time constraints. The second part will be more challenching, asking for not using multiplication or division or addition.

Part1: write a function that takes two integers as input and returns their gcd using the famous Euclidean algorithm.

Part 2: compute the Greatest Common Divisor of two integers without using multiplication addition or division operators (easy)


Yet Another one on Fibonacci serie - Knuth multiplication

Fibonacci serie - Knuth multiplication - Java

Leonardo Fibonacci,the Italian mathematician (from Pisa, 1170, one of the most important mathematician of all time) invented the famous serie defined recursively as follows:

    \[F(n) = F(n-1) + F(n-2)\]

    \[F(0) = 0\]

    \[F(1) = 1\]

We all know that as programmers once in life (at least) we will write some piece of code to compute the n_{th} element of the serie. Recursive implementation is straightforward but inefficient, due to the exponential number of recursive calls, that make "hard" to compute numbers bigger that 40. It is very easy to come up with some iterative or memoised implementation, that lowers the complexity to O(n) and allows to efficently compute "any" number of the serie.
In this post we will derive an implementation (JAVA) that have complexity O(log_2(n))1)using results from Knuth, Fibonacci multiplication,1988. Read On…

References   [ + ]

1. using results from Knuth, Fibonacci multiplication,1988

Project Euler - Problem 24

In this article we want to briefly discuss a possible solution to the project euler problem 24 which asks:

A permutation is an ordered arrangement of objects. For example, 3124 is one possible permutation of the digits 1, 2, 3 and 4. If all of the permutations are listed numerically or alphabetically, we call it lexicographic order. The lexicographic permutations of 0, 1 and 2 are:

012 021 102 120 201 210

What is the millionth lexicographic permutation of the digits 0, 1, 2, 3, 4, 5, 6, 7, 8 and 9?

One possible approach could obviously be to brute force all the lexicographical permutations and then take the millionth element out of the previously sorted collection of permutation (interpreted as number). The drawback of this solution is that unfortunately the number of permutations is more than exponential, is O(n!).

We want to show here a more mathematical/combinatorial approach. First of all start considering that how the permutations looks like.
First permutations are (thanks to Stefano Germano for typo correction) \{0,1,2,3,4,5,6,7,8,9\}, \{0,1,2,3,4,5,6,7,9,8\},\{0,1,2,3,4,5,6,8,7,9\},...,\{0,9,8,7,6,5,4,3,2,1\}.
A deeper look at the permutations reveals that are all the sequences starting with 0 and concatenated with a permutation of the ramaining 9 elements [1,2,3,4,5,6,7,8,9]. How many of them? 9!.
It means that the permutation we are looking for  has to be greater than the 9!=362880th permutation: \{0,9,8,7,6,5,4,3,2,1\}.
From this we can deduce that we can also skip all the sequences starting by 1 because they produces (again) 9! permutations and

2 \times 9! < 10^6. Moreover our millioth permutation is in the sequences starting with 2 because 3 \times 9! > 10^6 (hence is not in the permutations starting with 3). Read On…