Random numbers are everywhere in computing. They bring life and unpredictability to machines of otherwise relentless routine but sometimes we need to reign in their blissful chaos and exercise just a little bit more control. In this post I’ll take a look at a few ways of doing just that.

# Part 1: The Basics

In mathematics, probabilities are usually described as a proportion of one; so a 50% chance is 0.5, 100% is 1 and so on. This aligns neatly with random number generation in many languages: Math.random() in JavaScript, ActionScript and Java all produce a floating point number between 0 and 1. This can be used in a really simple way, lets say our program has to draw a grey monkey 30% of the time, the rest of the time she should be red:

```if(Math.random() < 0.3)
// draw grey monkey
else
// draw red monkey```

each time it is run our program will randomly chose to draw a grey or red monkey, overall the red monkey will be chosen 70% of the time and the grey the remaining 30%.

# Part 2: Many Outcomes

The above example is fine in simple cases but quickly becomes cumbersome when many outcomes are necessary. We could cater for red, yellow, blue, grey and green monkeys with multiple else if’s but it would be ugly, difficult to read and adjust the probabilities if needed. What we really want is a simple data structure which defines a series of outcomes and their individual probabilities:

```// array of prizes (value) with their specific probability of being awarded (p)
var outcomes = [
{ value: 'red',    p: 0.3 },
{ value: 'grey',   p: 0.2 },
{ value: 'yellow', p: 0.1 },
{ value: 'blue',   p: 0.05 },
{ value: 'green',  p: 0.35 }
];

getOutcome(outcomes);

function getOutcome(outcomes)
{
// n is the 'event' from which we'll determine an outcome
var n       = Math.random();
// p_total is the cumulation of all the potential outcomes we've considered
var p_total = 0;
var count   = 0;
var len     = outcomes.length;
while(count < len)
{
var outcome = outcomes[count];
p_total       += outcome.p;
if(n < p_total)
return outcome.value;

++count;
}
}
```

So our getOutcome function considers each outcome in turn until p_total (the cumulation of all probabilities we’ve considered so far) is greater than our random number, that signifies our hit so at that point the current value is returned. Let’s see it in action, here getOutcome is executed many times and the results outputted:

This is already pretty powerful, we have complete control over each outcome however we have to be careful that the sum of all the probabilities is 1; anything less means there’s a chance no outcome will occur, anything more could mean that certain outcomes cannot occur at all. If we want to adjust the balancing we’ll end up juggling lots of numbers less than one, it’s not very intuitive. Let’s look at how we can formulate the probabilities.

# Part 3: Quadratic Distribution

You might be familiar with the basic quadratc parabola of y = x2:

We’re going to use the area between the x-axis and the curve for different regions to represent the probability of our various outcomes.
The actual equation we’ll use is:

By changing the values of a, b and c we’ll be able to produce different types of curve and so give us probability distributions to suit many different situations. In the graph above, a=1, b=0 and c=0.

First we need to integrate the function. Fortunately (given the scope of integration) this is a simple one:

Using our integral we can find the area under the curve for any given region like so:

This gives us a total area for our outcomes. Next we divide this region into smaller portions, one for each outcome and calculate their individual areas:

Dividing each region by the total area gives its proportion; in other words the probability.

We can plug these probabilities straight back in to our getOutcome function. Now when we come to adjust the balancing we can simply experiment with different values for a, b and c until we find something suitable:

But wait!! Something isn’t quite right. Certain values of a, b and c produce some rather unusual results. Try inputting a=1, b=-3, c=-10 (lower bound=0, upper bound = 10) in the example above and you’ll see some negative probabilities; not good! The problem is that regions under the x-axis are considered to have negative area. This upsets all of our calculations; we’re not done with the maths just yet.

To calculate the true area of each region we need to separate it into subregions wherever the graph crosses the x-axis. The total area for the region is the sum of the magnitude of each of it’s parts. We need to find where the quadratic equation crosses the x-axis, i.e. where y=0, also known as finding the roots.

Not all quadratic curves cross the x-axis so we first determine if it does with the following formula:

if this produces a negative result our equation cannot be solved for zero; it doesn’t cross the x-axis. If it is positive (or zero) we can find a solution, and we do so with the following:

Most quadratic curves will have two solutions (if any), hence the two formulae above; this is a result of the ambiguity which arises when working with square roots.

Ok, that is all of the formulae we need 🙂 When it’s all put together we can produce a range of distributions and easily adjust the skewing of our probabilities by adjusting a number or two.

With the right values of a, b and c we can produce a broad range of distributions:

