-- ----------------------------------------------------------------------- -- Title: single_support_functions -- Last Mod: Tue Apr 2 14:33:40 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 single-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, single_primitive_functions; use math_constants, single_primitive_functions; package body single_support_functions is function reduced_logdel_2 (y: in single) return single 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 rational (even) 1,1 approximation to -- log_2( (1+y)/(1-y) )/y with relative error (omitting roundoff) < 1.91e-8. -- this delivers high relative accuracy for small y. -- Note: reduced_logdel_2( 0.0) = 0.0 exactly. -- p0: constant := -4.7845025; p1: constant := 1.2906108; q0: constant := -1.6581823; y2: constant single := y * y; begin return y * (p0 + p1 * y2) / (q0 + y2); end reduced_logdel_2; function reduced_log_2 (x: in single) return single 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 single; frac_part: out single; x: in single) 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: single := 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 := single( ex); frac_part := reduced_log_2( ma); end log_2_parts; function log_2 (x: in single) return single 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: single; begin -- log_2 log_2_parts( n, f, x); return n + f; end log_2; function reduced_two_to_the (x: in single) return single is -- -- return 2.0**x for x in the interval [-0.524, 0.524]. -- the coefficients are a Remez-optimized rational(3,2) approximation, -- with relative error (omitting roundoff) < 1e-8. -- The interval would be [-1/2, 1/2] but we allow some slop. -- p0: constant := 41.694187; p1: constant := 17.342325; p2: constant := 3.0047081; p3: constant := 0.23083161; q0: constant := p0; -- assures 2**0 = 1 q1: constant := -11.5578823; begin return (p0 + x * (p1 + x * (p2 + x * p3))) / (q0 + x * (q1 + x)); end reduced_two_to_the; function two_to_the (x: in single) return single is -- -- return 2.0**x, underflowing to zero for very negative x, -- and overflowing for very positive x. -- ix: single := floor( x + 0.5); frac: single := x - ix; begin return scale( reduced_two_to_the( frac), integer( ix)); end two_to_the; function reduced_sin (x: in single) return single 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) < 2e-8. -- this gives high relative accuracy for small x. -- -- scale factor for size of interval [0,pi/4] already reckoned in below a0: constant := 9.74495337094254e-01; a1: constant := 1.28892138844516e-01; a2: constant := -1.59024045807998e-01; a3: constant := -1.28509163870500e-02; a4: constant := 7.83587003622433e-03; a5: constant := 4.55929708079455e-04; pio8: constant := pi / 8.0; y: constant single := pio8 - x; begin return x * (a0 + y * (a1 + y * (a2 + y * (a3 + y * (a4 + y * a5))))); end reduced_sin; function reduced_cos (x: in single) return single 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) < 1.8e-9 -- -- scale factor for size of interval [0,pi/4] already reckoned in below a0: constant := 9.23879532410434e-01; a1: constant := 3.82683402033700e-01; a2: constant := -4.61939745322990e-01; a3: constant := -6.37789978530365e-02; a4: constant := 3.84943015536361e-02; a5: constant := 3.16859376215302e-03; a6: constant := -1.27611586733778e-03; pio8: constant := pi / 8.0; y: constant single := pio8 - x; begin return a0 + y * (a1 + y * (a2 + y * (a3 + y * (a4 + y * (a5 + y * a6))))); end reduced_cos; function reduced_arctan (x: in single) return single is -- -- returns arctan x in radians from [-pi/4,pi/4] for x in [-1,1]. -- coefficents from Abramowitz and Stegun p81, #4.4.49 . -- relative error < 2e-8. -- x2: constant single := x * x; -- a2: constant := -0.3333314528; a4: constant := 0.1999355085; a6: constant := -0.1420889944; a8: constant := 0.1065626393; a10: constant := -0.0752896400; a12: constant := 0.0429096138; a14: constant := -0.0161657367; a16: constant := 0.0028662257; begin return x * (1.0 + x2 * ( a2 + x2 * ( a4 + x2 * ( a6 + x2 * ( a8 + x2 * (a10 + x2 * (a12 + x2 * (a14 + x2 * a16)))))))); end reduced_arctan; end single_support_functions; -- $Header: s_support_functions_b.a,v 3.24 91/04/04 13:22:41 broman Exp $