File Coverage

blib/lib/Math/Cephes/Polynomial.pm

Criterion	Covered	Total	%
statement	299	359	83.2
branch	108	194	55.6
condition	10	27	37.0
subroutine	21	23	91.3
pod	1	20	5.0
total	439	623	70.4

package Math::Cephes::Polynomial;
use strict;
use warnings;
use vars qw(@EXPORT_OK $VERSION $MAXPOL $FMAXPOL $flag $fflag);
eval {require Math::Complex; import Math::Complex qw(Re Im)};
eval {local $^W=0; require Math::Fraction;};
$MAXPOL = 256;
$flag = 0;
$FMAXPOL = 256;
$fflag = 0;

require Exporter;
*import = \&Exporter::import;
@EXPORT_OK = qw(poly);
$VERSION = '0.5308';

require Math::Cephes;
require Math::Cephes::Fraction;
require Math::Cephes::Complex;

sub new {
    my ($caller, $arr) = @_;
    my $refer = ref($caller);
    my $class = $refer || $caller;
    die "Must supply data for the polynomial"
      unless ($refer or $arr);
    my ($type, $ref, $data, $n);
    if ($refer) {
      ($type, $ref, $n) =
	($caller->{type}, $caller->{ref}, $caller->{n});
      my $cdata = $caller->{data};
      if (ref($cdata) eq 'ARRAY') {
	$data = [ @$cdata ];
      }
      else {
	my ($f, $s) = ($type eq 'fract') ? ('n', 'd') : ('r', 'i');
	$data = {$f => [ @{$cdata->{$f}} ],
		 $s => [ @{$cdata->{$s}} ],
		};
      }
    }
    else {
      ($type, $ref, $data, $n) = get_data($arr);
    }
    bless { type => $type,
	    ref => $ref,
	    data => $data,
	    n => $n,
	}, $class;
}

sub poly {
  return Math::Cephes::Polynomial->new(shift);
}

sub coef {
    return $_[0]->{data};
}

sub get_data {
    my ($arr, $ref_in) = @_;
    die "Must supply an array reference" unless ref($arr) eq 'ARRAY';
    my $n = scalar @$arr - 1;
    my $ref = ref($arr->[0]);
    die "array data must be of type '$ref_in'"
	if (defined $ref_in and $ref_in ne $ref);
    my ($type, $data);
  SWITCH: {
      not $ref and do {
	  $type = 'scalar';
	  foreach (@$arr) {
	      die 'Found conflicting types in array data'
		  if ref($_);
	  }
	  $data = $arr;
	  set_max() unless $flag;
	  last SWITCH;
      };
      $ref eq 'Math::Cephes::Complex' and do {
	  $type = 'cmplx';
	  foreach (@$arr) {
	      die 'Found conflicting types in array data'
		  unless ref($_) eq $ref;
	      die "array data must be of type '$ref_in'"
		  if (defined $ref_in and $ref_in ne $ref);
	      push @{$data->{r}}, $_->r;
	      push @{$data->{i}}, $_->i;
	  }
	  set_max() unless $flag;
	  last SWITCH;
      };
      $ref eq 'Math::Complex' and do {
	  $type = 'cmplx';
	  foreach (@$arr) {
	      die 'Found conflicting types in array data'
		  unless ref($_) eq $ref;
	      die "array data must be of type '$ref_in'"
		  if (defined $ref_in and $ref_in ne $ref);
	      push @{$data->{r}}, Re($_);
	      push @{$data->{i}}, Im($_);
	  }
	  set_max() unless $flag;
	  last SWITCH;
      };
      $ref eq 'Math::Cephes::Fraction' and do {
	  $type = 'fract';
	  foreach (@$arr) {
	      die 'Found conflicting types in array data'
		  unless ref($_) eq $ref;
	      die "array data must be of type '$ref_in'"
		  if (defined $ref_in and $ref_in ne $ref);
	      my ($gcd, $n, $d) = Math::Cephes::euclid($_->n, $_->d);
	      push @{$data->{n}}, $n;
	      push @{$data->{d}}, $d;
	  }
	  set_fmax() unless $fflag;
	  last SWITCH;
      };
       $ref eq 'Math::Fraction' and do {
	  $type = 'fract';
	  foreach (@$arr) {
	      die 'Found conflicting types in array data'
		  unless ref($_) eq $ref;
	      die "array data must be of type '$ref_in'"
		  if (defined $ref_in and $ref_in ne $ref);
	      push @{$data->{n}}, $_->{frac}->[0];
	      push @{$data->{d}}, $_->{frac}->[1];
	  }
	  set_fmax() unless $fflag;
	  last SWITCH;
      };
      die "Unknown type '$ref' in array data";
  }
    return ($type, $ref, $data, $n);
}

