Jacob Salomon
and 1 contributors

NAME

Math::Yapp - Perl extension for working with Polynomials. Yes, I know there are *many!* Polynomial packages. And like them, I started it for (geeky) fun, then got obsessed with it as a learning experience. Enjoy!

AUTHOR

Jacob Salomon, jakesalomon@yahoo.com

COPYRIGHT AND LICENSE

Copyright (C) 2013, 2015 by Jacob Salomon

SYNOPSIS

use Math::Complex;

use Math::Yapp;

The "Math::Complex" is optional but a very strongly recommended, since solving a polynomial often yields complex numbers

Constructors

$yp = Math::Yapp->new(); # Degenerate polynomial - no terms

While the new() method is certainly available, your code will look cleaner if you use this form:

$yp = Yapp(parameter(s)); # Examples follow:

Real and complex coefficiencts

$yp = Yapp(1, -2.178, 3.14, -cplx(1,-3), 5); # Notice: Ascending in degree

Yields: 5x^4 +(-1+3i)x^3 +3.14x^2 -2.178x +1

Same as above, only passed as a reference to an array:

my @coef_list = (1, -2.178, 3.14, -cplx(1,-3), 5);

$yp = Yapp(\@coef_list);

Generated from a string

You can also generate the polynomial from a string: Note that the sign MUST precede each term with no intervening space; a space separating the sign from the term may mess up the matching pattern. Yes, the + sign is optional for the first term in the string.

$yp = Yapp("5x^4 +(-1+3i)x^3 +3.14x^2 -2.178x +1");

Note that the order of the terms is unimportant; they get sorted out automatically as I piece together the structure by exponent.

Copy Constructor

$yp1 = Yapp($yp); # Clones $yp to an identical copy w a different reference

Constructors by interpolation

$ypl = Yapp_interpolate(\@xvals, \@yvals); # Perform Lagrange interpolation

$yph = Yapp_interpolate(\@xvals, \@yvals, \@ypvals); # Hermite interpolation

Note that the interpolating constructors require array references; otherwise you run into the elementary error of array-mashing. The Lagrange form, just X and Y values, generates a polynomial that colocates at each of the points indicated in the X-Y pairs represented by the arrays. For The Hermite version the third array [reference] if for the desired derivative at each point.

Notes on Hermite Interpolation:

1. At this time, the author waives Hermite interpolation for more than the first derivative. :-)
2. The author has succssfully tested Hermite interpolation for up to 6 points in a 64-bit CygWin environment; at 7 points, the calculations seem to have run into some instability, possibly caused by rounding errors.

Limitations on interpolation

Ordinary Lagrange interpolation has been tested in 32- and 64-bit Cygwin environments for up to 7 points with accuracy to 11 decimal places. On the other hand, with Hermite interpolation with 1 derivative, testing began to fail at 5 points in the 32-bit and at 7 points in the 64-bit environment. And I'm not talking about missing by 6th decimal place; the errors became quite gross at the last points in the arrays. Hence, in the 64-bit environment, I limited the Hermite interpolation test to 6 points. Disappointing. (I really needed a Math::BigComplex module based on Math::BigFloat and demand 60-digit accuracy for some of this stuff!)

My plans for higher-derivative interpolation are hold for this. (Also for when I grok the "finite differences" algorithms.)

Constructing a Yapp polynomial from its roots

$yp = Yapp_by_roots(\@root_list); # Pass reference to an array of roots

$yp = Yapp_by_roots(1, 2, -4, cplx(3,4), cplx(3,-1)); # Pass a complete array to constructor

Arithmetic of Yapp Polynomials

Unary Operations

$yp2 = !$yp; # Change the signs of all coefficients

$yp2 = ~$yp; # Conjugate complex coefficients

Addition and Subtraction

$yp += 13; # Add a real value to the constant term

$yp += cplx(2, -5); # Add a complex number to the constant term

$yp += $yp3; # Add another polynomial to this one

$yp = $yp1 + $yp2; # Add two polynomials, term-by-term

Subtracting polynomials: Behaves pretty much like the adds so we are not inclucing all possible examples

$yp -= $yp3; # Subtract $yp3 from $yp in place

$yp = $yp1 - $yp2; # Subtract two polynomials, term-by-term

Multiplication and division:

$yp *= 42; # Multiply each coefficient by the same number

$yp = $yp1 * 42; # Multiply as above but return a new polynomial

$yp = 42 * $yp1; # (Same idea as above)

$yp *= $yp2; # In-place multiplication of a Yapp polynomial by another

