Bug #11183: Cumulative error on Complex::I ** 100000000000000000000000000000000

• Status: Closed
• Assignee: -
• Target version: -
• Backport: 2.0.0: UNKNOWN, 2.1: UNKNOWN, 2.2: UNKNOWN
• Tags: core

[ko1 (Koichi Sasada)]

p Complex::I ** 100000000000000000000000000000000

が、32bit 環境だと、

ruby 2.3.0dev (2015-05-21 trunk 50502) [i386-mswin32_110]
t.rb:1: warning: in a**b, b may be too big
(-0.08483712490438944+1.5651333429325577e+32i)

こんな感じで、変な値が出ます。

# ubuntu 64bit だと
(0.9943722416631883-0.10594264962575697i)

できれば、こっちも 1+0i にぴたっとしてるといいのでしょうが。

この問題に対する、うささんのコメント:

真面目にコードを読むと、Complex#** は引数(指数)がFixnumでない場合は複素数をいったん極座標形式に変換してからrとθの双方をそれぞれ計算し、最後にf_complex_polarを使ってComplexオブジェクトにしてるのですが、θはFloatなので当然に誤差が蓄積し、という感じです。まあ、誤差が蓄積しなかったとしても、f_complex_polarの中でせっかくθが0°・90°・180°・270°のケースを特別扱いしようとしてるのにdouble値を==で比較してるのでまったく救われそうにない、という感じではありますが。
θをn倍する計算に対して、nが整数である場合は360°は0°と同じという性質を利用して誤差を蓄積しない計算をすればだいぶマシになるとは思いますが、しかしめんどくさいっすね。

[usa (Usaku NAKAMURA)]

• File my_complex.diff my_complex.diff added

めんどくさいと言いつつ、手元には、添付のような「誤差を蓄積しない計算」を試してみたパッチがありますが(コンセプトを確認するためのもので実用ではない)、次は「45°とかは？」という話になることは明らかなのでびみょーです。
より賢い方法が求められていると思われます。

[JesseJohnson (Jesse Johnson)]

This behavior still exists in Ruby 3.2.2.

In order to fix this Complex#** would need to handle the special cases of self being 0+1i and 0-1i.

Correct:

irb(main):052:0> Complex(0, -1) ** 1000000000000000000
=> (1+0i)
irb(main):053:0> Complex(0, -1) ** 1000000000000000001
=> (0-1i)
irb(main):054:0> Complex(0, -1) ** 1000000000000000002
=> (-1+0i)
irb(main):055:0> Complex(0, -1) ** 1000000000000000003
=> (0+1i)
irb(main):056:0> Complex(0, -1) ** 1000000000000000004
=> (1+0i)
irb(main):057:0> Complex(0, 1) ** 1000000000000000000
=> (1+0i)
irb(main):058:0> Complex(0, 1) ** 1000000000000000001
=> (0+1i)
irb(main):059:0> Complex(0, 1) ** 1000000000000000002
=> (-1+0i)
irb(main):060:0> Complex(0, 1) ** 1000000000000000003
=> (0-1i)
irb(main):061:0> Complex(0, 1) ** 1000000000000000004
=> (1+0i)

Suboptimal:

irb(main):066:0> Complex(0, -1) ** 10000000000000000000
(irb):66: warning: in a**b, b may be too big
=> (-0.9125701512929312+0.40892018655135703i)
irb(main):067:0> Complex(0, -1) ** 10000000000000000001
(irb):67: warning: in a**b, b may be too big
=> (-0.9125701512929312+0.40892018655135703i)
irb(main):068:0> Complex(0, -1) ** 10000000000000000002
(irb):68: warning: in a**b, b may be too big
=> (-0.9125701512929312+0.40892018655135703i)
irb(main):069:0> Complex(0, -1) ** 10000000000000000003
(irb):69: warning: in a**b, b may be too big
=> (-0.9125701512929312+0.40892018655135703i)
irb(main):070:0> Complex(0, -1) ** 10000000000000000004
(irb):70: warning: in a**b, b may be too big
=> (-0.9125701512929312+0.40892018655135703i)
irb(main):071:0> Complex(0, 1) ** 10000000000000000000
(irb):71: warning: in a**b, b may be too big
=> (-0.9125701512929312-0.40892018655135703i)
irb(main):072:0> Complex(0, 1) ** 10000000000000000001
(irb):72: warning: in a**b, b may be too big
=> (-0.9125701512929312-0.40892018655135703i)
irb(main):073:0> Complex(0, 1) ** 10000000000000000002
(irb):73: warning: in a**b, b may be too big
=> (-0.9125701512929312-0.40892018655135703i)
irb(main):074:0> Complex(0, 1) ** 10000000000000000003
(irb):74: warning: in a**b, b may be too big
=> (-0.9125701512929312-0.40892018655135703i)
irb(main):075:0> Complex(0, 1) ** 10000000000000000004
(irb):75: warning: in a**b, b may be too big
=> (-0.9125701512929312-0.40892018655135703i)
irb(main):076:0> in):070:0> Complex(0, -1) ** 10000000000000000004
=> (-0.9125701512929312-0.40892018655135703i)

[kyanagi (Kouhei Yanagita)]

https://github.com/ruby/ruby/pull/8931

self の偏角が45度の整数倍であるような場合に特別扱いするようにしてみました。

1i ** (10 ** 100)         #=> 1+0i
1i ** (10 ** 100 + 1)     #=> 0+1i
(1+1i) ** (10 ** 100)     # warning: in a**b, b may be too big
                          #=> (Infinity+0.0i)
(1+1i) ** (10 ** 100 + 1) # warning: in a**b, b may be too big
                          #=> (Infinity+Infinity*i)

• self が 1, -1, 1i, -1i の場合
  • bignum 乗であっても正しい答えを返します。
• self が実数または純虚数で絶対値が1でない場合
• 偏角が 45 度の奇数倍である場合
  • 小さい数（大きい数）が0（Infinity）となりますが、その前提のもとで正しい答えを返します。

絶対値が 1 で偏角が 45 度の奇数倍である複素数は、実部・虚部（1/\sqrt{2}）をFloatで正確に表現することはできないため、
偏角が45度の奇数倍である複素数オブジェクトの n 乗は、絶対値が1からずれていくことになります。

なお、60 度など他の角度についてですが、偏角が 45 度の整数倍を除く有理数度である複素数は、Complex では正確には表現できません。（実部または虚部のいずれかが無理数になる）
従って、Complex オブジェクトの偏角は 45 度の整数倍であるか無理数であるかになるため、
偏角を n 倍することで n 乗を求める計算が誤差なしでできるのは、偏角が 45 度の整数倍の場合に限られることになります。

[kyanagi (Kouhei Yanagita)]

• Status changed from Open to Closed

Applied in changeset git|04eb40b633397d03e4cbce41418626f4fabdcb02.

---

[Bug #11183] Fix rb_complex_pow for special angles

Add a special treatment for when the argument of self is an
integral multiple of 45 degrees.

1i ** (10 ** 100) #=> 1+0i
1i ** (10 ** 100 + 1) #=> 0+1i
(1+1i) ** (10 ** 100) # warning: in ab, b may be too big
#=> (Infinity+0.0i)
(1+1i) ** (10 ** 100 + 1) # warning: in ab, b may be too big
#=> (Infinity+Infinity*i)