Mathematica

Counting permutations with roots

Posted on by John

My post from yesterday on permutation roots ends with a Mathematica code for finding the probability that a permutation of n elements has a kth root. This is done by finding the coefficient of xn in the generating function

\prod_{m=1}^\infty \exp_{\text{gcd}(m, k)} \left\frac{x^m}{m}\right)

I wanted to say more about this, and look at implementing the same code in SymPy. I was curious how well SymPy would do because I’ve noticed that LLMs often generate SymPy code since it’s an open source CAS.

Wilf [1] describes the infinite product above as the exponential generating function (egf) of f(n, k), the number of permutations of n objects that have a kth root. Since egfs have a n! term in the denominator, this is also the ordinary generating function (ogf) of the probability that a randomly chosen permutation on n objects has a kth root.

My first attempt at using Mathematica to probe the generating function was

expq[x_, q_] := MittagLefflerE[q, x^q]	 
p[n_, k_] :=  SeriesCoefficient[	 
    Product[expq[x^m/m, GCD[m, k]], {m, 1, Infinity}], {x, 0, n}]

This hung forever when I tried to use it on a small example. I realized, but apparently Mathematica did not, that Infinity could be replaced by n since terms higher than n do not contribute to the coefficient of xn. With that change, the code ran quickly.

This morning I tried converting the Mathematica code to Sympy; Claude did this in one shot. I also reproduced the table of f(n, k) values on page 150 of [1] to test the code. Since Wilf tabulated f(n, k), not f(n, k)/n!, I multiplied the results by n!.

Here is the output:

k = 2 [1, 1, 3, 12, 60, 270, 1890, 14280, 128520, 1096200]
k = 3 [1, 2, 4, 16, 80, 400, 2800, 22400, 181440, 1814400]
k = 4 [1, 1, 3, 12, 60, 270, 1890, 13020, 117180, 1039500]
k = 5 [1, 2, 6, 24, 96, 576, 4032, 32256, 290304, 2612736]
k = 6 [1, 1, 1, 4, 40, 190, 1330, 8680, 52920, 340200]
k = 7 [1, 2, 6, 24, 120, 720, 4320, 34560, 311040, 3110400]

and here is the SymPy code. I edited the main but the rest is verbatim from Claude.

from sympy import symbols, gcd, factorial, Rational, S

x = symbols('x')

def expq_coeffs(m, q, n):
    """
    Truncated (degree <= n) series coefficients of
        expq(x**m/m, q) = MittagLefflerE(q, (x**m/m)**q)
    Since q is a positive integer:
        E_q(y^q) = sum_j y^(q*j) / (q*j)!
    with y = x**m/m, so the term of degree m*q*j has coefficient
        1 / ( m**(q*j) * (q*j)! ).
    Returns a list c[0..n] of coefficients.
    """
    c = [S.Zero] * (n + 1)
    j = 0
    while m * q * j <= n:
        deg = m * q * j
        c[deg] += Rational(1, m**(q * j) * factorial(q * j))
        j += 1
    return c

def poly_mult_trunc(a, b, n):
    """Multiply two series (lists of coeffs, index = degree) truncated to degree n."""
    c = [S.Zero] * (n + 1)
    for i, ai in enumerate(a):
        if ai == 0:
            continue
        max_j = n - i
        for j2 in range(max_j + 1):
            bj = b[j2]
            if bj != 0:
                c[i + j2] += ai * bj
    return c

def p(n, k):
    """
    SymPy equivalent of:
        expq[x_, q_] := MittagLefflerE[q, x^q]
        p[n_, k_] := SeriesCoefficient[
            Product[expq[x^m/m, GCD[m, k]], {m, 1, n}], {x, 0, n}]
    """
    result = [S.Zero] * (n + 1)
    result[0] = S.One
    for m in range(1, n + 1):
        q = gcd(m, k)
        factor = expq_coeffs(m, q, n)
        result = poly_mult_trunc(result, factor, n)
    return result[n]

# example
if __name__ == "__main__":
    for k in range(2, 8):
        print("k =", k, [factorial(n)*p(n, k) for n in range(1,11)])

[1] Herbert Wilf. Generatingfunctionology. Available online here.

Closer look at an identity

Posted on by John

The previous post derived the identity

\cosh\Big( \operatorname{arccosh}(x) + \operatorname{arccosh}(y)\Big) = xy + \sqrt{x^2 -1} \sqrt{y^2 -1}

