#!/usr/bin/perl
# Copyright (c) 2008-2017 Martin Becker. All rights reserved.
# This package is free software; you can redistribute it and/or modify it
# under the same terms as Perl itself.
# Math::Polynomial usage example: calculating Legendre polynomials.
#
# Legendre polynomials are a special and well-known (to scientists,
# at least) kind of orthogonal polynomial series. This script generates
# the first few of them using a recursion formula and shows their
# orthogonality feature by calculating the related inner product of
# any two of them, yielding zero whenever two different polynomials
# are multiplied, and a positive value if a polynomial is multiplied
# by itself.
use strict;
use warnings;
use Math::Polynomial 1.000;
use Math::AnyNum;
my $max_degree = 5;
sub fmt_num {
my ($n, $d) = $_[0]->nude;
return 1 == $d? "$n": "$n/$d";
}
# adjust some printing options
Math::Polynomial->string_config({
fold_sign => 1,
prefix => q{},
suffix => q{},
convert_coeff => \&fmt_num,
});
# create p[0] = 1 and p[1] = x
# using arbitrary precision rational coefficients
my $one = Math::AnyNum->new('1');
my $p0 = Math::Polynomial->new($one);
my $p1 = $p0 << 1;
my @p = ($p0, $p1);
# recursion: (n+1)*p[n+1] = (2n+1)*x*p[n] - n*p[n-1]
foreach my $n (1..$max_degree-1) {
$p[$n+1] = ($p[$n] * $p1 * ($n+$n+1) - $p[$n-1] * $n) / ($n + 1);
}
# print polynomials
foreach my $n (0..$#p) {
print "P_$n = $p[$n]\n";
}
# demonstrate orthogonality
foreach my $n (0..$#p) {
foreach my $m (0..$n) {
my $s = ($p[$n] * $p[$m])->definite_integral(-$one, $one);
print " = ", fmt_num($s), "\n";
}
}
__END__