Using the random.org generator in Mathematica

An answer by the user Emre on the Mathematica StackExchange site introduced me to the random.org random number generation service by Mads Haahr. In the spirit of a previous post of mine, I wanted to be able to use the random.org service through the method plug-in framework for the random number generators of Mathematica. I have thus written a barebones set of routines for the purpose:

RandomOrg /: Random`InitializeGenerator[RandomOrg, ___] := RandomOrg["http://www.random.org/"]

RandomOrg[___]["GeneratesRealsQ"] := True;
RandomOrg[___]["GeneratesIntegersQ"] := True;

RandomOrg[___]["SeedGenerator"[seed_]] := RandomOrg["http://www.random.org/"]

RandomOrg[url_]["GenerateReals"[n_, {a_, b_}, prec_]] := {a + (b - a) Import[url <> "decimal-fractions/?col=1&format=plain&num=" <> ToString[n] <> "&dec=" <> ToString[Round[prec]] <> "&rnd=new", "List"], RandomOrg[url]}

RandomOrg[url_]["GenerateIntegers"[n_, {a_, b_}]] := {Import[url <> "integers/?col=1&base=10&format=plain" <> "&min=" <> ToString[a] <> "&max=" <> ToString[b] <> "&num=" <> ToString[n] <> "&rnd=new", "List"], RandomOrg[url]}

Here are a few examples of its use:

BlockRandom[SeedRandom[0, Method -> RandomOrg];
RandomReal[{0, 1}, 10, WorkingPrecision -> 20]]

BlockRandom[SeedRandom[0, Method -> RandomOrg];
RandomInteger[{1, 10}, 10]]

As a tiny reminder, the site has a daily quota in place that limits the amount of random numbers you can generate using their service in a single day. See their website for more details.

On making simulated landscapes

This entry is based on an idea by Paul Bourke for generating simulated landscapes through the application of the two-dimensional discrete Fourier transform to an array of uniformly-distributed random numbers. I was also inspired by some of the color gradients featured in GRASS, and decided to see if I could make a somewhat realistic color gradient for depicting a simulated landscape. Here is the color gradient that resulted from my many experiments:

topoFake[u_?NumericQ] :=
Blend[{{0, RGBColor[0, 1/10, 1/4]}, {1/5, RGBColor[0, 3/5, 3/5]}, {1/5, RGBColor[206/255, 197/255, 12/17]}, {1/4, RGBColor[39/85, 9/26, 5/44]}, {11/40, RGBColor[14/25, 27/40, 2/9]}, {23/80, RGBColor[9/14, 19/25, 12/85]}, {1/3, RGBColor[7/15, 29/42, 1/255]}, {3/8, RGBColor[3/10, 3/5, 5/98]}, {17/40, RGBColor[0, 38/85, 0]}, {19/40, RGBColor[0, 6/17, 0]}, {21/40, RGBColor[0, 3/10, 0]}, {23/40, RGBColor[0, 1/4, 0]}, {5/8, RGBColor[0, 1/8, 0]}, {5/8, RGBColor[39/85, 9/26, 5/44]}, {43/60, RGBColor[3/5, 25/54, 1/6]}, {97/120, RGBColor[19/23, 7/10, 1/3]}, {151/180, RGBColor[7/8, 3/4, 13/36]}, {313/360, RGBColor[39/44, 15/17, 2/5]}, {9/10, RGBColor[13/17, 9/11, 16/51]}, {9/10, RGBColor[1/2, 1/2, 1/2]}, {91/100, RGBColor[1/2, 1/2, 1/2]}, {91/100, RGBColor[1, 50/51, 50/51]}, {19/20, RGBColor[1, 50/51, 50/51]}, {19/20, RGBColor[3/4, 3/4, 3/4]}, {24/25, RGBColor[3/4, 3/4, 3/4]}, {24/25, RGBColor[1, 50/51, 50/51]}, {1, RGBColor[1, 50/51, 50/51]}}, u] /; 0

I wanted to try it out on a simpler function first, and decided to use the peaks() function of MATLAB as a test case:

Plot3D[3 (1 - x)^2 E^(-x^2 - (y + 1)^2) - 10 (x/5 - x^3 - y^5) E^(-x^2 - y^2) - 1/3 E^(-(x + 1)^2 - y^2), {x, -3, 3}, {y, -3, 3},
Axes -> None, BoundaryStyle -> None, Boxed -> False,
ColorFunction -> (topoFake[#3] &), Mesh -> False,
PerformanceGoal -> "Quality", PlotPoints -> 85, PlotRange -> All,
ViewPoint -> {-1.3, -2.4, 2.}]

Not too shabby, eh? Here’s what results if we try out the gradient on a modification of Bourke’s idea:

With[{n = 256},
BlockRandom[SeedRandom[20, Method -> "Legacy"];
ListPlot3D[Abs[InverseFourier[
Fourier[RandomReal[{0, 1}, {n, n}]] * Table[
Sech[(i/n - 0.5)/0.025] * Sech[(j/n - 0.5)/0.025],
{i, n}, {j, n}]]],
Axes -> None, Boxed -> False, BoundaryStyle -> None,
BoxRatios -> {1, 1, 1/10}, ColorFunction -> (topoFake[#3] &),
Mesh -> False, PerformanceGoal -> "Quality"]]]

The resulting archipelago looks a tad chaotic, but it’s a start. I’m still looking for a way to modify Bourke’s technique to produce bigger “islands”.

Faking marbles

Still on the matter of Perlin noise, I finally managed to figure out how to wrap Perlin noise on a sphere (that is, as a Texture[]). It was a bit inconvenient that it seems (at least, to me) that not one of the usual references on Perlin noise mentioned the explicit formulae needed for a proper spherical embedding of Perlin noise. Making use of the functions already defined in the doc page for Compile[], here are two Mathematica functions you can use for generating textures that can be wrapped on a sphere:

sphericalPerlin1 = With[{octaves = 4},
Compile[{{theta, _Real}, {phi, _Real}, {amplitude, _Real},
{frequency, _Real}, {gain, _Real}, {lacunarity, _Real}},
Module[{noiseVal = 0.0, x, y, z, sf, freq = frequency,
amp = amplitude},
sf = Sin[phi];
x = sf*Cos[theta] + 1; y = sf*Sin[theta] + 1; z = Cos[phi] + 1;
Do[
noiseVal += amp*signedNoise[x*freq, y*freq, z*freq];
freq *= lacunarity; amp *= gain,
{octaves}
];
noiseVal
], CompilationOptions -> {"InlineExternalDefinitions" -> True}
]]

sphericalPerlin2 = With[{octaves = 4},
Compile[{{theta, _Real}, {phi, _Real}, {amplitude, _Real},
{frequency, _Real}, {gain, _Real}, {lacunarity, _Real}},
Module[{noiseVal = 0.0, x, y, z, sf, freq = frequency,
amp = amplitude},
sf = Sin[phi];
x = sf*Cos[theta] + 1; y = sf*Sin[theta] + 1; z = Cos[phi] + 1;
Do[
noiseVal += amp*Abs[signedNoise[x*freq, y*freq, z*freq]];
freq *= lacunarity; amp *= gain,
{octaves}
];
noiseVal
], CompilationOptions -> {"InlineExternalDefinitions" -> True}
]]

As an example, here is how to use sphericalPerlin2[] along with DensityPlot[] to generate a texture:

DensityPlot[sphericalPerlin2[th, ph, 0.5, 2.0, 0.5, 2.0],
{th, 0, 2 Pi}, {ph, 0, Pi},
AspectRatio -> Automatic, Frame -> None, PerformanceGoal -> "Quality", PlotPoints -> 200, PlotRange -> All]

and here is the result, after adding appropriate coloring (with thanks to Rob Collyer for the color gradient):

In fact, this is the texture I used for upper leftmost “marble” in the first image.

On faking a speckled paint finish

I’ve been experimenting a lot with Perlin noise recently in Mathematica, thanks to the example implementation given in the documentation for Compile[].

(I must however note that the method implemented in the Mathematica help file is not actually the “classic” Perlin noise. In particular, the compiled fade[] routine is the Hermite quintic produced by InterpolatingPolynomial[{{{0}, 0, 0, 0}, {{1/2}, 1/2}, {{1}, 1, 0, 0}}, x], not the classical one that makes use of InterpolatingPolynomial[{{{0}, 0, 0}, {{1/2}, 1/2}, {{1}, 1, 0}}, x]. The Mathematica code also makes use of the centers of the edges of a cube as the gradients, stored in the grad variable within the function dot[], as proposed in Perlin’s update.)

The grid of swatches depicted above is just one sample of Perlin noise, colored differently with 49 of the gradient color schemes available in ColorData["Gradients"], and rendered with ReliefPlot[]. The choice of coloring scheme makes (at least to me) a heap of difference in the visual effect of Perlin noise. I like how the outputs look like samples of the fancy speckled paint I sometimes see in hardware stores.

I still haven’t fully digested the theory behind Perlin noise (though Gustavson’s notes seem to be one of the better treatments I’ve seen), but my experiments with the routines thus far have been quite encouraging.

On emulating the Texas Instruments random number generator

I recently needed to emulate the random number generator internally used by Texas Instruments calculators in Mathematica. Thanks to this United-TI forum post, I found that the algorithm being used is a combined multiplicative linear congruential generator due to Pierre L’Ecuyer. The particular generator used is listed on page 747 of the article; as noted there, the combined MLCG has a period of .

I figured this might be useful to other people, so I’m posting the routines needed to get RandomReal[] and related functions to emulate the rand() function (and ilk) on TI calculators:

TexasInstruments /: Random`InitializeGenerator[TexasInstruments, ___] := TexasInstruments[12345, 67890]