and said in a footnote that the identity holds at least for x > 1 and y > 1. That’s true, but let’s see why the footnote is necessary.

Let’s have Mathematica plot

\left| \cosh\Big( \operatorname{arccosh}(x) + \operatorname{arccosh}(y)\Big) - xy - \sqrt{x^2 -1} \sqrt{y^2 -1} \right|

The plot will be 0 where the identity above holds.

The plot is indeed flat for x > 1 and y > 1, and more, but not everywhere.

If we combine the two square roots

\left| \cosh\Big( \operatorname{arccosh}(x) + \operatorname{arccosh}(y)\Big) - xy - \sqrt{(x^2 -1) (y^2 -1)} \right|

and plot again we still get a valid identity for x > 1 and y > 1, but the plot changes.

This is because √a √b does not necessarily equal √(ab) when the arguments may be negative.

The square root function and the arccosh function are not naturally single-valued functions. They require branch cuts to force them to be single-valued, and the two functions require different branch cuts. I go into this in some detail here.

There is a way to reformulate our identity so that it holds everywhere. If we replace

\sqrt{z^2 - 1}

with

f(z) = \exp\left( \tfrac{1}{2}\left( \log(z - 1) + \log(z + 1) \right) \right)

which is equivalent for z > 1, the corresponding identity holds everywhere.

We can verify this with the following Mathematica code.

f[z_] := Exp[(1/2) (Log[z - 1 ] + Log[z + 1])]
FullSimplify[Cosh[ArcCosh[x] + ArcCosh[y]] - x y - f[x] f[y]]

This returns 0.

By contrast, the code

FullSimplify[
 Cosh[ArcCosh[x] + ArcCosh[y]] - x y - Sqrt[x^2 - 1] Sqrt[y^2 - 1]]

simply returns its input with no simplification, unless we add restrictions on x and y. The code

FullSimplify[
 Cosh[ArcCosh[x] + ArcCosh[y]] - x y - Sqrt[x^2 - 1] Sqrt[y^2 - 1], 
 Assumptions -> {x > -1 && y > -1}]

does return 0.

Related posts

sinh( arccosh(x) )

Posted on by John

I’ve written several posts about applying trig functions to inverse trig functions. I intended to write two posts, one about the three basic trig functions and one about their hyperbolic counterparts. But there’s more to explore here than I thought at first. For example, the mistakes that I made in the first post lead to a couple more posts discussing error detection and proofs.

I was curious about how Mathematica would handle these identities. Sometimes it doesn’t simplify expressions the way you expect, and for interesting reasons. It handled the circular functions as you might expect.

\renewcommand{\arraystretch}{2.2} \begin{array}{c|c|c|c} & \sin^{-1} & \cos^{-1} & \tan^{-1} \\ \hline \sin & x & \sqrt{1-x^{2}} & \dfrac{x}{\sqrt{1+x^2}} \\ \hline \cos & \sqrt{1-x^{2}} & x & \dfrac{1}{\sqrt{1 + x^2}} \\ \hline \tan & \dfrac{x}{\sqrt{1-x^{2}}} & \dfrac{\sqrt{1-x^{2}}}{x} & x \\ \end{array}

So, for example, if you enter Sin[ArcCos[x]] it returns √(1 − x²) as in the table above. Then I added an h on the end of all the function names to see whether it would reproduce the table of hyperbolic compositions.

\renewcommand{\arraystretch}{2.2} \begin{array}{c|c|c|c} & \sinh^{-1} & \cosh^{-1} & \tanh^{-1} \\ \hline \sinh & x & \sqrt{x^{2}-1} & \dfrac{x}{\sqrt{1-x^2}} \\ \hline \cosh & \sqrt{x^{2} + 1} & x & \dfrac{1}{\sqrt{1 - x^2}} \\ \hline \tanh & \dfrac{x}{\sqrt{x^{2}+1}} & \dfrac{\sqrt{x^{2}-1}}{x} & x \\ \end{array}

For the most part it did, but not entirely. The results were as expected except when applying sinh or cosh to arccosh. But Sinh[ArcCosh[x]] returns

\sqrt{\frac{x-1}{x+1}} (x+1)

and Tanh[ArcCosh[x]] returns

\frac{\sqrt{\frac{x-1}{x+1}} (x+1)}{x}

Why doesn’t Mathematica simplify as expected?

