Laguerre Polynomials

Code Golf Asked by golf69 on October 27, 2021

Laguerre polynomials are solutions to Laguerre’s equation, a second-order linear differential equation: $$xy”+(1-x)y’+ny=0$$. For a given value of n, the solution, y, is named $$L_n(x)$$.

The polynomials can be found without calculus using recursion:

$$L_0(x)=1$$

$$L_1(x)=1-x$$

$$L_{k+1}(x)=frac{(2k+1-x)L_k(x)-kL_{k-1}(x)}{k+1}$$

Summation can be used to the same end:

$$L_n(x)=sumlimits_{k=0}^{n}{nchoose k}frac{(-1)^k}{k!}x^k$$

$$L_n(x)=sumlimits_{i=0}^nprodlimits_{k=1}^ifrac{-(n-k+1)x}{k^2}$$

The first Laguerre polynomials are as follows:

Coefficients can be found here.

The Challenge

Given a nonnegative integer n and a real number x, find $$L_n(x)$$.

Rules

• This is so the shortest answer in bytes wins.

• Assume only valid input will be given.

• Error should be under one ten-thousandth (±0.0001) for the test cases.

Test Cases

Here, n is the first number and x is the second.

In: 1 2
Out: -1

In: 3 1.416
Out: -0.71360922

In: 4 8.6
Out: −7.63726667

In: 6 -2.1
Out: 91.86123261


C (gcc), 91 bytes

i;k;float f(n,x)float x;{float p,s=0;for(i=++n;k=i--;s+=p)for(p=1;--k;)p*=(k-n)*x/k/k;x=s;}


Try it online!

Straighforward implementation of polynomial expansion. Slightly golfed less

i;k;
float f(n,x)float x;{
float p,s=0;
for(i=++n;k=i--;s+=p)
for(p=1;--k;)
p*=(k-n)*x/k/k;
x=s;
}


Answered by ceilingcat on October 27, 2021

Fortran (GFortran), 69 68 bytes

read*,n,a
print*,sum([(product([((j-n-1)*a/j/j,j=1,i)]),i=0,n)])
end


-1 byte thanks to @ceilingcat

The program reads in an implicit integer n and real a. Summation and product operations are performed using arrays (initialized using implicit loops) with the intrinsics sum() and product().

Try it online!

Answered by Roninkoi on October 27, 2021

Pari/GP, 39 bytes

Using the formula $$L_n(x)=sum_{k=0}^n binom{n}{k}frac{(-1)^k}{k!} x^k$$.

l(n,x)=sum(k=0,n,n!*(-x)^k/(n-k)!/k!^2)


Try it online!

Pari/GP, 45 bytes

Using the generating function $$sum_{n=0}^infty x^n L_n(t)= frac{1}{1-x} e^{-xt/(1-x)}$$.

l(n,t)=Vec(exp(-x*t/(1-x)+O(x^n++))/(1-x))[n]


Try it online!

Answered by alephalpha on October 27, 2021

MATL, 5 bytes

_1iZh


Inputs are $$n$$, then $$x$$. Try it online! Or verify all test cases.

How it works

This uses the equivalence of Laguerre polynomials and the (confluent) hypergeometric function:

$$L_n(x) = {} _1F_1(-n,1,x)$$

_    % Implicit input: n. Negate
1    % Push 1
i    % Input: x
Zh   % Hypergeometric function. Implicit output


Answered by Luis Mendo on October 27, 2021

JavaScript (Node.js), 36 bytes

x=>(i=0,g=n=>n?1-x*n/++i/i*g(n-1):1)


Try it online!

Just convert the formula to this, and use recursive:

$$L_n(x) = sum_{i=0}^nprod_{k=1}^ifrac{-(n-k+1)x}{k^2}$$

Answered by tsh on October 27, 2021

Jelly, 11 bytes