sub as_string {
  my $self = shift;
  my ($type, $data, $n) =
    ($self->{type}, $self->{data}, $self->{n});
  my $d = shift || $n;
  $d = $n if $d > $n;
  my $string;
    for (my $j=0; $j<=$d; $j++) {
      my $coef;
    SWITCH: {
	$type eq 'fract' and do {
	      my $n = $data->{n}->[$j];
	      my $d = $data->{d}->[$j];
	      my $sgn = $n < 0 ? ' -' : ' +';
	      $coef = $sgn . ($j == 0? '(' : ' (') .
		abs($n) . '/' . abs($d) . ')';
	      last SWITCH;
	  };
	$type eq 'cmplx' and do {
	      my $re = $data->{r}->[$j];
	      my $im = $data->{i}->[$j];
	      my $sgn = $j == 0 ? ' ' :  ' + ';
	      $coef = $sgn . '(' . $re .
		( (int( $im / abs($im) ) == -1) ? '-' : '+' ) .
		  ( ($im < 0) ? abs($im) : $im) . 'I)';
	      last SWITCH;
	    };
	my $f = $data->[$j];
	my $sgn = $f < 0 ? ' -' : ' +';
	$coef = $j == 0 ? ' ' . $f :
	  $sgn . ' ' . abs($f);
      }
	$string .=  $coef . ($j > 0 ? "x^$j" : '');
    }
  return $string . "\n";
}