Near-linear distribution (a=0.1, b=0, c=-8, lower bound=9, upper bound=10)
Offset near-linear distribution (a=0.1, b=0, c=-6, lower bound=-10, upper bound=-9)
Exponential distribution (a=1, b=0, c=0, lower bound=0, upper bound=10)
Inverse exponential distribution (a=10, b=0, c=0, lower bound=-5, upper bound=0)
Approximation of a bell curve or normal distribution (a=-1, b=0, c=200, lower bound=-14, upper bound=14)

# Source

```
var assignQuadraticP = function(outcomes, a, b, c, lower_bound, upper_bound)
{
// our integrated quadratic function
var integral = function(x)
{
return 1 / 3 * a * Math.pow(x, 3) + 1 / 2 * b * Math.pow(x, 2) + c * x;
};

// this tests whether our quadraric formula crosses the x-axis within a given range
// returning the x coordinates of any intersects as an array. Note that the x-axis
// intersects have already been calculated and stored in the variable roots when this
// is invoked
var get_encapsulated_roots = function(lower_bound, upper_bound)
{
var contained = [];
var len       = roots.length;

for(var idx = 0; idx < len; ++idx)
{
var root = roots[idx];
if(root < upper_bound && root > lower_bound)
contained.push(root);
}

return contained;
};

// find the area under the quadratic curve for any given region. Any region with
// encapsulated roots is split into subregions and the magnitude of those subregions
// is summated
var find_area = function(lower_bound, upper_bound)
{
// find any roots within the region
var sub_regions = get_encapsulated_roots(lower_bound, upper_bound);
// sort encapsulated roots in ascending order
sub_regions.sort(function(a, b){return a - b;});

// if there are encapsulated roots...
if(sub_regions.length > 0)
{
sub_regions.unshift(lower_bound);
sub_regions.push(upper_bound);

var area_total = 0;
var len        = sub_regions.length;

// add up the total area from each subregion
for(var idx = 1; idx < len; ++idx)
area_total += find_area(sub_regions[idx - 1], sub_regions[idx]);

console.log('find area between: ', lower_bound, upper_bound, '=', area_total);
return area_total;
}
else
{
// if there are no encapsulated roots return the area between the boundaries
console.log('find area between: ', lower_bound, upper_bound, '=', integral(upper_bound), '-', integral(lower_bound), Math.abs(integral(upper_bound) - integral(lower_bound)));
return Math.abs(integral(upper_bound) - integral(lower_bound));
}
};

var roots    = [];
// test whether the curve crosses the x-axis
var exprss_1 = Math.pow(b, 2) - 4 * a * c;
// if it does then calculate the roots
{
var exprss_2 = 1 / 2 * a;

}

// find total area of our region
var area_total = find_area(lower_bound, upper_bound);
console.log(area_total);
var evaluated  = [];
var len        = outcomes.length;
var step_x     = (upper_bound - lower_bound) / len;

// for each possible outcome calculate its probability by finding the area of its assigned region
// as a proportion of the total area
for(var idx = 0; idx < len; ++idx)
{
evaluated.push({ value: outcomes[idx], p: find_area(lower_bound + (step_x * idx), lower_bound + (step_x * (idx + 1))) / area_total })
console.log(find_area(lower_bound + (step_x * idx), lower_bound + (step_x * (idx + 1))));
}

return evaluated;
};

var getOutcome = function(outcomes)
{
// n is the 'event' from which we'll determine an outcome
var n       = Math.random();
// p_total is the cumulation of all the potential outcomes we've considered
var p_total = 0;
var count   = 0;
var len     = outcomes.length;
while(count < len)
{
var outcome = outcomes[count];
p_total    += outcome.p;
if(n < p_total)
return outcome.value;

++count;
}
};

// define out possible outcomes, the order in which they're defined will dictate the
// order in which they're allocated a region when we start integrating
var outcomes = [ 'red', 'grey', 'yellow', 'blue', 'green' ];
var probabilities = assignQuadraticP(outcomes, a, b, c, low, high);
getOutcome(probabilities);
```

I love a bit of efficient low-power gadgetry so it was only a matter of time before I came up with a decent excuse to buy my slice of Raspberry Pi. Interest in these devices has been so strong that in the short time they’ve been around the community has already become huge; getting answers to any problems has been a doddle. With not much work I have mine running as a headless server backing up various machines with rsync over the network before it automatically publishes the backup logs to the web.

Aside from the environmental advantages of running such a low power backup device as this there is one other very good reason why only a Raspberry Pi would do for this job: for added data security it will live inside of a fire-proof safe. As good as this safe is at protecting it’s contents from the flames, it is just as good at holding on to any internal heat. Pretty much any other device would see internal temperatures creep higher and higher until problematic but our little hero remains cool and collected, producing heat no faster than it can be dissipated through the thick insulated walls of the safe.

A few final steps then to ensure all is well: I have the hard drive spin down when idle to reduce heat and wear and have the Pi log it’s core temperature at regular intervals, publishing the logs to the web as it goes. Secure, maintenance-free backup for around half the power of an energy saving light bulb.

