Luna Flow
日本語

float_backend チュートリアル

このチュートリアルでは Complex[Double] の解析関数を使って計算します。直交形式と極形式の相互変換、根と対数の計算、2 次方程式の解法、三角関数とその逆関数の評価を行い、分岐切断がどこにあるかを学びます。式とその数値的な扱いは float_backend の設計 にあります。

クイックスタート

moon add Luna-Flow/luna-complex@0.2.0
import {
  "Luna-Flow/luna-complex" @complex,
  "Luna-Flow/luna-complex/float_backend" @fb,
}
fn main {
  let z = @complex.Complex::new(-3.0, 4.0)
  println("|z|     = \{@fb.abs(z)}")
  println("sqrt(z) = \{@fb.sqrt(z)}")
  println("exp(z)  = \{@fb.exp(z)}")
}
|z|     = 5
sqrt(z) = 1 + 2i
exp(z)  = -0.032542999640154786 + -0.03767897757486585i

日常的なタスク

極形式

fn main {
  let z = @complex.Complex::new(1.0, 1.0)
  let r = @fb.abs(z)
  let theta = @fb.arg(z)
  println("r = \{r}, theta = \{theta}")
  println("back: \{@fb.polar(r, theta)}")
}
r = 1.4142135623730951, theta = 0.7853981633974483
back: 1.0000000000000002 + 1i

2 次方程式を解く

az2+bz+caz^2 + bz + c の根は (−b±b2−4ac)/(2a)(-b \pm \sqrt{b^2 - 4ac})/(2a) です。sqrt と div を使えば、この式は複素係数や負の判別式に対しても機能します:

fn main {
  // z^2 + 2z + 5 = 0
  let a = @complex.Complex::new(1.0, 0.0)
  let b = @complex.Complex::new(2.0, 0.0)
  let c = @complex.Complex::new(5.0, 0.0)
  let four = @complex.Complex::new(4.0, 0.0)
  let two_a = a + a
  let root = @fb.sqrt(b * b - four * a * c)
  println("z1 = \{@fb.div(-b + root, two_a)}")
  println("z2 = \{@fb.div(-b - root, two_a)}")
}
z1 = -1 + 2i
z2 = -1 + -2i

対数とべき乗

fn main {
  let z = @complex.Complex::new(0.0, 1.0)
  println("log(i)   = \{@fb.log(z)}")
  println("i^i      = \{@fb.pow(z, z)}")
  println("i^2      = \{@fb.pow_real(z, 2.0)}")
  println("log10(1000) = \{@fb.log_10(@complex.Complex::new(1000.0, 0.0))}")
}
log(i)   = 0 + 1.5707963267948966i
i^i      = 0.20787957635076193 + 0i
i^2      = -1 + 0i
log10(1000) = 2.9999999999999996 + 0i

ii=ei⋅iπ/2=e−π/2i^i = e^{i \cdot i\pi/2} = e^{-\pi/2} は実数です。整数の指数には厳密な二進累乗法を使うため、i2i^2 はちょうど −1-1 になります。

三角関数とその逆関数

fn main {
  let z = @complex.Complex::new(0.5, 0.5)
  let s = @fb.sin(z)
  println("sin(z)       = \{s}")
  println("asin(sin(z)) = \{@fb.asin(s)}")
  println("asin(2)      = \{@fb.asin_real(2.0)}")
}
sin(z)       = 0.5406126857131534 + 0.4573041531842493i
asin(sin(z)) = 0.5 + 0.4999999999999999i
asin(2)      = 1.5707963267948966 + 1.3169578969248166i

asin_real は Double を受け取り、[−1,1][-1, 1] の外側では複素数の結果を返します。

大きな引数

式はスケーリングされているため、非常に大きい入力や非常に小さい入力でもオーバーフローしません:

fn main {
  let huge = @complex.Complex::new(1.0e300, 1.0e300)
  println("|huge|     = \{@fb.abs(huge)}")
  println("log|huge|  = \{@fb.abs_log(huge)}")
  println("tan(1+50i) = \{@fb.tan(@complex.Complex::new(1.0, 50.0))}")
  let q = @fb.div(@complex.Complex::new(1.0, 0.0), @complex.Complex::new(1.0e-300, 1.0e-300))
  println("1/(tiny)   = \{q}")
}
|huge|     = 1.4142135623730952e+300
log|huge|  = 691.1221014884936
tan(1+50i) = 6.765311025183565e-44 + 1i
1/(tiny)   = 4.9999999999999995e+299 + -4.9999999999999995e+299i

さらに進んで

分岐切断

多価関数は分岐切断をまたぐと値が跳びます。sqrt の分岐切断は負の実軸上にあり、上側と下側から近づくと符号が逆の根が得られます:

fn main {
  println(@fb.sqrt(@complex.Complex::new(-4.0, 1.0e-12)))
  println(@fb.sqrt(@complex.Complex::new(-4.0, -1.0e-12)))
  println(@fb.sqrt(@complex.Complex::new(-4.0, 0.0)))
}
2.5e-13 + 2i
2.5e-13 + -2i
0 + 2i

分岐切断上そのものでは根は +2i+2i です。各関数の分岐切断は 設計 に一覧があります。

スカラーについてジェネリックなヘルパーを書く

能力トレイトを使うと、ヘルパーが実スカラーに何を必要とするかを表明できます。例えば、特殊値を調べるだけの関数は FloatingSpecialValues を要求し、Float と Double の両方で動作します:

fn[T : @fb.FloatingSpecialValues] finite_parts(re : T, im : T) -> Bool {
  !(@fb.FloatingSpecialValues::is_nan(re) || @fb.FloatingSpecialValues::is_inf(re) ||
  @fb.FloatingSpecialValues::is_nan(im) || @fb.FloatingSpecialValues::is_inf(im))
}

fn main {
  let z = @fb.exp(@complex.Complex::new(800.0, 1.0))
  println(finite_parts(z.re, z.im))
  println(finite_parts((1.0 : Float), (2.0 : Float)))
}
false
true

よくある落とし穴

  • 負の実軸。 そこでは arg と log は現在 π\pi ではなく 2π2\pi を返し(log(-1) は 2πi2\pi i)、整数でない指数による負の実数のべき乗はその角度を受け継ぎます。acos、acos_real、asec_real、acosh_real も負の入力に対して影響を受けます。float_backend API の警告を参照してください。
  • 逆数系の関数は極で中断する。 cot(0)、csc(0)、asec(0) などは、無限大を返さずにプログラムを停止させます。
  • 極端な値には / ではなく @fb.div を使う。 コアの演算子は c2+d2c^2 + d^2 を計算し、それがゼロにアンダーフローすると中断します。
  • exp は早くオーバーフローする。 Re⁡z>709.78\operatorname{Re} z > 709.78 に対する exp(z) は無限大または NaN の部分を持ちます。
  • Complex[Double] のみ。 解析関数は Complex[Float] を受け付けません。

次のステップ