📄 complex.pm
字号:
my ($z1, $z2, $inverted) = @_; if ($inverted) { return 1 if $z1 == 0 || $z2 == 1; return 0 if $z2 == 0 && Re($z1) > 0; } else { return 1 if $z2 == 0 || $z1 == 1; return 0 if $z1 == 0 && Re($z2) > 0; } my $w = $inverted ? &exp($z1 * &log($z2)) : &exp($z2 * &log($z1)); # If both arguments cartesian, return cartesian, else polar. return $z1->{c_dirty} == 0 && (not ref $z2 or $z2->{c_dirty} == 0) ? cplx(@{$w->_cartesian}) : $w;}## (_spaceship)## Computes z1 <=> z2.# Sorts on the real part first, then on the imaginary part. Thus 2-4i < 3+8i.#sub _spaceship { my ($z1, $z2, $inverted) = @_; my ($re1, $im1) = ref $z1 ? @{$z1->_cartesian} : ($z1, 0); my ($re2, $im2) = ref $z2 ? @{$z2->_cartesian} : ($z2, 0); my $sgn = $inverted ? -1 : 1; return $sgn * ($re1 <=> $re2) if $re1 != $re2; return $sgn * ($im1 <=> $im2);}## (_numeq)## Computes z1 == z2.## (Required in addition to _spaceship() because of NaNs.)sub _numeq { my ($z1, $z2, $inverted) = @_; my ($re1, $im1) = ref $z1 ? @{$z1->_cartesian} : ($z1, 0); my ($re2, $im2) = ref $z2 ? @{$z2->_cartesian} : ($z2, 0); return $re1 == $re2 && $im1 == $im2 ? 1 : 0;}## (_negate)## Computes -z.#sub _negate { my ($z) = @_; if ($z->{c_dirty}) { my ($r, $t) = @{$z->_polar}; $t = ($t <= 0) ? $t + pi : $t - pi; return (ref $z)->emake($r, $t); } my ($re, $im) = @{$z->_cartesian}; return (ref $z)->make(-$re, -$im);}## (_conjugate)## Compute complex's _conjugate.#sub _conjugate { my ($z) = @_; if ($z->{c_dirty}) { my ($r, $t) = @{$z->_polar}; return (ref $z)->emake($r, -$t); } my ($re, $im) = @{$z->_cartesian}; return (ref $z)->make($re, -$im);}## (abs)## Compute or set complex's norm (rho).#sub abs { my ($z, $rho) = @_; unless (ref $z) { if (@_ == 2) { $_[0] = $_[1]; } else { return CORE::abs($z); } } if (defined $rho) { $z->{'polar'} = [ $rho, ${$z->_polar}[1] ]; $z->{p_dirty} = 0; $z->{c_dirty} = 1; return $rho; } else { return ${$z->_polar}[0]; }}sub _theta { my $theta = $_[0]; if ($$theta > pi()) { $$theta -= pi2 } elsif ($$theta <= -pi()) { $$theta += pi2 }}## arg## Compute or set complex's argument (theta).#sub arg { my ($z, $theta) = @_; return $z unless ref $z; if (defined $theta) { _theta(\$theta); $z->{'polar'} = [ ${$z->_polar}[0], $theta ]; $z->{p_dirty} = 0; $z->{c_dirty} = 1; } else { $theta = ${$z->_polar}[1]; _theta(\$theta); } return $theta;}## (sqrt)## Compute sqrt(z).## It is quite tempting to use wantarray here so that in list context# sqrt() would return the two solutions. This, however, would# break things like## print "sqrt(z) = ", sqrt($z), "\n";## The two values would be printed side by side without no intervening# whitespace, quite confusing.# Therefore if you want the two solutions use the root().#sub sqrt { my ($z) = @_; my ($re, $im) = ref $z ? @{$z->_cartesian} : ($z, 0); return $re < 0 ? cplx(0, CORE::sqrt(-$re)) : CORE::sqrt($re) if $im == 0; my ($r, $t) = @{$z->_polar}; return (ref $z)->emake(CORE::sqrt($r), $t/2);}## cbrt## Compute cbrt(z) (cubic root).## Why are we not returning three values? The same answer as for sqrt().#sub cbrt { my ($z) = @_; return $z < 0 ? -CORE::exp(CORE::log(-$z)/3) : ($z > 0 ? CORE::exp(CORE::log($z)/3): 0) unless ref $z; my ($r, $t) = @{$z->_polar}; return 0 if $r == 0; return (ref $z)->emake(CORE::exp(CORE::log($r)/3), $t/3);}## _rootbad## Die on bad root.#sub _rootbad { my $mess = "Root '$_[0]' illegal, root rank must be positive integer.\n"; my @up = caller(1); $mess .= "Died at $up[1] line $up[2].\n"; die $mess;}## root## Computes all nth root for z, returning an array whose size is n.# `n' must be a positive integer.## The roots are given by (for k = 0..n-1):## z^(1/n) = r^(1/n) (cos ((t+2 k pi)/n) + i sin ((t+2 k pi)/n))#sub root { my ($z, $n, $k) = @_; _rootbad($n) if ($n < 1 or int($n) != $n); my ($r, $t) = ref $z ? @{$z->_polar} : (CORE::abs($z), $z >= 0 ? 0 : pi); my $theta_inc = pi2 / $n; my $rho = $r ** (1/$n); my $cartesian = ref $z && $z->{c_dirty} == 0; if (@_ == 2) { my @root; for (my $i = 0, my $theta = $t / $n; $i < $n; $i++, $theta += $theta_inc) { my $w = cplxe($rho, $theta); # Yes, $cartesian is loop invariant. push @root, $cartesian ? cplx(@{$w->_cartesian}) : $w; } return @root; } elsif (@_ == 3) { my $w = cplxe($rho, $t / $n + $k * $theta_inc); return $cartesian ? cplx(@{$w->_cartesian}) : $w; }}## Re## Return or set Re(z).#sub Re { my ($z, $Re) = @_; return $z unless ref $z; if (defined $Re) { $z->{'cartesian'} = [ $Re, ${$z->_cartesian}[1] ]; $z->{c_dirty} = 0; $z->{p_dirty} = 1; } else { return ${$z->_cartesian}[0]; }}## Im## Return or set Im(z).#sub Im { my ($z, $Im) = @_; return 0 unless ref $z; if (defined $Im) { $z->{'cartesian'} = [ ${$z->_cartesian}[0], $Im ]; $z->{c_dirty} = 0; $z->{p_dirty} = 1; } else { return ${$z->_cartesian}[1]; }}## rho## Return or set rho(w).#sub rho { Math::Complex::abs(@_);}## theta## Return or set theta(w).#sub theta { Math::Complex::arg(@_);}## (exp)## Computes exp(z).#sub exp { my ($z) = @_; my ($x, $y) = @{$z->_cartesian}; return (ref $z)->emake(CORE::exp($x), $y);}## _logofzero## Die on logarithm of zero.#sub _logofzero { my $mess = "$_[0]: Logarithm of zero.\n"; if (defined $_[1]) { $mess .= "(Because in the definition of $_[0], the argument "; $mess .= "$_[1] " unless ($_[1] eq '0'); $mess .= "is 0)\n"; } my @up = caller(1); $mess .= "Died at $up[1] line $up[2].\n"; die $mess;}## (log)## Compute log(z).#sub log { my ($z) = @_; unless (ref $z) { _logofzero("log") if $z == 0; return $z > 0 ? CORE::log($z) : cplx(CORE::log(-$z), pi); } my ($r, $t) = @{$z->_polar}; _logofzero("log") if $r == 0; if ($t > pi()) { $t -= pi2 } elsif ($t <= -pi()) { $t += pi2 } return (ref $z)->make(CORE::log($r), $t);}## ln## Alias for log().#sub ln { Math::Complex::log(@_) }## log10## Compute log10(z).#sub log10 { return Math::Complex::log($_[0]) * _uplog10;}## logn## Compute logn(z,n) = log(z) / log(n)#sub logn { my ($z, $n) = @_; $z = cplx($z, 0) unless ref $z; my $logn = $LOGN{$n}; $logn = $LOGN{$n} = CORE::log($n) unless defined $logn; # Cache log(n) return &log($z) / $logn;}## (cos)## Compute cos(z) = (exp(iz) + exp(-iz))/2.#sub cos { my ($z) = @_; return CORE::cos($z) unless ref $z; my ($x, $y) = @{$z->_cartesian}; my $ey = CORE::exp($y); my $sx = CORE::sin($x); my $cx = CORE::cos($x); my $ey_1 = $ey ? 1 / $ey : $Inf; return (ref $z)->make($cx * ($ey + $ey_1)/2, $sx * ($ey_1 - $ey)/2);}## (sin)## Compute sin(z) = (exp(iz) - exp(-iz))/2.#sub sin { my ($z) = @_; return CORE::sin($z) unless ref $z; my ($x, $y) = @{$z->_cartesian}; my $ey = CORE::exp($y); my $sx = CORE::sin($x); my $cx = CORE::cos($x); my $ey_1 = $ey ? 1 / $ey : $Inf; return (ref $z)->make($sx * ($ey + $ey_1)/2, $cx * ($ey - $ey_1)/2);}## tan## Compute tan(z) = sin(z) / cos(z).#sub tan { my ($z) = @_; my $cz = &cos($z); _divbyzero "tan($z)", "cos($z)" if $cz == 0; return &sin($z) / $cz;}## sec## Computes the secant sec(z) = 1 / cos(z).#sub sec { my ($z) = @_; my $cz = &cos($z); _divbyzero "sec($z)", "cos($z)" if ($cz == 0); return 1 / $cz;}## csc## Computes the cosecant csc(z) = 1 / sin(z).#sub csc { my ($z) = @_; my $sz = &sin($z); _divbyzero "csc($z)", "sin($z)" if ($sz == 0); return 1 / $sz;}## cosec## Alias for csc().#sub cosec { Math::Complex::csc(@_) }## cot## Computes cot(z) = cos(z) / sin(z).#sub cot { my ($z) = @_; my $sz = &sin($z); _divbyzero "cot($z)", "sin($z)" if ($sz == 0); return &cos($z) / $sz;}## cotan## Alias for cot().#sub cotan { Math::Complex::cot(@_) }## acos## Computes the arc cosine acos(z) = -i log(z + sqrt(z*z-1)).#sub acos { my $z = $_[0]; return CORE::atan2(CORE::sqrt(1-$z*$z), $z) if (! ref $z) && CORE::abs($z) <= 1; $z = cplx($z, 0) unless ref $z; my ($x, $y) = @{$z->_cartesian}; return 0 if $x == 1 && $y == 0; my $t1 = CORE::sqrt(($x+1)*($x+1) + $y*$y); my $t2 = CORE::sqrt(($x-1)*($x-1) + $y*$y); my $alpha = ($t1 + $t2)/2; my $beta = ($t1 - $t2)/2; $alpha = 1 if $alpha < 1; if ($beta > 1) { $beta = 1 } elsif ($beta < -1) { $beta = -1 } my $u = CORE::atan2(CORE::sqrt(1-$beta*$beta), $beta); my $v = CORE::log($alpha + CORE::sqrt($alpha*$alpha-1)); $v = -$v if $y > 0 || ($y == 0 && $x < -1); return (ref $z)->make($u, $v);}## asin## Computes the arc sine asin(z) = -i log(iz + sqrt(1-z*z)).#sub asin { my $z = $_[0]; return CORE::atan2($z, CORE::sqrt(1-$z*$z)) if (! ref $z) && CORE::abs($z) <= 1; $z = cplx($z, 0) unless ref $z; my ($x, $y) = @{$z->_cartesian}; return 0 if $x == 0 && $y == 0; my $t1 = CORE::sqrt(($x+1)*($x+1) + $y*$y); my $t2 = CORE::sqrt(($x-1)*($x-1) + $y*$y); my $alpha = ($t1 + $t2)/2; my $beta = ($t1 - $t2)/2; $alpha = 1 if $alpha < 1; if ($beta > 1) { $beta = 1 } elsif ($beta < -1) { $beta = -1 } my $u = CORE::atan2($beta, CORE::sqrt(1-$beta*$beta)); my $v = -CORE::log($alpha + CORE::sqrt($alpha*$alpha-1)); $v = -$v if $y > 0 || ($y == 0 && $x < -1); return (ref $z)->make($u, $v);}## atan## Computes the arc tangent atan(z) = i/2 log((i+z) / (i-z)).
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -