A short note on the generalized “smoothstep” function.

August 29, 2017

The Wikipedia page on Renderman’s “smoothstep” function,


features an entry on higher-order generalizations of this function, where each succeeding member has more of its higher order derivatives being set to zero at the endpoints 0 and 1.

Among other things, a formula for these generalized functions is presented, which I find rather unwieldy. I have thus written this short entry to present a much conceptually simpler formula.

One can use e.g. confluent divided differences to derive the following (IMNSHO cleaner) formula for n-th order smoothstep \mathrm{S}_n(x):

\displaystyle \mathrm{S}_n(x)=x^{n+1}\sum\limits_{k=0}^n\binom{n+k}{k}(1-x)^k

Here is a Mathematica demonstration of this identity for the first 21 members:

And @@ Table[InterpolatingPolynomial[{PadRight[{{0}, 0}, n + 2],
                                      PadRight[{{1}, 1}, n + 2]}, x] ==
             x^(n + 1) Sum[Binomial[n + k, k] (1 - x)^k, {k, 0, n}]
             // Simplify, {n, 0, 20}]

which should yield True.

In implementations for other languages (e.g. C++, JavaScript), one would want to use a suitable modification of Horner’s rule to evaluate the polynomial factor expressed in terms of 1-x, as well as exploiting the recurrence


for generating the coefficients.


A little personal update

July 14, 2017

Some people may have been wondering about my increasingly frequent absences on the Internet. I am writing this in an attempt to assuage them.

Read the rest of this entry »

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

February 1, 2017

In this post, I will present an extension of the method previously posted in Quick-and-dirty numerical inversion of the Laplace transform.

The method, due to Piessens, is effectively a Gauss-Kronrod extension of the complex Gaussian quadrature method by Salzer. Recall that the Gauss-Kronrod idea is to use two different quadrature abscissa/weight pairs, where one’s abscissas is a subset of the other, to be able to estimate the accuracy of the higher-order method.

Here is the implementation of the extended version of NInverseLaplaceTransform[]:

SalzerRuleData[n_Integer, prec_:MachinePrecision] :=
      Transpose[Map[{#, (-1)^(n + 1) # ((2 n - 1)/
                        HypergeometricPFQ[{1 - n, n - 1}, {}, #])^2/n} &,
                    \[FormalX] /.
                    NSolve[HypergeometricPFQ[{-n, n}, {}, \[FormalX]],
                           \[FormalX], prec]]]
piessensQ[n_Integer, x_] := Module[{cof},
        cof = LinearSolve[SparseArray[{{j_, k_} /; j <= k :> 1/((k - j)!
                                      Binomial[2 n + k - j, n + k - j])},
                                      {n + 1, n + 1}],
                          -Table[1/(j! Binomial[2 n + j, n + j]),
                                 {j, n + 1, 1, -1}]];
        Expand[FromDigits[Prepend[Reverse[cof], 1], x]]]
SalzerPiessensRuleData[n_Integer, prec_: MachinePrecision] :=
      Module[{c, pk, pl, sg, sk, sl},
             sl = \[FormalX] /.
             NSolve[HypergeometricPFQ[{-n, n}, {}, \[FormalX]], \[FormalX], prec];
             pl = \[FormalX] /.
             NSolve[piessensQ[n, \[FormalX]], \[FormalX], prec];
             sk = HypergeometricPFQ[{1 - n, n - 1}, {}, sl];
             sg = (-1)^(n + 1) sl ((2 n - 1)/sk)^2/n; c = n! Binomial[2 n, n];
             sk = sg - (2 n - 1) sl/(c n sk piessensQ[n, sl]);
             pk = 1/(c pl HypergeometricPFQ[{-n, n}, {}, pl]
                     (D[piessensQ[n, \[FormalX]], \[FormalX]] /.
                      \[FormalX] -> pl));
             {Join[sl, pl], Join[sk, pk]}]
Options[NInverseLaplaceTransform] = {Method -> Automatic,
        "Points" -> 6, WorkingPrecision -> MachinePrecision};

NInverseLaplaceTransform[f_, s_, t_, opts : OptionsPattern[]] :=
        Module[{met, xa, wt},
               met = Switch[OptionValue[Method],
                            "SalzerPiessens" | Automatic, SalzerPiessensRuleData,
                            "Salzer", SalzerRuleData, _, Return[$Failed]];
               {xa, wt} = met[OptionValue["Points"],
               Chop[(wt.Function[s, f][1/xa/t])/t]]

For error estimation purposes, one could then generate the approximant corresponding to both Salzer and Salzer-Piessens for the same setting of "Points" and look at their difference.

Using the Bessel function example from the previous post,

f1[t_] = NInverseLaplaceTransform[1/Sqrt[1 + s^2], s, t, 
          "Points" -> 9, WorkingPrecision -> 30];
f2[t_] = NInverseLaplaceTransform[1/Sqrt[1 + s^2], s, t, 
          Method -> "Salzer", "Points" -> 9, WorkingPrecision -> 30];
GraphicsGrid[{{Plot[{BesselJ[0, t], Re[f1[t]]}, {t, 0, 19}, 
                    Frame -> True, ImageSize -> Medium, 
                    PlotLegends -> {{TraditionalForm[BesselJ[0, t]],
                    PlotStyle -> {AbsoluteThickness[4], Automatic}],
              {Plot[Re[f1[t] - f2[t]], {t, 0, 19}, Frame -> True, 
                    PlotLabel -> Style["Difference between Salzer-Piessens
                                        and Salzer methods", Small],
                    PlotRange -> All], 
               Plot[Re[f1[t] - BesselJ[0, t]], {t, 0, 19}, Frame -> True, 
                    PlotLabel -> Style["Difference between Salzer-Piessens and
                                        true transform", Small],
                    PlotRange -> All]}}]


A compact Hilbert curve routine

January 30, 2017

Here is a relatively compact Mathematica routine for generating the b-th iterate of a Hilbert curve in [0,1]^n. The algorithm is due to Skilling:

HilbertCurve =
     Compile[{{x, _Integer}, {b, _Integer}, {n, _Integer}},
             Module[{t = BitXor[x, Quotient[x, 2]], p = 2, k, q, r, xx},
                    xx = Total[Table[BitAnd[Quotient[t, 2^(r (n - 1) + k)], 2^r],
                                     {r, b - 1, 0, -1}, {k, n - 1, 0, -1}]];
                    Do[q = p - 1;
                       Do[t = xx[[k]];
                          If[BitAnd[t, p] != 0,
                             xx[[1]] = xx[[1]] ~BitXor~ q,
                             t = BitAnd[t ~BitXor~ xx[[1]], q];
                             xx[[{1, k}]] = BitXor[xx[[{1, k}]], t]],
                          {k, n, 2, -1}];
                       If[BitAnd[xx[[1]], p] != 0,
                          xx[[1]] = xx[[1]] ~BitXor~ q];
                       p *= 2, {r, b - 1}];
                    2 xx/(p - 1) - 1], RuntimeAttributes -> {Listable}]

I’ve seen a variety of routines for generating the Hilbert curve before, but none were as compact or as general as this one. If you want a version that generates exact rational coordinates, it is easy to take out the requisite parts from the Compile[] function.

Here are the 2D and 3D versions of the curve:

With[{p = 6}, 
     Graphics[{ColorData[97, 1],
               Line[HilbertCurve[Range[0, 2^(2 p) - 1], p, 2]]}, Frame -> True]]


With[{p = 3}, 
     Graphics3D[Tube[HilbertCurve[Range[0, 2^(3 p) - 1], p, 3]], Axes -> True]]


Here’s a perspective projection of a four-dimensional Hilbert curve:

With[{p = 3, k = 2, f = 4, d = 3/2, r = 1/8}, 
     Graphics3D[{Directive[EdgeForm[], CapForm[None], GrayLevel[1/5], 
                           Glow[ColorData["Legacy", "SteelBlue"]],
                           Specularity[ColorData["Legacy", "Gold"], 16]], 
                 Tube[(f Delete[#, k])/(d - Extract[#, k]) & /@ 
                      HilbertCurve[Range[0, 2^(4 p) - 1], p, 4], r]}, 
                Boxed -> False, Lighting -> "Neutral"]]


(This visualization is admittedly a bit messy-looking; if you have suggestions for making a nice display of the four-dimensional version, tell me about it in the comments!)

Some terms left unsimplified

January 6, 2017

I had been slowly going through some OEIS entries in an attempt to add or optimize the Mathematica code found there. One thing that struck me with respect to the sequences that are related to the elliptic integrals/elliptic functions was that they were often left in the form that was spit out by a CAS. I had railed about the deficiencies of current software for generating simple expressions for elliptic integrals many times before, so this was but a confluence with the fact that not many people these days have experience in the analytical manipulation of these functions, as classical and important as they are.

I certainly have my work cut out for me…

Yes, I’m still here…

December 21, 2014

…and to those guys who showed concern and worry during my long-ish absence, I’m fine. Two important things, though: the hard drive containing a lot of my unpublished research crashed, including a good amount that I have not been able to back up; and, I have been quite busy with a (relatively) new job over the past few months, details of which I cannot discuss here.

In short, I’m fine, and I hope you guys will still be here when I’m actually back. Happy Holidays!

~ Jan

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