TexasInstruments[___]["GeneratesRealsQ"] := True;
TexasInstruments[___]["GeneratesIntegersQ"] := True;

TexasInstruments[___]["SeedGenerator"[seed_]] := If[seed == 0, TexasInstruments[12345, 67890],
TexasInstruments[Mod[40014 seed, 2147483563],
Mod[seed, 2147483399]]];

TexasInstruments[s1_, s2_]["GenerateReals"[n_, {a_, b_}, prec_]] :=
Module[{p = s1, q = s2, temp},
{a + (b - a) Table[p = Mod[40014 p, 2147483563]; q = Mod[40692 q, 2147483399]; temp = (p - q)/2147483563; If[temp < 0, temp += 1]; temp, {n}], TexasInstruments[p, q]}
]

TexasInstruments[s1_, s2_]["GenerateIntegers"[n_, {a_, b_}]] :=
Module[{p = s1, q = s2, temp},
{a + Floor[(b - a + 1) Table[p = Mod[40014 p, 2147483563]; q = Mod[40692 q, 2147483399]; temp = (p - q)/2147483563; If[temp < 0, temp += 1]; temp, {n}]], TexasInstruments[p, q]}
]

This is based on the PRNG method plug-in framework in Mathematica that is described here. For instance, here is the Mathematica equivalent of executing 0→rand:rand(10) (for Zilog Z80 calculators like the TI-83 Plus) or RandSeed 0:rand() (for Motorola 68k calculators like the TI-89):

(* 0→rand:rand(10) *)
BlockRandom[SeedRandom[0, Method -> TexasInstruments];
RandomReal[{0, 1}, 10, WorkingPrecision -> 20]]

Here’s corresponding code for randInt() and randM() on the Z80 calculators:

(* 0->rand:randInt(1,10,10) *)
BlockRandom[SeedRandom[0, Method -> TexasInstruments];
RandomInteger[{1, 10}, 10]]

(* 0->rand:randM(3,3) *)
BlockRandom[SeedRandom[0, Method -> TexasInstruments];
Reverse[Reverse /@ RandomInteger[{-9, 9}, {3, 3}]]]

Due to differences in precision, there is bound to be some discrepancy in the least significant figures of the results from the calculators and the results from the Mathematica emulation; still, I believe that this is a serviceable fake.