lib/flt/math.rb in flt-1.1.2 vs lib/flt/math.rb in flt-1.2.0
- old
+ new
@@ -1,12 +1,7 @@
require 'flt/dec_num'
-# TODO:
-# * convert arguments as Context#_convert does, to accept non DecNum arguments
-# * Better precision (currently is within about 2 ulps of the correct result)
-# * Speed optimization
-
module Flt
class DecNum
module Math
extend Flt # to access constructor methods DecNum
@@ -16,11 +11,10 @@
# Trigonometry
HALF = DecNum('0.5')
# Pi
- $no_cache = false
@@pi_cache = nil # truncated pi digits as a string
@@pi_cache_digits = 0
PI_MARGIN = 10
def pi(round_digits=nil)
@@ -66,11 +60,14 @@
local_context.precision = decimals if decimals
DecNum(1).exp
end
end
- # Cosine of angle in radians
+ # Cosine of an angle given in the units specified by DecNum.context.angle:
+ # * :rad for radians
+ # * :deg for degrees
+ # * :grad for gradians
def cos(x)
x = x.abs
rev_sign = false
s = nil
DecNum.context do |local_context|
@@ -80,10 +77,11 @@
if k % 2 == 0
x = pi_2 - x
else
rev_sign = !rev_sign
end
+ x = to_rad(x)
i, lasts, fact, num = 1, 0, 1, DecNum(x)
s = num
x2 = -x*x
while s != lasts
lasts = s
@@ -94,21 +92,24 @@
end
end
return rev_sign ? (-s) : (+s)
end
- # Sine of angle in radians
+ # Sine of an angle given in the units specified by DecNum.context.angle:
+ # * :rad for radians
+ # * :deg for degrees
+ # * :grad for gradians
def sin(x)
sign = x.sign
s = nil
DecNum.context do |local_context|
local_context.precision += 3 # extra digits for intermediate steps
x = x.abs if sign<0
x,k,pi_2 = reduce_angle2(x,2)
sign = -sign if k>1
x = pi_2 - x if k % 2 == 1
- x = +x
+ x = to_rad(x)
i, lasts, fact, num = 1, 0, 1, DecNum(x)
s = num
x2 = -x*x
while s != lasts
lasts = s
@@ -119,31 +120,35 @@
end
end
return (+s).copy_sign(sign)
end
+ # Tangent of an angle
def tan(x)
+DecNum.context do |local_context|
local_context.precision += 2 # extra digits for intermediate steps
s,c = sin(x), cos(x)
s/c
end
end
- # Arc-tangent.
+ # Arc-tangent in units specified by DecNum.context.angle
def atan(x)
s = nil
+ conversion = true
DecNum.context do |local_context|
local_context.precision += 2
if x == 0
return DecNum.zero
elsif x.abs > 1
if x.infinite?
- s = (pi*HALF).copy_sign(x)
+ s = (quarter_cycle).copy_sign(x)
+ conversion = false
break
else
- c = (pi*HALF).copy_sign(x)
+ # c = (quarter_cycle).copy_sign(x)
+ c = (HALF*pi).copy_sign(x)
x = 1 / x
end
end
local_context.precision += 2
x_squared = x ** 2
@@ -159,49 +164,49 @@
end
if c && c!= 0
s = c - s
end
end
- return +s
+ return conversion ? rad_to(s) : +s
end
def atan2(y, x)
abs_y = y.abs
abs_x = x.abs
y_is_real = !x.infinite?
if x != 0
if y_is_real
a = y!=0 ? atan(y / x) : DecNum.zero
- a += pi.copy_sign(y) if x < 0
+ a += half_cycle.copy_sign(y) if x < 0
return a
elsif abs_y == abs_x
x = DecNum(1).copy_sign(x)
y = DecNum(1).copy_sign(y)
- return pi * (2 - x) / (4 * y)
+ return half_cycle * (2 - x) / (4 * y)
end
end
if y != 0
return atan(DecNum.infinity(y.sign))
elsif x < 0
- return pi.copy_sign(x)
+ return half_cycle.copy_sign(x)
else
return DecNum.zero
end
end
def asin(x)
x = +x
return DecNum.context.exception(Num::InvalidOperation, 'asin needs -1 <= x <= 1') if x.abs > 1
if x == -1
- return -pi*HALF
+ return -quarter_cycle
elsif x == 0
return DecNum.zero
elsif x == 1
- return pi*HALF
+ return quarter_cycle
end
DecNum.context do |local_context|
local_context.precision += 3
x = x/(1-x*x).sqrt
@@ -210,33 +215,25 @@
+x
end
def acos(x)
- # We can compute acos(x) = pi/2 - asin(x)
- # but we must take care with x near 1, where that formula would cause loss of precision
-
return DecNum.context.exception(Num::InvalidOperation, 'acos needs -1 <= x <= 2') if x.abs > 1
if x == -1
- return pi
+ return half_cycle
elsif x == 0
- return +DecNum.context(:precision=>DecNum.context.precision+3){pi*HALF}
+ return quarter_cycle
elsif x == 1
return DecNum.zero
end
- # some identities:
- # acos(x) = pi/2 - asin(x) # (but this losses accuracy near x=+1)
- # acos(x) = pi/2 - atan(x/(1-x*x).sqrt) # this too
- # acos(x) = asin((1-x*x).sqrt) for x>=0; for x<=0 acos(x) = pi/2 - asin((1-x*x).sqrt)
-
if x < HALF
DecNum.context do |local_context|
local_context.precision += 3
x = x/(1-x*x).sqrt
- x = pi*HALF - atan(x)
+ x = quarter_cycle - atan(x)
end
else
# valid for x>=0
DecNum.context do |local_context|
local_context.precision += 3
@@ -253,11 +250,11 @@
local_context.precision += 3
(x*x + y*y).sqrt
end
end
- # TODO: degrees mode or radians/degrees conversion
+ # TODO: add angular units to context; add support for degrees
def pi2(decimals=nil)
decimals ||= DecNum.context.precision
DecNum.context(:precision=>decimals) do
pi(decimals)*2
@@ -296,38 +293,11 @@
end
sum
end
def modtwopi(x)
- return +DecNum.context(:precision=>DecNum.context.precision*3){x.modulo(pi2)}
- # This seems to be slower and less accurate:
- # prec = DecNum.context.precision
- # pi_2 = pi2(prec*2)
- # return x if x < pi_2
- # ex = x.fractional_exponent
- # DecNum.context do |local_context|
- # # x.modulo(pi_2)
- # local_context.precision *= 2
- # if ex > prec
- # # consider exponent separately
- # fd = nil
- # excess = ex - prec
- # x = x.scaleb(prec-ex)
- # # now obtain 2*prec digits from inv2pi after the initial excess digits
- # digits = nil
- # inv_2pi = inv2pi(local_context.precision+excess)
- # DecNum.context do |extended_context|
- # extended_context.precision += excess
- # digits = (inv2pi.scaleb(excess)).fraction_part
- # end
- # x *= digits*pi_2
- # end
- # # compute the fractional part of the division by 2pi
- # inv_2pi ||= inv2pi
- # x = pi_2*((x*inv2pi).fraction_part)
- # end
- # +x
+ return +DecNum.context(:precision=>DecNum.context.precision*3){x.modulo(one_cycle)}
end
# Reduce angle to [0,2Pi)
def reduce_angle(a)
modtwopi(a)
@@ -336,13 +306,124 @@
# Reduce angle to [0,Pi/k0) (result is not rounded to precision)
def reduce_angle2(a,k0=nil) # divisor of pi or nil for pi*2
# we could reduce first to pi*2 to avoid the mod k0 operation
k,r,divisor = DecNum.context do
DecNum.context.precision *= 3
- m = k0.nil? ? pi2 : pi/k0
+ m = k0.nil? ? one_cycle : half_cycle/k0
a.divmod(m)+[m]
end
[r, k.modulo(k0*2).to_i, divisor]
+ end
+
+ def one_cycle
+ case DecNum.context.angle
+ when :rad
+ pi2
+ when :deg
+ DecNum(360)
+ when :grad
+ DecNum(400)
+ end
+ end
+
+ def half_cycle
+ case DecNum.context.angle
+ when :rad
+ pi
+ when :deg
+ DecNum(180)
+ when :grad
+ DecNum(200)
+ end
+ end
+
+ def quarter_cycle
+ case DecNum.context.angle
+ when :rad
+ HALF*pi
+ when :deg
+ DecNum(90)
+ when :grad
+ DecNum(100)
+ end
+ end
+
+ def to_rad(x)
+ case DecNum.context.angle
+ when :rad
+ +x
+ else
+ +DecNum.context(:precision=>DecNum.context.precision+3){x*pi/half_cycle}
+ end
+ end
+
+ def to_deg(x)
+ case DecNum.context.angle
+ when :deg
+ +x
+ else
+ +DecNum.context(:precision=>DecNum.context.precision+3){x*DecNum(180)/half_cycle}
+ end
+ end
+
+ def to_grad(x)
+ case DecNum.context.angle
+ when :deg
+ +x
+ else
+ +DecNum.context(:precision=>DecNum.context.precision+3){x*DecNum(200)/half_cycle}
+ end
+ end
+
+ def to_angle(angular_units, x)
+ return +x if angular_units == DecNum.context.angle
+ case angular_units
+ when :rad
+ to_rad(x)
+ when :deg
+ to_deg(x)
+ when :grad
+ to_grad(x)
+ end
+ end
+
+ def rad_to(x)
+ case DecNum.context.angle
+ when :rad
+ +x
+ else
+ +DecNum.context(:precision=>DecNum.context.precision+3){x*half_cycle/pi}
+ end
+ end
+
+ def deg_to(x)
+ case DecNum.context.angle
+ when :deg
+ +x
+ else
+ +DecNum.context(:precision=>DecNum.context.precision+3){x*half_cycle/DecNum(180)}
+ end
+ end
+
+ def grad_to(x)
+ case DecNum.context.angle
+ when :grad
+ +x
+ else
+ +DecNum.context(:precision=>DecNum.context.precision+3){x*half_cycle/DecNum(200)}
+ end
+ end
+
+ def angle_to(x, angular_units)
+ return +x if angular_units == DecNum.context.angle
+ case angular_units
+ when :rad
+ rad_to(x)
+ when :deg
+ deg_to(x)
+ when :grad
+ grad_to(x)
+ end
end
#end
end # Math
\ No newline at end of file