$yp = $yp1 * $yp2; # Same as above, but not in-place

$yp /= 10; # Divide all coefficients by a number

Division by a polynomial is not defined in this package, although when I evaluate a polynomials at, say, x = 3, it is equivalent to dividing by the small polynomial (x - 3). Hence, for Yapp_eval (described later), you have the option of getting the quotient besides the evaluation

Inner Product

$fnum = $yapp1 . $yapp2 # Notice the dot-operator in place of *.

The inner product depends on the integral method, described later. This is the inner product for Legendre orhogonality, that is: Integral[-1,+1]($Y1 * $Y2)dx for real coefficients.

For complex coefficients, it is: Integral[-1,+1]($Y1 * ~$Y2)dz

Other inner product alternatives, like Laguerre, Hermite and others, are on the back burner but in separate modules yet to be written (at this time) like Math::Yapp::Laguerre and other, that inherit all basic methods but differ in the the inner product operator/method.

Documented methods and functions:

Format a polynomial for printing

Ysprint() formats a Yapp object into a string, suitable for printing: Example: printf("My yapp is: %s\n", $yp1->Ysprint());

By default, Ysprint formats the polynmial as follows:

  • Starting from the high-order exponent, working its way down

  • All coefficients are displayed with 3 decimal places

  • Zero coefficients are skipped.

This can be controlled by the following functions, which affect module-global variables:

  • Yapp_start_high(0); # Setting FALSE: Start from low-degree terms

  • Yapp_decimals(2); # Set number of decimal places to display

  • Yapp_print0(1); # Ysprint shall display 0-value coefficients

In all three cases, the newly setvalue is returned. Oh, and if you call it without a parameter, it will just return the currently set value.

You can also override the default precision by supplying a parameter to Ysprint;

printf("My yapp is: %s\n", $yp1->Ysprint(6));

The above example will print the coefficients with 6 decimal places, ignoring the default set by Yapp_decimals(),

Setting and retrieving the display variable

By default, when formatting the polynomial for display, the "variable" will be the letter "X". You can change this to any other string using the variable() method:

$yp->variable("Xy"); # Sets the "variable" to "Xy"

If you just want to know what variable is set for display, call variable() with no parameter:

my $varname = $yp->variable(); # Returns that variable string

Query individual coefficients

my $nterm = $yp->coefficient(3); # Retrieve the X^3 coefficient

Retrieve the nth degree coefficient:

my $high_expon = $yp->degree(); # Returns the degree of the polynomial

Evaluating a polynomial at specified values

my $realvar = $yp->Yapp_eval($areal);

The above plugs the parameter into the polynomial and return the value of the polynomial at that point. This works identically when plugging in a complex number but you should be prepared to have a complex number returned:

my $cplxvar = $yp->Yapp_eval($acplx);

(Of course, if the polynomial has complex coefficients and you plug in a real number, you are obviously at high risk of getting back a complex number.)

When you evaluate a polynomial at a certain value, say 3, it is like dividing the polynomial by (X - 3); the returned value is the "remainder". (You surely learned this in high school.) Now, what of the quotient? Well, just ask and you shall receive:

my ($eval, $quotient) = $yp->Yapp_eval($areal);

Of course, this quotient is of reduced degree.

Reduce the roots of a polynomial by a number

my $ypr = $yp->Yapp_reduce(3);

This applies Horner's method of successive synthetic division to the polyomial.

Negate the roots of a polynomial

$nr_yapp = $yp->Yapp_negate_roots();

This produces a polynomial whose roots are the negatives of the roots of the original polynomial.

Derivatives and Integrals

my $dyp = $yp->Yapp_derivative(n);

This returns [a ref to] another polynomial that is the nth derivative of this one. A couple of little notes on this:

  • If n is 0, it just returns the original Yapp reference and you've accomplished very little.

  • If n is omitted, it assumes a default value of 1:

my $i_ref = $yp->Yapp_integral();

This returns a reference to a polynomial whose derivative is the given Yapp. We insert a 0 for the arbitrary constant.

my $i_val = $yp->Yapp_integral($low, $high);

This calculates the value of the definite integral of the given Yapp in the given interval.

Solution Set for Polynomials

my @solutions = $yp->Yapp_solve(); # Solve for all real and complex roots

Inner Product space of Yapp polynomials

The inner-product operator

For these simplest of orthogonal polynomials, the inner product is the integral of ($Y1 * ~($Y2)) over the interval [-1, 1].

$cplx_num = $Y1 . $Y2; # Inner product of two Yapps