Ah Newton; does anyone else remember back before the 1900’s when we thought we knew almost everything about the Universe? No, I’m not that old either. I’m perhaps being a bit pedantic here in stipulating Newtonian gravity but I make the distinction for two important reasons:

1. I know far too little of quantum gravity as to even dare utter it’s name
2. we don’t yet have a complete picture of what gravity is, so making a true computer model is somewhat fanciful

Anyway titles aside, this post is just to show how a little lunchtime Maths can produce something really cool.

# The Equations

This is the star of the show the force (of gravity) exerted on a body is equal to the gravitational constant multiplied by the mass of the bodies divided by the square of the distance between them. The gravitational constant by the way is just a number, in our code we can tweak this to determine the strength of gravity.
This gives us a force, but we need to know the resultant acceleration from this force, enter the next equation

That is to say force is a result of mass times acceleration.

Now we have two representations of force we can substitute one equation into the other
For body1:

when rearranged to give the acceleration the m1‘s cancel each other out; we’ve just demonstrated why acceleration due to gravity on Earth is not affected by weight:

And for body 2 the equivalent is true:

We’re almost there, we just need some bodies and some way to determine their mass. For the demo I’ll be using planet-like bodies (aka circles) of varying sizes. I’ll assume they are of constant density which will make calculating the mass easy. The volume of a sphere:

This can just be multiplied by density to give the mass; the specific value for density can be tweaked in the code:

# Gas Cloud

In this demo the bodies merge when they collide. Given a few minutes you’ll find one large body dominating the whole system. This is (vaguely) like a cloud of stellar gas; individual particles gradually coalescing before collapsing to form a star.

We can see a few of the trademarks of gravity like orbits and gravitational catapulting, pretty cool 🙂

# An Extrasolar System

The previous example is a rather chaotic, with some careful tweaking it’s possible to produce something more stable:

# Notes

The code is written for clarity rather than performance or accuracy. The collision detection for example should be performed more frequently for fast moving objects to prevent tunnelling, this manifests itself as bodies failing to merge, especially around larger objects.

# Source

