laguerre.pas

maths lib with source

unit laguerre;
interface
uses Math, Ap, Sysutils;

function LaguerreCalculate(const n : Integer; const x : Double):Double;
function LaguerreSum(const C : TReal1DArray;
     const n : Integer;
     const x : Double):Double;
procedure LaguerreCoefficients(const N : Integer; var C : TReal1DArray);

implementation

(*************************************************************************
Calculation of the value of the Laguerre polynomial.

Parameters:
    n   -   degree, n>=0
    x   -   argument

Result:
    the value of the Laguerre polynomial Ln at x
*************************************************************************)
function LaguerreCalculate(const n : Integer; const x : Double):Double;
var
    a : Double;
    b : Double;
    i : Double;
begin
    Result := 1;
    a := 1;
    b := 1-x;
    if n=1 then
    begin
        Result := b;
    end;
    i := 2;
    while i<=n do
    begin
        Result := ((2*i-1-x)*b-(i-1)*a)/i;
        a := b;
        b := Result;
        i := i+1;
    end;
end;


(*************************************************************************
Summation of Laguerre polynomials using Clenshaw抯 recurrence formula.

This routine calculates c[0]*L0(x) + c[1]*L1(x) + ... + c[N]*LN(x)

Parameters:
    n   -   degree, n>=0
    x   -   argument

Result:
    the value of the Laguerre polynomial at x
*************************************************************************)
function LaguerreSum(const C : TReal1DArray;
     const n : Integer;
     const x : Double):Double;
var
    b1 : Double;
    b2 : Double;
    i : Integer;
begin
    b1 := 0;
    b2 := 0;
    i:=n;
    while i>=0 do
    begin
        Result := (2*i+1-x)*b1/(i+1)-(i+1)*b2/(i+2)+C[i];
        b2 := b1;
        b1 := Result;
        Dec(i);
    end;
end;


(*************************************************************************
Representation of Ln as C[0] + C[1]*X + ... + C[N]*X^N

Input parameters:
    N   -   polynomial degree, n>=0

Output parameters:
    C   -   coefficients
*************************************************************************)
procedure LaguerreCoefficients(const N : Integer; var C : TReal1DArray);
var
    i : Integer;
begin
    SetLength(C, N+1);
    C[0] := 1;
    i:=0;
    while i<=n-1 do
    begin
        C[i+1] := -C[i]*(n-i)/(i+1)/(i+1);
        Inc(i);
    end;
end;


end.