Mercurial > hg > Others > Rakudo
view src/core.c/Complex.pm6 @ 0:c341f82e7ad7 default tip
Rakudo branch in cr.ie.u-ryukyu.ac.jp
author | Shinji KONO <kono@ie.u-ryukyu.ac.jp> |
---|---|
date | Thu, 26 Dec 2019 16:50:27 +0900 |
parents | |
children |
line wrap: on
line source
my class X::Numeric::Real { ... }; my class Complex is Cool does Numeric { has num $.re; has num $.im; method !SET-SELF(Num() \re, Num() \im) { $!re = re; $!im = im; self } proto method new(|) {*} multi method new() { self.new: 0, 0 } multi method new(Real \re, Real \im) { nqp::create(self)!SET-SELF(re, im) } multi method WHICH(Complex:D: --> ValueObjAt:D) { nqp::box_s( nqp::concat( nqp::if( nqp::eqaddr(self.WHAT,Complex), 'Complex|', nqp::concat(nqp::unbox_s(self.^name), '|') ), nqp::concat($!re, nqp::concat('|', $!im)) ), ValueObjAt ) } method reals(Complex:D:) { (self.re, self.im); } method isNaN(Complex:D:) { self.re.isNaN || self.im.isNaN; } method coerce-to-real(Complex:D: $exception-target) { $!im ≅ 0e0 ?? $!re !! Failure.new(X::Numeric::Real.new(target => $exception-target, reason => "imaginary part not zero", source => self)) } multi method Real(Complex:D:) { self.coerce-to-real(Real); } # should probably be eventually supplied by role Numeric method Num(Complex:D:) { self.coerce-to-real(Num).Num; } method Int(Complex:D:) { self.coerce-to-real(Int).Int; } method Rat(Complex:D: $epsilon?) { self.coerce-to-real(Rat).Rat( |($epsilon // Empty) ); } method FatRat(Complex:D: $epsilon?) { self.coerce-to-real(FatRat).FatRat( |($epsilon // Empty) ); } multi method Bool(Complex:D:) { $!re != 0e0 || $!im != 0e0; } method Complex() { self } multi method Str(Complex:D:) { nqp::concat( $!re, nqp::concat( nqp::if(nqp::iseq_i( # we could have negative zero, so stringify nqp::ord(my $im := nqp::p6box_s($!im)),45),'','+'), nqp::concat( $im, nqp::if(nqp::isnanorinf($!im),'\\i','i') ) ) ) } multi method perl(Complex:D:) { '<' ~ self.Str ~ '>'; } method conj(Complex:D:) { Complex.new($.re, -$.im); } method abs(Complex $x:) { nqp::p6box_n(nqp::sqrt_n( nqp::add_n( nqp::mul_n($!re, $!re), nqp::mul_n($!im, $!im), ) )) } method polar() { $.abs, $!im.atan2($!re); } multi method log(Complex:D:) { my Num ($mag, $angle) = self.polar; Complex.new($mag.log, $angle); } method cis(Complex:D:) { self.cos + self.sin*Complex.new(0,1) } method sqrt(Complex:D:) { my Num $abs = self.abs; my Num $re = (($abs + self.re)/2).sqrt; my Num $im = (($abs - self.re)/2).sqrt; Complex.new($re, self.im < 0 ?? -$im !! $im); } multi method exp(Complex:D:) { my Num $mag = $!re.exp; Complex.new($mag * $!im.cos, $mag * $!im.sin); } method roots(Complex:D: Int() $n) { return NaN if $n < 1; return self if $n == 1; for $!re, $!im { return NaN if $_ eq 'Inf' || $_ eq '-Inf' || $_ eq 'NaN'; } my ($mag, $angle) = self.polar; $mag **= 1e0 / $n; (^$n).map: { $mag.unpolar( ($angle + $_ * 2e0 * pi) / $n) }; } method sin(Complex:D:) { $!re.sin * $!im.cosh + ($!re.cos * $!im.sinh)i; } method asin(Complex:D:) { (Complex.new(0e0, -1e0) * log((self)i + sqrt(1e0 - self * self))); } method cos(Complex:D:) { $!re.cos * $!im.cosh - ($!re.sin * $!im.sinh)i; } method acos(Complex:D:) { (pi / 2e0) - self.asin; } method tan(Complex:D:) { self.sin / self.cos; } method atan(Complex:D:) { ((log(1e0 - (self)i) - log(1e0 + (self)i))i / 2e0); } method sec(Complex:D:) { 1e0 / self.cos; } method asec(Complex:D:) { (1e0 / self).acos; } method cosec(Complex:D:) { 1e0 / self.sin; } method acosec(Complex:D:) { (1e0 / self).asin; } method cotan(Complex:D:) { self.cos / self.sin; } method acotan(Complex:D:) { (1e0 / self).atan; } method sinh(Complex:D:) { -((Complex.new(0e0, 1e0) * self).sin)i; } method asinh(Complex:D:) { (self + sqrt(1e0 + self * self)).log; } method cosh(Complex:D:) { (Complex.new(0e0, 1e0) * self).cos; } method acosh(Complex:D:) { (self + sqrt(self * self - 1e0)).log; } method tanh(Complex:D:) { -((Complex.new(0e0, 1e0) * self).tan)i; } method atanh(Complex:D:) { (((1e0 + self) / (1e0 - self)).log / 2e0); } method sech(Complex:D:) { 1e0 / self.cosh; } method asech(Complex:D:) { (1e0 / self).acosh; } method cosech(Complex:D:) { 1e0 / self.sinh; } method acosech(Complex:D:) { (1e0 / self).asinh; } method cotanh(Complex:D:) { 1e0 / self.tanh; } method acotanh(Complex:D:) { (1e0 / self).atanh; } method floor(Complex:D:) { Complex.new( self.re.floor, self.im.floor ); } method ceiling(Complex:D:) { Complex.new( self.re.ceiling, self.im.ceiling ); } proto method round(|) {*} multi method round(Complex:D:) { Complex.new( self.re.round, self.im.round ); } multi method round(Complex:D: Real() $scale) { Complex.new( self.re.round($scale), self.im.round($scale) ); } method truncate(Complex:D:) { Complex.new( self.re.truncate, self.im.truncate ); } method narrow(Complex:D:) { self == 0e0 ?? 0 !! $!re == 0e0 ?? self !! $!im / $!re ≅ 0e0 ?? $!re.narrow !! self; } } multi sub prefix:<->(Complex:D \a --> Complex:D) { my $new := nqp::create(Complex); nqp::bindattr_n( $new, Complex, '$!re', nqp::neg_n( nqp::getattr_n(nqp::decont(a), Complex, '$!re') ) ); nqp::bindattr_n( $new, Complex, '$!im', nqp::neg_n( nqp::getattr_n(nqp::decont(a), Complex, '$!im') ) ); $new; } multi sub abs(Complex:D \a --> Num:D) { my num $re = nqp::getattr_n(nqp::decont(a), Complex, '$!re'); my num $im = nqp::getattr_n(nqp::decont(a), Complex, '$!im'); nqp::p6box_n(nqp::sqrt_n(nqp::add_n(nqp::mul_n($re, $re), nqp::mul_n($im, $im)))); } multi sub infix:<+>(Complex:D \a, Complex:D \b --> Complex:D) { my $new := nqp::create(Complex); nqp::bindattr_n( $new, Complex, '$!re', nqp::add_n( nqp::getattr_n(nqp::decont(a), Complex, '$!re'), nqp::getattr_n(nqp::decont(b), Complex, '$!re'), ) ); nqp::bindattr_n( $new, Complex, '$!im', nqp::add_n( nqp::getattr_n(nqp::decont(a), Complex, '$!im'), nqp::getattr_n(nqp::decont(b), Complex, '$!im'), ) ); $new; } multi sub infix:<+>(Complex:D \a, Num(Real) \b --> Complex:D) { my $new := nqp::create(Complex); nqp::bindattr_n( $new, Complex, '$!re', nqp::add_n( nqp::getattr_n(nqp::decont(a), Complex, '$!re'), nqp::unbox_n(b) ) ); nqp::bindattr_n($new, Complex, '$!im', nqp::getattr_n(nqp::decont(a), Complex, '$!im'), ); $new } multi sub infix:<+>(Num(Real) \a, Complex:D \b --> Complex:D) { my $new := nqp::create(Complex); nqp::bindattr_n($new, Complex, '$!re', nqp::add_n( nqp::unbox_n(a), nqp::getattr_n(nqp::decont(b), Complex, '$!re'), ) ); nqp::bindattr_n($new, Complex, '$!im', nqp::getattr_n(nqp::decont(b), Complex, '$!im'), ); $new; } multi sub infix:<->(Complex:D \a, Complex:D \b --> Complex:D) { my $new := nqp::create(Complex); nqp::bindattr_n( $new, Complex, '$!re', nqp::sub_n( nqp::getattr_n(nqp::decont(a), Complex, '$!re'), nqp::getattr_n(nqp::decont(b), Complex, '$!re'), ) ); nqp::bindattr_n($new, Complex, '$!im', nqp::sub_n( nqp::getattr_n(nqp::decont(a), Complex, '$!im'), nqp::getattr_n(nqp::decont(b), Complex, '$!im'), ) ); $new } multi sub infix:<->(Complex:D \a, Num(Real) \b --> Complex:D) { my $new := nqp::create(Complex); nqp::bindattr_n( $new, Complex, '$!re', nqp::sub_n( nqp::getattr_n(nqp::decont(a), Complex, '$!re'), b, ) ); nqp::bindattr_n($new, Complex, '$!im', nqp::getattr_n(nqp::decont(a), Complex, '$!im') ); $new } multi sub infix:<->(Num(Real) \a, Complex:D \b --> Complex:D) { my $new := nqp::create(Complex); nqp::bindattr_n( $new, Complex, '$!re', nqp::sub_n( a, nqp::getattr_n(nqp::decont(b), Complex, '$!re'), ) ); nqp::bindattr_n($new, Complex, '$!im', nqp::neg_n( nqp::getattr_n(nqp::decont(b), Complex, '$!im') ) ); $new } multi sub infix:<*>(Complex:D \a, Complex:D \b --> Complex:D) { my num $a_re = nqp::getattr_n(nqp::decont(a), Complex, '$!re'); my num $a_im = nqp::getattr_n(nqp::decont(a), Complex, '$!im'); my num $b_re = nqp::getattr_n(nqp::decont(b), Complex, '$!re'); my num $b_im = nqp::getattr_n(nqp::decont(b), Complex, '$!im'); my $new := nqp::create(Complex); nqp::bindattr_n($new, Complex, '$!re', nqp::sub_n(nqp::mul_n($a_re, $b_re), nqp::mul_n($a_im, $b_im)), ); nqp::bindattr_n($new, Complex, '$!im', nqp::add_n(nqp::mul_n($a_re, $b_im), nqp::mul_n($a_im, $b_re)), ); $new; } multi sub infix:<*>(Complex:D \a, Num(Real) \b --> Complex:D) { my $new := nqp::create(Complex); my num $b_num = b; nqp::bindattr_n($new, Complex, '$!re', nqp::mul_n( nqp::getattr_n(nqp::decont(a), Complex, '$!re'), $b_num, ) ); nqp::bindattr_n($new, Complex, '$!im', nqp::mul_n( nqp::getattr_n(nqp::decont(a), Complex, '$!im'), $b_num, ) ); $new } multi sub infix:<*>(Num(Real) \a, Complex:D \b --> Complex:D) { my $new := nqp::create(Complex); my num $a_num = a; nqp::bindattr_n($new, Complex, '$!re', nqp::mul_n( $a_num, nqp::getattr_n(nqp::decont(b), Complex, '$!re'), ) ); nqp::bindattr_n($new, Complex, '$!im', nqp::mul_n( $a_num, nqp::getattr_n(nqp::decont(b), Complex, '$!im'), ) ); $new } multi sub infix:</>(Complex:D \a, Complex:D \b --> Complex:D) { my num $a_re = nqp::getattr_n(nqp::decont(a), Complex, '$!re'); my num $a_im = nqp::getattr_n(nqp::decont(a), Complex, '$!im'); my num $b_re = nqp::getattr_n(nqp::decont(b), Complex, '$!re'); my num $b_im = nqp::getattr_n(nqp::decont(b), Complex, '$!im'); my num $d = nqp::add_n(nqp::mul_n($b_re, $b_re), nqp::mul_n($b_im, $b_im)); my $new := nqp::create(Complex); nqp::bindattr_n($new, Complex, '$!re', nqp::div_n( nqp::add_n(nqp::mul_n($a_re, $b_re), nqp::mul_n($a_im, $b_im)), $d, ) ); nqp::bindattr_n($new, Complex, '$!im', nqp::div_n( nqp::sub_n(nqp::mul_n($a_im, $b_re), nqp::mul_n($a_re, $b_im)), $d, ) ); $new; } multi sub infix:</>(Complex:D \a, Real \b --> Complex:D) { Complex.new(a.re / b, a.im / b); } multi sub infix:</>(Real \a, Complex:D \b --> Complex:D) { Complex.new(a, 0e0) / b; } multi sub infix:<**>(Complex:D \a, Complex:D \b --> Complex:D) { (a.re == 0e0 && a.im == 0e0) ?? ( b.re == 0e0 && b.im == 0e0 ?? Complex.new(1e0, 0e0) !! Complex.new(0e0, 0e0) ) !! (b * a.log).exp } multi sub infix:<**>(Num(Real) \a, Complex:D \b --> Complex:D) { a == 0e0 ?? ( b.re == 0e0 && b.im == 0e0 ?? Complex.new(1e0, 0e0) !! Complex.new(0e0, 0e0) ) !! (b * a.log).exp } multi sub infix:<**>(Complex:D \a, Num(Real) \b --> Complex:D) { b == 0e0 ?? Complex.new(1e0, 0e0) !! (b * a.log).exp } multi sub infix:<==>(Complex:D \a, Complex:D \b --> Bool:D) { a.re == b.re && a.im == b.im } multi sub infix:<==>(Complex:D \a, Num(Real) \b --> Bool:D) { a.re == b && a.im == 0e0 } multi sub infix:<==>(Num(Real) \a, Complex:D \b --> Bool:D) { a == b.re && 0e0 == b.im } multi sub infix:<===>(Complex:D \a, Complex:D \b --> Bool:D) { a.WHAT =:= b.WHAT && a.re === b.re && a.im === b.im } multi sub infix:<≅>(Complex:D \a, Complex:D \b --> Bool:D) { a.re ≅ b.re && a.im ≅ b.im || a <=> b =:= Same } multi sub infix:<≅>(Complex:D \a, Num(Real) \b --> Bool:D) { a ≅ b.Complex } multi sub infix:<≅>(Num(Real) \a, Complex:D \b --> Bool:D) { a.Complex ≅ b } # Meaningful only for sorting purposes, of course. # We delegate to Real::cmp rather than <=> because parts might be NaN. multi sub infix:<cmp>(Complex:D \a, Complex:D \b --> Order:D) { a.re cmp b.re || a.im cmp b.im } multi sub infix:<cmp>(Num(Real) \a, Complex:D \b --> Order:D) { a cmp b.re || 0 cmp b.im } multi sub infix:<cmp>(Complex:D \a, Num(Real) \b --> Order:D) { a.re cmp b || a.im cmp 0 } multi sub infix:«<=>»(Complex:D \a, Complex:D \b --> Order:D) { my $tolerance = a && b ?? (a.re.abs + b.re.abs) / 2 * $*TOLERANCE # Scale slop to average real parts. !! $*TOLERANCE; # Don't want tolerance 0 if either arg is 0. # Fail unless imaginary parts are relatively negligible, compared to real parts. infix:<≅>(a.im, 0e0, :$tolerance) && infix:<≅>(b.im, 0e0, :$tolerance) ?? a.re <=> b.re !! Failure.new(X::Numeric::Real.new(target => Real, reason => "Complex is not numerically orderable", source => "Complex")) } multi sub infix:«<=>»(Num(Real) \a, Complex:D \b --> Order:D) { a.Complex <=> b } multi sub infix:«<=>»(Complex:D \a, Num(Real) \b --> Order:D) { a <=> b.Complex } proto sub postfix:<i>($, *% --> Complex:D) is pure {*} multi sub postfix:<i>(Real \a --> Complex:D) { Complex.new(0e0, a); } multi sub postfix:<i>(Complex:D \a --> Complex:D) { Complex.new(-a.im, a.re) } multi sub postfix:<i>(Numeric \a --> Complex:D) { a * Complex.new(0e0, 1e0) } multi sub postfix:<i>(Cool \a --> Complex:D) { a.Numeric * Complex.new(0e0, 1e0) } constant i = Complex.new(0e0, 1e0); # vim: ft=perl6 expandtab sw=4