sub add {
    my ($self, $b) = @_;
    my ($atype, $aref, $adata, $na) =
	($self->{type}, $self->{ref}, $self->{data}, $self->{n});
    my ($btype, $bref, $bdata, $nb) =
	ref($b) eq 'Math::Cephes::Polynomial' ?
	    ($b->{type}, $b->{ref}, $b->{data}, $b->{n}) :
		get_data($b, $aref);
    my $c = [];
    my $nc;
  SWITCH: {
      $atype eq 'fract' and do {
	  $nc = $na > $nb ? $na: $nb;
	  my $cn = [split //, 0 x ($nc+1)];
	  my $cd = [split //, 0 x ($nc+1)];
	Math::Cephes::fpoladd_wrap($adata->{n}, $adata->{d}, $na,
				   $bdata->{n}, $bdata->{d}, $nb,
				   $cn, $cd, $nc);
	  for (my $i=0; $i<=$nc; $i++) {
	      my ($gcd, $n, $d) = Math::Cephes::euclid($cn->[$i], $cd->[$i]);
	      push @$c, ($aref eq 'Math::Fraction' ?
			 Math::Fraction->new($n, $d) :
		       Math::Cephes::Fraction->new($n, $d) );
	  }
	  last SWITCH;
      };
      $atype eq 'cmplx' and do {
	  $nc = $na > $nb ? $na: $nb;
	  my $cr = [split //, 0 x ($nc+1)];
	  my $ci = [split //, 0 x ($nc+1)];
	Math::Cephes::poladd($adata->{r}, $na,
			     $bdata->{r}, $nb, $cr);
	Math::Cephes::poladd($adata->{i}, $na,
			     $bdata->{i}, $nb, $ci);
	  for (my $i=0; $i<=$nc; $i++) {
	      push @$c, ($aref eq 'Math::Complex' ?
		       Math::Complex->make($cr->[$i], $ci->[$i]) :
		       Math::Cephes::Complex->new($cr->[$i], $ci->[$i]) );
	  }
	  last SWITCH;
      };
      $nc = $na > $nb ? $na + 1 : $nb + 1;
      $c = [split //, 0 x $nc];
    Math::Cephes::poladd($adata, $na, $bdata, $nb, $c);
  }
    return wantarray ? (Math::Cephes::Polynomial->new($c), $nc) :
      Math::Cephes::Polynomial->new($c);

}

sub sub {
    my ($self, $b) = @_;
    my ($atype, $aref, $adata, $na) =
	($self->{type}, $self->{ref}, $self->{data}, $self->{n});
    my ($btype, $bref, $bdata, $nb) =
	ref($b) eq 'Math::Cephes::Polynomial' ?
	    ($b->{type}, $b->{ref}, $b->{data}, $b->{n}) :
		get_data($b, $aref);
    my $c = [];
    my $nc;
  SWITCH: {
      $atype eq 'fract' and do {
	  $nc = $na > $nb ? $na: $nb;
	  my $cn = [split //, 0 x ($nc+1)];
	  my $cd = [split //, 0 x ($nc+1)];
	Math::Cephes::fpolsub_wrap($bdata->{n}, $bdata->{d}, $nb,
				   $adata->{n}, $adata->{d}, $na,
				   $cn, $cd, $nc);
	  for (my $i=0; $i<=$nc; $i++) {
	      my ($gcd, $n, $d) = Math::Cephes::euclid($cn->[$i], $cd->[$i]);
	      push @$c, ($aref eq 'Math::Fraction' ?
			 Math::Fraction->new($n, $d) :
		       Math::Cephes::Fraction->new($n, $d) );
	  }
	  last SWITCH;
      };
      $atype eq 'cmplx' and do {
	  $nc = $na > $nb ? $na: $nb;
	  my $cr = [split //, 0 x ($nc+1)];
	  my $ci = [split //, 0 x ($nc+1)];
	Math::Cephes::polsub($bdata->{r}, $nb,
			     $adata->{r}, $na, $cr);
	Math::Cephes::polsub($bdata->{i}, $nb,
			     $adata->{i}, $na, $ci);
	  for (my $i=0; $i<=$nc; $i++) {
	      push @$c, ($aref eq 'Math::Complex' ?
		       Math::Complex->make($cr->[$i], $ci->[$i]) :
		       Math::Cephes::Complex->new($cr->[$i], $ci->[$i]) );
	  }
	  last SWITCH;
      };
      $nc = $na > $nb ? $na + 1 : $nb + 1;
      $c = [split //, 0 x $nc];
    Math::Cephes::polsub($bdata, $nb, $adata, $na, $c);
  }
    return wantarray ? (Math::Cephes::Polynomial->new($c), $nc) :
      Math::Cephes::Polynomial->new($c);

}

sub mul {
    my ($self, $b) = @_;
    my ($atype, $aref, $adata, $na) =
	($self->{type}, $self->{ref}, $self->{data}, $self->{n});
    my ($btype, $bref, $bdata, $nb) =
	ref($b) eq 'Math::Cephes::Polynomial' ?
	    ($b->{type}, $b->{ref}, $b->{data}, $b->{n}) :
		get_data($b, $aref);
    my $c = [];
    my $nc;
  SWITCH: {
      $atype eq 'fract' and do {
	  $nc = $na + $nb;
	  my $cn = [split //, 0 x ($nc+1)];
	  my $cd = [split //, 1 x ($nc+1)];
	Math::Cephes::fpolmul_wrap($adata->{n}, $adata->{d}, $na,
				   $bdata->{n}, $bdata->{d}, $nb,
				   $cn, $cd, $nc);
	  for (my $i=0; $i<=$nc; $i++) {
	      my ($gcd, $n, $d) = Math::Cephes::euclid($cn->[$i], $cd->[$i]);
	      push @$c, ($aref eq 'Math::Fraction' ?
			 Math::Fraction->new($n, $d) :
		       Math::Cephes::Fraction->new($n, $d) );
	  }
	  last SWITCH;
      };
      $atype eq 'cmplx' and do {
	  my $dc = $na + $nb + 3;
	  my $cr = [split //, 0 x $dc];
	  my $ci = [split //, 0 x $dc];
	  $nc = Math::Cephes::cpmul_wrap($adata->{r}, $adata->{i}, $na+1,
					 $bdata->{r}, $bdata->{i}, $nb+1,
					 $cr, $ci, $dc);
	  $cr = [ @{$cr}[0..$nc] ];
	  $ci = [ @{$ci}[0..$nc] ];
	  for (my $i=0; $i<=$nc; $i++) {
	      push @$c, ($aref eq 'Math::Complex' ?
		       Math::Complex->make($cr->[$i], $ci->[$i]) :
		       Math::Cephes::Complex->new($cr->[$i], $ci->[$i]) );
	  }
	  last SWITCH;
      };
      $nc = $na + $nb + 1;
      $c = [split //, 0 x $nc];
    Math::Cephes::polmul($adata, $na, $bdata, $nb, $c);
  }
    return wantarray ? (Math::Cephes::Polynomial->new($c), $nc) :
      Math::Cephes::Polynomial->new($c);
}

sub div {
    my ($self, $b) = @_;
    my ($atype, $aref, $adata, $na) =
	($self->{type}, $self->{ref}, $self->{data}, $self->{n});
    my ($btype, $bref, $bdata, $nb) =
	ref($b) eq 'Math::Cephes::Polynomial' ?
	    ($b->{type}, $b->{ref}, $b->{data}, $b->{n}) :
		get_data($b, $aref);
    my $c = [];
    my $nc;
  SWITCH: {
      $atype eq 'fract' and do {
	  $nc = $MAXPOL;
	  my $cn = [split //, 0 x ($nc+1)];
	  my $cd = [split //, 0 x ($nc+1)];
	Math::Cephes::fpoldiv_wrap($adata->{n}, $adata->{d}, $na,
				   $bdata->{n}, $bdata->{d}, $nb,
				   $cn, $cd, $nc);
	  for (my $i=0; $i<=$nc; $i++) {
	      my ($gcd, $n, $d) = Math::Cephes::euclid($cn->[$i], $cd->[$i]);
	      push @$c, ($aref eq 'Math::Fraction' ?
			 Math::Fraction->new($n, $d) :
		       Math::Cephes::Fraction->new($n, $d) );
	  }
	  last SWITCH;
      };
      $atype eq 'cmplx' and do {
	  die "Cannot do complex division";
	  last SWITCH;
      };
      $nc = $MAXPOL;
      $c = [split //, 0 x ($nc+1)];
    Math::Cephes::poldiv($adata, $na, $bdata, $nb, $c);
  }
    return wantarray ? (Math::Cephes::Polynomial->new($c), $nc) :
      Math::Cephes::Polynomial->new($c);
}

sub clr {
  my $self = shift;
  my ($atype, $aref, $adata, $na) =
	($self->{type}, $self->{ref}, $self->{data}, $self->{n});
  set_max() unless $flag;
  my $n = shift || $na;
  $n = $na if $n > $na;
 SWITCH: {
      $atype eq 'fract' and do {
	  for (my $i=0; $i<=$n; $i++) {
	    $self->{data}->{n}->[$i] = 0;
	    $self->{data}->{d}->[$i] = 1;
	  }
	  last SWITCH;
      };
      $atype eq 'cmplx' and do {
	  for (my $i=0; $i<=$n; $i++) {
	    $self->{data}->{r}->[$i] = 0;
	    $self->{data}->{i}->[$i] = 0;
	  }
	  last SWITCH;
      };
      for (my $i=0; $i<=$n; $i++) {
	$self->{data}->[$i] = 0;
      }
   }
}

sub sbt {
  my ($self, $b) = @_;
  my ($atype, $aref, $adata, $na) =
    ($self->{type}, $self->{ref}, $self->{data}, $self->{n});
  my ($btype, $bref, $bdata, $nb) =
    ref($b) eq 'Math::Cephes::Polynomial' ?
      ($b->{type}, $b->{ref}, $b->{data}, $b->{n}) :
	get_data($b, $aref);
  set_max() unless $flag;
  my $c = [];
  my $nc;
 SWITCH: {
    $atype eq 'fract' and do {
      $nc = ($na+1)*($nb+1);
      my $cn = [split //, 0 x ($nc+1)];
      my $cd = [split //, 0 x ($nc+1)];
      Math::Cephes::fpolsbt_wrap($bdata->{n}, $bdata->{d}, $nb,
				 $adata->{n}, $adata->{d}, $na,
				 $cn, $cd, $nc);
      $nc = $na * $nb;
      for (my $i=0; $i<=$nc; $i++) {
	my ($gcd, $n, $d) = Math::Cephes::euclid($cn->[$i], $cd->[$i]);
	push @$c, ($aref eq 'Math::Fraction' ?
		   Math::Fraction->new($n, $d) :
		   Math::Cephes::Fraction->new($n, $d) );
      }
      last SWITCH;
    };
    $atype eq 'cmplx' and do {
      die "Cannot do complex substitution";
      last SWITCH;
    };
    $nc = ($na+1)*($nb+1);
    $c = [split //, 0 x $nc];
    Math::Cephes::polsbt($bdata, $nb, $adata, $na, $c);
    $nc = $na*$nb;
    $c = [@$c[0..$nc]];
  }
  return wantarray ? (Math::Cephes::Polynomial->new($c), $nc) :
    Math::Cephes::Polynomial->new($c);
}

sub set_max {
  Math::Cephes::polini($MAXPOL);
    $flag = 1;
}

sub set_fmax {
  Math::Cephes::fpolini($FMAXPOL);
    $fflag = 1;
}

sub eval {
  my $self = shift;
  my $x = 0 || shift;
  my ($atype, $aref, $adata, $na) =
    ($self->{type}, $self->{ref}, $self->{data}, $self->{n});
  my $y;
 SWITCH: {
    $atype eq 'fract' and do {
      my $xref = ref($x);
      $y = Math::Cephes::Fraction->new(0, 1);
    FRACT: {
	not $xref and do {
	  $x = Math::Cephes::Fraction->new($x, 1);
	  last FRACT;
	};
	$xref eq 'Math::Cephes::Fraction' and do {
	  last FRACT;
	};
	$xref eq 'Math::Fraction' and do {
	  $x = Math::Cephes::Fraction->new($x->{frac}->[0], $x->{frac}->[1]);
	  last FRACT;
	};
	die "Unknown data type '$xref' for x";
      }
      Math::Cephes::fpoleva_wrap($adata->{n}, $adata->{d}, $na, $x, $y);
      $y = Math::Fraction->new($y->n, $y->d) if $aref eq 'Math::Fraction';
      last SWITCH;
    };
    $atype eq 'cmplx' and do {
      my $r = Math::Cephes::poleva($adata->{r}, $na, $x);
      my $i = Math::Cephes::poleva($adata->{i}, $na, $x);
      $y = $aref eq 'Math::Complex' ?
	Math::Complex->make($r, $i) :
	    Math::Cephes::Complex->new($r, $i);
      last SWITCH;
    };
    $y = Math::Cephes::poleva($adata, $na, $x);
  }
  return $y;
}

sub fract_to_real {
    my $in = shift;
    my $a = [];
    my $n = scalar @{$in->{n}} - 1;
    for (my $i=0; $i<=$n; $i++) {
	push @$a, $in->{n}->[$i] / $in->{d}->[$i];
    }
    return $a;
}

sub atn {
    my ($self, $bin) = @_;
    my $type = $self->{type};
    die "Cannot take the atan of a complex polynomial"
	if $type eq 'cmplx';
    my ($a, $b);
    my ($atype, $aref, $adata, $na) =
	($self->{type}, $self->{ref}, $self->{data}, $self->{n});
    die "Cannot take the atan of a complex polynomial"
	if $atype eq 'cmplx';
    $a = $atype eq 'fract' ? fract_to_real($adata) : $adata;

    my ($btype, $bref, $bdata, $nb) =
	ref($bin) eq 'Math::Cephes::Polynomial' ?
	    ($bin->{type}, $bin->{ref}, $bin->{data}, $bin->{n}) :
		get_data($bin);

    die "Cannot take the atan of a complex polynomial"
	if $btype eq 'cmplx';
    $b = $btype eq 'fract' ? fract_to_real($bdata) : $bdata;

    my $c = [split //, 0 x ($MAXPOL+1)];
  Math::Cephes::polatn($a, $b, $c, 16);
    return Math::Cephes::Polynomial->new($c);
}

sub sqt {
  my $self = shift;
  my $type = $self->{type};
  die "Cannot take the sqrt of a complex polynomial"
    if $type eq 'cmplx';
  my $a = $type eq 'fract' ? fract_to_real($self->{data}) : $self->coef;
  my $b = [split //, 0 x ($MAXPOL+1)];
  Math::Cephes::polsqt($a, $b, 16);
  return Math::Cephes::Polynomial->new($b);
}

sub sin {
  my $self = shift;
  my $type = $self->{type};
  die "Cannot take the sin of a complex polynomial"
    if $type eq 'cmplx';
  my $a = $type eq 'fract' ? fract_to_real($self->{data}) : $self->coef;
  my $b = [split //, 0 x ($MAXPOL+1)];
  Math::Cephes::polsin($a, $b, 16);
  return Math::Cephes::Polynomial->new($b);
}

sub cos {
  my $self = shift;
  my $type = $self->{type};
  die "Cannot take the cos of a complex polynomial"
    if $type eq 'cmplx';
  my $a = $type eq 'fract' ? fract_to_real($self->{data}) : $self->coef;
  my $b = [split //, 0 x ($MAXPOL+1)];
  Math::Cephes::polcos($a, $b, 16);
  return Math::Cephes::Polynomial->new($b);
}

sub rts {
  my $self = shift;
  my ($atype, $aref, $adata, $na) =
    ($self->{type}, $self->{ref}, $self->{data}, $self->{n});
  my ($a, $b, $ret);
  my $cof = [split //, 0 x ($na+1)];
  my $r = [split //, 0 x ($na+1)];
  my $i = [split //, 0 x ($na+1)];
 SWITCH: {
    $atype eq 'fract' and do {
      $adata = fract_to_real($adata);
      $ret = Math::Cephes::polrt_wrap($adata, $cof, $na, $r, $i);
      for (my $j=0; $j<$na; $j++) {
	push @$b,
	  Math::Cephes::Complex->new($r->[$j], $i->[$j]);
      }
      last SWITCH;
    };
    $atype eq 'cmplx' and do {
      die "Cannot do complex root finding";
      last SWITCH;
    };
    $ret = Math::Cephes::polrt_wrap($adata, $cof, $na, $r, $i);
    for (my $j=0; $j<$na; $j++) {
      push @$b,
	Math::Cephes::Complex->new($r->[$j], $i->[$j]);
    }
  }
  return wantarray ? ($ret, $b) : $b;
}

1;

__END__