version 1.1, 2013/10/29 21:01:21
|
version 1.2, 2019/11/10 23:52:07
|
Line 4
|
Line 4
|
# -- Jarkko Hietaniemi Since Mar 1997 |
# -- Jarkko Hietaniemi Since Mar 1997 |
# -- Daniel S. Lewart Since Sep 1997 |
# -- Daniel S. Lewart Since Sep 1997 |
# |
# |
# -- Stuart Raeburn: renamed package as LCMathComplex Oct 2013 |
# -- Stuart Raeburn: renamed package (rev. 1.55) as LCMathComplex Oct 2013 |
|
# renamed package (rev. 1.59_01) as LCMathComplex Nov 2019 |
# with minor changes to allow use in Safe Space |
# with minor changes to allow use in Safe Space |
# |
# |
|
|
package LONCAPA::LCMathComplex; |
package LONCAPA::LCMathComplex; |
|
|
use vars qw($VERSION @ISA @EXPORT @EXPORT_OK %EXPORT_TAGS $Inf $ExpInf); |
{ use 5.006; } |
|
use strict; |
|
|
|
our $VERSION = 1.59_01; |
|
|
$VERSION = 1.55; |
our ($Inf, $ExpInf); |
|
our ($vax_float, $has_inf, $has_nan); |
|
|
BEGIN { |
BEGIN { |
my %DBL_MAX = |
$vax_float = (pack("d",1) =~ /^[\x80\x10]\x40/); |
|
$has_inf = !$vax_float; |
|
$has_nan = !$vax_float; |
|
|
|
unless ($has_inf) { |
|
# For example in vax, there is no Inf, |
|
# and just mentioning the DBL_MAX (1.70141183460469229e+38) |
|
# causes SIGFPE. |
|
|
|
# These are pretty useless without a real infinity, |
|
# but setting them makes for less warnings about their |
|
# undefined values. |
|
$Inf = "Inf"; |
|
$ExpInf = "Inf"; |
|
return; |
|
} |
|
|
|
my %DBL_MAX = # These are IEEE 754 maxima. |
( |
( |
4 => '1.70141183460469229e+38', |
4 => '1.70141183460469229e+38', |
8 => '1.7976931348623157e+308', |
8 => '1.7976931348623157e+308', |
Line 25 BEGIN {
|
Line 47 BEGIN {
|
12 => '1.1897314953572317650857593266280070162E+4932', |
12 => '1.1897314953572317650857593266280070162E+4932', |
16 => '1.1897314953572317650857593266280070162E+4932', |
16 => '1.1897314953572317650857593266280070162E+4932', |
); |
); |
|
|
my $nvsize = 8; |
my $nvsize = 8; |
die "LONCAPA::LCMathComplex: Could not figure out nvsize\n" |
die "Math::Complex: Could not figure out nvsize\n" |
unless defined $nvsize; |
unless defined $nvsize; |
die "LONCAPA::LCMathComplex: Cannot not figure out max nv (nvsize = $nvsize)\n" |
die "LONCAPA::LCMathComplex: Cannot not figure out max nv (nvsize = $nvsize)\n" |
unless defined $DBL_MAX{$nvsize}; |
unless defined $DBL_MAX{$nvsize}; |
Line 37 BEGIN {
|
Line 60 BEGIN {
|
if ($^O eq 'unicosmk') { |
if ($^O eq 'unicosmk') { |
$Inf = $DBL_MAX; |
$Inf = $DBL_MAX; |
} else { |
} else { |
local $SIG{FPE} = { }; |
local $SIG{FPE} = sub { }; |
local $!; |
local $!; |
# We do want an arithmetic overflow, Inf INF inf Infinity. |
# We do want an arithmetic overflow, Inf INF inf Infinity. |
for my $t ( |
for my $t ( |
Line 56 BEGIN {
|
Line 79 BEGIN {
|
$Inf = $i; |
$Inf = $i; |
last; |
last; |
} |
} |
} |
} |
$Inf = $DBL_MAX unless defined $Inf; # Oh well, close enough. |
$Inf = $DBL_MAX unless defined $Inf; # Oh well, close enough. |
die "LONCAPA::LCMathComplex: Could not get Infinity" |
die "LONCAPA::LCMathComplex: Could not get Infinity" |
unless $Inf > $BIGGER_THAN_THIS; |
unless $Inf > $BIGGER_THAN_THIS; |
$ExpInf = exp(99999); |
$ExpInf = eval 'exp(99999)'; |
} |
} |
# print "# On this machine, Inf = '$Inf'\n"; |
# print "# On this machine, Inf = '$Inf'\n"; |
} |
} |
|
|
use strict; |
use warnings; |
|
no warnings 'syntax'; # To avoid the (_) warnings. |
|
|
my $i; |
my $i; |
my %LOGN; |
my %LOGN; |
Line 76 my $gre = qr'\s*([\+\-]?(?:(?:(?:\d+(?:_
|
Line 100 my $gre = qr'\s*([\+\-]?(?:(?:(?:\d+(?:_
|
|
|
require Exporter; |
require Exporter; |
|
|
@ISA = qw(Exporter); |
our @ISA = qw(Exporter); |
|
|
my @trig = qw( |
my @trig = qw( |
pi |
pi |
Line 90 my @trig = qw(
|
Line 114 my @trig = qw(
|
acsch acosech asech acoth acotanh |
acsch acosech asech acoth acotanh |
); |
); |
|
|
@EXPORT = (qw( |
our @EXPORT = (qw( |
i Re Im rho theta arg |
i Re Im rho theta arg |
sqrt log ln |
sqrt log ln |
log10 logn cbrt root |
log10 logn cbrt root |
Line 101 my @trig = qw(
|
Line 125 my @trig = qw(
|
|
|
my @pi = qw(pi pi2 pi4 pip2 pip4 Inf); |
my @pi = qw(pi pi2 pi4 pip2 pip4 Inf); |
|
|
@EXPORT_OK = @pi; |
our @EXPORT_OK = @pi; |
|
|
%EXPORT_TAGS = ( |
our %EXPORT_TAGS = ( |
'trig' => [@trig], |
'trig' => [@trig], |
'pi' => [@pi], |
'pi' => [@pi], |
); |
); |
|
|
use overload |
use overload |
|
'=' => \&_copy, |
|
'+=' => \&_plus, |
'+' => \&_plus, |
'+' => \&_plus, |
|
'-=' => \&_minus, |
'-' => \&_minus, |
'-' => \&_minus, |
|
'*=' => \&_multiply, |
'*' => \&_multiply, |
'*' => \&_multiply, |
|
'/=' => \&_divide, |
'/' => \&_divide, |
'/' => \&_divide, |
|
'**=' => \&_power, |
'**' => \&_power, |
'**' => \&_power, |
'==' => \&_numeq, |
'==' => \&_numeq, |
'<=>' => \&_spaceship, |
'<=>' => \&_spaceship, |
Line 124 use overload
|
Line 154 use overload
|
'log' => \&log, |
'log' => \&log, |
'sin' => \&sin, |
'sin' => \&sin, |
'cos' => \&cos, |
'cos' => \&cos, |
'tan' => \&tan, |
|
'atan2' => \&atan2, |
'atan2' => \&atan2, |
'""' => \&_stringify; |
'""' => \&_stringify; |
|
|
Line 165 sub _make {
|
Line 194 sub _make {
|
|
|
if (defined $p) { |
if (defined $p) { |
$p =~ s/^\+//; |
$p =~ s/^\+//; |
$p =~ s/^(-?)inf$/"${1}9**9**9"/e; |
$p =~ s/^(-?)inf$/"${1}9**9**9"/e if $has_inf; |
$q =~ s/^\+//; |
$q =~ s/^\+//; |
$q =~ s/^(-?)inf$/"${1}9**9**9"/e; |
$q =~ s/^(-?)inf$/"${1}9**9**9"/e if $has_inf; |
} |
} |
|
|
return ($p, $q); |
return ($p, $q); |
Line 190 sub _emake {
|
Line 219 sub _emake {
|
if (defined $p) { |
if (defined $p) { |
$p =~ s/^\+//; |
$p =~ s/^\+//; |
$q =~ s/^\+//; |
$q =~ s/^\+//; |
$p =~ s/^(-?)inf$/"${1}9**9**9"/e; |
$p =~ s/^(-?)inf$/"${1}9**9**9"/e if $has_inf; |
$q =~ s/^(-?)inf$/"${1}9**9**9"/e; |
$q =~ s/^(-?)inf$/"${1}9**9**9"/e if $has_inf; |
} |
} |
|
|
return ($p, $q); |
return ($p, $q); |
} |
} |
|
|
|
sub _copy { |
|
my $self = shift; |
|
my $clone = {%$self}; |
|
if ($self->{'cartesian'}) { |
|
$clone->{'cartesian'} = [@{$self->{'cartesian'}}]; |
|
} |
|
if ($self->{'polar'}) { |
|
$clone->{'polar'} = [@{$self->{'polar'}}]; |
|
} |
|
bless $clone,__PACKAGE__; |
|
return $clone; |
|
} |
|
|
# |
# |
# ->make |
# ->make |
# |
# |
Line 610 sub _conjugate {
|
Line 652 sub _conjugate {
|
# Compute or set complex's norm (rho). |
# Compute or set complex's norm (rho). |
# |
# |
sub abs { |
sub abs { |
my ($z, $rho) = @_; |
my ($z, $rho) = @_ ? @_ : $_; |
unless (ref $z) { |
unless (ref $z) { |
if (@_ == 2) { |
if (@_ == 2) { |
$_[0] = $_[1]; |
$_[0] = $_[1]; |
Line 671 sub arg {
|
Line 713 sub arg {
|
# Therefore if you want the two solutions use the root(). |
# Therefore if you want the two solutions use the root(). |
# |
# |
sub sqrt { |
sub sqrt { |
my ($z) = @_; |
my ($z) = @_ ? $_[0] : $_; |
my ($re, $im) = ref $z ? @{$z->_cartesian} : ($z, 0); |
my ($re, $im) = ref $z ? @{$z->_cartesian} : ($z, 0); |
return $re < 0 ? cplx(0, CORE::sqrt(-$re)) : CORE::sqrt($re) |
return $re < 0 ? cplx(0, CORE::sqrt(-$re)) : CORE::sqrt($re) |
if $im == 0; |
if $im == 0; |
Line 804 sub theta {
|
Line 846 sub theta {
|
# Computes exp(z). |
# Computes exp(z). |
# |
# |
sub exp { |
sub exp { |
my ($z) = @_; |
my ($z) = @_ ? @_ : $_; |
my ($x, $y) = @{$z->_cartesian}; |
return CORE::exp($z) unless ref $z; |
return (ref $z)->emake(CORE::exp($x), $y); |
my ($x, $y) = @{$z->_cartesian}; |
|
return (ref $z)->emake(CORE::exp($x), $y); |
} |
} |
|
|
# |
# |
Line 836 sub _logofzero {
|
Line 879 sub _logofzero {
|
# Compute log(z). |
# Compute log(z). |
# |
# |
sub log { |
sub log { |
my ($z) = @_; |
my ($z) = @_ ? @_ : $_; |
unless (ref $z) { |
unless (ref $z) { |
_logofzero("log") if $z == 0; |
_logofzero("log") if $z == 0; |
return $z > 0 ? CORE::log($z) : cplx(CORE::log(-$z), pi); |
return $z > 0 ? CORE::log($z) : cplx(CORE::log(-$z), pi); |
Line 884 sub logn {
|
Line 927 sub logn {
|
# Compute cos(z) = (exp(iz) + exp(-iz))/2. |
# Compute cos(z) = (exp(iz) + exp(-iz))/2. |
# |
# |
sub cos { |
sub cos { |
my ($z) = @_; |
my ($z) = @_ ? @_ : $_; |
return CORE::cos($z) unless ref $z; |
return CORE::cos($z) unless ref $z; |
my ($x, $y) = @{$z->_cartesian}; |
my ($x, $y) = @{$z->_cartesian}; |
my $ey = CORE::exp($y); |
my $ey = CORE::exp($y); |
Line 901 sub cos {
|
Line 944 sub cos {
|
# Compute sin(z) = (exp(iz) - exp(-iz))/2. |
# Compute sin(z) = (exp(iz) - exp(-iz))/2. |
# |
# |
sub sin { |
sub sin { |
my ($z) = @_; |
my ($z) = @_ ? @_ : $_; |
return CORE::sin($z) unless ref $z; |
return CORE::sin($z) unless ref $z; |
my ($x, $y) = @{$z->_cartesian}; |
my ($x, $y) = @{$z->_cartesian}; |
my $ey = CORE::exp($y); |
my $ey = CORE::exp($y); |
Line 1518 sub _stringify_polar {
|
Line 1561 sub _stringify_polar {
|
|
|
if (defined $format) { |
if (defined $format) { |
$r = sprintf($format, $r); |
$r = sprintf($format, $r); |
$theta = sprintf($format, $theta) unless defined $theta; |
$theta = sprintf($format, $t) unless defined $theta; |
} else { |
} else { |
$theta = $t unless defined $theta; |
$theta = $t unless defined $theta; |
} |
} |
Line 1553 Derived from Math::Complex.
|
Line 1596 Derived from Math::Complex.
|
|
|
Modified for use in Safe Space in LON-CAPA by removing the dependency on |
Modified for use in Safe Space in LON-CAPA by removing the dependency on |
Config.pm introduced in rev. 1.51 (which contains calls that are disallowed |
Config.pm introduced in rev. 1.51 (which contains calls that are disallowed |
in Safe Space). |
in Safe Space). In addition, Scalar::Util::set_prototype() is not used for |
|
abs(), cos(), exp(), log(), sin(), and sqrt(), to avoid warnings in logs: |
|
of type: "Use of uninitialized value" (Config.pm line 62). |
|
|
In this LON-CAPA-specific version, the following code changes were made. |
In this LON-CAPA-specific version, the following code changes were made. |
|
|
15,16d17 |
15,16d17 |
< use Config; |
< use Config; |
< |
< |
29,31c30 |
49,51c50 |
< my $nvsize = $Config{nvsize} || |
< my $nvsize = $Config{nvsize} || |
< ($Config{uselongdouble} && $Config{longdblsize}) || |
< ($Config{uselongdouble} && $Config{longdblsize}) || |
< $Config{doublesize}; |
< $Config{doublesize}; |
--- |
--- |
> my $nvsize = 8; |
> my $nvsize = 8; |
|
|
|
91,92d89 |
|
< use Scalar::Util qw(set_prototype); |
|
< |
|
96,109d92 |
|
< BEGIN { |
|
< # For certain functions that we override, in 5.10 or better |
|
< # we can set a smarter prototype that will handle the lexical $_ |
|
< # (also a 5.10+ feature). |
|
< if ($] >= 5.010000) { |
|
< set_prototype \&abs, '_'; |
|
< set_prototype \&cos, '_'; |
|
< set_prototype \&exp, '_'; |
|
< set_prototype \&log, '_'; |
|
< set_prototype \&sin, '_'; |
|
< set_prototype \&sqrt, '_'; |
|
< } |
|
< } |
|
|
Note: the value assigned to $nvsize is 8 by default. |
Note: the value assigned to $nvsize is 8 by default. |
|
|
Whenever ./UPDATE is run to install or update LON-CAPA, the code which |
Whenever ./UPDATE is run to install or update LON-CAPA, the code which |
Line 2060 messages like the following
|
Line 2123 messages like the following
|
|
|
LONCAPA::LCMathComplex::make: Cannot take real part of ... |
LONCAPA::LCMathComplex::make: Cannot take real part of ... |
LONCAPA::LCMathComplex::make: Cannot take real part of ... |
LONCAPA::LCMathComplex::make: Cannot take real part of ... |
LONCAPA::LCMathComplex:emake: Cannot take rho of ... |
LONCAPA::LCMathComplex::emake: Cannot take rho of ... |
LONCAPA::LCMathComplex::emake: Cannot take theta of ... |
LONCAPA::LCMathComplex::emake: Cannot take theta of ... |
|
|
=head1 BUGS |
=head1 BUGS |
Line 2084 L<Math::Trig>
|
Line 2147 L<Math::Trig>
|
|
|
=head1 AUTHORS |
=head1 AUTHORS |
|
|
Daniel S. Lewart <F<lewart!at!uiuc.edu>> |
Daniel S. Lewart <F<lewart!at!uiuc.edu>>, |
Jarkko Hietaniemi <F<jhi!at!iki.fi>> |
Jarkko Hietaniemi <F<jhi!at!iki.fi>>, |
Raphael Manfredi <F<Raphael_Manfredi!at!pobox.com>> |
Raphael Manfredi <F<Raphael_Manfredi!at!pobox.com>>, |
|
Zefram <zefram@fysh.org> |
|
|
=head1 LICENSE |
=head1 LICENSE |
|
|