home *** CD-ROM | disk | FTP | other *** search
/ PC Professionell 2004 December / PCpro_2004_12.ISO / files / webserver / tsw / TSW_3.4.0.exe / Apache2 / perl / Polynomial.pm < prev    next >
Encoding:
Perl POD Document  |  2002-09-12  |  23.0 KB  |  868 lines
  1. package Math::Cephes::Polynomial;
  2. use strict;
  3. use vars qw(@EXPORT_OK $VERSION $MAXPOL $FMAXPOL $flag $fflag);
  4. eval {require Math::Complex; import Math::Complex qw(Re Im)};
  5. eval {local $^W=0; require Math::Fraction;};
  6. $MAXPOL = 256;
  7. $flag = 0;
  8. $FMAXPOL = 256;
  9. $fflag = 0;
  10.  
  11. require Exporter;
  12. *import = \&Exporter::import;
  13. @EXPORT_OK = qw(poly);
  14. $VERSION = '0.36';
  15.  
  16. require Math::Cephes;
  17. require Math::Cephes::Fraction;
  18. require Math::Cephes::Complex;
  19.  
  20. sub new {
  21.     my ($caller, $arr) = @_;
  22.     my $refer = ref($caller);
  23.     my $class = $refer || $caller;
  24.     die "Must supply data for the polynomial" 
  25.       unless ($refer or $arr);
  26.     my ($type, $ref, $data, $n);
  27.     if ($refer) {
  28.       ($type, $ref, $n) =
  29.     ($caller->{type}, $caller->{ref}, $caller->{n});
  30.       my $cdata = $caller->{data};
  31.       if (ref($cdata) eq 'ARRAY') {
  32.     $data = [ @$cdata ];
  33.       }
  34.       else {
  35.     my ($f, $s) = ($type eq 'fract') ? ('n', 'd') : ('r', 'i');
  36.     $data = {$f => [ @{$cdata->{$f}} ],
  37.          $s => [ @{$cdata->{$s}} ], 
  38.         };
  39.       }
  40.     }
  41.     else {
  42.       ($type, $ref, $data, $n) = get_data($arr);
  43.     }
  44.     bless { type => $type,
  45.         ref => $ref,
  46.         data => $data,
  47.         n => $n,
  48.     }, $class;
  49. }
  50.  
  51. sub poly {
  52.   return Math::Cephes::Polynomial->new(shift);
  53. }
  54.  
  55. sub coef {
  56.     return $_[0]->{data};
  57. }
  58.  
  59. sub get_data {
  60.     my ($arr, $ref_in) = @_;
  61.     die "Must supply an array reference" unless ref($arr) eq 'ARRAY';
  62.     my $n = scalar @$arr - 1;
  63.     my $ref = ref($arr->[0]);
  64.     die "array data must be of type '$ref_in'"
  65.     if (defined $ref_in and $ref_in ne $ref);
  66.     my ($type, $data);
  67.   SWITCH: {
  68.       not $ref and do {
  69.       $type = 'scalar';
  70.       foreach (@$arr) {
  71.           die 'Found conflicting types in array data'
  72.           if ref($_);
  73.       }
  74.       $data = $arr;
  75.       set_max() unless $flag;
  76.       last SWITCH;
  77.       };
  78.       $ref eq 'Math::Cephes::Complex' and do {
  79.       $type = 'cmplx';
  80.       foreach (@$arr) {
  81.           die 'Found conflicting types in array data'
  82.           unless ref($_) eq $ref;
  83.           die "array data must be of type '$ref_in'"
  84.           if (defined $ref_in and $ref_in ne $ref);
  85.           push @{$data->{r}}, $_->r;
  86.           push @{$data->{i}}, $_->i;
  87.       }
  88.       set_max() unless $flag;
  89.       last SWITCH;
  90.       };
  91.       $ref eq 'Math::Complex' and do {
  92.       $type = 'cmplx';
  93.       foreach (@$arr) {
  94.           die 'Found conflicting types in array data'
  95.           unless ref($_) eq $ref;
  96.           die "array data must be of type '$ref_in'"
  97.           if (defined $ref_in and $ref_in ne $ref);
  98.           push @{$data->{r}}, Re($_);
  99.           push @{$data->{i}}, Im($_);
  100.       }
  101.       set_max() unless $flag;
  102.       last SWITCH;
  103.       };
  104.       $ref eq 'Math::Cephes::Fraction' and do {
  105.       $type = 'fract';
  106.       foreach (@$arr) {
  107.           die 'Found conflicting types in array data'
  108.           unless ref($_) eq $ref;
  109.           die "array data must be of type '$ref_in'"
  110.           if (defined $ref_in and $ref_in ne $ref);
  111.           my ($gcd, $n, $d) = Math::Cephes::euclid($_->n, $_->d);
  112.           push @{$data->{n}}, $n;
  113.           push @{$data->{d}}, $d;
  114.       }
  115.       set_fmax() unless $fflag;
  116.       last SWITCH;
  117.       };
  118.        $ref eq 'Math::Fraction' and do {
  119.       $type = 'fract';
  120.       foreach (@$arr) {
  121.           die 'Found conflicting types in array data'
  122.           unless ref($_) eq $ref;
  123.           die "array data must be of type '$ref_in'"
  124.           if (defined $ref_in and $ref_in ne $ref);
  125.           push @{$data->{n}}, $_->{frac}->[0];
  126.           push @{$data->{d}}, $_->{frac}->[1];
  127.       }
  128.       set_fmax() unless $fflag;
  129.       last SWITCH;
  130.       };
  131.       die "Unknown type '$ref' in array data";
  132.   }
  133.     return ($type, $ref, $data, $n);
  134. }
  135.     
  136. sub as_string {
  137.   my $self = shift;
  138.   my ($type, $data, $n) =
  139.     ($self->{type}, $self->{data}, $self->{n});
  140.   my $d = shift || $n;
  141.   $d = $n if $d > $n;
  142.   my $string;
  143.     for (my $j=0; $j<=$d; $j++) {
  144.       my $coef;
  145.     SWITCH: {
  146.     $type eq 'fract' and do {
  147.           my $n = $data->{n}->[$j];
  148.           my $d = $data->{d}->[$j];
  149.           my $sgn = $n < 0 ? ' -' : ' +';
  150.           $coef = $sgn . ($j == 0? '(' : ' (') . 
  151.         abs($n) . '/' . abs($d) . ')';
  152.           last SWITCH;
  153.       };
  154.     $type eq 'cmplx' and do {
  155.           my $re = $data->{r}->[$j];
  156.           my $im = $data->{i}->[$j];
  157.           my $sgn = $j == 0 ? ' ' :  ' + ';
  158.           $coef = $sgn . '(' . $re . 
  159.         ( (int( $im / abs($im) ) == -1) ? '-' : '+' ) . 
  160.           ( ($im < 0) ? abs($im) : $im) . 'I)';
  161.           last SWITCH;
  162.         };
  163.     my $f = $data->[$j];
  164.     my $sgn = $f < 0 ? ' -' : ' +';
  165.     $coef = $j == 0 ? ' ' . $f :
  166.       $sgn . ' ' . abs($f);
  167.       }
  168.     $string .=  $coef . ($j > 0 ? "x^$j" : '');
  169.     }
  170.   return $string . "\n";
  171. }
  172.  
  173. sub add {
  174.     my ($self, $b) = @_;
  175.     my ($atype, $aref, $adata, $na) =
  176.     ($self->{type}, $self->{ref}, $self->{data}, $self->{n});
  177.     my ($btype, $bref, $bdata, $nb) =
  178.     ref($b) eq 'Math::Cephes::Polynomial' ?
  179.         ($b->{type}, $b->{ref}, $b->{data}, $b->{n}) :
  180.         get_data($b, $aref);
  181.     my $c = [];
  182.     my $nc;
  183.   SWITCH: {
  184.       $atype eq 'fract' and do {
  185.       $nc = $na > $nb ? $na: $nb;
  186.       my $cn = [split //, 0 x ($nc+1)];
  187.       my $cd = [split //, 0 x ($nc+1)];
  188.     Math::Cephes::fpoladd_wrap($adata->{n}, $adata->{d}, $na, 
  189.                    $bdata->{n}, $bdata->{d}, $nb, 
  190.                    $cn, $cd, $nc);
  191.       for (my $i=0; $i<=$nc; $i++) {
  192.           my ($gcd, $n, $d) = Math::Cephes::euclid($cn->[$i], $cd->[$i]);
  193.           push @$c, ($aref eq 'Math::Fraction' ? 
  194.              Math::Fraction->new($n, $d) :
  195.                Math::Cephes::Fraction->new($n, $d) );
  196.       }
  197.       last SWITCH;
  198.       };
  199.       $atype eq 'cmplx' and do {
  200.       $nc = $na > $nb ? $na: $nb;
  201.       my $cr = [split //, 0 x ($nc+1)];
  202.       my $ci = [split //, 0 x ($nc+1)];
  203.     Math::Cephes::poladd($adata->{r}, $na, 
  204.                  $bdata->{r}, $nb, $cr);
  205.     Math::Cephes::poladd($adata->{i}, $na, 
  206.                  $bdata->{i}, $nb, $ci);
  207.       for (my $i=0; $i<=$nc; $i++) {
  208.           push @$c, ($aref eq 'Math::Complex' ?
  209.                Math::Complex->make($cr->[$i], $ci->[$i]) :
  210.                Math::Cephes::Complex->new($cr->[$i], $ci->[$i]) );
  211.       }
  212.       last SWITCH;
  213.       };
  214.       $nc = $na > $nb ? $na + 1 : $nb + 1;
  215.       $c = [split //, 0 x $nc];
  216.     Math::Cephes::poladd($adata, $na, $bdata, $nb, $c);
  217.   }
  218.     return wantarray ? (Math::Cephes::Polynomial->new($c), $nc) : 
  219.       Math::Cephes::Polynomial->new($c);
  220.     
  221. }
  222.  
  223. sub sub {
  224.     my ($self, $b) = @_;
  225.     my ($atype, $aref, $adata, $na) =
  226.     ($self->{type}, $self->{ref}, $self->{data}, $self->{n});
  227.     my ($btype, $bref, $bdata, $nb) =
  228.     ref($b) eq 'Math::Cephes::Polynomial' ?
  229.         ($b->{type}, $b->{ref}, $b->{data}, $b->{n}) :
  230.         get_data($b, $aref);
  231.     my $c = [];
  232.     my $nc;
  233.   SWITCH: {
  234.       $atype eq 'fract' and do {
  235.       $nc = $na > $nb ? $na: $nb;
  236.       my $cn = [split //, 0 x ($nc+1)];
  237.       my $cd = [split //, 0 x ($nc+1)];
  238.     Math::Cephes::fpolsub_wrap($bdata->{n}, $bdata->{d}, $nb, 
  239.                    $adata->{n}, $adata->{d}, $na, 
  240.                    $cn, $cd, $nc);
  241.       for (my $i=0; $i<=$nc; $i++) {
  242.           my ($gcd, $n, $d) = Math::Cephes::euclid($cn->[$i], $cd->[$i]);
  243.           push @$c, ($aref eq 'Math::Fraction' ? 
  244.              Math::Fraction->new($n, $d) :
  245.                Math::Cephes::Fraction->new($n, $d) );
  246.       }
  247.       last SWITCH;
  248.       };
  249.       $atype eq 'cmplx' and do {
  250.       $nc = $na > $nb ? $na: $nb;
  251.       my $cr = [split //, 0 x ($nc+1)];
  252.       my $ci = [split //, 0 x ($nc+1)];
  253.     Math::Cephes::polsub($bdata->{r}, $nb, 
  254.                  $adata->{r}, $na, $cr);
  255.     Math::Cephes::polsub($bdata->{i}, $nb, 
  256.                  $adata->{i}, $na, $ci);
  257.       for (my $i=0; $i<=$nc; $i++) {
  258.           push @$c, ($aref eq 'Math::Complex' ?
  259.                Math::Complex->make($cr->[$i], $ci->[$i]) :
  260.                Math::Cephes::Complex->new($cr->[$i], $ci->[$i]) );
  261.       }
  262.       last SWITCH;
  263.       };
  264.       $nc = $na > $nb ? $na + 1 : $nb + 1;
  265.       $c = [split //, 0 x $nc];
  266.     Math::Cephes::polsub($bdata, $nb, $adata, $na, $c);
  267.   }
  268.     return wantarray ? (Math::Cephes::Polynomial->new($c), $nc) : 
  269.       Math::Cephes::Polynomial->new($c);
  270.     
  271. }
  272.  
  273. sub mul {
  274.     my ($self, $b) = @_;
  275.     my ($atype, $aref, $adata, $na) =
  276.     ($self->{type}, $self->{ref}, $self->{data}, $self->{n});
  277.     my ($btype, $bref, $bdata, $nb) =
  278.     ref($b) eq 'Math::Cephes::Polynomial' ?
  279.         ($b->{type}, $b->{ref}, $b->{data}, $b->{n}) :
  280.         get_data($b, $aref);
  281.     my $c = [];
  282.     my $nc;
  283.   SWITCH: {
  284.       $atype eq 'fract' and do {
  285.       $nc = $na + $nb;
  286.       my $cn = [split //, 0 x ($nc+1)];
  287.       my $cd = [split //, 1 x ($nc+1)];
  288.     Math::Cephes::fpolmul_wrap($adata->{n}, $adata->{d}, $na, 
  289.                    $bdata->{n}, $bdata->{d}, $nb, 
  290.                    $cn, $cd, $nc);
  291.       for (my $i=0; $i<=$nc; $i++) {
  292.           my ($gcd, $n, $d) = Math::Cephes::euclid($cn->[$i], $cd->[$i]);
  293.           push @$c, ($aref eq 'Math::Fraction' ? 
  294.              Math::Fraction->new($n, $d) :
  295.                Math::Cephes::Fraction->new($n, $d) );
  296.       }
  297.       last SWITCH;
  298.       };
  299.       $atype eq 'cmplx' and do {
  300.       my $dc = $na + $nb + 3;
  301.       my $cr = [split //, 0 x $dc];
  302.       my $ci = [split //, 0 x $dc];
  303.       $nc = Math::Cephes::cpmul_wrap($adata->{r}, $adata->{i}, $na+1, 
  304.                      $bdata->{r}, $bdata->{i}, $nb+1, 
  305.                      $cr, $ci, $dc);
  306.       $cr = [ @{$cr}[0..$nc] ];
  307.       $ci = [ @{$ci}[0..$nc] ];
  308.       for (my $i=0; $i<=$nc; $i++) {
  309.           push @$c, ($aref eq 'Math::Complex' ?
  310.                Math::Complex->make($cr->[$i], $ci->[$i]) :
  311.                Math::Cephes::Complex->new($cr->[$i], $ci->[$i]) );
  312.       }
  313.       last SWITCH;
  314.       };
  315.       $nc = $na > $nb ? $na + 3 : $nb + 3;
  316.       $c = [split //, 0 x $nc];
  317.     Math::Cephes::polmul($adata, $na, $bdata, $nb, $c);
  318.   }
  319.     return wantarray ? (Math::Cephes::Polynomial->new($c), $nc) : 
  320.       Math::Cephes::Polynomial->new($c);
  321. }
  322.  
  323. sub div {
  324.     my ($self, $b) = @_;
  325.     my ($atype, $aref, $adata, $na) =
  326.     ($self->{type}, $self->{ref}, $self->{data}, $self->{n});
  327.     my ($btype, $bref, $bdata, $nb) =
  328.     ref($b) eq 'Math::Cephes::Polynomial' ?
  329.         ($b->{type}, $b->{ref}, $b->{data}, $b->{n}) :
  330.         get_data($b, $aref);
  331.     my $c = [];
  332.     my $nc;
  333.   SWITCH: {
  334.       $atype eq 'fract' and do {
  335.       $nc = $MAXPOL;
  336.       my $cn = [split //, 0 x ($nc+1)];
  337.       my $cd = [split //, 0 x ($nc+1)];
  338.     Math::Cephes::fpoldiv_wrap($adata->{n}, $adata->{d}, $na, 
  339.                    $bdata->{n}, $bdata->{d}, $nb, 
  340.                    $cn, $cd, $nc);
  341.       for (my $i=0; $i<=$nc; $i++) {
  342.           my ($gcd, $n, $d) = Math::Cephes::euclid($cn->[$i], $cd->[$i]);
  343.           push @$c, ($aref eq 'Math::Fraction' ? 
  344.              Math::Fraction->new($n, $d) :
  345.                Math::Cephes::Fraction->new($n, $d) );
  346.       }
  347.       last SWITCH;
  348.       };
  349.       $atype eq 'cmplx' and do {
  350.       die "Cannot do complex division";
  351.       last SWITCH;
  352.       };
  353.       $nc = $MAXPOL;
  354.       $c = [split //, 0 x ($nc+1)];
  355.     Math::Cephes::poldiv($adata, $na, $bdata, $nb, $c);
  356.   }
  357.     return wantarray ? (Math::Cephes::Polynomial->new($c), $nc) : 
  358.       Math::Cephes::Polynomial->new($c);
  359. }
  360.  
  361. sub clr {
  362.   my $self = shift;
  363.   my ($atype, $aref, $adata, $na) =
  364.     ($self->{type}, $self->{ref}, $self->{data}, $self->{n});
  365.   set_max() unless $flag;
  366.   my $n = shift || $na;
  367.   $n = $na if $n > $na;
  368.  SWITCH: {
  369.       $atype eq 'fract' and do {
  370.       for (my $i=0; $i<=$n; $i++) {
  371.         $self->{data}->{n}->[$i] = 0;
  372.         $self->{data}->{d}->[$i] = 1;
  373.       }
  374.       last SWITCH;
  375.       };
  376.       $atype eq 'cmplx' and do {
  377.       for (my $i=0; $i<=$n; $i++) {
  378.         $self->{data}->{r}->[$i] = 0;
  379.         $self->{data}->{i}->[$i] = 0;
  380.       }
  381.       last SWITCH;
  382.       };
  383.       for (my $i=0; $i<=$n; $i++) {
  384.     $self->{data}->[$i] = 0;
  385.       }
  386.    }
  387. }
  388.  
  389. sub sbt {
  390.   my ($self, $b) = @_;
  391.   my ($atype, $aref, $adata, $na) =
  392.     ($self->{type}, $self->{ref}, $self->{data}, $self->{n});
  393.   my ($btype, $bref, $bdata, $nb) =
  394.     ref($b) eq 'Math::Cephes::Polynomial' ?
  395.       ($b->{type}, $b->{ref}, $b->{data}, $b->{n}) :
  396.     get_data($b, $aref);
  397.   set_max() unless $flag;
  398.   my $c = [];
  399.   my $nc;
  400.  SWITCH: {
  401.     $atype eq 'fract' and do {
  402.       $nc = ($na+1)*($nb+1);
  403.       my $cn = [split //, 0 x ($nc+1)];
  404.       my $cd = [split //, 0 x ($nc+1)];
  405.       Math::Cephes::fpolsbt_wrap($bdata->{n}, $bdata->{d}, $nb,
  406.                  $adata->{n}, $adata->{d}, $na,
  407.                  $cn, $cd, $nc);
  408.       $nc = $na * $nb;
  409.       for (my $i=0; $i<=$nc; $i++) {
  410.     my ($gcd, $n, $d) = Math::Cephes::euclid($cn->[$i], $cd->[$i]);
  411.     push @$c, ($aref eq 'Math::Fraction' ? 
  412.            Math::Fraction->new($n, $d) :
  413.            Math::Cephes::Fraction->new($n, $d) );
  414.       }
  415.       last SWITCH;
  416.     };
  417.     $atype eq 'cmplx' and do {
  418.       die "Cannot do complex substitution";
  419.       last SWITCH;
  420.     };
  421.     $nc = ($na+1)*($nb+1);
  422.     $c = [split //, 0 x $nc];
  423.     Math::Cephes::polsbt($bdata, $nb, $adata, $na, $c);
  424.     $nc = $na*$nb;
  425.     $c = [@$c[0..$nc]];
  426.   }
  427.   return wantarray ? (Math::Cephes::Polynomial->new($c), $nc) : 
  428.     Math::Cephes::Polynomial->new($c);
  429. }
  430.  
  431. sub set_max {
  432.   Math::Cephes::polini($MAXPOL);
  433.     $flag = 1;
  434. }
  435.  
  436. sub set_fmax {
  437.   Math::Cephes::fpolini($FMAXPOL);
  438.     $fflag = 1;
  439. }
  440.  
  441. sub eval {
  442.   my $self = shift;
  443.   my $x = 0 || shift;
  444.   my ($atype, $aref, $adata, $na) =
  445.     ($self->{type}, $self->{ref}, $self->{data}, $self->{n});
  446.   my $y;
  447.  SWITCH: {
  448.     $atype eq 'fract' and do {
  449.       my $xref = ref($x);
  450.       $y = Math::Cephes::Fraction->new(0, 1);
  451.     FRACT: {
  452.     not $xref and do {
  453.       $x = Math::Cephes::Fraction->new($x, 1);
  454.       last FRACT;
  455.     };
  456.     $xref eq 'Math::Cephes::Fraction' and do {
  457.       last FRACT;
  458.     };
  459.     $xref eq 'Math::Fraction' and do {
  460.       $x = Math::Cephes::Fraction->new($x->{frac}->[0], $x->{frac}->[1]);
  461.       last FRACT;
  462.     };
  463.     die "Unknown data type '$xref' for x";
  464.       }
  465.       Math::Cephes::fpoleva_wrap($adata->{n}, $adata->{d}, $na, $x, $y);
  466.       $y = Math::Fraction->new($y->n, $y->d) if $aref eq 'Math::Fraction';
  467.       last SWITCH;
  468.     };
  469.     $atype eq 'cmplx' and do {
  470.       my $r = Math::Cephes::poleva($adata->{r}, $na, $x);
  471.       my $i = Math::Cephes::poleva($adata->{i}, $na, $x);
  472.       $y = $aref eq 'Math::Complex' ?
  473.     Math::Complex->make($r, $i) :
  474.         Math::Cephes::Complex->new($r, $i);
  475.       last SWITCH;
  476.     };
  477.     $y = Math::Cephes::poleva($adata, $na, $x);
  478.   }
  479.   return $y;
  480. }
  481.  
  482. sub fract_to_real {
  483.     my $in = shift;
  484.     my $a = [];
  485.     my $n = scalar @{$in->{n}} - 1;
  486.     for (my $i=0; $i<=$n; $i++) {
  487.     push @$a, $in->{n}->[$i] / $in->{d}->[$i];
  488.     }
  489.     return $a;
  490. }
  491.  
  492. sub atn {
  493.     my ($self, $bin) = @_;
  494.     my $type = $self->{type};
  495.     die "Cannot take the atan of a complex polynomial"
  496.     if $type eq 'cmplx';
  497.     my ($a, $b);
  498.     my ($atype, $aref, $adata, $na) =
  499.     ($self->{type}, $self->{ref}, $self->{data}, $self->{n});
  500.     die "Cannot take the atan of a complex polynomial"
  501.     if $atype eq 'cmplx';
  502.     $a = $atype eq 'fract' ? fract_to_real($adata) : $adata;
  503.  
  504.     my ($btype, $bref, $bdata, $nb) =
  505.     ref($bin) eq 'Math::Cephes::Polynomial' ?
  506.         ($bin->{type}, $bin->{ref}, $bin->{data}, $bin->{n}) :
  507.         get_data($bin);
  508.     
  509.     die "Cannot take the atan of a complex polynomial"
  510.     if $btype eq 'cmplx';
  511.     $b = $btype eq 'fract' ? fract_to_real($bdata) : $bdata;
  512.  
  513.     my $c = [split //, 0 x ($MAXPOL+1)];
  514.   Math::Cephes::polatn($a, $b, $c, 16);
  515.     return Math::Cephes::Polynomial->new($c);
  516. }
  517.  
  518. sub sqt {
  519.   my $self = shift;
  520.   my $type = $self->{type};
  521.   die "Cannot take the sqrt of a complex polynomial"
  522.     if $type eq 'cmplx';
  523.   my $a = $type eq 'fract' ? fract_to_real($self->{data}) : $self->coef;
  524.   my $b = [split //, 0 x ($MAXPOL+1)];
  525.   Math::Cephes::polsqt($a, $b, 16);
  526.   return Math::Cephes::Polynomial->new($b);
  527. }
  528.  
  529. sub sin {
  530.   my $self = shift;
  531.   my $type = $self->{type};
  532.   die "Cannot take the sin of a complex polynomial"
  533.     if $type eq 'cmplx';
  534.   my $a = $type eq 'fract' ? fract_to_real($self->{data}) : $self->coef;
  535.   my $b = [split //, 0 x ($MAXPOL+1)];
  536.   Math::Cephes::polsin($a, $b, 16);
  537.   return Math::Cephes::Polynomial->new($b);
  538. }
  539.  
  540. sub cos {
  541.   my $self = shift;
  542.   my $type = $self->{type};
  543.   die "Cannot take the cos of a complex polynomial"
  544.     if $type eq 'cmplx';
  545.   my $a = $type eq 'fract' ? fract_to_real($self->{data}) : $self->coef;
  546.   my $b = [split //, 0 x ($MAXPOL+1)];
  547.   Math::Cephes::polcos($a, $b, 16);
  548.   return Math::Cephes::Polynomial->new($b);
  549. }
  550.  
  551. sub rts {
  552.   my $self = shift;
  553.   my ($atype, $aref, $adata, $na) =
  554.     ($self->{type}, $self->{ref}, $self->{data}, $self->{n});
  555.   my ($a, $b, $ret);
  556.   my $cof = [split //, 0 x ($na+1)];
  557.   my $r = [split //, 0 x ($na+1)];
  558.   my $i = [split //, 0 x ($na+1)];
  559.  SWITCH: {
  560.     $atype eq 'fract' and do {
  561.       $adata = fract_to_real($adata);
  562.       $ret = Math::Cephes::polrt_wrap($adata, $cof, $na, $r, $i);
  563.       for (my $j=0; $j<$na; $j++) {
  564.     push @$b,
  565.       Math::Cephes::Complex->new($r->[$j], $i->[$j]);
  566.       }
  567.       last SWITCH;
  568.     };
  569.     $atype eq 'cmplx' and do {
  570.       die "Cannot do complex root finding";
  571.       last SWITCH;
  572.     };
  573.     $ret = Math::Cephes::polrt_wrap($adata, $cof, $na, $r, $i);
  574.     for (my $j=0; $j<$na; $j++) {
  575.       push @$b,
  576.     Math::Cephes::Complex->new($r->[$j], $i->[$j]);
  577.     }
  578.   }
  579.   return wantarray ? ($ret, $b) : $b;
  580. }
  581.  
  582. 1;
  583.  
  584. __END__
  585.  
  586. =head1 NAME
  587.  
  588. Math::Cephes::Polynomial - Perl interface to the cephes math polynomial routines
  589.  
  590. =head1 SYNOPSIS
  591.  
  592.   use Math::Cephes::Polynomial qw(poly);
  593.   # 'poly' is a shortcut for Math::Cephes::Polynomial->new
  594.  
  595.   require Math::Cephes::Fraction; # if coefficients are fractions
  596.   require Math::Cephes::Complex;  # if coefficients are complex
  597.  
  598.   my $a = poly([1, 2, 3]);           # a(x) = 1 + 2x + 3x^2
  599.   my $b = poly([4, 5, 6, 7];         # b(x) = 4 + 5x + 6x^2 + 7x^3
  600.   my $c = $a->add($b);               # c(x) = 5 + 7x + 9x^2 + 7x^3
  601.   my $cc = $c->coef;
  602.   for (my $i=0; $i<4; $i++) {
  603.      print "term $i: $cc->[$i]\n";
  604.   }
  605.   my $x = 2;
  606.   my $r = $c->eval($x);
  607.   print "At x=$x, c(x) is $r\n";
  608.  
  609.   my $u1 = Math::Cephes::Complex->new(2,1);
  610.   my $u2 = Math::Cephes::Complex->new(1,-3);
  611.   my $v1 = Math::Cephes::Complex->new(1,3);
  612.   my $v2 = Math::Cephes::Complex->new(2,4);
  613.   my $z1 = Math::Cephes::Polynomial->new([$u1, $u2]);
  614.   my $z2 = Math::Cephes::Polynomial->new([$v1, $v2]);
  615.   my $z3 = $z1->add($z2);
  616.   my $z3c = $z3->coef;
  617.   for (my $i=0; $i<2; $i++) {
  618.      print "term $i: real=$z3c->{r}->[$i], imag=$z3c->{i}->[$i]\n";
  619.   }
  620.   $r = $z3->eval($x);
  621.   print "At x=$x, z3(x) has real=", $r->r, " and imag=", $r->i, "\n";
  622.  
  623.   my $a1 = Math::Cephes::Fraction->new(1,2);
  624.   my $a2 = Math::Cephes::Fraction->new(2,1);
  625.   my $b1 = Math::Cephes::Fraction->new(1,2);
  626.   my $b2 = Math::Cephes::Fraction->new(2,2);
  627.   my $f1 = Math::Cephes::Polynomial->new([$a1, $a2]);
  628.   my $f2 = Math::Cephes::Polynomial->new([$b1, $b2]);
  629.   my $f3 = $f1->add($f2);
  630.   my $f3c = $f3->coef;
  631.   for (my $i=0; $i<2; $i++) {
  632.      print "term $i: num=$f3c->{n}->[$i], den=$f3c->{d}->[$i]\n";
  633.   }
  634.   $r = $f3->eval($x);
  635.   print "At x=$x, f3(x) has num=", $r->n, " and den=", $r->d, "\n";
  636.   $r = $f3->eval($a1);
  637.   print "At x=", $a1->n, "/", $a1->d, 
  638.       ", f3(x) has num=", $r->n, " and den=", $r->d, "\n";
  639.  
  640. =head1 DESCRIPTION
  641.  
  642. This module is a layer on top of the basic routines in the
  643. cephes math library to handle polynomials. In the following,
  644. a Math::Cephes::Polynomial object is created as
  645.  
  646.   my $p = Math::Cephes::Polynomial->new($arr_ref);
  647.  
  648. where C<$arr_ref> is a reference to an array which can
  649. consist of one of
  650.  
  651. =over
  652.  
  653. =item * floating point numbers, for polynomials with floating
  654. point coefficients,
  655.  
  656. =item * I<Math::Cephes::Fraction> or I<Math::Fraction> objects,
  657. for polynomials with fractional coefficients,
  658.  
  659. =item * I<Math::Cephes::Complex> or I<Math::Complex> objects,
  660. for polynomials with complex coefficients,
  661.  
  662. =back
  663.  
  664. The maximum degree of the polynomials handled is set by default
  665. to 256 - this can be changed by setting I<$Math::Cephes::Polynomial::MAXPOL>.
  666.  
  667. A copy of a I<Math::Cephes::Polynomial> object may be done as
  668.  
  669.   my $p_copy = $p->new();
  670.  
  671. and a string representation of the polynomial may be gotten through
  672.  
  673.   print $p->as_string;
  674.  
  675. =head2 Methods
  676.  
  677. The following methods are available.
  678.  
  679. =over 4
  680.  
  681. =item I<coef>: get coefficients of the polynomial
  682.  
  683.  SYNOPSIS:
  684.  
  685.  my $c = $p->coef;
  686.  
  687.  DESCRIPTION:
  688.  
  689. This returns an array reference containing the coefficients of 
  690. the polynomial. 
  691.  
  692. =item I<clr>: set a polynomial identically equal to zero
  693.  
  694.  SYNOPSIS:
  695.  
  696.  $p->clr($n);
  697.  
  698.  DESCRIPTION:
  699.  
  700. This sets the coefficients of the polynomial identically to 0, 
  701. up to $p->[$n]. If $n is omitted, all elements are set to 0. 
  702.  
  703. =item I<add>: add two polynomials
  704.  
  705.  SYNOPSIS:
  706.  
  707.  $c = $a->add($b);
  708.  
  709.  DESCRIPTION:
  710.  
  711. This sets $c equal to $a + $b.
  712.  
  713. =item I<sub>: subtract two polynomials
  714.  
  715.  SYNOPSIS:
  716.  
  717.  $c = $a->sub($b);
  718.  
  719.  DESCRIPTION:
  720.  
  721. This sets $c equal to $a - $b.
  722.  
  723. =item I<mul>: multiply two polynomials
  724.  
  725.  SYNOPSIS:
  726.  
  727.  $c = $a->mul($b);
  728.  
  729.  DESCRIPTION:
  730.  
  731. This sets $c equal to $a * $b.
  732.  
  733. =item I<div>: divide two polynomials
  734.  
  735.  SYNOPSIS:
  736.  
  737.  $c = $a->div($b);
  738.  
  739.  DESCRIPTION:
  740.  
  741. This sets $c equal to $a / $b, expanded by a Taylor series. Accuracy
  742. is approximately equal to the degree of the polynomial, with an
  743. internal limit of about 16.
  744.  
  745. =item I<sbt>: change of variables
  746.  
  747.  SYNOPSIS:
  748.  
  749.  $c = $a->sbt($b);
  750.  
  751.  DESCRIPTION:
  752.  
  753. If a(x) and b(x) are polynomials, then
  754.  
  755.      c(x) = a(b(x))
  756.  
  757. is a polynomial found by substituting b(x) for x in a(x). This method is
  758. not available for polynomials with complex coefficients.
  759.  
  760. =item I<eval>: evaluate a polynomial
  761.  
  762.  SYNOPSIS:
  763.  
  764.  $s = $a->eval($x);
  765.  
  766.  DESCRIPTION:
  767.  
  768. This evaluates the polynomial at the value $x. The returned value
  769. is of the same type as that used to represent the coefficients of
  770. the polynomial.
  771.  
  772. =item I<sqt>: square root of a polynomial
  773.  
  774.  SYNOPSIS:
  775.  
  776.  $b = $a->sqt();
  777.  
  778.  DESCRIPTION:
  779.  
  780. This finds the square root of a polynomial, evaluated by a 
  781. Taylor expansion. Accuracy is approximately equal to the 
  782. degree of the polynomial, with an internal limit of about 16.
  783. This method is not available for polynomials with complex coefficients.
  784.  
  785. =item I<sin>: sine of a polynomial
  786.  
  787.  SYNOPSIS:
  788.  
  789.  $b = $a->sin();
  790.  
  791.  DESCRIPTION:
  792.  
  793. This finds the sine of a polynomial, evaluated by a 
  794. Taylor expansion. Accuracy is approximately equal to the 
  795. degree of the polynomial, with an internal limit of about 16.
  796. This method is not available for polynomials with complex coefficients.
  797.  
  798. =item I<cos>: cosine of a polynomial
  799.  
  800.  SYNOPSIS:
  801.  
  802.  $b = $a->cos();
  803.  
  804.  DESCRIPTION:
  805.  
  806. This finds the cosine of a polynomial, evaluated by a 
  807. Taylor expansion. Accuracy is approximately equal to the 
  808. degree of the polynomial, with an internal limit of about 16.
  809. This method is not available for polynomials with complex coefficients.
  810.  
  811. =item I<atn>: arctangent of the ratio of two polynomials
  812.  
  813.  SYNOPSIS:
  814.  
  815.  $c = $a->atn($b);
  816.  
  817.  DESCRIPTION:
  818.  
  819. This finds the arctangent of the ratio $a / $b of two polynomial, 
  820. evaluated by a Taylor expansion. Accuracy is approximately equal to the 
  821. degree of the polynomial, with an internal limit of about 16.
  822. This method is not available for polynomials with complex coefficients.
  823.  
  824. =item I<rts>: roots of a polynomial
  825.  
  826.  SYNOPSIS:
  827.  
  828.   my $w = Math::Cephes::Polynomial->new([-2, 0, -1, 0, 1]);
  829.   my ($flag, $r) = $w->rts();
  830.   for (my $i=0; $i<4; $i++) {
  831.     print "Root $i has real=", $r->[$i]->r, " and imag=", $r->[$i]->i, "\n";
  832.   }
  833.  
  834.  DESCRIPTION:
  835.  
  836. This finds the roots of a polynomial. I<$flag>, if non-zero, indicates
  837. a failure of some kind. I<$roots> in an array reference of
  838. I<Math::Cephes::Complex> objects holding the
  839. real and complex values of the roots found.
  840. This method is not available for polynomials with complex coefficients.
  841.  
  842.  ACCURACY:
  843.  
  844. Termination depends on evaluation of the polynomial at
  845. the trial values of the roots.  The values of multiple roots
  846. or of roots that are nearly equal may have poor relative
  847. accuracy after the first root in the neighborhood has been
  848. found.
  849.  
  850. =back
  851.  
  852. =head1 BUGS
  853.  
  854. Please report any to Randy Kobes <randy@theoryx5.uwinnipeg.ca>
  855.  
  856. =head1 COPYRIGHT
  857.  
  858. The C code for the Cephes Math Library is
  859. Copyright 1984, 1987, 1989, 2002 by Stephen L. Moshier, 
  860. and is available at http://www.netlib.org/cephes/.
  861. Direct inquiries to 30 Frost Street, Cambridge, MA 02140.
  862.  
  863. The perl interface is copyright 2000, 2002 by Randy Kobes.
  864. This library is free software; you can redistribute it and/or
  865. modify it under the same terms as Perl itself.
  866.  
  867. =cut
  868.