-- ----------------------------------------------------------------------- -- Title: double_support_functions -- Last Mod: Tue Apr 2 08:03:18 1991 -- Author: Vincent Broman -- Copyright 1990 Vincent Broman -- Permission granted to copy, modify, or compile this software for -- one's own use, provided that this copyright notice is preserved intact. -- Permission granted to distribute compiled binary copies of this -- software which are linked in with some other application. -- Permission granted to distribute other copies of this software, -- provided that (1) any copy which is not source code, i.e. not in the -- form in which the software is usually maintained, must be accompanied -- by a copy of the source code from which it was compiled, and (2) the -- one distributing it must refrain from imposing on the recipient -- further restrictions on the distribution of this software. -- -- Description: -- -- a double-precision-only implementation of functions -- supporting an implementation of generic_elementary_functions. -- -- Exceptions: No checking done. An accidental overflow could be caused -- by an argument way outside the spec domain of a function. -- ----------------------------------------------------------------------- -- with math_constants, double_primitive_functions; use math_constants, double_primitive_functions; package body double_support_functions is function reduced_logdel_2 (y: in double) return double is -- -- return the base-2 logarithm of (1+y)/(1-y), -- y in the interval [ -3+sqrt(8) .. 3-sqrt(8) ]. -- the coefficients are from a chebychev economized polynomial approximation -- to log_2( (1+y)/(1-y) )/y with relative error (omitting roundoff) < 9.7e-18. -- this delivers high relative accuracy for small x. -- Note: reduced_logdel_2( 0.0) = 0.0 exactly. -- p0: constant := 2.8853900817779268; p1: constant := 0.96179669392598977; p2: constant := 0.57707801634548050; p3: constant := 0.41219858584899499; p4: constant := 0.32059853412766626; p5: constant := 0.26233439170687922; p6: constant := 0.22091211514218184; p7: constant := 0.21366836734234657; -- y2: double := y * y; begin return y * (p0 + y2 * (p1 + y2 * (p2 + y2 * (p3 + y2 * (p4 + y2 * (p5 + y2 * (p6 + y2 * p7))))))); end reduced_logdel_2; function reduced_log_2 (x: in double) return double is -- -- return the base-2 logarithm of x, x in the interval [sqrt(1/2),sqrt(2)]. -- begin return reduced_logdel_2( (x - 1.0) / (x + 1.0)); end reduced_log_2; procedure log_2_parts (int_part: out double; frac_part: out double; x: in double) is -- -- return the integer and fractional part of the -- base-2 logarithm of x. -- int_part is a floating pt integer, and abs( frac_part) <= 1/2. -- ex: integer := exponent( x); ma: double := mantissa( x); sqrt_half: constant := sqrt_2 / 2.0; begin if ma < sqrt_half then ma := 2.0 * ma; ex := ex - 1; end if; -- sqrt_half <= ma < sqrt_2 -- and x = ma * 2 ** ex int_part := double( ex); frac_part := reduced_log_2( ma); end log_2_parts; function log_2 (x: in double) return double is -- -- return the base-2 logarithm of x, assuming x > 0.0 -- -- the value is exact for x an integer power of 2. -- n, f: double; begin -- log_2 log_2_parts( n, f, x); return n + f; end log_2; function reduced_two_to_the (x: in double) return double is -- -- return 2.0**x for x in the interval [-0.524, 0.524]. -- the coefficients are a Remez-optimized rational(8,2) approximation, -- with relative error (omitting roundoff) < 7.9e-18 . -- The interval would be [-1/2, 1/2] but we allow some slop. -- p0: constant := 187.371286015885939; p1: constant := 103.902608044422697; p2: constant := 28.0082502747857420; p3: constant := 4.8535553356237245; p4: constant := 0.60076449452675390; p5: constant := 0.0555224326757751590; p6: constant := 0.00384831074260226161; p7: constant := 0.000190479152083616768; p8: constant := 0.0000054921973342160052; q0: constant := p0; -- assures 2**0 = 1 q1: constant := -25.9732705753797810; begin return (p0 + x * (p1 + x * (p2 + x * (p3 + x * (p4 + x * (p5 + x * (p6 + x * (p7 + x * p8)))))))) / (q0 + x * (q1 + x)); end reduced_two_to_the; function two_to_the (x: in double) return double is -- -- return 2.0**x, underflowing to zero for very negative x, -- and overflowing for very positive x. -- ix: double := floor( x + 0.5); frac: double := x - ix; begin return scale( reduced_two_to_the( frac), integer( ix)); end two_to_the; function reduced_sin (x: in double) return double is -- -- return sin x for x in the interval [0,pi/4]. -- the coefficients are a chebychev-economized taylor series for -- sin((1-x)*pi/8)/((1-x)*pi/8) on [-1,1], -- with absolute error (omitting roundoff) < 9.9e-19. -- this gives high relative accuracy for small x. -- -- scale factor for size of interval [0,pi/4] already reckoned in below a0: constant := 9.74495358404432644e-01; a1: constant := 1.28892142233166662e-01; a2: constant := -1.59026534208421952e-01; a3: constant := -1.28510922444164585e-02; a4: constant := 7.87893558955578277e-03; a5: constant := 4.58212136760670637e-04; a6: constant := -1.86638149794387119e-04; a7: constant := -8.47651507694942824e-06; a8: constant := 2.58375120256637880e-06; a9: constant := 9.62537066878746960e-08; a10: constant := -2.33664460675007227e-08; a11: constant := -7.37964897324995069e-10; pio8: constant := pi / 8.0; y: constant double := pio8 - x; begin return x * (a0 + y * (a1 + y * (a2 + y * (a3 + y * (a4 + y * (a5 + y * (a6 + y * (a7 + y * (a8 + y * (a9 + y * (a10 + y * a11))))))))))); end reduced_sin; function reduced_cos (x: in double) return double is -- -- return cos x for x in the interval [0,pi/4]. -- the coefficients are a chebychev-economized taylor series for -- cos((1-x)*pi/8) on [-1,1], -- with absolute error (omitting roundoff) < 8.1e-20 -- -- scale factor for size of interval [0,pi/4] already reckoned in below a0: constant := 9.23879532511286756e-01; a1: constant := 3.82683432365089769e-01; a2: constant := -4.61939766255643376e-01; a3: constant := -6.37805720608478201e-02; a4: constant := 3.84949805213034383e-02; a5: constant := 3.18902860301776349e-03; a6: constant := -1.28316601736991906e-03; a7: constant := -7.59292519052434937e-05; a8: constant := 2.29136787544634377e-05; a9: constant := 1.05456702512338564e-06; a10: constant := -2.54595221365290357e-07; a11: constant := -9.55627162087695516e-09; a12: constant := 1.92304873377856603e-09; pio8: constant := pi / 8.0; y: constant double := pio8 - x; begin return a0 + y * (a1 + y * (a2 + y * (a3 + y * (a4 + y * (a5 + y * (a6 + y * (a7 + y * (a8 + y * (a9 + y * (a10 + y * (a11 + y * a12))))))))))); end reduced_cos; function reduced_series_arctan (x: in double) return double is -- -- returns arctan x in radians for x in [-tan pi/20, tan pi/20]. -- relative error < 6.5e-20 plus roundoff. -- the coefficients come from a remez-optimized rational[3,3] approximation -- of arctan(x)/x, minimimizing relative error. -- x2: constant double := x * x; p0: constant := 1.24554717745473876e+01; p1: constant := 1.59002191199817663e+01; p2: constant := 4.87692830958320933e+00; p3: constant := 2.07772731164829364e-01; q0: constant := p0; q1: constant := 2.00520430448308923e+01; q2: constant := 9.06984830295275812e+00; -- q3 would be 1 begin return x * ((p0 + x2 * (p1 + x2 * (p2 + x2 * p3))) / (q0 + x2 * (q1 + x2 * (q2 + x2)))); end reduced_series_arctan; function reduced_arctan (x: in double) return double is -- -- returns arctan x in radians from [-pi/4,pi/4] for x in [-1,1]. -- the reduction formula below comes from the angle-sum formula for tangent. -- the interval [0,1] is broken twice at tan(pi/20) and tan(3 pi/20). -- ax: double := abs( x); aat: double; tan_pio20: constant := 0.15838444032453629384; tan_3pio20: constant := 0.50952544949442881051; -- the angles following are accurate arctangents of the rounded -- 51 bit approximations to the tangent of 2*pi/20 and 4*pi/20. tan_2pio20: constant := 1463308222879852.0 / 2.0 ** 52; -- 51 bits d_2pio20: constant := 0.31415926535897938744; tan_4pio20: constant := 1636028329196881.0 / 2.0 ** 51; -- 51 bits d_4pio20: constant := 0.62831853071795872905; begin -- reduced_arctan if ax <= tan_3pio20 then if ax <= tan_pio20 then aat := reduced_series_arctan( ax); else aat := d_2pio20 + reduced_series_arctan( (ax - tan_2pio20) / (1.0 + ax * tan_2pio20)); end if; else aat := d_4pio20 + reduced_series_arctan( (ax - tan_4pio20) / (1.0 + ax * tan_4pio20)); end if; return copy_sign( aat, x); end reduced_arctan; end double_support_functions; -- $Header: d_support_functions_b.a,v 3.24 91/04/04 13:22:36 broman Exp $