File Coverage

blib/lib/Math/Polynomial/Chebyshev2.pm
Criterion Covered Total %
statement 52 52 100.0
branch 12 14 85.7
condition n/a
subroutine 8 8 100.0
pod 2 2 100.0
total 74 76 97.3


line stmt bran cond sub pod time code
1             # -*- coding: utf-8-unix -*-
2              
3             package Math::Polynomial::Chebyshev2;
4              
5 2     2   57709 use 5.008;
  2         11  
6 2     2   513 use utf8;
  2         14  
  2         9  
7 2     2   56 use strict;
  2         4  
  2         32  
8 2     2   9 use warnings;
  2         3  
  2         53  
9              
10 2     2   10 use Carp qw< croak >;
  2         2  
  2         71  
11 2     2   605 use Math::Polynomial;
  2         8031  
  2         779  
12              
13             our $VERSION = '0.02';
14             our @ISA = qw< Math::Polynomial >;
15              
16             =pod
17              
18             =encoding UTF-8
19              
20             =head1 NAME
21              
22             Math::Polynomial::Chebyshev2 - Chebyshev polynomials of the second kind
23              
24             =head1 SYNOPSIS
25              
26             use Math::Polynomial::Chebyshev2;
27              
28             # create a Chebyshev polynomial of the second kind of order 7
29             my $p = Math::Polynomial::Chebyshev2 -> chebyshev2(7);
30              
31             # get the location of all roots
32             my @xn = $p -> roots();
33              
34             # use higher accuracy
35             use Math::BigFloat;
36             Math::BigFloat -> accuracy(60);
37             my $n_mbf = Math::BigFloat -> new(7);
38             my $p_mbf = Math::Polynomial::Chebyshev2 -> chebyshev2($n_mbf);
39              
40             =head1 DESCRIPTION
41              
42             This package extends Math::Polynomial, so each instance polynomial created by
43             this module is a subclass of Math::Polynomial.
44              
45             The Chebyshev polynomials of the second kind are orthogonal with respect to the
46             weight function sqrt(1-x^2).
47              
48             The first Chebyshev polynomials of the second kind are
49              
50             U₀(x) = 1
51             U₁(x) = 2 x
52             U₂(x) = 4 x^2 - 1
53             U₃(x) = 8 x^3 - 4 x
54             U₄(x) = 16 x^4 - 12 x^2 + 1
55             U₅(x) = 32 x^5 - 32 x^3 + 6 x
56             U₆(x) = 64 x^6 - 80 x^4 + 24 x^2 - 1
57             U₇(x) = 128 x^7 - 192 x^5 + 80 x^3 - 8 x
58             U₈(x) = 256 x^8 - 448 x^6 + 240 x^4 - 40 x^2 + 1
59             U₉(x) = 512 x^9 - 1024 x^7 + 672 x^5 - 160 x^3 + 10 x
60              
61             =head1 CLASS METHODS
62              
63             =head2 Constructors
64              
65             =over 4
66              
67             =item I
68              
69             Cchebyshev2($n)> creates a new polynomial of
70             order C<$n>, where C<$n> is a non-negative integer.
71              
72             # create a Chebyshev polynomial of the second kind of order 7
73             $p = Math::Polynomial::Chebyshev2 -> chebyshev2(7);
74              
75             # do the same, but with Math::BigFloat coefficients
76             use Math::BigFloat;
77             $n = Math::BigFloat -> new(7);
78             $p = Math::Polynomial::Chebyshev2 -> chebyshev2($n);
79              
80             # do the same, but with Math::Complex coefficients
81             use Math::Complex;
82             $n = Math::Complex -> new(7);
83             $p = Math::Polynomial::Chebyshev2 -> chebyshev2($n);
84              
85             =cut
86              
87             sub chebyshev2 {
88 10     10 1 34993 my $class = shift;
89 10         15 my $n = shift;
90              
91 10 50       26 croak "order must be an integer" unless $n == int $n;
92              
93 10         14 my $zero = $n - $n;
94 10         17 my $one = $n ** 0;
95 10         14 my $two = $one + $one;
96              
97 10         16 my $c = [];
98 10 100       26 if ($n == 0) {
    100          
99 1         2 $c = [ $one ];
100             } elsif ($n == 1) {
101 1         3 $c = [ $zero, $two ];
102             } else {
103 8         13 my $a = [ $one ];
104 8         13 my $b = [ $zero, $two ];
105              
106 8         18 for (my $i = 2 ; $i <= $n ; ++$i) {
107 36         48 $c = [ $zero, map { $two * $_ } @$b ];
  156         209  
108              
109 36         67 for (my $i = 0 ; $i <= $#$a ; ++$i) {
110 120         191 $c -> [$i] -= $a -> [$i];
111             }
112              
113 36         47 $a = $b;
114 36         63 $b = $c;
115             }
116             }
117              
118 10         32 return $class -> new(@$c);
119             }
120              
121             =pod
122              
123             =back
124              
125             =head2 Property Accessors
126              
127             =over 4
128              
129             =item I
130              
131             C<$p-Eroots()> returns the location of all root of C<$p>. All roots are
132             located in the open interval (-1,1).
133              
134             # get the roots of a polynomial
135             @x = $p -> roots();
136              
137             =cut
138              
139             sub roots {
140 10     10 1 50960 my $self = shift;
141 10 50       26 croak 'array context required' unless wantarray;
142              
143 10         25 my $n = $self -> degree();
144              
145             # Quick exit for the simple case N = 0.
146              
147 10 100       53 return () if $n == 0;
148              
149             # Quick exit for the simple case N = 1.
150              
151 9         29 my $zero = $self -> coeff_zero();
152 9 100       35 return $zero if $n == 1;
153              
154             # The general case, when N > 0.
155              
156 8         17 my $one = $self -> coeff_one();
157 8         36 my $pi = atan2 $zero, -$one;
158 8         17 my $c = $pi / ($n + 1);
159              
160 8         11 my @x = ();
161              
162             # First compute all roots in the open interval (0,1).
163              
164 8         21 @x = map { cos($c * $_) } 1 .. int($n / 2);
  20         60  
165              
166             # Now create an array with all extremas on the closed interval [-1,1].
167              
168 8 100       15 @x = (map({ -$_ } @x),
  20         45  
169             ($n % 2 ? $zero : ()),
170             reverse(@x));
171              
172 8         26 return @x;
173             }
174              
175             =pod
176              
177             =back
178              
179             =head1 BUGS
180              
181             Please report any bugs through the web interface at
182             L
183             (requires login). We will be notified, and then you'll automatically be
184             notified of progress on your bug as I make changes.
185              
186             =head1 SUPPORT
187              
188             You can find documentation for this module with the perldoc command.
189              
190             perldoc Math::Polynomial::Chebyshev
191              
192             You can also look for information at:
193              
194             =over 4
195              
196             =item * GitHub Source Repository
197              
198             L
199              
200             =item * RT: CPAN's request tracker
201              
202             L
203              
204             =item * CPAN Ratings
205              
206             L
207              
208             =item * MetaCPAN
209              
210             L
211              
212             =item * CPAN Testers Matrix
213              
214             L
215              
216             =back
217              
218             =head1 SEE ALSO
219              
220             =over
221              
222             =item *
223              
224             The Perl module L.
225              
226             =item *
227              
228             The Wikipedia page L.
229              
230             =back
231              
232             =head1 LICENSE AND COPYRIGHT
233              
234             Copyright (c) 2020 Peter John Acklam.
235              
236             This program is free software; you may redistribute it and/or modify it under
237             the same terms as Perl itself.
238              
239             =head1 AUTHOR
240              
241             Peter John Acklam Epjacklam (at) gmail.comE.
242              
243             =cut
244              
245             1;