$cplx_num = $Y1->Yapp_innerProd($Y2); # Alternative form of above

$cplx_num = Yapp_innerProd($Y1, $Y2); # Another alternative form

But it's really intended to be invoked by the first form, with the "dot" operator.

The norm function

The simplest norm function in any inner-product space is simply the square root of the inner product of a vector with itself. This will always return a real number, even with complex coefficients.

my $norm = $Yp->Yapp_Norm(); # Returns the (Legendre) norm of a polynomial

Orthognality:

my $perp = $Y1->Yapp_Orthogonal($Y2); # True/False: Is $Y2 orthogonal to $Y1? my $perp = Yapp_Orthogonal($Y1, $Y2); # Either form is just fine

Note: With complex coefficients, $Y1 . $Y2 and $Y2 . $Y1 are conjugates.

Missing functions:

  • A Gram-Schmidt orthogonalization process

  • Corrolary to above: A Gram-Schmidt orthogonalization process

  • Generation of a sequence of the first N orthogonal Yapp polynomials using the recursion relation common to all classes orthogonal polynomials.

I've rushed (ahem ;-) this out before getting those done because I found a nasty bug in the conjugation operator function.

DESCRIPTION

Man, if that synopsis don't say it all, what can I possibly add? :-)

OK, as mentioned above, this is a fun project. The plan, not necessarily all implemented at the first release, is to provide many kinds of operations on the polynomials that intimidated us in high school and even college (if you took PDE or Numerical Analysis). The usual high-school functions are ordinary arithmetic on these algebraic expressions, as well as plugging a value in there to get the result. (Horner's synthetic division saves a lot of work.)

Addendum: A Note on Inner products and orthogonal polynomials;

At this stage, the plan is to provide the inner product algorithms for various classes of inner-product spaces but in separate modules, for example: Math::Yapp:Tchebycheff or Math::Yapp::Laguerre. These will use the overloaded "dot" operator.

EXPORT

This package export the following functions by default. (They are few enough to be an unlikey source of name-space pollution, IMHO)

  • Yapp(): The constructor that is NOT a method, so you don't have to type Math::Yapp::new().

  • Yapp_interpolate(): The constructor by Lagrange and Hermite interpolation

  • Yapp_by_roots(): Construct a polynomial by its roots

  • Yapp_decimals(): Sets or retrives a global setting for how many decimal places to display with each %f-formatted variable.

  • Yapp_print0(): By default, printing a Yapp will skip any term whose coefficient is 0. This exported function sets an internal global flag to either display 0-coefficient terms (1) or to skip them (0)

  • Yapp_start_high(): By default, printing a Yapp will start from the highest coefficient, the way we are accustomed to writing a polynomial. For some testing purposes and, I imagine, some other applications, it may be neater to start printing from the low-degree coefficients. This functions sets an intenal global flag to print high-to-low (default: 1) or low-to-high (0).

Bugs

Note that the current release of Math::Yapp uses the default floating-point library of its host system. It was developed in a Cygwin environment running under Windows-7. I discovered some limitations to the 64-bit FP operations when solving polynomials of degree higher thatn 8 or using Hermite interpolation of more than 6 points. I have researched Math::MPC a bit and hope to use that in a future release of this module. However, I encountered some errors when trying to compile the required MPC, MPFR and GMP C libraries. So that plan is going on the back burner.

There is some sloppy code in method Ysprint() which produces the correct output but I would really rather figure out why I needed to add said sloppy code. The bug this covers up is that the first degree term should display without exponent, not as X^1. Similarly, the constant term should display with neither variable nor exponent, rather than as X^0. ANd it does, but only due to the afterthought corrections.

It might be argued that my failure to include higher-order derivatives in the Hermite interpolation scheme is a bug. Perhaps by the time I publish the next release of this module I will have an understanding of Newton's Method of Divided Differences. (Sorry, no promises on that account.)

I have been unable to get unary negation and conjugation operators (! and ~) to work in-place. That is: While I can happily set $Y1 = ~ $Y1, I cannot simply say ~$Y1.

The usual disclaimers :-)

This library is free software; you can redistribute it and/or modify it under the same terms as Perl itself, either Perl version 5.14.2 or, at your option, any later version of Perl 5 you may have available.

Acknowledements

Thanks to John Altom, formerly of CCM Consulting Services, for his help in some basic (but arguably forgettable) details of the EXPORT behavior.

An old debt of gratitude to [the presumably late] Professor Stanley Preisler, who taught Numerical Analysis at Polytechnic University in the late '70s. The course was quite over my head but some stuff remained, as you can see.