[100] 010 Chapter 10 Polynomial, Tutorial by www.msharpmath.com
[100] 010 revised on 2012.12.03 / cemmathThe Simple is the Best
Chapter 10 Polynomial
10-1 Class ‘poly’
10-2 Operations for Polynomials
10-3 Evaluation of Polynomials
10-4 Subscripts
10-5 Piecewise Continuous Polynomials
10-6 Member Functions of Polynomial
10-7 Sequence and Series
10-8 Special Polynomials
10-9 Coefficient Functions
10-10 Summary
■ class functions.
■ ascending-order polynomial.
■ descending-order polynomial.
■ 'poly' by known roots.
■ Newton polynomial.
■ unary operations.
■ binary operations
■ compound assignments.
■ compound operator for deconvolution.
■ piecewise continuous polynomials.
■ ‘double’-returning.
■ ‘poly’-returning.
■ ‘matrix’-returning.
■ special polynomials.
Polynomials are one of the most important data types in science and engineering. To accommodate this, Cemmath employs class ‘poly’ to handle a number of valuable operations of polynomials. Due to a structure similar to matrix, conversion from matrix to polynomial and vice versa can be easily performed just by assignment operation.
Section 10-1 Class ‘poly’
In order to treat polynomials efficiently, class ‘poly’ is innovated in Cemmath. This contains a large number of member functions to handle various operations of polynomials.
■ Syntax of poly. In Cemmath, polynomials are expressed by default in ascending order of power, i.e.
px=a0+a1x+a2x2+a3x3+ …+anxn+ …
For practical reason, the maximum degree of polynomials is set to be 100, and a serious error will be encountered when this limit is exceeded.
The class ‘poly’ is used to create a polynomial, and the syntax is
poly(a0,a1,a2, …, an); // list of coefficients
poly( [a0,a1,a2,…,an] ); // a single matrix
poly( A ); // A is a matrix
[a0,a1,a2,…,an] .apoly
A .apoly
The zero polynomial is created simply by
poly() or .[]
#> poly();
ans = poly( 0 )
= 0
#> poly(1,2,3,4,5);
ans = poly( 1 2 3 4 5 )
= 5x^4 +4x^3 +3x^2 +2x +1
#> poly( [ 1,2,3,4,5 ] );
ans = poly( 1 2 3 4 5 )
= 5x^4 +4x^3 +3x^2 +2x +1
#> [ 1,2,3,4,5 ].apoly;
ans = poly( 1 2 3 4 5 )
= 5x^4 +4x^3 +3x^2 +2x +1
■ Descending-order polynomial. When a descending-order in power is preferred, either of the followings
[an, …,a2,a1,a0].dpoly; // general order, n >= 4
.[an, …,a2,a1,a0]; // general order, n >= 4
.[ a ] // constant polynomial
.[ a,b ] // 1st-order poly, line y = ax+b
.[ a,b,c ] // 2nd-order poly, parabola y = ax^2 + bx + c
.[ a,b,c,d ] // 3rd-order poly, cubic y = ax^3 + bx^2 + cx +d
poly(an, …,a2,a1,a0).rev;
poly(A).rev;
can be used ('rev' is an abbreviation of reverse).
%> descending order by Matrix
#> .[] ; // zero polynomial not a 0 x 0 matrix
ans = poly( 0 )
= 0
#> [1,2] .dpoly; // 1st-order polynomial, line
y = x + 2, root = -2
#> .[1,2] ; // a leading dot is equal to [].dpoly
y = x + 2, root = -2
#> .[1,2,3] ; // 2nd-order poly, parabola
y = x^2 + 2x + 3
roots
[ -1 - i 1.414 ]
[ -1 + i 1.414 ]
vertex at ( -1, 2 )
y_minimum = 2
y-axis intersect ( 0, 3 )
symmetry at x = -1
directrix at y = 1.75
focus at ( -1, 2.25 )
#> .[1,2,3,4] ; // 3rd-order poly, cubic
y = x^3 + 2x^2 + 3x + 4
roots
[ -1.651 ]
[ -0.1747 - i 1.547 ]
[ -0.1747 + i 1.547 ]
no maxima, no minima
inflection pt. ( -0.666667, 2.59259 )
y-axis intersect ( 0, 4 )
#> .[1,2,3,4,5] ; // no information if order >= 4
ans = poly( 5 4 3 2 1 )
= x^4 +2x^3 +3x^2 +4x +5
#> (1:10).dpoly; // a leading dot is only for .[] grammar
ans = poly( 10 9 8 7 6 5 4 3 2 1 )
= x^9 +2x^8 +3x^7 +4x^6 +5x^5
+6x^4 +7x^3 +8x^2 +9x +10
■ 'poly' by known roots. When roots of a polynomial are given, a polynomial can be obtained by
[ x1,x2, ... , xn ] .roots // (x-x1)(x-x2) .. (x-xn),
.roots(x1,x2, ... , xn) // (x-x1)(x-x2) .. (x-xn)
.roots( A ) // A is a matrix including complex numbers
#> [ 0,1,2,3 ] .roots ;
ans = poly( 0 -6 11 -6 1 )
= x^4 -6x^3 +11x^2 -6x
#> A = [ 1+1i; 1-1i ];
A =
[ 1 + i 1 ]
[ 1 - i 1 ]
#> .roots(A);
ans = poly( 2 -2 1 )
= x^2 -2x +2
#> poly(A) ; // be careful ! This disregards complex parts
ans = poly( 1 1 )
= x +1
■ Newton polynomial. A Newton polynomial is of the form
px=a0+a1x-x1+a2x-x1x-x2
+a3x-x1x-x2(x-x3)+ …
poly(a0,a1,...) .newton( [x0,x1,x2, ...] )
This special polynomial occurs in many applications including interpolation and curve-fitting. For example
%> 1 + 2(x-3) + 3(x-3)(x-4) + 4(x-3)(x-4)(x-5)
#> poly(1,2,3,4) .newton( [3,4,5] );
ans = poly( -209 169 -45 4 )
= 4x^3 -45x^2 +169x -209
#> 1 + 2*[3].roots + 3*[3,4].roots + 4*[3,4,5].roots ;
ans = poly( -209 169 -45 4 )
= 4x^3 -45x^2 +169x -209
#> 1 + 2*.[1,-3] + 3*.[1,-3]*.[1,-4] + 4*.[1,-3]*.[1,-4]*.[1,-5] ;
ans = poly( -209 169 -45 4 )
= 4x^3 -45x^2 +169x -209
■ Constant polynomial. Assignment by a constant results in a constant polynomial.
#> p = (1:5).dpoly;
p = poly( 5 4 3 2 1 )
= x^4 +2x^3 +3x^2 +4x +5
#> p = 2 ;
p = poly( 2 )
= 2
#> (1:5).dpoly - (1:5).dpoly ;
ans = poly( 0 )
= 0
■ Declaring variables as polynomials. Variables can be declared as polynomials by .[] or poly().
%> poly a,b,c,d; // can be alternatively written as
#> a = b = c = d = .[] ;; // declaring variables as polynomials
#> a; b; c; d;
a = poly( 0 )
= 0
b = poly( 0 )
= 0
c = poly( 0 )
= 0
d = poly( 0 )
= 0
■ Plot of polynomial. Since a polynomial is of course a kind of function, member function ‘plot’ draws a polynomial over an interval a≤x≤b by the following syntax
.plot[n=101,g=1](a,b)
where preceding data must a polynomial, either named or unnamed. It is believed here that readers are familiar with the span operation ‘[n=101,g=1](a,b)’. An example is
%> plot of polynomial
#> p = .[ 1, -7, 14, 5 ];
#> p.plot(0,5); // equivalent to plot.x(0,5) ( p(x) )
The resulting plot is shown in Figure 1.
Figure 1 Plot of a polynomial
■ Displaying Polynomials. Consider two polynomials
px=2+3x+4x2+x3
qx=x+2x2+3x3+…+19x19+20x20
and create them by the following Cemmath commands
%> default display
#> p = poly(2,3,4,1);
#> q = poly(0:20);
The results shown on the screen are
p = poly( 2 3 4 1 )
= x^3 +4x^2 +3x +2
q = poly( 0 1 2 3 4 5 6 7 8 9 10
11 12 13 14 15 16 17 18 19 20 )
= 20x^20 +19x^19 +18x^18 +17x^17 +16x^16
+15x^15 +14x^14 +13x^13 +12x^12 +11x^11
+10x^10 +9x^9 +8x^8 +7x^7 +6x^6
+5x^5 +4x^4 +3x^3 +2x^2 +x
The maximum column length displayed on the screen is set by default to be 10 (excluding a constant coefficient a0). The maximum column length, the delimiters (i.e. a pair of parenthesis) and number format for polynomial can be changed, for example, by
%> display
#> poly.format("%5.1f");
#> poly.left (" "); // left delimiter blank
#> poly.right(" "); // right delimiter blank
#> poly.maxcol(5); // maximum column length
#> p = poly(2,3,4,1);
#> q = poly(0:20);
Then, blanks replace both ‘(’ and ‘)’, until delimiters are changed again. Also, the number format as well as the maximum column length will be kept until redefined. Thus, newly displayed results on the screen are
ans = 5
p = 2.0 3.0 4.0 1.0
q = 0.0 1.0 2.0 3.0 4.0 5.0
6.0 7.0 8.0 9.0 10.0
11.0 12.0 13.0 14.0 15.0
16.0 17.0 18.0 19.0 20.0
Section 10-2 Operations for Polynomials
To explain the syntax for polynomial operations, representative polynomials p(x), q(x) and r(x) are used.
■ Unary Operations. The syntax for unary operations of polynomials is explained
|p| // termwise absolute
-p // unary minus
p' // derivative, d(p(x))/dx
p~ // integration
p.' = p` // p(x)-p(x-1), finite-difference, p.' = p.diff = p`
p.~ // p(1)+p(2)+...+p(n), cumulative sum p.~ = p.cumsum
The derivative (denoted by prime) and integration (denoted by tilde) can be repeated depending on the order. Several examples can be seen below.
%> unary operations
#> p = .[1,-4,-3,2]; // polynomial not a matrix due to a leading dot
p = poly( 2 -3 -4 1 )
= x^3 -4x^2 -3x +2
#> |p|; ans = poly( 2 3 4 1 )
#> -p; ans = poly( -2 3 4 -1 )
#> p' ; ans = poly( -3 -8 3 )
#> p'' ; ans = poly( -8 6 )
#> p''' ; ans = poly( 6 )
#> p~ ; ans = poly( 0 2 -1.5 -1.33333 0.25 )
#> p~~ ; ans = poly( 0 0 1 -0.5 -0.333333 0.05 )
#> p~~~ ; ans = poly( 0 0 0 0.333333 -0.125 -0.0666667 0.00833333 )
#> .[1,0,0,0] ' ; // dp(x)/dx = d(x^3)/dx
ans = poly( 0 0 3 )
= 3x^2
#> .[ 1,0,0,0 ].' ; // p(x)-p(x-1) = x^3 - (x-1)^3
ans = poly( 1 -3 3 )
= 3x^2 -3x +1
#> .[1,0,0,0] ~ ; // int_0^x x^3 dx, integration constant is zero
ans = poly( 0 0 0 0 0.25 )
= 0.25x^4
#> .[ 1,0,0,0 ].~ ; // 1^3 + 2^3 + 3^3 + .. + x^3 = [1/2(x)(x+1)]^2
ans = poly( 0 0 0.25 0.5 0.25 )
= 0.25x^4 +0.5x^3 +0.25x^2
It is noteworthy that the integration constants are always set to be zero.
■ Binary Operations with Scalars. Polynomials can operate with ‘double’ through binary operations
p + d, d + p // addition
p - d, d - p // subtraction
p * d, d * p // multiplication
p / d, d \ p // right and left divisions
p ^ n // right power only
p .^ d, d .^ p // element-by-element power
p %% d // synthetic division
#> p = .[ 1,1 ]; // x + 1
p = poly( 1 1 )
= x +1
#> for.i(1,6) p^i ; // Pascal's triangle
ans = poly( 1 1 )
= x +1
ans = poly( 1 2 1 )
= x^2 +2x +1
ans = poly( 1 3 3 1 )
= x^3 +3x^2 +3x +1
ans = poly( 1 4 6 4 1 )
= x^4 +4x^3 +6x^2 +4x +1
ans = poly( 1 5 10 10 5 1 )
= x^5 +5x^4 +10x^3 +10x^2 +5x +1
■ Synthetic Division. A synthetic division is defined as
px=a0+a1x+a2x2+a3x3+ …+anxn
px %% c=b0+b1(x-c)+b2(x-c)2+b3(x-c)3+ …+bn(x-c)n
the coefficients are determined by comparing two expressions such that
px=a0+a1x+a2x2+ …+anxn
=b0+b1(x-c)+b2(x-c)2+…+bn(x-c)n
For example,
x3+4x2+3x+2=x+33-5x+32+6x+3+2
can be confirmed by
%> synthetic division
#> p = .[1,4,3,2];
p = poly( 2 3 4 1 )
= x^3 +4x^2 +3x +2
#> p %% -3;
ans = poly( 2 6 -5 1 )
= x^3 -5x^2 +6x +2
#> ans ( .[1,3] ); // ans = p %% -3 ; from the above command
ans = poly( 2 3 4 1 )
= x^3 +4x^2 +3x +2
■ Binary Operations. Binary operations between polynomials are also presented
p + q // addition
p - q // subtraction
p * q // multiplication
p ¥ q // quotient (left division)
p / q // quotient (right division)
p % q // remainder
p %% q // synthetic division
p == q // conditional equality
p != q // conditional inequality
#> p = .[3,4,5]; // 3x^2 +4x + 5
p = poly( 5 4 3 )
= 3x^2 +4x +5
#> q = .[1,2]; // x +2
q = poly( 2 1 )
= x +2
#> p + q ; // (3x^2 +4x +5) + (x +2)
ans = poly( 7 5 3 )
= 3x^2 +5x +7
#> p - q ; // (3x^2 +4x +5) - (x +2)
ans = poly( 3 3 3 )
= 3x^2 +3x +3
#> p * q ; // (3x^2 +4x +5) * (x +2)
ans = poly( 10 13 10 3 )
= 3x^3 +10x^2 +13x +10
#> p / q ; // (3x^2 +4x +5) = (x +2)(3x-2) + 9
ans = poly( -2 3 )
= 3x -2
#> q \ p ; // (3x^2 +4x +5) = (x +2)(3x-2) + 9
ans = poly( -2 3 )
= 3x -2
#> p % q ; // (3x^2 +4x +5) = (x +2)(3x-2) + 9
ans = poly( 9 )
= 9
#> p == q ; // (3x^2 +4x +5) == (x+2)
ans = 0
#> p != q ; // (3x^2 +4x +5) != (x+2)
ans = 1
In Cemmath, deconvolution of polynomials a=bq+r can be found efficiently by using a compound operator ‘/%’, similar to the class ‘double’. The syntax is
(q,r) = a /% b;
This can be confirmed by
%> deconvolution
#> a = .[1,4,3,2]; b = .[1,1,3];
a = poly( 2 3 4 1 )
= x^3 +4x^2 +3x +2
b = poly( 3 1 1 )
= x^2 +x +3
#> (q,r) = a /% b; // a = bq + r, x^3+4x^2+3x+2 = (x^2+x+3)(x+3) -3x-7
q = poly( 3 1 )
= x +3
r = poly( -7 -3 )
= -3x -7
The synthetic division by higher degrees can be also treated by an operator ‘%%’. For example
x5+8x4+6x3+3x2-10x+2x2+x+12=-6x+4x2+x+12+-9x-8x2+x+1+(x+6)
This has a nature similar to synthetic division
x5+8x4+6x3+3x2-10x+2
=x+6x2+x+12+-9x-8x2+x+1+(-6x+4)
An operator ‘%%’ for this case can be applied
%> synthetic division (degree >= 2)
#> p = .[ 1,8,6,3,-10,2 ];
p = poly( 2 -10 3 6 8 1 )
= x^5 +8x^4 +6x^3 +3x^2 -10x +2
#> q = .[1,1,1]; // equal to q = poly(1,1,1);
q = poly( 1 1 1 )
= x^2 +x +1
#> A = p %% q; // matrix for ascending poly
A =
[ 4 -6 ]
[ -8 -9 ]