三角関数群の実装
from マルチサイクル RISC-V CPU を作成したい
レイトレの実行に必要な sin cos atan の3つの三角関数を実装する。
マクローリン展開
sin, cos, atan の各関数は、マクローリン展開した以下の式から求めることができる。
三角関数
min_caml_cos (float -> float)
cos x のマクローリン展開
https://gyazo.com/21f2d4fee58157f01376454cfcf1ce91
min_caml_sin (float -> float) → 三角関数群の実装
sin x のマクローリン展開
https://gyazo.com/ff74e6fcfd06ef5cac694db486f9005c
min_caml_atan (float -> float)
Arctanのマクローリン展開
https://gyazo.com/66bb8342f18b07c1b12e0b088b2a63d3
Rubyでsin関数の練習
マクローリン展開を使って、Rubyでsin関数を実装。
sin x のマクローリン展開
https://gyazo.com/ff74e6fcfd06ef5cac694db486f9005c
code:sin.rb
def my_sin(x)
x2 = x * x
x3 = x2 * x
x5 = x3 * x2
x7 = x5 * x2
x9 = x7 * x2
x11 = x9 * x2
x13 = x11 * x2
x - x3 / 6 + x5 / 120 - x7 / 5040 + x9 / 362880 - x11 / 39916800 + x13 / 6227020800
end
# Math.sin と比較（CSVへ出力）
open("sin.csv", "w") do |f|
f.puts "x,my_sin,Math.sin"
-12.step(12, 0.1) do |i|
x = i.to_f
f.puts "#{x},#{my_sin(x)},#{Math.sin(x)}"
end
end
x の絶対値が大きくなるほど誤差が増える。sin関数は同じ波が繰り返されるので、x の値を 「-3.14 から 3.14」の範囲に帰着させてあげたい
https://gyazo.com/62bbff8c067880f8deacd142801547ca
code:sin.rb
def my_sin(x)
x = adjustx(x)
x2 = x * x
x3 = x2 * x
x5 = x3 * x2
x7 = x5 * x2
x9 = x7 * x2
x11 = x9 * x2
x13 = x11 * x2
x - x3 / 6 + x5 / 120 - x7 / 5040 + x9 / 362880 - x11 / 39916800 + x13 / 6227020800
end
# x が -PI から PI までの範囲に収まるようにして返す
# 例)
# x = 3 → 3
# x = -3 → -3
# x = 3.2 → -3.083185307179586
# x = -3.2 → 3.083185307179586
# x = 6.5 → 0.641592653589793
# x = -6.5 → -0.641592653589793
def adjustx(x)
if x > Math::PI
adjustx(x - 2 * Math::PI)
elsif x < -Math::PI
adjustx(x + 2 * Math::PI)
else
x
end
end
# Math.sin と比較（CSVへ出力）
open("sin.csv", "w") do |f|
f.puts "x,my_sin,Math.sin"
-30.step(30, 0.1) do |i|
x = i.to_f
f.puts "#{x},#{my_sin(x)},#{Math.sin(x)}"
end
end
https://gyazo.com/4a56b29eadf9e3fd6119d5240b3513a5
計算量を減らしてみた。9次まで減らしても大丈夫そう。
code:sin.rb
def my_sin(x)
x = adjustx(x)
x2 = x * x
x3 = x2 * x
x5 = x3 * x2
x7 = x5 * x2
x9 = x7 * x2
x - x3 / 6 + x5 / 120 - x7 / 5040 + x9 / 362880
end
# x が -PI から PI までの範囲に収まるようにして返す
# 例)
# x = 3 → 3
# x = -3 → -3
# x = 3.2 → -3.083185307179586
# x = -3.2 → 3.083185307179586
# x = 6.5 → 0.641592653589793
# x = -6.5 → -0.641592653589793
def adjustx(x)
if x > Math::PI
adjustx(x - 2 * Math::PI)
elsif x < -Math::PI
adjustx(x + 2 * Math::PI)
else
x
end
end
# Math.sin と比較（CSVへ出力）
open("sin.csv", "w") do |f|
f.puts "x,my_sin,Math.sin"
-15.step(15, 0.1) do |i|
x = i.to_f
f.puts "#{x},#{my_sin(x)},#{Math.sin(x)}"
end
end
OCamlで実装
code:sin.ml
(* sin 関数 *)
let rec sin x =
(* x が -PI から PI の間に収まるよう調整する *)
let rec adjust x =
if x > 3.141592653589793 then
adjust (x -. 6.283185307179586)
else if x < -3.141592653589793 then
adjust (x +. 6.283185307179586)
else
x
in
let x = adjust x in
let x2 = x *. x in
let x3 = x *. x2 in
let x5 = x3 *. x2 in
let x7 = x5 *. x2 in
let x9 = x7 *. x2 in
x -. x3 /. 6.0 +. x5 /. 120.0 -. x7 /. 5040.0 +. x9 /. 362880.0
in
(* 度数法からラジアンへ変換 *)
let rec deg_to_rad deg =
deg *. 3.141592653589793 /. 180.0
in
(* sin(0°) * 1000 *)
print_int (int_of_float (sin (deg_to_rad 0.0) *. 10000.0));
print_newline ();
(* sin(30°) * 1000 *)
print_int (int_of_float (sin (deg_to_rad 30.0) *. 10000.0));
print_newline ();
(* sin(45°) * 1000 *)
print_int (int_of_float (sin (deg_to_rad 45.0) *. 10000.0));
print_newline ();
(* sin(60°) * 1000 *)
print_int (int_of_float (sin (deg_to_rad 60.0) *. 10000.0));
print_newline ();
(* sin(90°) * 1000 *)
print_int (int_of_float (sin (deg_to_rad 90.0) *. 10000.0));
print_newline ();
(* アジャストのテスト *)
print_int (int_of_float (sin (deg_to_rad 0.0) *. 10000.0)); (* => 0 *)
print_newline ();
print_int (int_of_float (sin (deg_to_rad 720.0) *. 10000.0)); (* => 0 *)
print_newline ()
cos 関数
Rubyで cos 関数。adjustx は sin と同じものが使える。
code:cos.rb
def my_cos(x)
x = adjustx(x)
x2 = x * x
x4 = x2 * x2
x6 = x4 * x2
x8 = x6 * x2
x10 = x8 * x2
1 - x2 / 2 + x4 / 24 - x6 / 720 + x8 / 40320 - x10 / 3628800
end
# x が -PI から PI までの範囲に収まるようにして返す
# 例)
# x = 3 → 3
# x = -3 → -3
# x = 3.2 → -3.083185307179586
# x = -3.2 → 3.083185307179586
# x = 6.5 → 0.641592653589793
# x = -6.5 → -0.641592653589793
def adjustx(x)
if x > Math::PI
adjustx(x - 2 * Math::PI)
elsif x < -Math::PI
adjustx(x + 2 * Math::PI)
else
x
end
end
# Math.cos と比較（CSVへ出力）
open("cos.csv", "w") do |f|
f.puts "x,my_cos,Math.cos"
-10.step(10, 0.1) do |i|
x = i.to_f
f.puts "#{x},#{my_cos(x)},#{Math.cos(x)}"
end
end
atan 関数
（うまく実装できなかったが、呼ばれてる雰囲気もなかったので、常に1.0を返してお茶を濁した）
code:atan.rb
# マクローリン展開でatan(x)を求める
def my_atan(x)
x2 = x * x
x3 = x2 * x
x5 = x3 * x2
x7 = x5 * x2
x9 = x7 * x2
x11 = x9 * x2
x13 = x11 * x2
x15 = x13 * x2
x17 = x15 * x2
x19 = x17 * x2
x21 = x19 * x2
x23 = x21 * x2
x25 = x23 * x2
x27 = x25 * x2
x29 = x27 * x2
x31 = x29 * x2
x33 = x31 * x2
x35 = x33 * x2
x37 = x35 * x2
x39 = x37 * x2
x41 = x39 * x2
x43 = x41 * x2
x45 = x43 * x2
x47 = x45 * x2
x49 = x47 * x2
x51 = x49 * x2
x53 = x51 * x2
x55 = x53 * x2
x57 = x55 * x2
x59 = x57 * x2
x61 = x59 * x2
x63 = x61 * x2
x65 = x63 * x2
x67 = x65 * x2
x69 = x67 * x2
x71 = x69 * x2
x73 = x71 * x2
x75 = x73 * x2
x77 = x75 * x2
x79 = x77 * x2
x81 = x79 * x2
x83 = x81 * x2
x85 = x83 * x2
x87 = x85 * x2
x89 = x87 * x2
x91 = x89 * x2
x93 = x91 * x2
x95 = x93 * x2
x97 = x95 * x2
x99 = x97 * x2
x - x3 / 3 + x5 / 5 - x7 / 7 + x9 / 9 - x11 / 11 + x13 / 13 - x15 / 15 + x17 / 17 - x19 / 19 + x21 / 21 - x23 / 23 + x25 / 25 - x27 / 27 + x29 / 29 - x31 / 31 + x33 / 33 - x35 / 35 + x37 / 37 - x39 / 39 + x41 / 41 - x43 / 43 + x45 / 45 - x47 / 47 + x49 / 49 - x51 / 51 + x53 / 53 - x55 / 55 + x57 / 57 - x59 / 59 + x61 / 61 - x63 / 63 + x65 / 65 - x67 / 67 + x69 / 69 - x71 / 71 + x73 / 73 - x75 / 75 + x77 / 77 - x79 / 79 + x81 / 81 - x83 / 83 + x85 / 85 - x87 / 87 + x89 / 89 - x91 / 91 + x93 / 93 - x95 / 95 + x97 / 97 - x99 / 99
end
# Math.atan と比較（CSVへ出力）
open("atan.csv", "w") do |f|
f.puts "x,my_atan,Math.atan"
# (-1.5).step(1.5, 0.1) do |i|
(-1.1).step(1.1, 0.1) do |i|
x = i.to_f
f.puts "#{x},#{my_atan(x)},#{Math.atan(x)}"
# f.puts "#{x},#{Math.atan(x)},#{Math.atan(x)}"
end
end
atan は繰り返しが無いから、いい感じに収束させれない...どうしたものか...
→ atanに渡される値を確認
atan 呼ばれてなさそうだから、適当な値を返そうかな... → 常に1.0を返すようにしたけど、問題なかった
https://gyazo.com/9d00185e027318dc972ad2aaed2c2c88
三角関数表
https://www.math.s.chiba-u.ac.jp/~yasuda/sysKOU/cit-H20/trig-table.pdf
https://gyazo.com/52b74a9326601b73bf8c257da459a253
「角度の正規化」の有無
x の絶対値が大きくなるほど誤差が増える。sin関数は同じ波が繰り返されるので、x の値を 「-3.14 から 3.14」の範囲に帰着させてあげると良い。これを「角度の正規化」と呼ぶらしい。
正規化なしだとこんな感じで、
https://gyazo.com/5d963233fea7b639fdbafec188337d87
正規化ありだとこんな感じ。
https://gyazo.com/7064e89ae3400677c3375a51e1ab0ccf
正規化なしのコードと
code:ruby
def my_sin(x)
x2 = x * x
x3 = x2 * x
x5 = x3 * x2
x7 = x5 * x2
x9 = x7 * x2
x - x3 / 6 + x5 / 120 - x7 / 5040 + x9 / 362880
end
正規化ありのコード
code:ruby
def my_sin(x)
x = degree_normalize(x)
x2 = x * x
x3 = x2 * x
x5 = x3 * x2
x7 = x5 * x2
x9 = x7 * x2
x - x3 / 6 + x5 / 120 - x7 / 5040 + x9 / 362880
end
# x が -PI から PI までの範囲に収まるようにして返す（角度の正規化）
# 例)
# x = 3 → 3
# x = -3 → -3
# x = 3.2 → -3.083185307179586
# x = -3.2 → 3.083185307179586
# x = 6.5 → 0.641592653589793
# x = -6.5 → -0.641592653589793
def degree_normalize(x)
if x > Math::PI
degree_normalize(x - 2 * Math::PI)
elsif x < -Math::PI
degree_normalize(x + 2 * Math::PI)
else
x
end
end