Roll simultaneously a hundred 100-sided dice, and add the resulting values. The set of possible outcomes is, of course, $\{100, 101, \dotsc, 10000\}.$ Note that there is only one way to obtain the outcome 100—all dice showing a 1. There is also only one way to obtain the outcome 10,000—all dice showing 100. For the outcome 101, there are exactly 100 different ways to obtain it: all dice except one show a 1, and one die shows a 2. The follow-up question is:

Write a (fast) script that computes the different ways in which we can obtain each of the possible outcomes.

The obvious first strategy is to come up with a function of three variables $(v=\text{value}, n=\text{number of dice}, s=\text{sides per die})$ that expresses the number of different ways to obtain the value $v$ rolling $n$ $s-$sided dice. We may re-write this function iteratively in the following fashion:

 $f(v,n,s)$ $= \displaystyle{\sum_{k=1}^s} f(v-k, n-1, s) = \displaystyle{\sum_{k_1=1}^s \sum_{k_2=1}^s} f\big( v-k_1-k_2, n-2,s)$ $= \displaystyle{\sum_{k_1=1}^s \sum_{k_2=1}^s \dotsb \sum_{k_{d-1}}^s} f(v-k_1-\dotsb-k_{d-1}, 1, s\big)$

with the convention that $f(v,1,s)=0$ provided $v>s$ or $v<1,$ and of course $f(v,1,s)=1$ for $v \in \{1,\dotsc,s\}.$

A simple code (in ruby) that takes advantage of this strategy could be written as follows:

#!/usr/bin/ruby

idice=ARGV[0].to_i
isides=ARGV[1].to_i

def vwnd(value,dice,sides)
if value <  1
return 0
elsif dice==1 && value > sides
return 0
else
collect=0
for step in 1...sides+1
collect=collect+vwnd(value-step,dice-1,sides)
end
return collect
end
end

for value in idice...idice+(idice*isides-idice)/2+1
stuff = vwnd(value,idice,isides)
print value, " and ", idice*isides-value+idice, " with ", isides, "-sided ",
idice, " dice: ", stuff, "\n"
end


But this is not the best way to attack this problem: this script chokes even with relative small choices of number of dice and sides per die. A different approach to the functional relationship above offers a better chance. The focus is on the second expression:

 $f(v,n,s) =$ $\displaystyle{\sum_{k_1=1}^s \sum_{k_2=1}^s} f\big( v-k_1-k_2, n-2,s)$ $=$ $f(v-2,n-2,s) + f(v-3,n-2,s) + \dotsb + f(v-1-s,n-2,s)$ $+$ $f(v-3,n-2,s) + f(v-4,n-2,s) + \dotsb + f(v-2-s,n-2,s)$ $\vdots$ $+$ $f(v-s-1,n-2,s) + f(v-s-2,n-2,s) + \dotsb + f(v-2s,n-2,s)$

The term $f(v-2,n-2,s)$ appears exactly once; the term $f(v-3,n-2,s)$ appears exactly twice; the term $f(v-4,n-2,s)$ appears exactly three times. You start seeing now how this can work in our advantage, right? Let me make it more clear with a small example:

• 1 dice with 6 faces each: There are six possible outcomes, and each of them can be obtained in a single way.

$\begin{array}{r}\text{outcomes:}\\ \text{ways to obtain it:}\end{array} \begin{array} {|c|c|c|c|c|c|}\hline 1&2&3&4&5&6\\ \hline 1&1&1&1&1&1\\ \hline\end{array}$

• 2 dice with 6 faces each: There are eleven possible outcomes, and each of them can be obtained in different ways:

$\begin{array}{r}\text{outcomes:}\\ \text{ways to obtain it:}\end{array} \begin{array}{|c|c|c|c|c|c|c|c|c|c|c|} \hline 2&3&4&5&6&7&8&9&10&11&12\\ \hline1&2&3&4&5&6&5&4&3&2&1\\ \hline\end{array}$

A strategy to obtain this result from the previous is by collecting the previous array of "ways to obtain the outcomes"—in this case, the vector $\big(1,1,1,1,1,1\big)$—and creating the following matrix of size $(6+5) \times 6$ from it:

$\left(\begin{array}{ccccccccccc} 1&1&1&1&1&1&0&0&0&0&0\\ 0&1&1&1&1&1&1&0&0&0&0\\ 0&0&1&1&1&1&1&1&0&0&0\\ 0&0&0&1&1&1&1&1&1&0&0\\ 0&0&0&0&1&1&1&1&1&1&0\\ 0&0&0&0&0&1&1&1&1&1&1\end{array}\right)$

Add all the elements on each column, to obtain the desired result. Do you see now how to iterate this process?
• 3 dice with 6 faces each: We start with the previous result—the array $\big( 1, 2, 3, 4, 5, 6, 5, 4, 3, 2, 1\big)$—and add the elements of the columns of the matrix of size $(11+5) \times 6$ below:

$\left( \begin{array}{cccccccccccccccc} 1&2&3&4&5&6&5&4&3&2&1&0&0&0&0&0\\ 0&1&2&3&4&5&6&5&4&3&2&1&0&0&0&0\\ 0&0&1&2&3&4&5&6&5&4&3&2&1&0&0&0\\ 0&0&0&1&2&3&4&5&6&5&4&3&2&1&0&0\\ 0&0&0&0&1&2&3&4&5&6&5&4&3&2&1&0\\ 0&0&0&0&0&1&2&3&4&5&6&5&4&3&2&1\\ \end{array}\right)$

The result is, in this case, $\big( 1, 3, 6, 10, 15, 21, 25, 27, 27, 25, 21, 15, 10, 6, 3, 1\big).$

The method that we propose is then as follows:

To obtain a vector that indicates the possible ways to compute the sums of all outcomes of rolling $n$ $s-$sided dice by iteration of the above procedure. Start with $m=1$ die, until reaching $m=n$ dice. The initial vector has all ones and length $s$, and each successive $k-$th step of this algorithm will provide us with a new vector of length $s+(k-1)(s-1).$ For the last step $(k=n)$ we have, of course, $s+(n-1)(s-1) = ns-n+1$ possible outcomes.

A ruby script that uses this idea could look like this:

#!/usr/bin/ruby

dice = ARGV[0].to_i
sides= ARGV[1].to_i
original_permutation=Array.new(sides) {1}

def helper(pattern,target,step,size)
if step==1
return target
else
target+=[0]
pattern.each_index{ |x| target[x+1-step+size] += pattern[x] }
return helper(pattern,target,step-1,size)
end
end

def count_permutations(dice,sides,ary)
if dice==1
return ary
else
return puts count_permutations(dice-1,sides,helper(ary,ary,sides,sides))
end
end

puts count_permutations(dice,sides,original_permutation)


## Miscellaneous

### Generalizations

It is a great exercise to code this fun script in a language more apt to scientific computing: python, sage, matlab, C/C++, Ocaml, and of course look for tricks to simplify and speed it up. A compiled version of this algorithm should not take much time to obtain the desired result.  You can even go a little bit farther and code a new script that will compute the number of different ways to come up with the possible outcomes of rolling $n$ $s-$sided dice, where the faces of each die belong in a given set $\{F_1, F_2, \dotsc, F_s\}.$ The values $F_k$ are allowed to be equal, or non-consecutive… anything, really.

### Helper: more than a mere helper

The small recursive function that computes the key step of the algorithm is also the key to other similar counting problems. We called it helper for lack of a better word, and you can use it to aid you in other related tasks. For instance, to compute all the binomial coefficients $\binom{k}{n}$ for any $n \in \mathbb{N}$ and any $0 \leq k \leq n,$ you could issue:

first_bc = Array.new(1) {1}
def binomial_coeffs(n,ary)
if n==0
return ary
else
return binomial_coeffs(n-1,helper(ary,ary,2,2))
end
end

n=ARGV[0].to_i
puts binomial_coeffs(n,first_bc)


### Solution to our problem

The script presented above computes the required different ways to come up with the $10000-100+1=9901$ possible outcomes of rolling a hundred 100-sided dice… in about a minute! (mileage may vary depending on your computer, of course)