Why didn’t Sinh[ ArcCosh[x] ] just return √(x² − 1)? The expression it returned is equivalent to this: just square the (x + 1) term, bring it inside the radical, and simplify. That line of reasoning is correct for some values of x but not for others. For example, Sinh[ArcCosh[2]] returns −√3 but √(2² − 1) = √3. The expression Mathematica returns for Sinh[ArcCosh[x]] correctly evaluates to −√3.

Defining ArcCosh

To understand what’s going on, we have to look closer at what arccosh(x) means. You might say it is a function that returns the number whose hyperbolic cosine equals x. But cosh is an even function: cosh(−x) = cosh(x), so we can’t say the value. OK, so we define arccosh(x) to be the positive number whose hyperbolic cosine equals x. That works for real values of x that are at least 1. But what do we mean by, for example, arccosh(1/2)? There is no real number y such that cosh(y) = 1/2.

To rigorously define inverse hyperbolic cosine, we need to make a branch cut. We cannot define arccosh as an analytic function over the entire complex plane. But if we remove (−∞, 1], we can. We define arccosh(x) for real x > 1 to be the positive real number y such that cosh(y) = x, and define it for the rest of the complex plane (with our branch cut (−∞, 1] removed) by analytic continuation.

If we look up ArcCosh in Mathematica’s documentation, it says “ArcCosh[z] has a branch cut discontinuity in the complex z plane running from −∞ to +1.” But what about values of x that lie on the branch cut? For example, we looked at ArcCosh[-2] above. We can extend arccosh to the entire complex plane, but we cannot extend it as an analytic function.

So how do we define arccosh(x) for x in (−∞, 1]? We could define it to be the limit of arccosh(z) as z approaches x for values of z not on the branch cut. But we have to make a choice: do we approach x from above or from below? That is, we can define arccosh(x) for real x ≤ 1 by

\text{arccosh}(x) = \lim_{\varepsilon \to 0^+} \text{arccosh}(x + \varepsilon i)

or by

\text{arccosh}(x) = \lim_{\varepsilon \to 0^-} \text{arccosh}(x + \varepsilon i)

but we have to make a choice because the two limits are not the same. For example, ArcCosh[-2 + 0.001 I] returns 1.31696 + 3.14102 I but ArcCosh[-2 - 0.001 I] returns 1.31696 - 3.14102 I. By convention, we choose the limit from above.

Defining square root

Where did we go wrong when we assumed Mathematica’s expression for sinh(arccosh(x))

\sqrt{\frac{x-1}{x+1}} (x+1)

could be simplified to √(x² − 1)? We implicitly assumed √(x + 1)² = (x + 1). And that’s true, if x ≥ − 1, but not for smaller x. Just as we have be careful about how we define arccosh, we have to be careful about how we define square root.

The process of defining the square root function for all complex numbers is analogous to the process of defining arccosh. First, we define square root to be what we expect for positive real numbers. Then we make a branch cut, in this case (−∞, 0]. Then we define it by analytic continuation for all values not on the cut. Then finally, we define it along the cut by continuity, taking the limit from above.

Once we’ve defined arccosh and square root carefully, we can see that the expressions Mathematica returns for sinh(arccosh(x)) and tanh(arccosh(x)) are correct for all complex inputs, while the simpler expressions in the table above implicitly assume we’re working with values of x for which arccosh(x) is real.

Making assumptions explicit

If we are only concerned with values of x ≥ − 1 we can tell Mathematica this, and it will simplify expressions accordingly. If we ask it for

    Simplify[Sinh[ArcCosh[x]], Assumptions -> {x >= -1}]

it will return √(x² − 1).

Related posts

Solving H_n = 100

Posted on by John

The previous post includes code for solving the equation

Hn = m

i.e. finding the value of n for which the nth harmonic number is the closest to m. It works well for small values of m. It works for large m in the sense that the solution is very close to m, but it’s not necessarily the best solution.

For example, set m = 100. The code returns

n = 15092688622113830917200248731913020965388288

and indeed for that value of n,

Hn − 100 ≈ 3 × 10−15

and that’s as much as we could hope for with IEEE 754 floats.

The approximation

n = exp(m −γ)

is very good for large values of m. Using Mathematica we can find the exact value of n.

f[n_] := Log[n] + EulerGamma + 1/(2 n) - 1/(12 n^2)
n = Floor[Exp[100 - EulerGamma]];
N[f[n], 50]
100.00000000000000000000000000000000000000000000900
N[f[n - 1], 50]
99.999999999999999999999999999999999999999999942747

So

n = 15092688622113788323693563264538101449859497

A similar process can find the solution to

Hn = 1000

