1 // Copyright 2010 The Go Authors. All rights reserved.
2 // Use of this source code is governed by a BSD-style
3 // license that can be found in the LICENSE file.
4 5 package cmplx
6 7 import "math"
8 9 // The original C code, the long comment, and the constants
10 // below are from http://netlib.sandia.gov/cephes/c9x-complex/clog.c.
11 // The go code is a simplified version of the original C.
12 //
13 // Cephes Math Library Release 2.8: June, 2000
14 // Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
15 //
16 // The readme file at http://netlib.sandia.gov/cephes/ says:
17 // Some software in this archive may be from the book _Methods and
18 // Programs for Mathematical Functions_ (Prentice-Hall or Simon & Schuster
19 // International, 1989) or from the Cephes Mathematical Library, a
20 // commercial product. In either event, it is copyrighted by the author.
21 // What you see here may be used freely but it comes with no support or
22 // guarantee.
23 //
24 // The two known misprints in the book are repaired here in the
25 // source listings for the gamma function and the incomplete beta
26 // integral.
27 //
28 // Stephen L. Moshier
29 // moshier@na-net.ornl.gov
30 31 // Complex power function
32 //
33 // DESCRIPTION:
34 //
35 // Raises complex A to the complex Zth power.
36 // Definition is per AMS55 # 4.2.8,
37 // analytically equivalent to cpow(a,z) = cexp(z clog(a)).
38 //
39 // ACCURACY:
40 //
41 // Relative error:
42 // arithmetic domain # trials peak rms
43 // IEEE -10,+10 30000 9.4e-15 1.5e-15
44 45 // Pow returns x**y, the base-x exponential of y.
46 // For generalized compatibility with [math.Pow]:
47 //
48 // Pow(0, ±0) returns 1+0i
49 // Pow(0, c) for real(c)<0 returns Inf+0i if imag(c) is zero, otherwise Inf+Inf i.
50 func Pow(x, y complex128) complex128 {
51 if x == 0 { // Guaranteed also true for x == -0.
52 if IsNaN(y) {
53 return NaN()
54 }
55 r, i := real(y), imag(y)
56 switch {
57 case r == 0:
58 return 1
59 case r < 0:
60 if i == 0 {
61 return complex(math.Inf(1), 0)
62 }
63 return Inf()
64 case r > 0:
65 return 0
66 }
67 panic("not reached")
68 }
69 modulus := Abs(x)
70 if modulus == 0 {
71 return complex(0, 0)
72 }
73 r := math.Pow(modulus, real(y))
74 arg := Phase(x)
75 theta := real(y) * arg
76 if imag(y) != 0 {
77 r *= math.Exp(-imag(y) * arg)
78 theta += imag(y) * math.Log(modulus)
79 }
80 s, c := math.Sincos(theta)
81 return complex(r*c, r*s)
82 }
83