cos.cr

require "benchmark"

def fake_cos(x : Float64) : Float64
  1.0
end

def taylor_cos(x : Float64) : Float64
  x = x.abs
  kf, x = x.divmod(Math::PI)
  k = kf.to_i
  sign = 1.0-(k&1)*2
  x2 = x*x
  sign * (1 + 
         x2*(1.0/2)*(-1 +
         x2*(1.0/(3*4))*(1 +
         x2*(1.0/(5*6))*(-1 +
         x2*(1.0/(7*8))*(1 +
         x2*(1.0/(9*10))*(-1 +
         x2*(1.0/(11*12))*(1 +
         x2*(1.0/(13*14))*(-1))))))))
end

F4 = 4.16666666666666666e-2  # 1.0/(2*3*4)
F5 = 8.33333333333333333e-3  # 1.0/(2*3*4*5)
F8 = 2.48015873015873015e-5  # 1.0/(2*3*4*5*6*7*8)
F9 = 2.7557319223985893e-6   # 1.0/(2*3*4*5*6*7*8*9)
F12 = 2.08767569878681e-9    # 1.0/(2*3*4*5*6*7*8*9*10*11*12)
F13 = 1.6059043836821613e-10 # 1.0/(2*3*4*5*6*7*8*9*10*11*12*13)
M2 = 0.5                     # 1.0/2
M3 = 1.66666666666666666e-1  # 1.0/(2*3)
M6 = 3.33333333333333333e-2  # 1.0/(5*6)
M7 = 2.3809523809523808e-2   # 1.0/(6*7)
M10 = 1.1111111111111111e-2  # 1.0/(9*10)
M11 = 9.0909090909090909e-3  # 1.0/(10*11)
M14 = 5.4945054945054945e-3  # 1.0/(13*14)
M15 = 4.7619047619047619e-3  # 1.0/(14*15)

def taylor_cos2(x : Float64) : Float64
  x = x.abs
  kf, x = x.divmod(Math::PI)
  k = kf.to_i
  sign = 1.0-(k&1)*2
  x2 = x*x
  x4 = x2*x2
  x8 = x4*x4
  sign * (1 - x2*M2 + 
         (1 - x2*M6)*x4*F4 +
         (1 - x2*M10)*x8*F8 +
         (1 - x2*M14)*x8*x4*F12)
end

def taylor_cos3(x : Float64) : Float64
  x = x.abs
  kf, x = x.divmod(Math::PI/2)
  k = kf.to_i
  sign = 1.0-((k+1)&2)
  if k&1 == 1
    x = Math::PI/2 - x
  end

  if x < Math::PI/4
    x2 = x*x
    x4 = x2*x2
    co = (1 - x2 * M2 + ((1 - x2 * M6) * F4 + x4*F8)*x4)
    sign * co
  else
    x = Math::PI/2 - x
    x2 = x*x
    x4 = x2*x2
    si = x * (1 - x2 * M3 + ((1 - x2 * M7) * F5 + x4*F9)*x4)
    sign * si
  end
end

def sum_cos
  sum = 1.0
  10000.times{|i| sum += yield i/1000.0 }
  sum
end

20.times do |i|
  c = Math.cos(i/2.0)
  p [i/2.0, c, ((taylor_cos2(i/2.0)-c)/c).abs, ((taylor_cos3(i/2.0)-c)/c).abs]
end

#exit

total = 0.0
total1 = 0.0
total2 = 0.0
total3 = 0.0
Benchmark.ips do |x|
  #x.report("fake") { total += sum_cos{|f| fake_cos(f)} }
  x.report("std") { total += sum_cos{|f| Math.cos(f)} }
  x.report("cos") { total1 += sum_cos{|f| taylor_cos(f)} }
  x.report("cos2") { total2 += sum_cos{|f| taylor_cos2(f)} }
  x.report("cos3") { total3 += sum_cos{|f| taylor_cos3(f)} }
end
p [total, total1, total2, total3]