is

n = 110611511026604935641074705584421138393028001852577373936470952377218354575172401275457597579044729873152469512963401398362087144972181770571895264066114088968182356842977823764462179821981744448731785408629116321919957856034605877855212667092287520105386027668843119590555646814038787297694678647529533718769401069269427475868793531944696435696745559289326610132208504257721469829210704462876574915362273129090049477919400226313586033

For this calculation you’ll need to increase the precision from 50 digits to something like 500 digits, something more than 435 because n is a 435-digit number.

In case you’re wondering whether my function for computing harmonic numbers is accurate enough, it’s actually overkill, with error O(1/120n4).

Time to factor big integers Python and Mathematica

Posted on by John

This post will look at the time required to factor n − 1 each of the following prime numbers in Python (SymPy) and Mathematica. The next post will explain why I wanted to factor these numbers.

p = 2254 + 4707489544292117082687961190295928833
q = 2254 + 4707489545178046908921067385359695873
r = 2254 + 45560315531419706090280762371685220353
s = 2254 + 45560315531506369815346746415080538113

Here are the timing results.

    |   |   Python | Mathematica |
    |---+----------+-------------|
    | p |    0.913 |       0.616 |
    | q |    0.003 |       0.002 |
    | r |  582.107 |      14.915 |
    | s | 1065.925 |      20.763 |

This is hardly a carefully designed benchmark, but it’s enough to suggest Mathematica can be a couple orders of magnitude faster than Python.

Here are the factorizations.

p − 1 = 234 × 3 × 4322432633228119 × 129942003317277863333406104563609448670518081918257
q − 1 = 233 × 3 × 5179 × 216901160674121772178243990852639108850176422522235334586122689
r − 1 = 232 × 32 × 463 × 539204044132271846773 × 8999194758858563409123804352480028797519453
s − 1 = 232 × 32 × 1709 × 24859 × 1690502597179744445941507 × 10427374428728808478656897599072717

Powers of 3 + √2

Posted on by John

Yesterday’s post looked at the distribution of powers of x mod 1. For almost all x > 1 the distribution is uniform in the limit. But there are exceptions, and the post raised the question of whether 3 + √2 is an exception.

A plot made it look like 3 + √2 is an exception, but that turned out to be a numerical problem.

A higher precision calculation showed that the zeros on the right end of the plot were erroneous.

So this raises the question of how to calculate (3 + √2)n accurately for large n. The way I created the second plot was to use bc to numerically calculate the powers of 3 + √2. In this post, I’ll look at using Mathematica to calculate the powers symbolically.

For all positive integers n,

(3 + √2)n = an + bn√2

where an and bn are positive integers. We want to compute the a and b values.

If you ask Mathematica to compute (3 + √2)n it will simply echo the expression. But if you use the Expand function it will give you want. For example

    Expand[(3 + Sqrt[2])^10]

returns

    1404491 + 993054 √2

We can use the Coefficient function to split a + b √2 into a and b.

    parts[n_] := 
        Module[{x = (3 + Sqrt[2])^n}, 
            {Coefficient[x, Sqrt[2], 0], Coefficient[x, Sqrt[2], 1]}]

Now parts[10] returns the pair {1404491, 993054}.

Here’s something interesting. If we set

(3 + √2)n = an + bn√2

as above, then the two halves of the expression on the right are asymptotically equal. That is, as n goes to infinity, the ratio

an /  bn√2

converges to 1.

We can see this by defining

    ratio[n_] := 
        Module[ {a = Part[ parts[n], 1], b = Part[parts[n], 2]}, 
        N[a / (b Sqrt[2])]]

and evaluating ratio at increasing values of n. ratio[12] returns 1.00001 and ratio[13] returns 1, not that the ratio is exactly 1, but it is as close to 1 as a floating point number can represent.

This seems to be true more generally, as we can investigate with the following function.

    ratio2[p_, q_, r_, n_] := 
        Module[{x = (p + q Sqrt[r])^n}, 
            N[Coefficient[x, Sqrt[r], 0]/(Coefficient[x, Sqrt[r], 1] Sqrt[r])]]

When r is a prime and

(p + q√r)n = an + bn√r

then it seems that the ratio an / bn √r converges to 1 as n goes to infinity. For example, ratio2[3, 5, 11, 40] returns 1, meaning that the two halves of the expression for (3 + 5√11)n are asymptotically equal.

I don’t know whether the suggested result is true, or how to prove it if it is true. Feels like a result from algebraic number theory, which is not something I know much about.

