## Using the random.org generator in Mathematica

May 15, 2012

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 /: RandomInitializeGenerator[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.

## A short note on Costa’s minimal surface

May 14, 2012

The other day, I finally managed to simplify Alfred Gray’s parametric equations for Costa’s minimal surface. I might edit this post later, with details on how to manipulate the Weierstrass elliptic functions that show up in the equations, but enjoy these for now:

$\displaystyle x=\frac{\pi u}{2}-\Re\left(\frac{\zeta(u+iv;g_2,0)}{2}-\frac{\pi\,\wp(u+iv;g_2,0)}{\wp^\prime(u+iv;g_2,0)}\right) \\ y=\frac{\pi v}{2}+\Im\left(\frac{\zeta(u+iv;g_2,0)}{2}+\frac{\pi\,\wp(u+iv;g_2,0)}{\wp^\prime(u+iv;g_2,0)}\right) \\ z=\sqrt{\frac{\pi}{8}}\log\left|\frac1{\frac12+\frac{\wp(u+iv;g_2,0)}{\sqrt{g_2}}}-1\right|$

Here, $\wp, \wp^\prime$ and $\zeta$ are respectively the Weierstrass elliptic function, its derivative, and the Weierstrass zeta function, with invariants $g_2=\left(\frac12\mathrm{B}\left(\frac14,\frac14\right)\right)^4=\frac{\Gamma(1/4)^8}{16\pi^2}$ and $0$, and $\mathrm{B}(x,y)$ and $\Gamma(x)$ are the usual beta and gamma functions. The invariants are the ones corresponding to the semi-periods $\omega_1=\frac12$ and $\omega_3=\frac{i}{2}$. The parameter ranges are $0 < u < 1$ and $0 < v < 1$.

I was very delighted at my result, that I decided to make a stylized plot of the surface in Mathematica, using Perlin’s simplex noise for the coloring:

## FNV hashing in Mathematica

May 4, 2012

I recently needed to do Fowler/Noll/Vo hashing in Mathematica. I figured other people might be have a need for this hash function as well, so here’s my Mathematica implementation:

(* FNV hash parameters: bit length, FNV prime, FNV offset *) $FNVData = Append[#, Function[{bi, pr}, Fold[BitXor[Mod[pr #1, 2^bi], #2] &, 0, ToCharacterCode["chongo /\\../\\"]]] @@ #] & /@ Transpose[{2^Range[5, 10], 2^(8 Quotient[2^Range[5, 10] + 5, 12]) + 2^8 + Map[FromDigits[#, 16] &, {"93", "B3", "3B", "63", "57", "8D"}]}];  FNV[dat_String, bi : _Integer : 32, shift : (0 | 1) : 1] := With[{k = IntegerExponent[bi, 2] - 4}, Fold[BitXor[Mod[$FNVData[[k, 2]] #1, 2^bi], #2] &, shift \$FNVData[[k, 3]], ToCharacterCode[dat]]] /; 32 <= bi <= 1024 && BitAnd[bi, bi - 1] == 0 

I made sure to provide for both shift-1 and shift-0 versions in the implementation (though as mentioned on the FNV page, FNV-1 is what should be preferred). The FNV page has a few suggestions to extend FNV hashing to bit lengths that are not a power of two, but I did not bother with them. (Maybe somebody else can try!) It should not be too hard to modify the code to do FNV-1a hashing.

## Quick-and-dirty numerical inversion of the Laplace transform

May 4, 2012

I present here a set of Mathematica routines implementing a quick-and-dirty method for numerically evaluating the inverse Laplace transform. As it is well-known, the problem is a numerically ill-conditioned one, and thus checking the sensibility of your results is very much important here. The value of the following method, however, compared to more elaborate transform inversion methods, is that if one wants a quick estimate for arguments not too far from the imaginary axis, the approximation derived from this method does not take too much effort to evaluate.

The method, due to Herbert Salzer, is a complex-valued $n$-point Gaussian quadrature, where the nodes are the (complex!) roots of (a special case of) the Bessel polynomial. Due to this, if you want to use this method, the function you wish to transform should be able to evaluate complex arguments.

The method below does not give an error estimate; what you can try for verifying the accuracy of your results is to evaluate for an initial number of points (the default is ten points in the routine below), and then increase the number of points taken, and then compare the results you have. (You might also wish to use a sufficiently high value of the option WorkingPrecision.)

Here are the the Mathematica routines:

SalzerRuleData[n_Integer, prec_:MachinePrecision] := Block[{x}, Transpose[Map[{#, (-1)^(n + 1) # ((2 n - 1)/HypergeometricPFQ[{1 - n, n - 1}, {}, #])^2/n} &, x /. NSolve[HypergeometricPFQ[{-n, n}, {}, x], x, prec]]]]

Options[NInverseLaplaceTransform] = {Points -> 10, WorkingPrecision -> MachinePrecision}; NInverseLaplaceTransform[f_, s_, t_, opts:OptionsPattern[]] := Module[{xa, wt}, {xa, wt} = SalzerRuleData[OptionValue[Points], OptionValue[WorkingPrecision]]; Chop[(wt.Function[s, f][1/xa/t])/t]]

As an example, we try the inverse Laplace transform of the function $\frac1{(1+s)^2}$, which is $t\exp(-t)$:
 f[t_] = NInverseLaplaceTransform[(1 + s)^(-2), s, t];
 With[{t = N[2]}, {t Exp[-t], f[t], t Exp[-t] - f[t]}] {0.270671, 0.270671 + 5.68434*10^-14 I, 2.91759*10^-8 - 5.68434*10^-14 I}

(The tiny imaginary part is an artifact of the evaluation; you can use Chop[] if you wish.)

As another example, here is a comparison of the approximate inverse transform of $\frac1{\sqrt{1+s^2}}$ with the true inverse transform, $J_0(t)$ (where $J_0(t)$ is the Bessel function of the first kind):

f[t_] = NInverseLaplaceTransform[1/Sqrt[1 + s^2], s, t, WorkingPrecision -> 25];
 GraphicsRow[{Plot[{f[t], BesselJ[0, t]}, {t, 0, 15}, Frame -> True], Plot[Re[f[t] - BesselJ[0, t]], {t, 0, 15}, Frame -> True, PlotRange -> All]}]

The plot on the left shows that the true function and the approximation are almost indistinguishable within the range plotted. The plot on the right shows the difference between the approximate and the true function.

## Bunnies, teapots, and noise

March 16, 2012

I should be able to have regular Internet access soon. For now, I’d like to share some of the experiments I did with coloring 3D objects with the help of Ken Perlin’s simplex noise in Mathematica. (Relatedly, I’ve found that the implementation of improved Perlin noise supplied in the Mathematica help file has a few subtle flaws, which I was able to fix. That would be for another blog entry, though.)

Here are some of the results of coloring the Utah teapot with simplex noise:

Here are some results of coloring the Stanford bunny with simplex noise:

I never would have thought that Mathematica would ever be able to do graphics as fancy as this. I’ll only note that none of these used the Texture[] construct; both models in fact used simplex noise as solid textures (i.e. an RGBColor[] associated with each point in the model).

I might write a (series of) more descriptive post(s) on these things much later; for now, these pictures, entirely generated in Mathematica (no post-processing at all!) will have to suffice.

## On making simulated landscapes

March 1, 2012

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

February 29, 2012

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.