\$ time ./justin 100 100
1 100 5050 171700 4421275 91962520 1609344100 24370067800 325949656825 39113958
81900 42634215112710 426342151127100 3943664897925675 33976189889821200 2742363
89824985400 2084196562669889040 14980162794189827475 102217581419177646300 6644
14279224654700950 4126362365711013405900 24551856075980529765105 14029632043417
4455800600 771629762387959506903300
...
102217581419177646300 14980162794189827475 2084196562669889040 2742363898249854
00 33976189889821200 3943664897925675 426342151127100 42634215112710 3911395881
900 325949656825 24370067800 1609344100 91962520 4421275 171700 5050 100 1

real	1m16.462s
user	1m10.895s
sys	0m1.863s


## Fast scripts for computing multiple dice results

While looking for ideas to implement voronoi in sage, we stumbled upon a beautiful paper written by a group of japanese computer graphic professionals from the universities of Hokkaido and Tokyo: A Method for Creating Mosaic Images Using Voronoi Diagrams. The first step of their algorithm is simple yet brilliant: Start with any given image, and superimpose an hexagonal tiling of the plane. By a clever approximation scheme, modify the tiling to become a voronoi diagram that adaptively minimizes some approximation error. As a consequence, the resulting voronoi diagram is somehow adapted to the desired contours of the original image.
 (Fig. 1) (Fig. 2) (Fig. 3) (Fig. 4)
In a second step, they manually adjust the Voronoi image interactively by moving, adding, or deleting sites. They also take the liberty of adding visual effects by hand: emphasizing the outlines and color variations in each Voronoi region, so they look like actual pieces of stained glass (Fig. 4).
The following is a simple script in sage to initialize the mosaic. It loads the original image into the array img (Fig. 1), and creates the hexagonal voronoi image newimg (Fig. 2). Each cell on this image has a unique color, and this is obviously chosen so to minimize the $L_2$ error by constants:
from numpy import *
import scipy.ndimage
import scipy.misc
img=img[:,:,0:3]
newimg=zeros(img.shape)

# The parameter *scale* modifies the size of the initial hexagonal cells.
scale=sqrt(3.0)*4.0

# Creation of the initial hexagonal tiling.
# Centers are stored in *lattice*, cells in *voronoi*
lattice=[]
A=zeros(newimg[:,:,0].shape)
length=newimg[:,:,0].shape[0]
width=newimg[:,:,0].shape[1]

for k in range(floor(length/scale/sqrt(3)-0.5)+1):
for j in range(floor(width/scale/3-0.5)+1):
xcoord=Integer(max(0,min(floor(scale*sqrt(3)*(2*k+1)/2),length-1)))
ycoord=Integer(max(0,min(floor(3*scale*(2*j+1)/2),width-1)))
A[xcoord,ycoord]=1
lattice.append([xcoord,ycoord,1.0])

for k in range(floor(length/scale/sqrt(3)+1)):
for j in range(floor(width/scale/3+1)):
xcoord=Integer(max(0,min(floor(scale*k*sqrt(3)),length-1)))
ycoord=Integer(max(0,min(floor(3*scale*j),width-1)))
A[xcoord,ycoord]=1
lattice.append([xcoord,ycoord,1.0])
s=[[0,1,0],[1,1,1],[0,1,0]]
segmentation,segments=scipy.ndimage.label(A,s)
L1,L2=scipy.ndimage.distance_transform_edt(segmentation==0,
return_distances=False,
return_indices=True)
voronoi=segmentation[L1,L2]

# Resulting image is stored in the array *newimg*
rimg=img[:,:,0]
gimg=img[:,:,1]
bimg=img[:,:,2]

for k in range(1,voronoi.max()+1):
base=sum(voronoi==k)
rmean = sum(multiply(voronoi==k,rimg))/base
gmean = sum(multiply(voronoi==k,gimg))/base
bmean = sum(multiply(voronoi==k,bimg))/base
newimg[voronoi==k,0]=rmean
newimg[voronoi==k,1]=gmean
newimg[voronoi==k,2]=bmean

