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!
Jacob Salomon, jakesalomon@yahoo.com
Copyright (C) 2013, 2015 by Jacob Salomon
use Math::Complex;
use Math::Yapp;
The "Math::Complex" is optional but a very strongly recommended, since solving a polynomial often yields complex numbers
$yp = Math::Yapp->new(); # Degenerate polynomial - no terms
$yp = Math::Yapp->new();
While the new() method is certainly available, your code will look cleaner if you use this form:
$yp = Yapp(parameter(s)); # Examples follow:
$yp = Yapp(parameter(s));
$yp = Yapp(1, -2.178, 3.14, -cplx(1,-3), 5); # Notice: Ascending in degree
$yp = Yapp(1, -2.178, 3.14, -cplx(1,-3), 5);
Yields: 5x^4 +(-1+3i)x^3 +3.14x^2 -2.178x +1
my @coef_list = (1, -2.178, 3.14, -cplx(1,-3), 5);
$yp = Yapp(\@coef_list);
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.
$yp1 = Yapp($yp); # Clones $yp to an identical copy w a different reference
$yp1 = Yapp($yp);
$ypl = Yapp_interpolate(\@xvals, \@yvals); # Perform Lagrange interpolation
$ypl = Yapp_interpolate(\@xvals, \@yvals);
$yph = Yapp_interpolate(\@xvals, \@yvals, \@ypvals); # Hermite interpolation
$yph = Yapp_interpolate(\@xvals, \@yvals, \@ypvals);
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:
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.)
$yp = Yapp_by_roots(\@root_list); # Pass reference to an array of roots
$yp = Yapp_by_roots(\@root_list);
$yp = Yapp_by_roots(1, 2, -4, cplx(3,4), cplx(3,-1)); # Pass a complete array to constructor
$yp = Yapp_by_roots(1, 2, -4, cplx(3,4), cplx(3,-1));
$yp2 = !$yp; # Change the signs of all coefficients
$yp2 = !$yp;
$yp2 = ~$yp; # Conjugate complex coefficients
$yp2 = ~$yp;
$yp += 13; # Add a real value to the constant term
$yp += 13;
$yp += cplx(2, -5); # Add a complex number to the constant term
$yp += cplx(2, -5);
$yp += $yp3; # Add another polynomial to this one
$yp += $yp3;
$yp = $yp1 + $yp2; # Add two polynomials, term-by-term
$yp = $yp1 + $yp2;
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 -= $yp3;
$yp = $yp1 - $yp2; # Subtract two polynomials, term-by-term
$yp = $yp1 - $yp2;
$yp *= 42; # Multiply each coefficient by the same number
$yp *= 42;
$yp = $yp1 * 42; # Multiply as above but return a new polynomial
$yp = $yp1 * 42;
$yp = 42 * $yp1; # (Same idea as above)
$yp = 42 * $yp1;
$yp *= $yp2; # In-place multiplication of a Yapp polynomial by another
$yp *= $yp2;
$yp = $yp1 * $yp2; # Same as above, but not in-place
$yp = $yp1 * $yp2;
$yp /= 10; # Divide all coefficients by a number
$yp /= 10;
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
$fnum = $yapp1 . $yapp2 # Notice the dot-operator in place of *.
$fnum = $yapp1 . $yapp2
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.
Ysprint() formats a Yapp object into a string, suitable for printing: Example: printf("My yapp is: %s\n", $yp1->Ysprint());
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(),
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
my $nterm = $yp->coefficient(3); # Retrieve the X^3 coefficient
my $high_expon = $yp->degree(); # Returns the degree of the polynomial
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.
my $ypr = $yp->Yapp_reduce(3);
This applies Horner's method of successive synthetic division to the polyomial.
$nr_yapp = $yp->Yapp_negate_roots();
This produces a polynomial whose roots are the negatives of the roots of the original polynomial.
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.
my @solutions = $yp->Yapp_solve(); # Solve for all real and complex roots
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 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
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.
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.
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.)
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.
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).
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.
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.
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.
To install Math::Yapp, copy and paste the appropriate command in to your terminal.
cpanm
cpanm Math::Yapp
CPAN shell
perl -MCPAN -e shell install Math::Yapp
For more information on module installation, please visit the detailed CPAN module installation guide.