Update: An anonymous person on X suggested a clever and simple proof. Observe that

\begin{align*} a_n &= \frac{(3+\sqrt{2})^n + (3-\sqrt{2})^n}{2} \\ b_n\sqrt{2} &= \frac{(3 + \sqrt{2})^n - (3-\sqrt{2})^n}{2} \end{align*}

In this form it’s clear that the ratio an / bn √2 converges to 1, and the proof can be generalized to cover more.

How Mathematica Draws a Dragonfly

Posted on by John

Mathematica includes code to draw various whimsical images. For example, if you enter the following command in Mathematica

    Entity["PopularCurve", "DragonflyCurve"][
        EntityProperty["PopularCurve", "Image"]]

you get an image of a dragonfly.

It draws such images with Fourier series. You can tell by asking for the parameterization of the curve. If you enter

    Entity["PopularCurve", "DragonflyCurve"][
        EntityProperty["PopularCurve", "ParametricEquations"]]

you’ll get the following, after some rearrangement.

    Function[t, {
        (7714/27) Sin[47/20 - t] + (1527/37) Sin[16/5 - 2t] + 
        (3202/39) Sin[108/41 - 3t] + … + 2/9 Sin[15/19 - 81 t], 
        (9406/37) Sin[29/7 - t] + (3591/53) Sin[28/13 - 2t] + 
        (1111/20) Sin[9/23 - 3t] + … -(3/29) Sin[8/23 + 81 t]
        }]

The function is a parameterized curve, taking t to (x(t), y(t)) where x(t) and y(t) are Fourier series including frequencies up to sin(81t). Each of the sine terms has a phase shift that could be eliminated by expressing sin(φ + ωt) as a linear combination of sin(ωt) and cos(ωt).

Presumably somebody drew the dragonfly, say in Adobe Illustrator or Inkscape, then did a Fourier transform of a sampling of the curve.

To make sure Mathematica wasn’t doing anything behind the scenes that I wasn’t aware of, I reproduced the dragonfly curve by porting the Mathematica code to Python.

The number of Fourier components needed to draw an image depends on the smoothness and complexity of the image. The curve for π, for example, the highest frequency component is sin(32t).

The triceratops curve is more complicated and Mathematica uses frequencies up to sin(188t).

 

n-queens in 3D

Posted on by John

It’s possible to place n queens on an n × n chess board so that no queen is attacking any other queen, provided n > 3. A queen is able to move any number of squares in any horizontal, vertical, or diagonal direction. Another way to put it is that the queen can move in any of 8 = 3² − 1 directions, in the direction of any cell in a 3 × 3 grid centered on the queen.

What if we extend chess to three dimensions? Now we have an n × n × n cube. Now a queen is able to move in 26 = 3³ − 1 directions, in the direction of any cell in a 3 × 3 × 3 cube centered on the queen.

Klarner [1] proved that it is possible to place n² queens in an n × n × n cube so that no queen is attacking any other, provided gcd(n, 210) = 1, i.e. provided the smallest prime factor of n is larger than 7. The condition gcd(n, 210) = 1 is sufficient, and it is conjectured to be necessary as well [2].

Klarner constructed a solution as follows: Place the queens on

(i, j, (3i + 5j) mod n)

for i and j running from 1 to n.

One way to visualize the queens in 3D is to draw a 2D grid where each cell contains the vertical coordinate of the corresponding queen. This grid will be a Latin square, and so the 3D queen placement problem is also known as the Latin queen problem.

Here’s a visualization of Klarner’s example.

Note that if you only pay attention to a particular number, you have a solution to the 11 queens problem in a square. That’s because every slice of a 3D solution is a solution to the 2D problem.

I played around with plotting the points in three dimensions. Here’s one view.

A slight rotation of the cube gives a substantially different perspective. This would be better as an animation, available here.

Mathematica code

Here’s the code that made the images above.

Flat grid:

    Grid[
        Table[
            Mod[3 i + 5 j, 11], 
            {i, 0, 10}, 
            {j, 0, 10}
        ], 
        Frame -> All
    ]

Static 3D view:

     ListPointPlot3D[
         Table[
             {i, j, Mod[3 i + 5 j, 11]}, 
             {i, 0, 10}, 
             {j, 0, 10}
          ], 
          BoxRatios -> {1, 1, 1}, 
          PlotStyle -> {PointSize[0.02]}
      ]

Animation:

    Animate[
        ListPointPlot3D[
            Table[
                {i, j, Mod[3 i + 5 j, 11]}, 
                {i, 0, 10}, 
                {j, 0, 10}
            ], 
        BoxRatios -> {1, 1, 1}, 
        PlotStyle -> {PointSize[0.02]}, 
        ViewPoint -> {2 Cos[t], 2 Sin[t], 1}
        ], 
        {t, 0, Pi, 0.05}, 
        AnimationRunning -> True, 
        AnimationRate -> 4
    ]

    Export["klerner_animation.gif", %]

Related posts

[1] D.A. Klarner, Queen squares, J. Recreational Math. 12 (3) (1979/1980) 177–178.

[2] Jordan Bell and Brett Stevens. A survey of known results and research areas for n-queens. Discrete Mathematics 309 (2009) 1–31

Duplicating a hand-drawn contour plot

Posted on by John

Like the previous post, this post also seeks to approximately reproduce a hand-drawn plot. This time the goal is reproduce figure 7.3 from A&S page 298.

This plot is a visualizing of the function of a complex variable

w(z) = exp(−z²) erfc(− iz)

where erfc is the complementary error function.

A&S calls the graph above an “altitude chart.” This could be a little misleading since it’s the overlay of two plots. One plot is the absolute value of w(z), which could well be called altitude. But it also contains a plot of the phase. To put it another way, if we denote the values of the function in polar form r exp(iθ) the altitude chart is an overlay of a plot of r and a plot of θ.

We begin by defining

f[z_] := Exp[-z^2] Erfc[-I z]

The following code reproduces the lines of constant phase fairly well.

ContourPlot[Arg[f[x + I y]], {x, 0.1, 3.1}, {y, -2.5, 3}, 
    Contours -> 20, ContourShading -> None, AspectRatio -> 1.6]

The lines of constant absolute value take a little more effort to reproduce. If we let Mathematica pick where to put the contour lines, they will not be distributed the same way they were in A&S.

ContourPlot[Abs[f[x + I y]], {x, 0, 3.1}, {y, -2.6, 3}, 
    Contours -> 20, ContourShading -> None, AspectRatio -> 1.6]

We can duplicated the spacing in the original plot by providing Mathematica a list of contour values rather than number of contour values.

ContourPlot[Abs[f[x + I y]], {x, 0, 3.1}, {y, -2.6, 3}, 
    Contours -> {0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 2, 3, 4, 5, 10, 100}, 
    ContourShading -> None, AspectRatio -> 1.6]

(For reasons I don’t understand, Mathematica does not draw the contours corresponding to w = 10 and w = 100.)

When I overlay the phase and absolute value plots with the Show command I get a plot approximately reproducing the original.

Related posts

Reproducing a hand-drawn plot

Posted on by John

The plots in old (i.e. pre-computer) math books are impressive. These plots took a lot of effort to produce, and so they weren’t created lightly. Consequently they tend to be aesthetically and mathematically interesting. A few weeks ago I recreated a plot from A&S using Mathematica, and today I’d like to do the same for the plot below from a different source [1].

Here is my approximate reproduction.

I’ll give the mathematical details below for those who are interested.

The plot shows “normalized associated Legendre functions of the first kind.” There are two families of Legendre polynomials, denoted P and Q; we’re interested in the former. “Associated” means polynomials that are derived the Legendre polynomials by taking derivatives. Normalized means the polynomials are multiplied by constants so that their squares integrate to 1 over [−1, 1].

Mathematica has a function LegendreP[n, m, x] that implements the associated Legendre polynomials Pnm(x). I didn’t see that Mathematica has a function for the normalized version of these functions, so I rolled by own.

f[n_, m_, x_] := (-1)^n LegendreP[n, m, x] 
        Sqrt[(2 n + 1) Factorial[n - m]/(2 Factorial[n + m])]

I added the alternating sign term up front after discovering that apparently the original plot used a different convention for defining Pnm than Mathematica uses.

I make my plot by stacking the plots created by

Plot[Table[f[n, n, x], {n, 1, 8}],  {x, 0, 1}]

and

Plot[Table[f[n + 4, n, x], {n, 0, 4}],  {x, 0, 1}]

The original plot shows P4(x). I used the fact that this equals P40(x) to simplify the code. I also left out the flat line plotting P0 because I thought that looked better.

Related post: Duplicating a Hankel function plot from A&S.

[1] Tables Of Functions With Formulae And Curves by Fritz Emde, published in 1945. Available on Archive.org.