```<html>
<style>
canvas { border: solid 1px #999; }
</style>
<script type="text/javascript">
document.newtGDemo1 =
{
FPS        : 30,
WIDTH      : 500,
HEIGHT     : 300,
NUM_BODIES : 40,

init: function()
{
this.cvs        = document.createElement('canvas');
this.cvs.width  = this.WIDTH;
this.cvs.height = this.HEIGHT;
this.ctx        = this.cvs.getContext('2d');

document.getElementById('newtGDemo1').appendChild(this.cvs);
},

start: function()
{
if(this.updateTmr)
clearInterval(this.updateTmr);

// this is short for body directory; an array of references to our body instances
this.bodyDir = []

// configure Body class
document.Body.prototype.CTX      = this.ctx;
document.Body.prototype.WIDTH    = this.WIDTH;
document.Body.prototype.HEIGHT   = this.HEIGHT;
document.Body.prototype.BODY_DIR = this.bodyDir;
document.Body.prototype.PREVENT_WANDER   = true;
document.Body.prototype.RESPAWN_ON_MERGE = true;

for(var count = 0; count < this.NUM_BODIES; ++count)
new document.Body();

this.lastUpdate = (new Date()).getTime();
this.updateTmr  = setInterval(this.evUpdate, 1000 / this.FPS);
},

stop: function()
{
clearInterval(this.updateTmr);
},

reset: function()
{
var len = this.bodyDir.length;
for(var idx = 0; idx < len; ++idx)
this.bodyDir[idx].destroy();
},

evUpdate: function()
{
document.newtGDemo1.update.apply(document.newtGDemo1);
},

update: function()
{
if(this.bodyDir !== document.Body.prototype.BODY_DIR)
{
clearInterval(this.updateTmr);
return;
}

// clear the canvas
this.cvs.width = this.WIDTH;
this.ctx.fillStyle = document.Body.prototype.COLOUR;

// update each body
var dt = (new Date()).getTime() - this.lastUpdate;
if(dt > 100)
dt = 100;

var count = 0;
while(count < this.bodyDir.length)
{
var body = this.bodyDir[count];
body.update(dt);
++count;
}

this.lastUpdate = (new Date()).getTime();
}
}

document.Body = function(def)
{
if (!def)
def = {}

// take definition values or default
this.density = def.density === undefined ? this.DEFAULT_DENSITY : def.density;
this.r       = def.r === undefined ? this.rndWithinRange(this.MAX_R, this.MIN_R) : def.r;
this.x       = def.x === undefined ? Math.random() * this.WIDTH : def.x;
this.y       = def.y === undefined ? Math.random() * this.HEIGHT : def.y;
this.xVel    = def.xVel === undefined ? this.rndWithinRange(this.MAX_VEL, -this.MAX_VEL) : def.xVel;
this.yVel    = def.yVel === undefined ? this.rndWithinRange(this.MAX_VEL, -this.MAX_VEL) : def.yVel;

this.mass        = this.getMass();
this.TOTAL_MASS += this.mass;

// register instance with body directory
this.BODY_DIR.push(this);
}

document
.Body.prototype =
{
DEFAULT_DENSITY  : 3,
COLOUR           : '#ffffff',
G                : 0.00005,
BODY_DIR         : [],
CTX              : null,
CIRCLE_RAD       : Math.PI * 2,
MAX_VEL          : 0.05,
MIN_R            : 0.4,
MAX_R            : 3,
TOTAL_MASS       : 0,
PREVENT_WANDER   : true,
RESPAWN_ON_MERGE : true,

destroy: function()
{
var len = this.BODY_DIR.length;
for(var idx = 0; idx < len; ++idx)
{
if(this.BODY_DIR[idx] === this)
{
this.BODY_DIR.splice(idx, 1);
break;
}
}
},

recycle: function()
{
// randomize properties and place off stage
this.r    = this.rndWithinRange(this.MAX_R, this.MIN_R);
this.x    = Math.random() < 0.5 ? Math.random() * this.WIDTH + this.WIDTH : -Math.random() * this.WIDTH;
this.y    = Math.random() < 0.5 ? Math.random() * this.HEIGHT + this.HEIGHT : -Math.random() * this.HEIGHT;
this.xVel = this.rndWithinRange(this.MAX_VEL, -this.MAX_VEL);
this.yVel = this.rndWithinRange(this.MAX_VEL, -this.MAX_VEL);

this.mass        = this.getMass();
this.TOTAL_MASS += this.mass;
},

rndWithinRange: function(max, min)
{
var range = max - min;
return min + Math.random() * range;
},

getMass: function()
{
// calculate mass based on density and radius
return this.density * Math.PI * Math.pow(this.r, 3) * 4 / 3;
},

{
// calculate radius based on mass and density
return Math.pow((3 * this.mass) / (4 * this.density * Math.PI), 1 / 3);
},

// kill off a body and merge it's properties into self
merge: function(body)
{
// combine mass
new_mass  = this.mass + body.mass;
// combine velocities, the proportion of a body's velocity used is dictated by it's mass as a proportion of combined mass
this.xVel = (this.xVel * this.mass / new_mass) + (body.xVel * body.mass / new_mass);
this.yVel = (this.yVel * this.mass / new_mass) + (body.yVel * body.mass / new_mass);

// ensure larger body dictates merged coordinates
if(body.mass > this.mass)
{
this.x = body.x;
this.y = body.y;
}

// density tends to increase with size; ensure densest body dictates merged density
if(body.density > this.density)
this.density = body.density;

// set mass and recalculate radius using merged properties
this.mass = new_mass;

if(this.RESPAWN_ON_MERGE)
body.recycle();
else
body.destroy();
},

// takes our derived equation and applies the acceleration due to gravity with the chosen body (m2)
// dt or delta time is the time elapsed since the last update
resolveVelocity: function(dt, body)
{
if(body === this)
return;

var dx    = body.x - this.x;
var dy    = body.y - this.y;
var r_sqr = Math.pow(dx, 2) + Math.pow(dy, 2);
var r     = Math.sqrt(r_sqr);

// if bodies get too close (allowing some overlap) then merge them
if(r < (this.r + body.r) / 2)
{
this.merge(body);
return;
}

// this is our derived equation
var accel = this.G * body.mass / r_sqr;

// multiplying acceleration by time gives us the change in velocity. Do this for both horizontal and vertical components
this.xVel += dt * accel * dx / r;
this.yVel += dt * accel * dy / r;
},

update: function(dt)
{
// apply gravitational acceleration caused by each individual body
for(var idx = 0; idx < this.BODY_DIR.length; ++idx)
this.resolveVelocity(dt, this.BODY_DIR[idx]);

// this just decelerates a body if it is off screen and moving further off screen, prevents a body from wandering off forever
if(this.PREVENT_WANDER)
{
if((this.x < 0 && this.xVel < 0) || (this.x > this.WIDTH && this.xVel > 0))
this.xVel *= 0.9;

if((this.y < 0 && this.yVel < 0) || (this.y > this.HEIGHT && this.yVel > 0))
this.yVel *= 0.9;
}

// move!
this.x += this.xVel * dt;
this.y += this.yVel * dt;

// draw
this.CTX.beginPath();
this.CTX.arc(this.x, this.y, this.r, 0, this.CIRCLE_RAD, false);
this.CTX.fill();
}
}
</script>