In the second step, each cell of the voronoi diagram is allowed to change its center to the eight immediate neighborhood pixels. Among the nine possible cells arising (counting the original as well), we keep the one that offers the smallest error. We repeat this second step as many times as necessary, updating the diagram each time accordingly (see Fig. 3). A naïve script to realize this task in sage could go like this—note that the algorithm is designed with pedagogical purposes in mind, instead of seeking an optimal code that employes the usual trickery to gain speed and minimize memory allocation:
for k in range(len(lattice)):
center  = lattice[k]
xcenter = center[0]
ycenter = center[1]
ecenter = center[2]
errors = zeros((3,3))
for xoffset in range(-1,2):
for yoffset in range(-1,2):
testxcenter=xcenter+xoffset
testycenter=ycenter+yoffset
if testxcenter>-1
and testxcenter < length
and testxcenter > -1
and testycenter < width:
testlattice=lattice
testlattice[k] = [xcenter+xoffset,ycenter+yoffset,ecenter]
testA=zeros(A.shape)

for location in testlattice:
testA[location[0],location[1]]=1
testsegmentation,testsegments=scipy.ndimage.label(testA,s)
L1,L2=scipy.ndimage.distance_transform_edt(
testsegmentation==0,
return_distances=False,
return_indices=True)
testvoronoi=testsegmentation[L1,L2]
newimg=zeros(img.shape)
cell = testvoronoi[xcenter,ycenter]
base=sum(testvoronoi==cell)
rmean=sum(multiply(testvoronoi==cell,rimg))/base
gmean=sum(multiply(testvoronoi==cell,gimg))/base
bmean=sum(multiply(testvoronoi==cell,bimg))/base
newimg[testvoronoi==cell,0]=rmean
newimg[testvoronoi==cell,1]=gmean
newimg[testvoronoi==cell,2]=bmean
error = 0.0
for step in range(img.shape[2]):
error+=sum(img[testvoronoi==cell,step]
-newimg[testvoronoi==cell,step])^2
errors[xoffset+1,yoffset+1] = sqrt(error)
minerror = errors.min()
minerrorposition = errors.argmin()
xshift = minerrorposition//3 -1
yshift = minerrorposition%3 -1
newxcenter = xcenter+xshift
newycenter = ycenter+yshift
if newxcenter > -1
and newxcenter < length
and newycenter > -1
and newycenter < width
and minerror < ecenter:
lattice[k]=[newxcenter,newycenter,minerror]

It is very desirable to design a third step that will actually accomplish the task of imposing the real contours of the image, and enhances them accordingly. In a future post we will present some ideas in this direction, based on two techniques: "constrained Voronoi diagrams", and "optimal subdivision of cells by linear regression."

## Geostationary Satellite Finder

Recall the First Spherical Law of Cosines:

Given a unit sphere, a spherical triangle on the surface of the sphere is defined by the great circles connecting three points $u$, $v$, and $w$ on the sphere. If the lengths of these three sides are $a$ (from $u$ to $v),$ $b$ (from $u$ to $w),$ and $c$ (from $v$ to $w),$ and the angle of the corner opposite $c$ is $C,$ then $$\cos c = \cos a \cos b + \sin a \sin b \cos C$$

In any decent device and for most computer languages, this formula should give well-conditioned results down to distances as small as around three feet, and thus can be used to compute an accurate geodetic distance between two given points in the surface of the Earth (well, ok, assuming the Earth is a perfect sphere). The geodetic form of the law of cosines is rearranged from the canonical one so that the latitude can be used directly, rather than the colatitude, and reads as follows: Given points $p_1$ and $p_2$ with positions $(lat_1, long_1)$ and $(lat_2, long_2)$ respectively, the distance $d$ between the two points is given by the following formula.

$$\cos\displaystyle{\frac{d}{R_\oplus}}=\sin(lat_1)\sin(lat_2) + \cos(lat_1)\cos(lat_2)\cos(long_2-long_1),$$

where $R_\oplus=3,959$ is the radius of the Earth in miles (well, ok, the average radius of the Earth…)

A nice application of this formula can be used for geolocation purposes, and we recently had the pleasure to assist Thumb Mobile to write such functionality for one of their clients.

Go to www.lizardsthicket.com in your mobile device, and click on "Find a Location." This fires up the location services of your browser. When you accept, your latitude and longitude are tracked. After a fast, reliable and resource-efficient algorithm, the page offers the location of the restaurant from the Lizard's chain that is closest to you. Simple, right?

## Geolocation

Tizona Scientific Solutions LLC

230 Corley Woods Drive
Lexington, SC 29072

voice: (803) 386-1822
info@tizona-sci.com

Thumb-mobile