use
vars
qw(@EXPORT_OK $VERSION $MAXPOL $FMAXPOL $flag $fflag)
;
$MAXPOL
= 256;
$flag
= 0;
$FMAXPOL
= 256;
$fflag
= 0;
*import
= \
&Exporter::import
;
@EXPORT_OK
=
qw(poly)
;
$VERSION
=
'0.36'
;
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
?
$na
+ 3 :
$nb
+ 3;
$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;