cŻ÷Ż!$ƲṚḅN}  A dyadic Link accepting $$n$$ on the left and $$x$$ on the right which yields $$L_n(x)$$. Try it online! How? This makes the observation that $$L_n(x)=sumlimits_{k=0}^{n}{nchoose k}frac{(-1)^k}{k!}x^k=sumlimits_{k=0}^{n}{(-x)^k}frac{nchoose k}{k!}$$ which is the evaluation of a base $$-x$$ number with n+1 digits of the form $$frac{nchoose k}{k!}$$. cŻ÷Ż!$ƲṚḅN} - Link: n, x
Ż          -   zero-range (n) -> [0, 1, 2, ..., n]
c           -   (n) binomial (that) -> [nC0, nC1, nC2, ..., nCn]
$- last two links as a monad - g(n): Ż - zero-range (n) -> [0, 1, 2, ..., n] ! - factorial (that) -> [0!, 1!, 2!, ..., n!] ÷ - division -> [nC0÷0!, nC1÷1!, nC2÷2!, ..., nCn÷n!] Ṛ - reverse -> [nCn÷n!, ..., nC2÷2!, nC1÷1!, nC0÷0!] } - use the chain's right argument for: N - negate -> -x ḅ - convert from base (-x) -> -xⁿnCn÷n!+...+-x²nC2÷2!+-x¹nC1÷1!+-x°nC0÷0!  Answered by Jonathan Allan on October 27, 2021 APL (Dyalog Unicode), 16 bytes 1⊥⍨0,⎕×(-÷⌽×⌽)⍳⎕  Try it online! A full program that takes n and x from two separate lines of stdin. How it works 1⊥⍨0,⎕×(-÷⌽×⌽)⍳⎕ ⍳⎕ ⍝ Take n and generate 1..n (-÷⌽×⌽) ⍝ Compute i÷(n+1-i)^2 for i←1..n 0,⎕× ⍝ Multiply x to each and prepend 0, call it B 1⊥⍨ ⍝ Convert all ones from base B to single number  The mixed base conversion looks like this: 1..n: ... n-3 n-2 n-1 1 B: 0 ... (n-3)x/4^2 (n-2)x/3^2 (n-1)x/2^2 nx digits: 1 ... 1 1 1 1 digit values: x^n/n! ... (nC3 x^3/3!) (nC2 x^2/2!) (nC1 x^1/1!) (nC0 x^0/0!)  It is essentially a fancy way to write the sum of product scan over 1, nx, (n-1)x/2^2, (n-2)x/3^2, .... This happens to be shorter than a more straightforward -x-base conversion (evaluating a polynomial at -x): APL (Dyalog Unicode), 18 bytes (-⎕)⊥⌽1,(!÷⍨⊢!≢)⍳⎕  Try it online! Answered by Bubbler on October 27, 2021 J, 37 20 bytes -5 thanks to @Bubbler Calculates the polynomial adapted from the summation formula and uses J's p. operator to calculate that polynomial with a given x. (p.-)~i.((!]/)%!)@,]  Try it online! J, 45 byte Alternative Recursive function. 1:-@.[~ ::((>:@]%~($:*[-~1+2*])-]*(\$:<:))<:)


Try it online!

How it works

We define a hook (fg), which is x f (g n). f is (p.-)~ so it will be evaluated as ((i.((!]/)%!)@,]) n) p. (- x).

(p.-)~i.((!]/)%!)@,]
i.         @,] enumerate 3 -> 0 1 2, append 3 -> 0 1 2 3, …
(!]/)       3 over i
%      divided by
!     !i
-                 negate x
p.                  apply -x to the polynomial expressed in J as
1 3 1.5 0.166667, so 1-3(-x)+1.5(-x)^2+0.16(-x)^3


Answered by xash on October 27, 2021

Japt-x, 28 27 26 bytes

