``````#!/usr/bin/perl

# 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";
}

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 "<P_\$n, P_\$m> = ", fmt_num(\$s), "\n";
}
}

__END__
``````