1. my class X::Numeric::Real { ... };
  2. my class Complex is Cool does Numeric {
  3. has num $.re;
  4. has num $.im;
  5. method !SET-SELF(Num() \re, Num() \im) {
  6. $!re = re;
  7. $!im = im;
  8. self
  9. }
  10. proto method new(|) { * }
  11. multi method new() { self.new: 0, 0 }
  12. multi method new(Real \re, Real \im) { nqp::create(self)!SET-SELF(re, im) }
  13. multi method WHICH(Complex:D:) {
  14. nqp::box_s(
  15. nqp::concat(
  16. nqp::if(
  17. nqp::eqaddr(self.WHAT,Complex),
  18. 'Complex|',
  19. nqp::concat(nqp::unbox_s(self.^name), '|')
  20. ),
  21. nqp::concat($!re, nqp::concat('|', $!im))
  22. ),
  23. ObjAt
  24. )
  25. }
  26. method reals(Complex:D:) {
  27. (self.re, self.im);
  28. }
  29. method isNaN(Complex:D:) {
  30. self.re.isNaN || self.im.isNaN;
  31. }
  32. method coerce-to-real(Complex:D: $exception-target) {
  33. $!im ≅ 0e0
  34. ?? $!re
  35. !! Failure.new(X::Numeric::Real.new(target => $exception-target, reason => "imaginary part not zero", source => self))
  36. }
  37. multi method Real(Complex:D:) { self.coerce-to-real(Real); }
  38. # should probably be eventually supplied by role Numeric
  39. method Num(Complex:D:) { self.coerce-to-real(Num).Num; }
  40. method Int(Complex:D:) { self.coerce-to-real(Int).Int; }
  41. method Rat(Complex:D: $epsilon?) {
  42. self.coerce-to-real(Rat).Rat( |($epsilon // Empty) );
  43. }
  44. method FatRat(Complex:D: $epsilon?) {
  45. self.coerce-to-real(FatRat).FatRat( |($epsilon // Empty) );
  46. }
  47. multi method Bool(Complex:D:) {
  48. $!re != 0e0 || $!im != 0e0;
  49. }
  50. method Complex() { self }
  51. multi method Str(Complex:D:) {
  52. nqp::concat(
  53. $!re,
  54. nqp::concat(
  55. nqp::if(nqp::iseq_i(nqp::ord($!im),45),'','+'),
  56. nqp::concat(
  57. $!im,
  58. nqp::if(nqp::isnanorinf($!im),'\\i','i')
  59. )
  60. )
  61. )
  62. }
  63. multi method perl(Complex:D:) {
  64. '<' ~ self.Str ~ '>';
  65. }
  66. method conj(Complex:D:) {
  67. Complex.new($.re, -$.im);
  68. }
  69. method abs(Complex $x:) {
  70. nqp::p6box_n(nqp::sqrt_n(
  71. nqp::add_n(
  72. nqp::mul_n($!re, $!re),
  73. nqp::mul_n($!im, $!im),
  74. )
  75. ))
  76. }
  77. method polar() {
  78. $.abs, $!im.atan2($!re);
  79. }
  80. multi method log(Complex:D:) {
  81. my Num ($mag, $angle) = self.polar;
  82. Complex.new($mag.log, $angle);
  83. }
  84. method sqrt(Complex:D:) {
  85. my Num $abs = self.abs;
  86. my Num $re = (($abs + self.re)/2).sqrt;
  87. my Num $im = (($abs - self.re)/2).sqrt;
  88. Complex.new($re, self.im < 0 ?? -$im !! $im);
  89. }
  90. multi method exp(Complex:D:) {
  91. my Num $mag = $!re.exp;
  92. Complex.new($mag * $!im.cos, $mag * $!im.sin);
  93. }
  94. method roots(Complex:D: Int() $n) {
  95. return NaN if $n < 1;
  96. return self if $n == 1;
  97. for $!re, $!im {
  98. return NaN if $_ eq 'Inf' || $_ eq '-Inf' || $_ eq 'NaN';
  99. }
  100. my ($mag, $angle) = self.polar;
  101. $mag **= 1e0 / $n;
  102. (^$n).map: { $mag.unpolar( ($angle + $_ * 2e0 * pi) / $n) };
  103. }
  104. method sin(Complex:D:) {
  105. $!re.sin * $!im.cosh + ($!re.cos * $!im.sinh)i;
  106. }
  107. method asin(Complex:D:) {
  108. (Complex.new(0e0, -1e0) * log((self)i + sqrt(1e0 - self * self)));
  109. }
  110. method cos(Complex:D:) {
  111. $!re.cos * $!im.cosh - ($!re.sin * $!im.sinh)i;
  112. }
  113. method acos(Complex:D:) {
  114. (pi / 2e0) - self.asin;
  115. }
  116. method tan(Complex:D:) {
  117. self.sin / self.cos;
  118. }
  119. method atan(Complex:D:) {
  120. ((log(1e0 - (self)i) - log(1e0 + (self)i))i / 2e0);
  121. }
  122. method sec(Complex:D:) {
  123. 1e0 / self.cos;
  124. }
  125. method asec(Complex:D:) {
  126. (1e0 / self).acos;
  127. }
  128. method cosec(Complex:D:) {
  129. 1e0 / self.sin;
  130. }
  131. method acosec(Complex:D:) {
  132. (1e0 / self).asin;
  133. }
  134. method cotan(Complex:D:) {
  135. self.cos / self.sin;
  136. }
  137. method acotan(Complex:D:) {
  138. (1e0 / self).atan;
  139. }
  140. method sinh(Complex:D:) {
  141. -((Complex.new(0e0, 1e0) * self).sin)i;
  142. }
  143. method asinh(Complex:D:) {
  144. (self + sqrt(1e0 + self * self)).log;
  145. }
  146. method cosh(Complex:D:) {
  147. (Complex.new(0e0, 1e0) * self).cos;
  148. }
  149. method acosh(Complex:D:) {
  150. (self + sqrt(self * self - 1e0)).log;
  151. }
  152. method tanh(Complex:D:) {
  153. -((Complex.new(0e0, 1e0) * self).tan)i;
  154. }
  155. method atanh(Complex:D:) {
  156. (((1e0 + self) / (1e0 - self)).log / 2e0);
  157. }
  158. method sech(Complex:D:) {
  159. 1e0 / self.cosh;
  160. }
  161. method asech(Complex:D:) {
  162. (1e0 / self).acosh;
  163. }
  164. method cosech(Complex:D:) {
  165. 1e0 / self.sinh;
  166. }
  167. method acosech(Complex:D:) {
  168. (1e0 / self).asinh;
  169. }
  170. method cotanh(Complex:D:) {
  171. 1e0 / self.tanh;
  172. }
  173. method acotanh(Complex:D:) {
  174. (1e0 / self).atanh;
  175. }
  176. method floor(Complex:D:) {
  177. Complex.new( self.re.floor, self.im.floor );
  178. }
  179. method ceiling(Complex:D:) {
  180. Complex.new( self.re.ceiling, self.im.ceiling );
  181. }
  182. proto method round(|) {*}
  183. multi method round(Complex:D:) {
  184. Complex.new( self.re.round, self.im.round );
  185. }
  186. multi method round(Complex:D: Real() $scale) {
  187. Complex.new( self.re.round($scale), self.im.round($scale) );
  188. }
  189. method truncate(Complex:D:) {
  190. Complex.new( self.re.truncate, self.im.truncate );
  191. }
  192. method narrow(Complex:D:) {
  193. self == 0e0 ?? 0 !!
  194. $!re == 0e0 ?? self !!
  195. $!im / $!re ≅ 0e0
  196. ?? $!re.narrow
  197. !! self;
  198. }
  199. }
  200. multi sub prefix:<->(Complex:D \a --> Complex:D) {
  201. my $new := nqp::create(Complex);
  202. nqp::bindattr_n( $new, Complex, '$!re',
  203. nqp::neg_n(
  204. nqp::getattr_n(nqp::decont(a), Complex, '$!re')
  205. )
  206. );
  207. nqp::bindattr_n( $new, Complex, '$!im',
  208. nqp::neg_n(
  209. nqp::getattr_n(nqp::decont(a), Complex, '$!im')
  210. )
  211. );
  212. $new;
  213. }
  214. multi sub abs(Complex:D \a --> Num:D) {
  215. my num $re = nqp::getattr_n(nqp::decont(a), Complex, '$!re');
  216. my num $im = nqp::getattr_n(nqp::decont(a), Complex, '$!im');
  217. nqp::p6box_n(nqp::sqrt_n(nqp::add_n(nqp::mul_n($re, $re), nqp::mul_n($im, $im))));
  218. }
  219. multi sub infix:<+>(Complex:D \a, Complex:D \b --> Complex:D) {
  220. my $new := nqp::create(Complex);
  221. nqp::bindattr_n( $new, Complex, '$!re',
  222. nqp::add_n(
  223. nqp::getattr_n(nqp::decont(a), Complex, '$!re'),
  224. nqp::getattr_n(nqp::decont(b), Complex, '$!re'),
  225. )
  226. );
  227. nqp::bindattr_n( $new, Complex, '$!im',
  228. nqp::add_n(
  229. nqp::getattr_n(nqp::decont(a), Complex, '$!im'),
  230. nqp::getattr_n(nqp::decont(b), Complex, '$!im'),
  231. )
  232. );
  233. $new;
  234. }
  235. multi sub infix:<+>(Complex:D \a, Num(Real) \b --> Complex:D) {
  236. my $new := nqp::create(Complex);
  237. nqp::bindattr_n( $new, Complex, '$!re',
  238. nqp::add_n(
  239. nqp::getattr_n(nqp::decont(a), Complex, '$!re'),
  240. nqp::unbox_n(b)
  241. )
  242. );
  243. nqp::bindattr_n($new, Complex, '$!im',
  244. nqp::getattr_n(nqp::decont(a), Complex, '$!im'),
  245. );
  246. $new
  247. }
  248. multi sub infix:<+>(Num(Real) \a, Complex:D \b --> Complex:D) {
  249. my $new := nqp::create(Complex);
  250. nqp::bindattr_n($new, Complex, '$!re',
  251. nqp::add_n(
  252. nqp::unbox_n(a),
  253. nqp::getattr_n(nqp::decont(b), Complex, '$!re'),
  254. )
  255. );
  256. nqp::bindattr_n($new, Complex, '$!im',
  257. nqp::getattr_n(nqp::decont(b), Complex, '$!im'),
  258. );
  259. $new;
  260. }
  261. multi sub infix:<->(Complex:D \a, Complex:D \b --> Complex:D) {
  262. my $new := nqp::create(Complex);
  263. nqp::bindattr_n( $new, Complex, '$!re',
  264. nqp::sub_n(
  265. nqp::getattr_n(nqp::decont(a), Complex, '$!re'),
  266. nqp::getattr_n(nqp::decont(b), Complex, '$!re'),
  267. )
  268. );
  269. nqp::bindattr_n($new, Complex, '$!im',
  270. nqp::sub_n(
  271. nqp::getattr_n(nqp::decont(a), Complex, '$!im'),
  272. nqp::getattr_n(nqp::decont(b), Complex, '$!im'),
  273. )
  274. );
  275. $new
  276. }
  277. multi sub infix:<->(Complex:D \a, Num(Real) \b --> Complex:D) {
  278. my $new := nqp::create(Complex);
  279. nqp::bindattr_n( $new, Complex, '$!re',
  280. nqp::sub_n(
  281. nqp::getattr_n(nqp::decont(a), Complex, '$!re'),
  282. b,
  283. )
  284. );
  285. nqp::bindattr_n($new, Complex, '$!im',
  286. nqp::getattr_n(nqp::decont(a), Complex, '$!im')
  287. );
  288. $new
  289. }
  290. multi sub infix:<->(Num(Real) \a, Complex:D \b --> Complex:D) {
  291. my $new := nqp::create(Complex);
  292. nqp::bindattr_n( $new, Complex, '$!re',
  293. nqp::sub_n(
  294. a,
  295. nqp::getattr_n(nqp::decont(b), Complex, '$!re'),
  296. )
  297. );
  298. nqp::bindattr_n($new, Complex, '$!im',
  299. nqp::neg_n(
  300. nqp::getattr_n(nqp::decont(b), Complex, '$!im')
  301. )
  302. );
  303. $new
  304. }
  305. multi sub infix:<*>(Complex:D \a, Complex:D \b --> Complex:D) {
  306. my num $a_re = nqp::getattr_n(nqp::decont(a), Complex, '$!re');
  307. my num $a_im = nqp::getattr_n(nqp::decont(a), Complex, '$!im');
  308. my num $b_re = nqp::getattr_n(nqp::decont(b), Complex, '$!re');
  309. my num $b_im = nqp::getattr_n(nqp::decont(b), Complex, '$!im');
  310. my $new := nqp::create(Complex);
  311. nqp::bindattr_n($new, Complex, '$!re',
  312. nqp::sub_n(nqp::mul_n($a_re, $b_re), nqp::mul_n($a_im, $b_im)),
  313. );
  314. nqp::bindattr_n($new, Complex, '$!im',
  315. nqp::add_n(nqp::mul_n($a_re, $b_im), nqp::mul_n($a_im, $b_re)),
  316. );
  317. $new;
  318. }
  319. multi sub infix:<*>(Complex:D \a, Num(Real) \b --> Complex:D) {
  320. my $new := nqp::create(Complex);
  321. my num $b_num = b;
  322. nqp::bindattr_n($new, Complex, '$!re',
  323. nqp::mul_n(
  324. nqp::getattr_n(nqp::decont(a), Complex, '$!re'),
  325. $b_num,
  326. )
  327. );
  328. nqp::bindattr_n($new, Complex, '$!im',
  329. nqp::mul_n(
  330. nqp::getattr_n(nqp::decont(a), Complex, '$!im'),
  331. $b_num,
  332. )
  333. );
  334. $new
  335. }
  336. multi sub infix:<*>(Num(Real) \a, Complex:D \b --> Complex:D) {
  337. my $new := nqp::create(Complex);
  338. my num $a_num = a;
  339. nqp::bindattr_n($new, Complex, '$!re',
  340. nqp::mul_n(
  341. $a_num,
  342. nqp::getattr_n(nqp::decont(b), Complex, '$!re'),
  343. )
  344. );
  345. nqp::bindattr_n($new, Complex, '$!im',
  346. nqp::mul_n(
  347. $a_num,
  348. nqp::getattr_n(nqp::decont(b), Complex, '$!im'),
  349. )
  350. );
  351. $new
  352. }
  353. multi sub infix:</>(Complex:D \a, Complex:D \b --> Complex:D) {
  354. my num $a_re = nqp::getattr_n(nqp::decont(a), Complex, '$!re');
  355. my num $a_im = nqp::getattr_n(nqp::decont(a), Complex, '$!im');
  356. my num $b_re = nqp::getattr_n(nqp::decont(b), Complex, '$!re');
  357. my num $b_im = nqp::getattr_n(nqp::decont(b), Complex, '$!im');
  358. my num $d = nqp::add_n(nqp::mul_n($b_re, $b_re), nqp::mul_n($b_im, $b_im));
  359. my $new := nqp::create(Complex);
  360. nqp::bindattr_n($new, Complex, '$!re',
  361. nqp::div_n(
  362. nqp::add_n(nqp::mul_n($a_re, $b_re), nqp::mul_n($a_im, $b_im)),
  363. $d,
  364. )
  365. );
  366. nqp::bindattr_n($new, Complex, '$!im',
  367. nqp::div_n(
  368. nqp::sub_n(nqp::mul_n($a_im, $b_re), nqp::mul_n($a_re, $b_im)),
  369. $d,
  370. )
  371. );
  372. $new;
  373. }
  374. multi sub infix:</>(Complex:D \a, Real \b --> Complex:D) {
  375. Complex.new(a.re / b, a.im / b);
  376. }
  377. multi sub infix:</>(Real \a, Complex:D \b --> Complex:D) {
  378. Complex.new(a, 0e0) / b;
  379. }
  380. multi sub infix:<**>(Complex:D \a, Complex:D \b --> Complex:D) {
  381. (a.re == 0e0 && a.im == 0e0)
  382. ?? ( b.re == 0e0 && b.im == 0e0
  383. ?? Complex.new(1e0, 0e0)
  384. !! Complex.new(0e0, 0e0)
  385. )
  386. !! (b * a.log).exp
  387. }
  388. multi sub infix:<**>(Num(Real) \a, Complex:D \b --> Complex:D) {
  389. a == 0e0
  390. ?? ( b.re == 0e0 && b.im == 0e0
  391. ?? Complex.new(1e0, 0e0)
  392. !! Complex.new(0e0, 0e0)
  393. )
  394. !! (b * a.log).exp
  395. }
  396. multi sub infix:<**>(Complex:D \a, Num(Real) \b --> Complex:D) {
  397. b == 0e0 ?? Complex.new(1e0, 0e0) !! (b * a.log).exp
  398. }
  399. multi sub infix:<==>(Complex:D \a, Complex:D \b --> Bool:D) { a.re == b.re && a.im == b.im }
  400. multi sub infix:<==>(Complex:D \a, Num(Real) \b --> Bool:D) { a.re == b && a.im == 0e0 }
  401. multi sub infix:<==>(Num(Real) \a, Complex:D \b --> Bool:D) { a == b.re && 0e0 == b.im }
  402. multi sub infix:<===>(Complex:D \a, Complex:D \b --> Bool:D) {
  403. a.WHAT =:= b.WHAT && a.re === b.re && a.im === b.im
  404. }
  405. multi sub infix:<≅>(Complex:D \a, Complex:D \b --> Bool:D) { a.re ≅ b.re && a.im ≅ b.im || a <=> b =:= Same }
  406. multi sub infix:<≅>(Complex:D \a, Num(Real) \b --> Bool:D) { a ≅ b.Complex }
  407. multi sub infix:<≅>(Num(Real) \a, Complex:D \b --> Bool:D) { a.Complex ≅ b }
  408. # Meaningful only for sorting purposes, of course.
  409. # We delegate to Real::cmp rather than <=> because parts might be NaN.
  410. multi sub infix:<cmp>(Complex:D \a, Complex:D \b --> Order:D) { a.re cmp b.re || a.im cmp b.im }
  411. multi sub infix:<cmp>(Num(Real) \a, Complex:D \b --> Order:D) { a cmp b.re || 0 cmp b.im }
  412. multi sub infix:<cmp>(Complex:D \a, Num(Real) \b --> Order:D) { a.re cmp b || a.im cmp 0 }
  413. multi sub infix:«<=>»(Complex:D \a, Complex:D \b --> Order:D) {
  414. my $tolerance = a && b
  415. ?? (a.re.abs + b.re.abs) / 2 * $*TOLERANCE # Scale slop to average real parts.
  416. !! $*TOLERANCE; # Don't want tolerance 0 if either arg is 0.
  417. # Fail unless imaginary parts are relatively negligible, compared to real parts.
  418. infix:<≅>(a.im, 0e0, :$tolerance) && infix:<≅>(b.im, 0e0, :$tolerance)
  419. ?? a.re <=> b.re
  420. !! Failure.new(X::Numeric::Real.new(target => Real, reason => "Complex is not numerically orderable", source => "Complex"))
  421. }
  422. multi sub infix:«<=>»(Num(Real) \a, Complex:D \b --> Order:D) { a.Complex <=> b }
  423. multi sub infix:«<=>»(Complex:D \a, Num(Real) \b --> Order:D) { a <=> b.Complex }
  424. proto sub postfix:<i>(\a --> Complex:D) is pure { * }
  425. multi sub postfix:<i>(Real \a --> Complex:D) { Complex.new(0e0, a); }
  426. multi sub postfix:<i>(Complex:D \a --> Complex:D) { Complex.new(-a.im, a.re) }
  427. multi sub postfix:<i>(Numeric \a --> Complex:D) { a * Complex.new(0e0, 1e0) }
  428. multi sub postfix:<i>(Cool \a --> Complex:D) { a.Numeric * Complex.new(0e0, 1e0) }
  429. constant i = Complex.new(0e0, 1e0);