ò@l *VpX /Xl ²*JpX /(U-X l


Try it

Japt, 30 29 28 bytes

ò x@l *VpX /Xl ²*JpX /(U-X l


Try it

Explanation

ò x@l *VpX /Xl ²*JpX /(U-X l
ò                               // Create a array [0, 1, ..., U]
x                             // sum the array after mapping through
@                            // Function(X)
l                           //    U!
*VpX                      //    times V ** X
/Xl ²                //    divided by X! ** 2
*JpX            //    times (-1) ** X
/(U-X l    //    divided by (U - X)!

• U is the first input
• V is the second input
• ** represents exponentiation
• ! represents factorial

Answered by Mukundan314 on October 27, 2021

05AB1E, 16 bytes

1λèN·<I-₁*N<₂*-N/


Try it online. (No test suite for all test cases at once, since there seems to be a bug in the recursive environment..)

Explanation:

 λ                # Create a recursive environment
è               # to get the 0-based n'th value afterwards
# (where n is the first implicit input)
# (which will be output implicitly as result in the end)
1                 # Starting with a(-1)=0 and a(0)=1,
# and for every other a(N), we'll:
#  (implicitly push a(N-1))
N·             #  Push N doubled
<            #  Decrease it by 1
I-          #  Decrease it by the second input x
*         #  Multiply it by the implicit a(N-1)
N<       #  Push N-1
₂*     #  Multiply it by a(N-2)
-    #  Decrease the a(N-1)*(2N-1-x) by this (N-1)*a(N-2)
N/  #  And divide it by N: (a(N-1)*(2N-1-x)-(N-1)*a(N-2))/N


Answered by Kevin Cruijssen on October 27, 2021

Charcoal, 29 bytes

⊞υ¹ＦＮ⊞υ×⌈υＬυＩ↨Ｅυ∕⌈υ×ιＸ§⮌υκ²±Ｎ


Try it online! Link is to verbose version of code. Uses a slightly modified version of the summation given in the question. Explanation:

⊞υ¹ＦＮ⊞υ×⌈υＬυ


Calculate the factorials from $$0!$$ to $$n!$$.

Ｉ↨Ｅυ∕⌈υ×ιＸ§⮌υκ²±Ｎ


For each index $$i$$ from $$0$$ to $$n$$ calculate $$frac{n!}{i!(n-i)!^2}$$ and then perform base conversion from base $$-x$$ which multiplies each term by $$(-1)^{n-i}x^{n-i}$$ and takes the sum.

If we set $$k=n-i$$ we see that we calculate $$sumlimits_{k=0}^{n}{frac{n!(-1)^k}{(n-k)!k!^2}x^k}=sumlimits_{k=0}^{n}{nchoose k}frac{(-1)^k}{k!}x^k$$ as required.

Answered by Neil on October 27, 2021

Python 3.8 (pre-release), 61 bytes

L=lambda k,x:k<1or[1-x,L(w:=k-1,x)*(k+w-x)-L(k-2,x)*w][k>1]/k


Try it online!

Answered by ovs on October 27, 2021

Python 2, 53 bytes

f=lambda n,x:n<1or((2*n-1-x)*f(n-1,x)-~-n*f(n-2,x))/n


Try it online!

Answered by xnor on October 27, 2021

JavaScript (ES6),  48 42  41 bytes

Expects (x)(n). May output true instead of 1.

x=>g=k=>k<1||((x-k---k)*g(k)+k*g(k-1))/~k


Try it online!

Answered by Arnauld on October 27, 2021

Wolfram Language (Mathematica), 9 bytes

LaguerreL


Try it online!

Answered by ZaMoC on October 27, 2021

Python 3.8 (pre-release), 66 bytes

L=lambda n,x:((2*n-1-x)*L(d:=n-1,x)-d*L(n-2,x))/n if n>1else 1-n*x


Try it online!

Direct implementation of the recursive algorithm, with one interesting part: L(1,x) and L(0,x) can be combined as L(n,x)=1-n*x.

Could save 2 bytes using L=lambda n,x:n>1and((2*n-1-x)*L(d:=n-1,x)-d*L(n-2,x))/n or 1-n*x, but L(n)` is not necessarily zero.

Answered by fireflame241 on October 27, 2021