complex pow in 1024bit

complex_pow.rb

def complex_fixed_point_1024bits_pow(x, n, bits: 1024)
  ans = base = 1 << bits
  x2 = Complex.rect(*x.rect.map{(_1.to_r * base).round})
  mult = ->(a, b) { Complex.rect(*(a * b).rect.map { _1 / base }) }
  n.digits(2).each do |a|
    ans = mult[ans, x2] if a == 1
    x2 = mult[x2, x2]
  end
  ans.fdiv base
end

p complex_fixed_point_1024bits_pow(Complex.polar(1,1), 1000000000000000000).abs

n = 100000000000000
x = Complex.polar(1, 1)
p abs: complex_fixed_point_1024bits_pow(x, n).abs
p actual: complex_fixed_point_1024bits_pow(x, n)
p fixnum: x ** n
p float: x ** n.to_f

# =>
# 33285375491.354782
# {abs: 1.0024257758410897}
# {actual: (-0.9810531507137839-0.20589354420468883i)}
# {fixnum: (-0.9822320247697417-0.2047636321393639i)}
# {float: (-0.9778282879685325-0.2094083074964523i)}