Power of Complex Numbers

[RainbowZephyr]

Feature Request

Hello, we use crystal in a few projects for mathematical analysis and we are quite satisfied with the performance. The only issue we face is raising a complex number to a power, for example $(1+2i)^5$, this feature is not available as far as I know.

Could this feature be implemented?

Thanks in advance.

[jzakiya]

The easiest way to do this is to convert the number: z = a + bi into polar form z = (a^2 + b^2)^(1/2) arctan (b/a) = ||z|| e^(i*theta)

then z^n = ||z||^n e^(ntheta)i, which you can convert back into complex form if you need to.

Of course, for rigorous accuracy, the angle should be - theta + 2k*pi

[erdian718]

Perhaps this feature will be added in the future, and in the current version, if you don't care too much about performance, you can do so:

def pow(x, y)
  Math.exp(y * Math.log(x))
end

[erdian718]

I benchmarked both algorithms with almost no difference in performance, but the first algorithm was simpler and more versatile.

require "complex"
require "benchmark"

a = Complex.new(1.1, 2.2)
b = -0.5

def pow1(a, b)
  Math.exp(b * Math.log(a))
end

def pow2(a : Complex, b) : Complex
  b = b.to_f64
  abs, phase = a.polar
  abs **= b
  phase *= b
  Complex.new(abs*Math.cos(phase), abs*Math.sin(phase))
end

Benchmark.ips do |x|
  x.report("pow1") { pow1(a, b) }
  x.report("pow2") { pow2(a, b) }
end

pow1  44.89M ( 22.28ns) (± 1.30%)  0.0B/op        fastest
pow2  44.41M ( 22.52ns) (± 1.36%)  0.0B/op   1.01× slower

I created a PR #13802

[jzakiya]

Nice work!

Comparing each with different inputs give me these results.

a, b = Complex.new(1.1, 2.2), -0.5
puts "z = #{a}, power = #{b}"
puts pow1(a,b)
puts pow2(a,b)
puts

a, b = Complex.new(5, -8), 7
puts "z = #{a}, power = #{b}"
puts pow1(a,b)
puts pow2(a,b)
puts

a, b = Complex.new(10.5, 7), 20.5
puts "z = #{a}, power = #{b}"
puts pow1(a,b)
puts pow2(a,b)
puts

a, b = Complex.new(100, 45), -4
puts "z = #{a}, power = #{b}"
puts pow1(a,b)
puts pow2(a,b)
puts
----------------------------------------------
z = 1.1 + 2.2i, power = -0.5
0.542391000989624 - 0.33521607380366547i
0.542391000989624 - 0.33521607380366547i

z = 5.0 - 8.0i, power = 7
4623084.999999997 - 4781047.999999989i
4623085.000000002 - 4781047.9999999935i

z = 10.5 + 7.0i, power = 20.5
3.2481564813785184e+22 - 1.8267792339489547e+22i
3.2481564813785234e+22 - 1.8267792339489576e+22i

z = 100.0 + 45.0i, power = -4
-8.321340786258395e-10 - 6.8653527489774366e-9i
-8.321340786258384e-10 - 6.865352748977428e-9i

[jzakiya]

Comparing more inputs show pow2 gives more accurate (exact) results with integer inputs.

a, b = Complex.new(100, 0), 3
puts "z = #{a}, power = #{b}"
puts pow1(a,b)
puts pow2(a,b)
puts

a, b = Complex.new(0, 100), 3
puts "z = #{a}, power = #{b}"
puts pow1(a,b)
puts pow2(a,b)
puts

a, b = Complex.new(1000000, 0), -3
puts "z = #{a}, power = #{b}"
puts pow1(a,b)
puts pow2(a,b)
puts

a, b = Complex.new(0, 1000000), -3
puts "z = #{a}, power = #{b}"
puts pow1(a,b)
puts pow2(a,b)
puts

a, b = Complex.new(1000000, 0), 0.5
puts "z = #{a}, power = #{b}"
puts pow1(a,b)
puts pow2(a,b)
puts

a, b = Complex.new(0, 1000000), 0.5
puts "z = #{a}, power = #{b}"
puts pow1(a,b)
puts pow2(a,b)
puts
---------------------------------
z = 100.0 + 0.0i, power = 3
1000000.0000000013 + 0.0i
1000000.0 + 0.0i

z = 0.0 + 100.0i, power = 3
-1.836970198721032e-10 - 1000000.0000000013i
-1.8369701987210297e-10 - 1000000.0i

z = 1000000.0 + 0.0i, power = -3
1.0000000000000014e-18 - 0.0i
1.0e-18 - 0.0i

z = 0.0 + 1000000.0i, power = -3
-1.8369701987210323e-34 + 1.0000000000000014e-18i
-1.8369701987210297e-34 + 1.0e-18i

z = 1000000.0 + 0.0i, power = 0.5
999.9999999999998 + 0.0i
1000.0 + 0.0i

z = 0.0 + 1000000.0i, power = 0.5
707.1067811865474 + 707.1067811865473i
707.1067811865476 + 707.1067811865474i

Though brevity of code is nice, it should not trump accuracy. Since there is no significance performance difference between the versions, I'd go with pow2 because it gives more accurate results with more types of inputs.

[erdian718]

Thanks for the reminder @jzakiya, pow2 does have better numerical stability, I improved the benchmark and I'll update the PR later.

require "complex"
require "benchmark"

def pow1(a, b)
  Math.exp(b * Math.log(a))
end

def pow2(a : Complex, b : Number) : Complex
  abs, phase = a.polar
  r = abs ** b
  t = phase * b
  Complex.new(r*Math.cos(t), r*Math.sin(t))
end

def pow2(a : Number, b : Complex)
  a = a.abs
  r = a ** b.real
  t = b.imag * Math.log(a)
  Complex.new(r*Math.cos(t), r*Math.sin(t))
end

def pow2(a : Complex, b : Complex)
  abs, phase = a.polar
  r = (abs ** b.real) / Math.exp(phase * b.imag)
  t = phase * b.real + Math.log(abs) * b.imag
  Complex.new(r*Math.cos(t), r*Math.sin(t))
end

puts "complex ** number"
a1 = Complex.new(1.1, 2.2)
b1 = -0.5
Benchmark.ips do |x|
  x.report("pow1") { pow1(a1, b1) }
  x.report("pow2") { pow2(a1, b1) }
end

puts "number ** complex"
a2 = -0.5
b2 = Complex.new(1.1, 2.2)
Benchmark.ips do |x|
  x.report("pow1") { pow1(a2, b2) }
  x.report("pow2") { pow2(a2, b2) }
end

puts "complex ** complex"
a3 = Complex.new(1.1, 2.2)
b3 = Complex.new(3.3, 4.4)
Benchmark.ips do |x|
  x.report("pow1") { pow1(a3, b3) }
  x.report("pow2") { pow2(a3, b3) }
end

complex ** number
pow1  45.37M ( 22.04ns) (± 1.03%)  0.0B/op   1.00× slower
pow2  45.37M ( 22.04ns) (± 1.03%)  0.0B/op        fastest
number ** complex
pow1 789.34M (  1.27ns) (± 2.20%)  0.0B/op   1.00× slower
pow2 790.19M (  1.27ns) (± 2.10%)  0.0B/op        fastest
complex ** complex
pow1  45.53M ( 21.97ns) (± 0.44%)  0.0B/op        fastest
pow2  45.39M ( 22.03ns) (± 1.20%)  0.0B/op   1.00× slower