/* ==================================================================== ** Implement the standard and inverse trigonometric and hyperbolic ** functions according to the recipes of William Kahan in unpublished ** notes ``Elementary Functions from Kernels'' (24-Oct-1985). The ** formula on the second page of that paper for sin(y) is incorrect; a ** correction is provided below. In his second formula for atan(x), the ** multiplier s is wrong on the first term: it must be applied to the ** entire formula. ** ** Usage: ** gcc -Wmissing-prototypes -Wall -pedantic -g -g3 \ ** [-DUSE_FLOAT | -DUSE_DOUBLE | -DUSE_LONG_DOUBLE ] \ ** [-DUSE_FMA] \ ** [-DUSE_HP_COS] \ ** [-DUSE_HP_SIN] \ ** [-DUSE_SIN_T_FORM] \ ** [-DUSE_PI_HALF_SUM | -DUSE_TWOSUM] \ ** wkkernels.c -lm && ./a.out [args] > foo.out ** ** If a working precision is not chosen, the default is -DUSE_DOUBLE. ** ** If USE_FMA is defined, the fused multiply-add function is used ** instead of the equivalent traditional expressions --- Kahan's paper ** predates the invention of the FMA. ** ** If USE_HP_COS is defined, the three largest terms in the cosine ** series are computed in the highest available precision, and the final ** sum is then rounded to working precision. ** ** If USE_HP_SIN is defined, the three largest terms in the sine series are ** computed in the highest available precision, and the final sum is ** then rounded to working precision. ** ** If USE_SIN_T_FORM is defined, for small arguments, sin(x) is computed ** from the T(x) function; otherwise, it is computed from Taylor series ** summation. ** ** If USE_TWOSUM is defined, PI/2 - w is computed with error ** compensation using the 2sum algorithm; otherwise, it is computed from ** (PI_HALF_HI - w) + PI_HALF_LO. This option overrides any definition ** of USE_PI_HALF_SUM. ** ** If USE_PI_HALF_SUM is defined, PI/2 - w is computed with error ** compensation from an extended precision version of PI/2. This ** option is ignored if USE_TWOSUM is also defined. ** ** If USE_TWOSUM or USE_PI_HALF_SUM is defined, PI/2 - w is computed ** EXACTLY; otherwise, for w in [0, 1/2], errors in that difference ** increase from 0.500 to 0.770 ulps. ** ** The default executable, a.out, accepts one to four optional ** arguments: ** ** (1) [integer] N (sample size); ** (2) [integer] 0 for no plot data, nonzero for plot data; ** (3) [string] suffix for plot filenames; ** (4) [string] one of FE_DOWNWARD, FE_TONEAREST, FE_TOWARDZERO, ** or FE_UPWARD (default if omitted, or invalid, ** is FE_TONEAREST). ** ** The plot data files are recorded in "/tmp/TESTNAME.dat.*.suffix.mode" ** files, for TESTNAME = wk_cos, wk_sin, ..., The numeric suffix after ** ".dat." is the floating-point precision in bits, and any requested ** suffix appears next, with a separating dot, and any nondefault ** rounding mode, with a separating dot, ends the filename. ** ** Summarize the largest errors with: ** ** % grep rmax foo.out | sort -k3g | tail ** ** Make data plots like this [point plots are the default, so ** the "using" clause is not needed]: ** ** % gnuplot ** gnuplot> set pointsize 0.25 ** gnuplot> set key font "large, 20" ** gnuplot> set xrange [-2:2] ** gnuplot> set yrange [0:2] ** gnuplot> plot "/tmp/wk_acos.dat.24" ** gnuplot> plot "/tmp/wk_acos.dat.24" using ($1):($2) ** gnuplot> set term x11 1 ** gnuplot> plot 0.5, 1, "/tmp/wk_acos.dat.24" using ($1):($2) ** gnuplot> set term x11 2 ** gnuplot> plot 0.5, 1, "/tmp/wk_asin.dat.24" using ($1):($2) ** gnuplot> set term x11 3 ** gnuplot> plot 0.5, 1, "/tmp/wk_atan.dat.24" using ($1):($2) ** ** The required kernels are atan(), copysign(), exp(), expm1(), ** fabs(), log(), log1p(), sqrt(), tan(), and possibly, sin() for ** small arguments. ** ** The atan(x) function is needed only for x <= sqrt(2) - 1 ~= 0.414, ** for which at most 47 terms are needed to sum its Taylor series in ** binary128 arithmetic, and fewer terms in shorter formats. ** ** The tan(x) function is only needed for |x| <= PI/8, which can be ** evaluated by a Taylor series with at most 57 terms for binary128, or ** 33 (binary80), 27 (binary64), and 13 (binary32). It does, however, ** need a compile-time table of coefficients, because they are ** complicated values depending on Bernouilli numbers, with no simple ** term recurrence. ** ** Kahan's functions assume that trigonometric argument reduction has ** already been performed: ** ** cos, sin, tan: arguments in [-PI/2, PI/2] ** ** I have added code in each to return distinct QNaNs for arguments ** outside that range, and set errno = EDOM. ** ** Kahan's functions ensure that any needed argument transformations are ** EXACT in base-2 arithmetic, and are designed to (mostly) guarantee ** monotonicity. Several functions need an argument reduction y = PI/2 ** - x, and that must be done with a double-length value of PI/2, split ** into exact high, and accurate low, parts. We supply that reduction ** with the function PI_HALF_minus_w(w). It has three different ** implementations, selectable at compile time by -DUSE_xxx options. ** ** Kahan assumes IEEE 754 default of round to nearest, ties to even, so ** tiny arguments reduce to the first term of the function's Taylor ** series. I changed all of those assignments to use the first two ** terms to accommodate other rounding modes. ** ** [25-Aug-2025] ** ================================================================== */ #include #include #include #include #include #include #include static const char * rounding_mode = "FE_TONEAREST"; static const char * suffix = (const char *)NULL; static FILE * fpdata = (FILE *)NULL; static int save_data = 0; static int N = 10240; /* NB: overridable from command line */ /* snoop loop counters */ static int __K_ACOS__ = 0; /* in wk_acos_series() */ static int __K_ACOS_ASIN__ = 0; /* in wk_acos_asin_series() */ static int __K_ASIN__ = 0; /* in wk_asin_series() */ static int __K_ATAN__ = 0; /* in wk_atan_series() */ static int __K_COS__ = 0; /* in wk_cos_series() */ static int __K_SIN__ = 0; /* in wk_sin_series() */ static int n_exact = 0; static int n_finite = 0; // extern double qnan(const char *payload); typedef long double hp_t; #define CONS(a,b) (a##b) #define HP(x) CONS(x, L) #if defined(USE_LONG_DOUBLE) typedef long double fp_t; static const fp_t EPS = LDBL_EPSILON; static const fp_t MAXNORMAL = LDBL_MAX; static const fp_t MINNORMAL = LDBL_MIN; static const fp_t MINSUBNORMAL = LDBL_TRUE_MIN; #define P LDBL_MANT_DIG #define FP(x) CONS(x, L) #elif defined(USE_FLOAT) typedef float fp_t; static const fp_t EPS = FLT_EPSILON; static const fp_t MAXNORMAL = FLT_MAX; static const fp_t MINNORMAL = FLT_MIN; static const fp_t MINSUBNORMAL = FLT_TRUE_MIN; #define P FLT_MANT_DIG #define FP(x) CONS(x, F) #else /* assume USE_DOUBLE */ typedef double fp_t; #define FP(x) (x) static const fp_t EPS = DBL_EPSILON; static const fp_t MAXNORMAL = DBL_MAX; static const fp_t MINNORMAL = DBL_MIN; static const fp_t MINSUBNORMAL = DBL_TRUE_MIN; #define P DBL_MANT_DIG #endif /* defined(USE_LONG_DOUBLE) */ static const fp_t __TINY__ = FP(1.0) / (fp_t)(1LL << ((P + 1)/2)); #define FP128(x) CONS(x, L) #define FP80(x) CONS(x, L) #define FP64(x) (x) #define FP32(x) CONS(x, F) static const fp_t PI = (long double)3.141592653589793238462643383279502884197L; /* ** Convention: UPPERCASE function names are ALWAYS the highest ** available precision, and are used for accuracy comparisons against ** other functions in this file. Their type, hp_t, is always the ** highest available precision. The type fp_t is intended to be that ** of the precision of functions to be tested. */ #define COPYSIGN copysignl #define E expm1l #define EXP expl #define EXPM1 expm1l #define FABS fabsl #define FMA fmal #define ISINF isinf /* generic macro: no type suffix */ #define ISNAN isnan /* generic macro: no type suffix */ #define L log1pl #define LOG logl #define QNAN qnanl #define SQRT sqrtl #define ACOS acosl #define ASIN asinl #define ATAN atanl #define ACOSH acoshl #define ASINH asinhl #define ATANH atanhl #define COS cosl #define SIN sinl #define TAN tanl #define COSH coshl #define SINH sinhl #define TANH tanhl long double wk_tan_series_128 (long double x); long double wk_tan_series_80 (long double x); double wk_tan_series_64 (double x); float wk_tan_series_32 (float x); float qnanf (const char * payload); double qnan (const char * payload); long double qnanl (const char * payload); hp_t hp_PI_HALF_minus_w_alt (hp_t w); hp_t hp_twosum (hp_t a, hp_t b, hp_t * err); fp_t PI_HALF_minus_w (fp_t w); fp_t PI_HALF_minus_w_alt (fp_t w); fp_t T (fp_t x); fp_t oneulp (void); fp_t randl (fp_t x); fp_t randlab (fp_t a, fp_t b); fp_t range_error (const char * payload); fp_t twosum (fp_t a, fp_t b, fp_t * err); fp_t unirand (void); fp_t wk_acos (fp_t x); fp_t wk_acos_asin_series (fp_t x); fp_t wk_acos_series (fp_t x); fp_t wk_asin (fp_t x); fp_t wk_asin_series (fp_t x); fp_t wk_atan (fp_t x); fp_t wk_atan_series (fp_t x); fp_t wk_acosh (fp_t x); fp_t wk_acosh_alt (fp_t x); fp_t wk_asinh (fp_t x); fp_t wk_asinh_alt (fp_t x); fp_t wk_atanh (fp_t x); fp_t wk_cos (fp_t x); fp_t wk_cos_series (fp_t x); fp_t wk_sin (fp_t x); fp_t wk_sin_series (fp_t x); fp_t wk_tan (fp_t x); fp_t wk_tan_series (fp_t x); fp_t wk_cosh (fp_t x); fp_t wk_sinh (fp_t x); fp_t wk_tanh (fp_t x); fp_t wk_fma (fp_t x, fp_t y, fp_t z); fp_t wk_sqrt (fp_t x); void begin (const char * name); void end (void); void report (const char * name, const char * suffix, fp_t lo, fp_t hi, int n, fp_t xmax, fp_t amax, fp_t rmax, fp_t rsum); void save (fp_t x, fp_t r, fp_t a, hp_t e); void show_options (void); void test_PI_HALF_minus_w (void); void test_PI_HALF_minus_w_alt (void); void test_signed_zero (void); void test_zero (const char * name, fp_t wk_any(fp_t x), int symmetric); void test_wk_any (const char * name, fp_t wk_any(fp_t x), hp_t any(hp_t x), fp_t lo, fp_t hi, int n); void test_wk_any_rand (const char * name, fp_t wk_any(fp_t x), hp_t any(hp_t x), fp_t lo, fp_t hi, int n); void test_wk_acos (void); void test_wk_acos_series (void); void test_wk_acosh (void); void test_wk_acosh_alt (void); void test_wk_asin (void); void test_wk_asin_series (void); void test_wk_asinh (void); void test_wk_asinh_alt (void); void test_wk_atan (void); void test_wk_atan_series (void); void test_wk_atanh (void); void test_wk_cos (void); void test_wk_cos_series (void); void test_wk_cosh (void); void test_wk_sin (void); void test_wk_sin_series (void); void test_wk_sinh (void); void test_wk_tan (void); void test_wk_tan_series_128 (void); void test_wk_tan_series_32 (void); void test_wk_tan_series_64 (void); void test_wk_tan_series_80 (void); void test_wk_tan_series (void); void test_wk_tanh (void); static const fp_t ZERO = FP(0.0); static const fp_t FOURTH = FP(0.25); static const fp_t HALF = FP(0.5); static const fp_t ONE = FP(1.0); static const fp_t TWO = FP(2.0); static const fp_t THREE = FP(3.0); static const fp_t FOUR = FP(4.0); static const fp_t SIX = FP(6.0); #define PRINTF (void)printf #define FMT_dat "%.9g" #define FMT_L_ulps "%.6Lg" #if P == 113 #define FMT_L_g "%.36Lg" #define FMT_L_a "%.29La" static fp_t PI_HALF_HI = FP(4078020416507594100416871540687068.0) / FP(2596148429267413814265248164610048.0); static fp_t PI_HALF_LO = FP(4.33585846996875529104874722961539082031431045e-35); #elif P == 64 #define FMT_L_g "%.21Lg" #define FMT_L_a "%.16La" static fp_t PI_HALF_HI = FP(14488038916154245684.0) / FP(9223372036854775808.0); static fp_t PI_HALF_LO = FP(8.33374291852087832820e-20); #elif P == 53 #define FMT_L_g "%.17Lg" #define FMT_L_a "%.13La" static fp_t PI_HALF_HI = FP(7074237752028440.0) / FP(4503599627370496.0); static fp_t PI_HALF_LO = FP(6.1232339957367660e-17); #elif P == 24 #define FMT_L_g "%.9Lg" #define FMT_L_a "%.6La" static fp_t PI_HALF_HI = FP(13176794.0) / FP(8388608.0); static fp_t PI_HALF_LO = FP(7.54979013e-08); #else #error "No split of PI/2 implemented for this P" #endif /* P == 113 */ #define PI_HALF (PI_HALF_HI + PI_HALF_LO) #define PI_FOURTH (FP(0.50) * (PI_HALF_HI + PI_HALF_LO)) #define PI_EIGHTH (FP(0.25) * (PI_HALF_HI + PI_HALF_LO)) void begin(const char * name) { /* action at start of block of tests of a single function */ char s[FILENAME_MAX]; const char * mode; if (rounding_mode == (const char *)NULL) mode = (const char *)NULL; else if (strcmp(rounding_mode, "FE_TONEAREST") == 0) mode = (const char *)NULL; else mode = rounding_mode; errno = 0; n_exact = 0; n_finite = 0; if (save_data) { if (fpdata != (FILE *)NULL) (void)fclose(fpdata); (void)snprintf(s, sizeof(s), "/tmp/%s.dat.%d%s%s%s%s", name, (int)P, ((suffix == (const char *)NULL) ? "" : "."), ((suffix == (const char *)NULL) ? "" : suffix), ((mode == (const char *)NULL) ? "" : "."), ((mode == (const char *)NULL) ? "" : mode)); fpdata = fopen(s, "w"); } } void end(void) { /* action at end of block of tests of a single function */ if (fpdata != (FILE *)NULL) { (void)fclose(fpdata); fpdata = (FILE *)NULL; } if (errno != 0) PRINTF("\nAt test end, errno = %d: %s\n", errno, strerror(errno)); if ((n_finite > 0) && (n_exact > 0)) { PRINTF("Total finite results = %d\n", n_finite); PRINTF("Total exact results = %d\n", n_exact); PRINTF("Correctly rounded results = %.2f%%\n", 100.0 * (double)n_exact / (double)n_finite); } } fp_t oneulp(void) { return (EPS); } fp_t PI_HALF_minus_w(fp_t w) { /* return accurate value of PI/2 - w */ fp_t result; #if defined(USE_TWOSUM) fp_t err1, err2, sum1, sum2; volatile fp_t err3; sum1 = twosum(-w, PI_HALF_LO, &err1); sum2 = twosum(PI_HALF_HI, sum1, &err2); err3 = err1 + err2; result = sum2 + err3; #elif defined(USE_PI_HALF_SUM) result = PI_HALF_minus_w_alt(w); #else result = (PI_HALF_HI - w) + PI_HALF_LO; #endif /* defined(USE_TWOSUM) */ return (result); } fp_t PI_HALF_minus_w_alt(fp_t w) { /* ** Compute PI/2 - w from a multiple-precision representation of PI/2 ** as a sum of EXACT 20- or 21-bit terms, summing from smallest to ** largest. The hexadecimal values of the terms are computed like ** this: ** ** % maple ** > Digits := 75: ** > c := evalf(Pi/2): ** > ** > scale := 2**20: ** > ** > for k from 0 to 6 do ** > h := trunc(scale * c): ** > printf("FP(0x%06xp%+d),\n", h, -(k + 1) *20): ** > c := scale * c - h: ** > end do: */ fp_t err, sum, t; int k; static const fp_t pi_half[] = { FP(0x1921fbp-20), FP(0x054442p-40), FP(0x0d1846p-60), FP(0x09898cp-80), FP(0x0c5170p-100), FP(0x01b839p-120), FP(0x0a2520p-140) }; sum = -w; err = ZERO; for (k = sizeof(pi_half) / sizeof(pi_half[0]) - 1; k >= 0; --k) { sum = twosum(sum, pi_half[k], &t); err += t; } return (sum + err); } hp_t hp_PI_HALF_minus_w_alt(hp_t w) { /* ** % maple ** > Digits := 75: ** > c := evalf(Pi/2): ** > ** > scale := 2**20: ** > ** > for k from 0 to 6 do ** > h := trunc(scale * c): ** > printf("FP(0x%06xp%+d),\n", h, -(k + 1) *20): ** > c := scale * c - h: ** > end do: */ hp_t err, sum, t; int k; static const hp_t pi_half[] = { HP(0x1921fbp-20), HP(0x054442p-40), HP(0x0d1846p-60), HP(0x09898cp-80), HP(0x0c5170p-100), HP(0x01b839p-120), HP(0x0a2520p-140) }; sum = -w; err = HP(0.0); for (k = sizeof(pi_half) / sizeof(pi_half[0]) - 1; k >= 0; --k) { sum = hp_twosum(sum, pi_half[k], &t); err += t; } return (sum + err); } hp_t hp_twosum(hp_t a, hp_t b, hp_t * err) { /* result + err = exact sum of a + b */ volatile hp_t s, v; if (FABS(b) > FABS(a)) { /* ensure |a| >= |b| */ s = b; b = a; a = s; } s = a + b; v = s - a; *err = b - v; return (s); } float qnanf(const char * payload) { volatile float x, y; x = y = ZERO; return (x / y); } double qnan(const char * payload) { volatile double x, y; x = y = ZERO; return (x / y); } long double qnanl(const char * payload) { volatile long double x, y; x = y = ZERO; return (x / y); } fp_t range_error(const char * payload) { errno = EDOM; return (QNAN(payload)); } fp_t twosum(fp_t a, fp_t b, fp_t * err) { /* result + err = exact sum of a + b */ volatile fp_t s, v; if (FABS(b) > FABS(a)) { /* ensure |a| >= |b| */ s = b; b = a; a = s; } s = a + b; v = s - a; *err = b - v; return (s); } fp_t unirand(void) { return ((fp_t)drand48()); } fp_t randl(fp_t x) { /* ** Return a pseudo random number logarithmically distributed over ** (1, exp(x)). Thus a*Randl(ln(b/a)) is logarithmically ** distributed in (a, b). */ return (EXP(x * unirand())); } fp_t randlab(fp_t a, fp_t b) { /* TO DO: test the three alternatives */ /* ** Return a pseudo random number logarithmically distributed over ** (a, b). */ #if 1 return (a * randl(LOG(b / a))); #elif 0 return (a * randl(LOG(FABS(b)) - LOG(FABS(a)))); #else static fp_t a_save = ZERO; static fp_t b_save = ZERO; static fp_t log_b_a = ZERO; if ( (a_save == ZERO) || (b_save == ZERO) || (a != a_save) || (b != b_save) ) { a_save = a; b_save = b; log_b_a = LOG(FABS(b)) - LOG(FABS(a)); } return (a * randl(log_b_a)); #endif } fp_t T(fp_t x) { return (TWO * wk_tan_series(HALF * x)); } fp_t wk_acos(fp_t x) { fp_t result; if (ISNAN(x) || (FABS(x) > ONE)) result = range_error("0xdead0001"); else result = TWO * wk_atan(wk_sqrt((ONE - x) / (ONE + x))); return (result); } fp_t wk_acos_asin_series(fp_t x) { /* ** % maple ** > taylor(arccos(x), x = 0, 8); ** Pi 3 5 7 9 ** ---- - x - 1/6 x - 3/40 x - 5/112 x + O(x ) ** 2 ** ** > tackx := proc(k, x) ** return ((2*k)! * x**(2*k + 1)) / ** (2**(2*k) * (k!)**2 * (2*k + 1)) ** end proc: ** ** > tackx(0,x); tackx(1,x); tackx(2,x); tackx(3,x); ** ** x ** ** 3 ** x ** ---- ** 6 ** ** 5 ** 3 x ** ---- ** 40 ** ** 7 ** 5 x ** ---- ** 112 ** ** > simplify(tackx(k, x) / tackx(k - 1, x)); ** 2 2 ** x (2 k - 1) ** ------------- ** 2 ** 4 k + 2 k ** ** > simplify(arccos(-x) + arccos(x)); ** Pi ** ** The series does not converge computationally for |x| > 1/2, but ** for that limit we find: ** ** > for k from 0 to 55 by 5 do ** printf("%2d\t%15.7g\n", k, tackx(k, 1/2)) ** end do: ** 0 0.5 ** 5 1.092391e-05 ** 10 4.000824e-09 ** 15 2.170048e-12 ** 20 1.390536e-15 ** 25 9.776508e-19 ** 30 7.292818e-22 ** 35 5.668289e-25 ** 40 4.540705e-28 ** 45 3.722558e-31 ** 50 3.108160e-34 ** 55 2.633930e-37 ** ** For the cutoff used in wk_acos_series() and wk_asin_series(), ** x = 183/256, we need 104 terms to binary128 convergence, 55 for ** binary80, 44 for binary64, and 17 for binary32. */ fp_t new_sum, result, sum, t0, t1, t2, tk, xx; int k, kkm1; __K_ACOS_ASIN__ = 0; if (ISNAN(x) || (FABS(x) > FP(31.0) / FP(32.0))) result = range_error("0xdead0002"); else { xx = x * x; t0 = x; t1 = t0 * xx / SIX; t2 = t1 * FP(9.0) * xx / FP(20.0); tk = t2; sum = COPYSIGN(ZERO, x); for (k = 3; k < 1055; ++k) { kkm1 = 2 * k - 1; tk *= (fp_t)(kkm1 * kkm1) * xx / (fp_t)(4 * k * k + 2 * k); new_sum = sum + tk; if (new_sum == sum) break; sum = new_sum; } sum += t2; sum += t1; sum += t0; result = sum; __K_ACOS_ASIN__ = k; } return (result); } fp_t wk_acos_series(fp_t x) { fp_t result, w; /* ** Series convergence is slow for x > 1/2, so use acos(x) = ** asin(sqrt(1 - x**2)) (if x > 0): a suitable compromise cutoff ** for series length is x = 180/256 = 0x1.68p-1 = 0.703_125, for ** which sqrt(1 - x**2) = 0.711_06... However, to avoid infinite ** recursion, we use a cutoff here of 183/256 = +0x1.6ep-1 = ** 0.714_843_75. */ __K_ACOS__ = 0; w = FABS(x); if (ISNAN(x) || (w > ONE)) result = range_error("0xdead0003"); else if (w > FP(+0x1.6ep-1)) { result = wk_asin_series(wk_sqrt(wk_fma(-x, x, ONE))); /* for x < 0, acos(x) = PI - asin(sqrt(1 - x**2)) */ if (x < ZERO) result = TWO * PI_HALF_minus_w(HALF * result); __K_ACOS__ = __K_ASIN__; } else { result = PI_HALF_minus_w(wk_acos_asin_series(x)); __K_ACOS__ = __K_ACOS_ASIN__; } return (result); } fp_t wk_acosh(fp_t x) { fp_t result; if (ISNAN(x) || (x < ONE)) result = range_error("0xdead0004"); else if (ISINF(TWO * x)) result = TWO * x; else result = L(wk_sqrt(x - ONE) * (wk_sqrt(x - ONE) + wk_sqrt(x + ONE))); return (result); } fp_t wk_acosh_alt(fp_t x) { fp_t result; if (ISNAN(x) || (x < ONE)) result = range_error("0xdead0005"); else if (x > ONE / __TINY__) result = LOG(x) + LOG(TWO); else if ( (FP(1.25) < x) && (x <= ONE / __TINY__) ) result = LOG(TWO * x - ONE / (x + wk_sqrt(wk_fma(x, x, -ONE)))); else { fp_t xm1; xm1 = x - ONE; result = L(xm1 + wk_sqrt(wk_fma(xm1, xm1, TWO * xm1))); } return (result); } fp_t wk_asin(fp_t x) { fp_t result, w; w = FABS(x); if (ISNAN(x) || (w > ONE)) result = range_error("0xdead0006"); else if (x == ZERO) result = x; else if (w < __TINY__) result = x + x * x * x / SIX; else if ( (__TINY__ <= w) && (w <= HALF) ) result = wk_atan(x / wk_sqrt(wk_fma(-w, w, ONE))); else { fp_t one_minus_w; one_minus_w = ONE - w; result = wk_atan(x / wk_sqrt(wk_fma(-one_minus_w, one_minus_w, TWO * one_minus_w))); } return (result); } fp_t wk_asin_series(fp_t x) { fp_t result, s, w; /* ** Series convergence is slow for x > 1/2, so use asin(x) = ** acos(sqrt(1 - x**2)): a suitable compromise cutoff for series ** length is x = 180/256 = 0x1.68p-1 = 0.703_125, for which ** sqrt(1 - x**2) = 0.711_06... */ w = FABS(x); s = COPYSIGN(ONE, x); if (ISNAN(x) || (w > ONE)) result = range_error("0xdead0007"); else if (w > FP(+0x1.68p-1)) { result = s * wk_acos_series(wk_sqrt(wk_fma(-x, x, ONE))); __K_ASIN__ = __K_ACOS__; } else { result = wk_acos_asin_series(x); __K_ASIN__ = __K_ACOS_ASIN__; } return (result); } fp_t wk_asinh(fp_t x) { fp_t result, s, w; w = FABS(x); s = COPYSIGN(ONE, x); if (ISNAN(x)) result = range_error("0xdead0008"); else if (ISINF(TWO * w)) result = QNAN("0xdead0009"); else if (x == ZERO) result = x; else if (w < __TINY__) result = x - x * x * x / SIX; else { fp_t w_inv; w_inv = ONE / w; result = s * L(w + w / (w_inv + wk_sqrt(wk_fma(w_inv, w_inv, ONE)))); } return (result); } fp_t wk_asinh_alt(fp_t x) { fp_t result, s, w; /* # TO DO: check: wk_sinh(3e20) vs sinh(3e20) */ w = FABS(x); s = COPYSIGN(ONE, x); if (ISNAN(x)) result = range_error("0xdead0010"); else if (w < __TINY__) result = x - x * x * x / SIX; else if (ISINF(TWO * w)) result = s * TWO * w; else if (w > FOUR / THREE) { if (w < ONE / __TINY__) result = s * LOG(TWO * w + ONE / (w + wk_sqrt(wk_fma(w, w, ONE)))); else result = s * (LOG(w) + LOG(TWO)); } else { fp_t w_inv; w_inv = ONE / w; result = s * L(w + w / (ONE / w + wk_sqrt(wk_fma(w_inv, w_inv, ONE)))); } return (result); } fp_t wk_atan(fp_t x) { fp_t result, s, w; w = FABS(x); s = COPYSIGN(ONE, x); /* # NB: first test added by NHFB because */ /* # > taylor(arctan(x), x = 0, 10); */ /* # 3 5 7 9 11 */ /* # x - 1/3 x + 1/5 x - 1/7 x + 1/9 x + O(x ) */ if (ISNAN(x)) result = range_error("0xdead0011"); else if (x == ZERO) result = x; else if (w < __TINY__) result = x - x * x * x / THREE; else if (w > ONE) result = s * PI_HALF_minus_w(wk_atan(ONE / w)); else if ( ((wk_sqrt(TWO) - ONE) < w) && (w < (wk_sqrt(TWO) + ONE)) ) { // This is wrong in Kahan's paper: // result = s * PI_FOURTH + wk_atan((x - s) / (x + s)); result = s * (PI_FOURTH + wk_atan((x - s) / (x + s))); } else /* x in [__TINY__, sqrt(2) - 1] ~= [__TINY__, 0.4143] */ /* Kahan does not handle this range */ result = wk_atan_series(x); return (result); } fp_t wk_atan_series(fp_t x) { fp_t new_sum, result, sum, t0, t1, t1_more, t2, tk, xx, xxx; int k; static const fp_t XMAX = FP(0.41421356237309504880168872420969808); /* sqrt(2) - 1 ~= 0.414213... */ /* > series(arctan(x), x = 0, 10); ** 3 5 7 9 11 ** x - 1/3 x + 1/5 x - 1/7 x + 1/9 x + O(x ) ** ** > tkx := proc(k, x) ** local m: ** m := 2*k + 1: ** return (-1)**k * x**m / m ** end proc: ** ** > tkx(0,x); ** ** x ** ** > tkx(1,x); ** ** 3 ** x ** - ---- ** 3 ** ** > tkx(2,x); ** ** 5 ** x ** ---- ** 5 ** ** > simplify(tkx(k,x) / tkx(k-1,x)); ** ** 2 ** x (1 - 2 k) ** ------------ ** 2 k + 1 ** ** > u128 := 2**(-112); ** ** > for k from 40 to 47 do ** printf("%3d\t%8.3f\n", ** k, tkx(k, xmax) / tkx(0, xmax) / u128) ** end do: ** 40 15.305 ** 41 -2.563 ** 42 0.429 ** 43 -0.072 ** 44 0.012 ** 45 -0.002 ** 46 0.000 ** 47 -0.000 */ if (ISNAN(x) || (FABS(x) > XMAX)) result = range_error("0xdead0012"); else if (x == ZERO) result = x; else { xx = x * x; xxx = x * xx; t0 = x; /* # Improve accuracy: x**3 / 3 = x**3 / 4 + x**3 / 12 */ /* t1 = -xxx / 3 */ t1 = -xxx * FOURTH; t1_more = -xxx / FP(12.0); /* with xmax = SQRT(2) - 1, */ /* evalf(ta(2, xmax) / ta(1, xmax)) -> -0.1029 */ /* so it is likely not useful to split 1/5 = 1/4 - 1/20, */ /* and we need t2 anyway for the loop recurrence */ t2 = xxx * xx / FP(5.0); tk = t2; sum = COPYSIGN(ZERO, x); /* ** for largest x used in this file, ** x = SQRT(2) - 1 ~= 0.414_213..., ** we need at most 42 terms in binary128 arithmetic */ for (k = 3; k < 50; ++k) { tk *= -xx * (fp_t)(2*k - 1) / (fp_t)(2*k + 1); new_sum = sum + tk; if (new_sum == sum) break; sum = new_sum; } sum += t2; sum += t1_more; sum += t1; sum += t0; result = sum; __K_ATAN__ = k; } return (result); } fp_t wk_atanh(fp_t x) { fp_t result, s, w; w = FABS(x); s = COPYSIGN(ONE, x); /* # > taylor(arctanh(x), x = 0, 10); */ /* # 3 5 7 9 11 */ /* # x + 1/3 x + 1/5 x + 1/7 x + 1/9 x + O(x ) */ if (ISNAN(x) || (w > ONE)) result = range_error("0xdead0013"); else if (w < __TINY__) result = x + x * x * x / THREE; else result = s * HALF * L(TWO * w / (ONE - w)); return (result); } fp_t wk_cos(fp_t x) { fp_t q, result, y, w; w = FABS(x); if (ISNAN(x) || (w > PI_HALF)) result = range_error("0xdead0014"); else if (w < __TINY__) result = wk_fma(-HALF * x, x, ONE); else if ( (PI_FOURTH <= w) && (w <= PI_HALF) ) result = wk_sin(PI_HALF_minus_w(w)); else /* must have x in (__TINY__, PI / 4) */ { #if 1 /* this computation is about 0.5 ulps worse than COS(x), but still acceptable */ y = x; q = T(y); q *= q; if (q < FP(4.0) / FP(15.0)) result = ONE - TWO / (ONE + FOUR / q); else result = FP(0.75) + (wk_fma(-TWO, q, ONE) + FOURTH * q) / (FOUR + q); #else result = COS(x); #endif } return (result); } fp_t wk_cos_series(fp_t x) { /* ** % maple ** > taylor(cos(x), x = 0, 10); ** 2 4 6 8 10 ** 1 - 1/2 x + 1/24 x - 1/720 x + 1/40320 x + O(x ) ** ** > tcos := proc(k, x) (-1)**k * x**(2*k) / (2*k)! end proc: ** ** > for k from 0 to 4 do printf("%2d\t%a\n", k, tcos(k,x)) end do: ** 0 1 ** 1 -1/2*x^2 ** 2 1/24*x^4 ** 3 -1/720*x^6 ** 4 1/40320*x^8 ** ** > k := 'k': x := 'x': simplify(tcos(k, x) / tcos(k - 1, x)); ** 2 ** x ** - ---------- ** 2 ** 4 k - 2 k ** ** > u := 2**(-23): ** ** > for k from 4 to 8 do ** printf("%2d\t%.7f ulps32\n", k, tcos(k,Pi/2) / u) ** end do: ** 4 7711.3141040 ulps32 ** 5 -211.4100546 ulps32 ** 6 3.9517682 ulps32 ** 7 -0.0535747 ulps32 ** 8 0.0005508 ulps32 ** ** > for k from 4 to 8 do ** printf("%2d\t%.7f ulps32\n", k, tcos(k,Pi/4) / u) ** end do: ** 4 30.1223207 ulps32 ** 5 -0.2064551 ulps32 ** 6 0.0009648 ulps32 ** 7 -0.0000033 ulps32 ** 8 0.0000000 ulps32 ** ** The worst-case term counts for the 4 IEEE 754 formats are 5, ** 9, 10 and 15 for |x| = PI/4. Because of the rapid ** convergence, we could allow |x| = Pi/2, for which we need 7, ** 11, 13, and 18 terms. However, if x > 1, there has to be ** subtraction loss in the series to reach a result in [-1, 1], ** so a smaller limit on |x| is advisable. ** ** The first TWO terms contribute more than 97.80% of the function ** value, and THREE terms, 99.96%, at x = Pi/4: ** ** % maple ** > tcos := proc(k,x) (-1)**k * x**(2*k + 1) / (2*k + 1)! end proc; ** > plot((tcos(0,x) + tcos(1,x)) / cos(x), x = -Pi/4 .. Pi/4); ** ** ** HHHHHHHHHHHH1H*HHHHHHHHHHHHH ** HHHHH + HHHHH ** HHH + HHH ** HH + HH ** HHH + HHH ** HH 0.995 + HH ** HH + HH ** H + H ** HH + HH ** HH 0.99 + HH ** H + H ** HH + HH ** H + H ** HH + HH ** H 0.985 + H ** HH + HH ** H + H ** H + H ** HH 0.98 + HH ** H + H ** *++-++-++-+-++-++-++-+-++-++-++++-++-++-+-++-++-++-+-++-++-++* ** -0.6 -0.4 -0.2 0 0.2 0.4 0.6 ** ** Thus, improvements in accuracy must come primarily from the sum ** 1 - x**2/2 + x**4 / 24. */ fp_t new_sum, result, sum, t0, t1, t2, tk, xx; int k; __K_SIN__ = 0; if (ISNAN(x) || (FABS(x) > PI_FOURTH)) result = range_error("0xdead0015"); else { xx = x * x; t0 = ONE; t1 = -xx / TWO; t2 = t1 * (-xx / FP(12.0)); tk = t2; sum = ZERO; for (k = 3; k < 20; ++k) { tk *= -xx / (fp_t)(2 * k * (2 * k - 1)); new_sum = sum + tk; if (new_sum == sum) break; sum = new_sum; } #if defined(USE_HP_COS) { hp_t hp_sum, hp_t0, hp_t1, hp_t2, hp_x, hp_xx; hp_x = (hp_t)x; hp_xx = hp_x * hp_x; hp_t0 = HP(1.0); hp_t1 = -hp_xx * HP(0.5); hp_t2 = hp_t1 * (-hp_xx / HP(12.0)); hp_sum = (hp_t)sum; hp_sum += hp_t2; hp_sum += hp_t1; hp_sum += hp_t0; result = (fp_t)hp_sum; } #else sum += t2; sum += t1; sum += t0; result = sum; #endif /* defined(USE_HP_COS) */ __K_COS__ = k; } return (result); } fp_t wk_cosh(fp_t x) { fp_t h, result, w; w = FABS(x); h = HALF * EXP(w); result = h + FOURTH / (h); return (result); } fp_t wk_fma(fp_t x, fp_t y, fp_t z) { fp_t result; #if defined(USE_FMA) #if P == 113 result = fmal(x, y, z); #elif P == 64 result = fmal(x, y, z); #elif P == 53 result = fma(x, y, z); #elif P == 24 result = fmaf(x, y, z); #else result = QNAN("0xdead0016"); #endif #else /* prevent compiler use of FMA */ volatile fp_t xy; xy = x * y; result = xy + z; #endif /* defined(USE_FMA) */ return (result); } fp_t wk_sin(fp_t x) { fp_t result, s, w; w = FABS(x); s = COPYSIGN(ONE, x); if (ISNAN(x) || (w > PI_HALF)) result = range_error("0xdead0017"); else if (x == ZERO) result = x; else if (w < __TINY__) result = x - x * x * x / SIX; else if ( (PI_FOURTH <= w) && (w <= PI_HALF) ) result = s * wk_cos(PI_HALF_minus_w(w)); else /* w in [__TINY__, PI/4) */ { #if defined(USE_SIN_T_FORM) fp_t t; /* formula used in Intel 8087 (Palmer & Morse book, p. 102 */ /* error in ulps is about twice that of series sum */ t = T(x); result = t / (ONE + FOURTH * t * t); #else /* more accurate */ result = wk_sin_series(x); #endif /* defined(USE_SIN_T_FORM) */ } return (result); } fp_t wk_sin_series(fp_t x) { /* ** % maple ** > taylor(sin(x), x = 0, 10); ** 3 5 7 9 11 ** x - 1/6 x + 1/120 x - 1/5040 x + 1/362880 x + O(x ) ** ** > simplify(tk(k,x) / tk(k - 1,x)); ** 2 ** x ** - ---------- ** 2 ** 4 k + 2 k ** ** > for k from 1 to 6 do ** printf("%2d\t%a\n", k, tk(k,x)/tk(k-1,x)) ** end do; ** 1 -1/6*x^2 ** 2 -1/20*x^2 ** 3 -1/42*x^2 ** 4 -1/72*x^2 ** 5 -1/110*x^2 ** 6 -1/156*x^2 ** ** > for k from 2 to 15 do printf("%2d\t%15.7g\n", ** k, (tk(k, Pi/4) / tk(0, Pi/4))) ** end do: ** 2 0.003170869 ** 3 -4.657027e-05 ** 4 3.989845e-07 ** 5 -2.237397e-09 ** 6 8.847045e-12 ** 7 -2.598715e-14 ** 8 5.893450e-17 ** 9 -1.062975e-19 ** 10 1.561183e-22 ** 11 -1.903194e-25 ** 12 1.956642e-28 ** 13 -1.719310e-31 ** 14 1.306104e-34 ** 15 -8.663127e-38 ** ** The worst-case term counts for the 4 IEEE 754 formats are 5, ** 8, 10 and 15 for x = |PI/4|. Because of the rapid ** convergence, we could allow x = |Pi/2|, for which we need 6, ** 10, 12, and 18 terms. However, if x > 1, there has to be ** subtraction loss in the series to reach a result in [-1, 1], ** so a smaller limit on |x| is advisable. ** ** The first TWO terms contribute more than 99.65% of the function ** value: ** ** % maple ** > tkx := proc(k,x) (-1)**k * x**(2*k + 1) / (2*k + 1)! end proc; ** > plot((tkx(0,x) + tkx(1,x)) / sin(x), x = -Pi/4 .. Pi/4); ** ** HHHHHHHHHHHHHH1H*HHHHHHHHHHHHHHH ** HHHHH + HHHHH ** HHH + HHH ** HHH 0.9995 + HHH ** HH + HH ** HHH 0.999 + HHH ** HH + HH ** HH + HH ** H 0.9985 + H ** HH + HH ** HH + HH ** HH 0.998 + HH ** H + H ** HH 0.9975 + HH ** H + H ** HH + HH ** H 0.997 + H ** HH + HH ** H 0.9965 + H ** -+-++--0.6+-++--0.4++-+--0.2-+-+-+0++-+-+-0.2+-+-++0.4+-++-+0.6+-++-+- ** ** Thus, improvements in accuracy must come primarily from the sum ** x - x**3 / 6. */ fp_t new_sum, result, sum, t0, t1, t2, tk, xx; int k; __K_SIN__ = 0; if (ISNAN(x) || (FABS(x) > PI_FOURTH)) result = range_error("0xdead0018"); else if (x == ZERO) result = x; else { xx = x * x; t0 = x; t1 = t0 * (-xx / SIX); t2 = t1 * (-xx / FP(20.0)); tk = t2; sum = COPYSIGN(ZERO, x); for (k = 3; k < 20; ++k) { tk *= -xx / (fp_t)(2 * k * (2 * k + 1)); new_sum = sum + tk; if (new_sum == sum) break; sum = new_sum; } #if defined(USE_HP_SIN) { hp_t hp_sum, hp_t0, hp_t1, hp_t2, hp_x, hp_xx; hp_x = (hp_t)x; hp_xx = hp_x * hp_x; hp_t0 = hp_x; hp_t1 = hp_t0 * (-hp_xx / HP(6.0)); hp_t2 = hp_t1 * (-hp_xx / HP(20.0)); hp_sum = (hp_t)sum; hp_sum += hp_t2; hp_sum += hp_t1; hp_sum += hp_t0; result = (fp_t)hp_sum; } #else sum += t2; sum += t1; sum += t0; result = sum; #endif /* defined(USE_HP_SIN) */ __K_SIN__ = k; } return (result); } fp_t wk_sinh(fp_t x) { fp_t e, result, s, w; w = FABS(x); s = COPYSIGN(ONE, x); if (ISNAN(x)) result = range_error("0xdead0019"); else if (w < __TINY__) result = x + x * x * x / SIX; else { e = EXPM1(w); result = s * (e + e / (ONE + e)) * HALF; } return (result); } fp_t wk_sqrt(fp_t x) { fp_t result; #if P == 113 result = sqrtl(x); #elif P == 64 result = sqrtl(x); #elif P == 53 result = sqrt(x); #elif P == 24 result = sqrtf(x); #else result = QNAN("0xdead0020"); #endif return (result); } fp_t wk_tan(fp_t x) { fp_t result, s, u, w; w = FABS(x); s = COPYSIGN(ONE, x); if (ISNAN(x) || (w > PI_HALF)) result = range_error("0xdead0021"); else if (w < __TINY__) result = x + x * x * x / THREE; else if (w < PI_EIGHTH) result = HALF * T(TWO * x); else if (FP(0.375) * PI < w) result = TWO * s / T(TWO * PI_HALF_minus_w(w)); else { /* u = T(2*w - PI_HALF) */ u = T(-PI_HALF_minus_w(TWO * w)); result = s * (TWO + u) / (TWO - u); } return (result); } fp_t wk_tan_series(fp_t x) { #if P == 113 return (wk_tan_series_128(x)); #elif P == 64 return (wk_tan_series_80(x)); #elif P == 53 return (wk_tan_series_64(x)); #elif P == 24 return (wk_tan_series_32(x)); #else return (QNAN("0xdead0022")); #endif } float wk_tan_series_32(float x) { /* ** > t14 := convert(taylor(tan(x), x = 0, 14), polynom); ** 3 5 17 7 62 9 1382 11 21844 13 ** x + 1/3 x + 2/15 x + --- x + ---- x + ------ x + ------- x ** 315 2835 155925 6081075 ** ** > Digits := 12: # because matula(24) == 9 ** ** > lprint(evalf(convert((t14 - x - x**3/3) / ((2/15) * x**5), ** horner))); ** [then prettyprint output to get s = ... below] ** ** > u32:= 2**(-23): ** ** > for k from 3 to 8 do ** printf("%3d\t%8.3f\n", k, tk(k, Pi/8) / tk(1,Pi/8) / u32) ** end do: ** 3 26599.176 ** 4 -1660.304 ** 5 103.755 ** 6 -6.485 ** 7 0.405 ** 8 -0.025 */ float result, s, sum, t1, t2, t3_7, xx, xxx, xxx_12; if (ISNAN(x) || (x > (float)PI_EIGHTH)) result = range_error("0xdead0023"); else { xx = x * x; s = FP32(1.0) + FP32(7.5) * (FP32(0.539682539683e-1) + (FP32(0.218694885362e-1) + (FP32(0.886323552990e-2) + FP32(0.359212803657e-2) * xx) * xx) * xx) * xx; xxx = xx * x; t3_7 = s * (FP32(2.0) * xx * xxx / FP32(15.0)); /* in t2, replace factor 1/3 by 1/4 + 1/12 to reduce sum error */ /* t2 = xxx / 3 */ xxx_12 = xxx / FP32(12.0); t2 = FP32(0.25) * xxx; t1 = x; sum = t3_7; sum += xxx_12; sum += t2; sum += t1; result = sum; } return (result); } double wk_tan_series_64(double x) { /* ** > t30 := convert(taylor(tan(x), x = 0, 30), polynom): ** ** > Digits := 21: # because matula(53) == 17 ** ** > lprint(evalf(convert((t30 - x - x**3/3) / ((2/15) * x**5), ** horner))); ** [then prettyprint output to get s = ... below] ** ** u64:= 2**(-53): ** ** Digits := 21: ** ** for k from 10 to 16 do ** printf("%3d\t%8.3f\n", k, tk(k, Pi/8)/tk(1,Pi/8)/u64) ** end do: ** 10 -106242.961 ** 11 6640.185 ** 12 -415.012 ** 13 25.938 ** 14 -1.621 ** 15 0.101 ** 16 -0.006 */ double result, s, sum, t1, t2, t3_15, xx, xxx, xxx_12; if (ISNAN(x) || (x > (double)PI_EIGHTH)) result = range_error("0xdead0024"); else { xx = x * x; s = FP64(1.0) + FP64(7.5) * (FP64(0.539682539682539682540e-1) + (FP64(0.218694885361552028219e-1) + (FP64(0.886323552990219656886e-2) + (FP64(0.359212803657248101693e-2) + (FP64(0.145583438705131826825e-2) + (FP64(0.590027440945585981378e-3) + (FP64(0.239129114243552481486e-3) + (FP64(0.969153795692945032560e-4) + (FP64(0.392783238833168340534e-4) + (FP64(0.159189050693289647407e-4) + (FP64(0.645168921565543076319e-5) + FP64(0.261477115129075455426e-5) * xx) * xx) * xx) * xx) * xx) * xx) * xx) * xx) * xx) * xx) * xx) * xx; xxx = xx * x; t3_15 = s * (TWO * xx * xxx / FP64(15.0)); /* in t2, replace factor 1/3 by 1/4 + 1/12 to reduce sum error */ /* t2 = xxx / 3 */ xxx_12 = xxx / FP64(12.0); t2 = FP64(0.25) * xxx; t1 = x; sum = t3_15; sum += xxx_12; sum += t2; sum += t1; result = sum; } return (result); } #if LDBL_MANT_DIG == 64 long double wk_tan_series_80(long double x) { /* ** > u80 := 2**(-63): ** ** > for k from 15 to 18 do printf("%3d\t%8.3f\n", k, ** tk(k, Pi/8)/tk(1,Pi/8)/u80) end do: ** 15 103.753 ** 16 -6.485 ** 17 0.405 ** 18 -0.025 ** ** > t34 := convert(taylor(tan(x), x = 0, 34), polynom): ** ** > Digits := 25: # because matula(64) == 21 ** ** > lprint(evalf(convert((t34 - x - x**3/3) / ((2/15) * x**5), ** horner))); ** [then prettyprint output to get s = ... below] */ long double result, s, sum, t1, t2, t3_17, xx, xxx, xxx_12; if (ISNAN(x) || (x > (long double)PI / FP80(8.0))) result = range_error("0xdead0025"); else { xx = x * x; s = FP80(1.0) + FP80(7.5) * (FP80(0.5396825396825396825396825e-1) + (FP80(0.2186948853615520282186949e-1) + (FP80(0.8863235529902196568863236e-2) + (FP80(0.3592128036572481016925461e-2) + (FP80(0.1455834387051318268249485e-2) + (FP80(0.5900274409455859813780760e-3) + (FP80(0.2391291142435524814857315e-3) + (FP80(0.9691537956929450325595875e-4) + (FP80(0.3927832388331683405337081e-4) + (FP80(0.1591890506932896474074428e-4) + (FP80(0.6451689215655430763190842e-5) + (FP80(0.2614771151290754554263594e-5) + (FP80(0.1059726832010465435091355e-5) + FP80(0.4294911078273805854820351e-6) * xx) * xx) * xx) * xx) * xx) * xx) * xx) * xx) * xx) * xx) * xx) * xx) * xx) * xx; xxx = xx * x; t3_17 = s * (TWO * xx * xxx / FP80(15.0)); /* in t2, replace factor 1/3 by 1/4 + 1/12 to reduce sum error */ /* t2 = xxx / FP80(3.0) */ xxx_12 = xxx / FP80(12.0); t2 = FP80(0.25) * xxx; t1 = x; sum = t3_17; sum += xxx_12; sum += t2; sum += t1; result = sum; } return (result); } #endif /* LDBL_MANT_DIG == 64 */ #if LDBL_MANT_DIG == 113 long double wk_tan_series_128(long double x) { /* ** > taylor(tan(x), x = 0, 9); ** 3 5 17 7 9 ** x + 1/3 x + 2/15 x + --- x + O(x ) ** 315 ** ** > u128 := 2**(-112): ** ** > B := proc(n) bernoulli(n) end proc: ** ** > tk := proc(k,x) (-1)**(k + 1) * ((2**(2 * k) * ** (2**(2 * k) - 1)) * B(2*k) / ** factorial(2*k)) * x**(2*k - 1) ** end proc: ** ** > for k from 0 to 4 do printf("%2d\t%a\n", k, tk(k,x)) end do: ** 0 0 ** 1 x ** 2 1/3*x^3 ** 3 2/15*x^5 ** 4 17/315*x^7 ** ** > for k from 26 to 30 do ** printf("%3d\t%8.3f\n", k, tk(k, Pi/8) / tk(1,Pi/8) / u128) ** end do: ** 27 207.506 ** 28 12.969 ** 29 0.811 ** 30 0.051 ** ** > t62 := convert(taylor(tan(x), x = 0, 62), polynom): ** ** > Digits := 40: # because matula(113) == 36 ** ** > lprint(evalf(convert((t62 - x - x**3/3) / ((2/15) * x**5), ** horner))); ** [then prettyprint output to get s = ... below] */ long double result, s, sum, t1, t2, t3_rest, xx, xxx, xxx_12; if (ISNAN(x) || (x > (long double)PI / FP128(8.0))) result = range_error("0xdead0026"); else { xx = x * x; s = FP128(1.) + FP128(7.5) * (FP128(0.5396825396825396825396825396825396825397e-1) + (FP128(0.2186948853615520282186948853615520282187e-1) + (FP128(0.8863235529902196568863235529902196568863e-2) + (FP128(0.3592128036572481016925461369905814350259e-2) + (FP128(0.1455834387051318268249485180702111919043e-2) + (FP128(0.5900274409455859813780759937000210887544e-3) + (FP128(0.2391291142435524814857314588851128093206e-3) + (FP128(0.9691537956929450325595875000388809377859e-4) + (FP128(0.3927832388331683405337080809312333509521e-4) + (FP128(0.1591890506932896474074427981657004305531e-4) + (FP128(0.6451689215655430763190842315302558044632e-5) + (FP128(0.2614771151290754554263594256409741097370e-5) + (FP128(0.1059726832010465435091355394124777619004e-5) + (FP128(0.4294911078273805854820351280397272242615e-6) + (FP128(0.1740661896357164777986229231330272921937e-6) + (FP128(0.7054636946400968325214518544480556569837e-7) + (FP128(0.2859136662305253908331094534141646943575e-7) + (FP128(0.1158764443279885220019793202909179039631e-7) + (FP128(0.4696295398230901628819043472467672382366e-8) + (FP128(0.1903336833931275921319706595462384034737e-8) + (FP128(0.7713933635359062290478933072162728548990e-9) + (FP128(0.3126339545892087045778120013440899955440e-9) + (FP128(0.1267057693030540106009628808122343862489e-9) + (FP128(0.5135191408039367740076049513054295176989e-10) + (FP128(0.2081214686770047419897227195060971086049e-10) + (FP128(0.8434845419094338293849148188558706453424e-11) + (FP128(0.3418514086811155814038680292485683222778e-11) + FP128(0.1385471574294846899390054465352388234911e-11) * xx) * xx) * xx) * xx) * xx) * xx) * xx) * xx) * xx) * xx) * xx) * xx) * xx) * xx) * xx) * xx) * xx) * xx) * xx) * xx) * xx) * xx) * xx ) * xx) * xx) * xx) * xx) * xx; xxx = xx * x; t3_rest = s * (TWO * xx * xxx / FP128(15.0)); /* in t2, replace factor 1/3 by 1/4 + 1/12 to reduce sum error */ /* t2 = xxx / 3 */ xxx_12 = xxx / FP128(12.0); t2 = FP128(0.25) * xxx; t1 = x; sum = t3_rest; sum += xxx_12; sum += t2; sum += t1; result = sum; } return (result); } #endif /* LDBL_MANT_DIG == 113 */ fp_t wk_tanh(fp_t x) { fp_t e, result, s, w; w = FABS(x); s = COPYSIGN(ONE, x); e = EXPM1(-TWO * w); if (ISNAN(x)) result = range_error("0xdead0027"); else if (x == ZERO) result = x; else if (w < __TINY__) result = x - x * x * x / THREE; else result = -s * e / (TWO + e); return (result); } /**********************************************************************/ void report(const char * name, const char * suffix, fp_t lo, fp_t hi, int n, fp_t xmax, fp_t amax, fp_t rmax, fp_t rsum) { PRINTF("------------------------------------" "------------------------------------\n"); PRINTF("%s() %s on [%.7Lg, %.7Lg] with %d samples\n", name, suffix, (long double)lo, (long double)hi, n); if (n > 0) { PRINTF("xmax = " FMT_L_g " = " FMT_L_a "\n", (long double)xmax, (long double)xmax); PRINTF("amax = " FMT_L_g " = " FMT_L_a "\n", (long double)amax, (long double)amax); PRINTF("rmax = " FMT_L_ulps " ulps\n", (long double)(rmax / oneulp())); PRINTF("ravg = " FMT_L_ulps " ulps\n", (long double)((rsum / (fp_t)n) / oneulp())); } } void save(fp_t x, fp_t r, fp_t a, hp_t e) { if (fpdata != (FILE *)NULL) (void)fprintf(fpdata, FMT_dat "\t" FMT_dat "\t" FMT_dat "\t" FMT_L_g "\n", (double)x, (double)r / (double)oneulp(), (double)a, (long double)e); } void show_options(void) { PRINTF("Compiler options:\n"); #if defined(USE_DECIMAL_FLOAT) PRINTF("\tUSE_DECIMAL_FLOAT\n"); #endif #if defined(USE_DECIMAL_DOUBLE) PRINTF("\tUSE_DECIMAL_DOUBLE\n"); #endif #if defined(USE_DECIMAL_LONG_DOUBLE) PRINTF("\tUSE_DECIMAL_LONG_DOUBLE\n"); #endif #if defined(USE_DOUBLE) PRINTF("\tUSE_DOUBLE\n"); #endif #if defined(USE_FLOAT) PRINTF("\tUSE_FLOAT\n"); #endif #if defined(USE_FMA) PRINTF("\tUSE_FMA\n"); #endif #if defined(USE_HP_COS) PRINTF("\tUSE_HP_COS\n"); #endif #if defined(USE_HP_SIN) PRINTF("\tUSE_HP_SIN\n"); #endif #if defined(USE_LONG_DOUBLE) PRINTF("\tUSE_LONG_DOUBLE\n"); #endif #if defined(USE_PI_HALF_SUM) PRINTF("\tUSE_PI_HALF_SUM\n"); #endif #if defined(USE_SIN_T_FORM) PRINTF("\tUSE_SIN_T_FORM\n"); #endif #if defined(USE_TWOSUM) PRINTF("\tUSE_TWOSUM\n"); #endif #if defined(USE_XP_TAN) PRINTF("\tUSE_XP_TAN\n"); #endif PRINTF("\n"); } void test_PI_HALF_minus_w(void) { /* ** NB: w = fabs(x) is always positive, so negative args are never required, ** However, we test arguments of both signs. */ begin("PI_HALF_minus_w"); test_wk_any("PI_HALF_minus_w", PI_HALF_minus_w, hp_PI_HALF_minus_w_alt, ZERO, PI_HALF, N); test_wk_any_rand("PI_HALF_minus_w", PI_HALF_minus_w, hp_PI_HALF_minus_w_alt, MINNORMAL, PI_HALF, N); test_wk_any("PI_HALF_minus_w", PI_HALF_minus_w, hp_PI_HALF_minus_w_alt, -PI_HALF, -ZERO, N); test_wk_any_rand("PI_HALF_minus_w", PI_HALF_minus_w, hp_PI_HALF_minus_w_alt, -PI_HALF, -MINNORMAL, N); end(); } void test_PI_HALF_minus_w_alt(void) { /* ** NB: w = fabs(x) is always positive, so negative args are never required, ** However, we test arguments of both signs. */ begin("PI_HALF_minus_w_alt"); test_wk_any("PI_HALF_minus_w_alt", PI_HALF_minus_w_alt, hp_PI_HALF_minus_w_alt, ZERO, PI_HALF, N); test_wk_any_rand("PI_HALF_minus_w_alt", PI_HALF_minus_w_alt, hp_PI_HALF_minus_w_alt, MINNORMAL, PI_HALF, N); test_wk_any("PI_HALF_minus_w_alt", PI_HALF_minus_w_alt, hp_PI_HALF_minus_w_alt, -PI_HALF, -ZERO, N); test_wk_any_rand("PI_HALF_minus_w_alt", PI_HALF_minus_w_alt, hp_PI_HALF_minus_w_alt, -PI_HALF, -MINNORMAL, N); end(); } void test_signed_zero(void) { PRINTF("Test of signed zero handling\n\n"); test_zero("wk_acos", wk_acos, 1); test_zero("wk_asin", wk_asin, 0); test_zero("wk_atan", wk_atan, 0); PRINTF("\n"); test_zero("wk_acos_series", wk_acos_series, 1); test_zero("wk_asin_series", wk_asin_series, 0); test_zero("wk_atan_series", wk_atan_series, 0); PRINTF("\n"); /* BB: acosh(0) is undefined: needs x >= 1 */ test_zero("wk_asinh", wk_asinh, 0); test_zero("wk_atanh", wk_atanh, 0); PRINTF("\n"); test_zero("wk_cos", wk_cos, 1); test_zero("wk_sin", wk_sin, 0); test_zero("wk_tan", wk_tan, 0); PRINTF("\n"); test_zero("wk_cos_series", wk_cos_series, 1); test_zero("wk_sin_series", wk_sin_series, 0); test_zero("wk_tan_series", wk_tan_series, 0); PRINTF("\n"); test_zero("wk_cosh", wk_cosh, 1); test_zero("wk_sinh", wk_sinh, 0); test_zero("wk_tanh", wk_tanh, 0); PRINTF("\n"); } void test_zero(const char * name, fp_t wk_any(fp_t x), int symmetric) { fp_t f_pos_zero, f_neg_zero, zero, neg_zero; int status; zero = ZERO; neg_zero = COPYSIGN(zero, -ONE); f_pos_zero = wk_any(zero); f_neg_zero = wk_any(neg_zero); if (symmetric) { status = (f_neg_zero == f_pos_zero) && (COPYSIGN(ONE, f_neg_zero) == ONE) && (COPYSIGN(ONE, f_pos_zero) == ONE); } else { status = (f_neg_zero == f_pos_zero) && (COPYSIGN(ONE, f_neg_zero) == -ONE) && (COPYSIGN(ONE, f_pos_zero) == ONE); } PRINTF("%-15s : %s\n", name, status ? "PASS" : "\tFAILED"); } void test_wk_any(const char * name, fp_t wk_any(fp_t x), hp_t any(hp_t x), fp_t lo, fp_t hi, int n) { fp_t amax, dx, rmax, rsum, x, xmax; int k, n_inf, n_nan, n_test; dx = (hi - lo) / (fp_t)n; if ((lo + dx) == lo) return; /* guard against dx too small */ amax = -MAXNORMAL; xmax = lo; rmax = -MAXNORMAL; rsum = ZERO; k = 0; n_inf = 0; n_nan = 0; n_test = 0; for (x = lo; x < hi; x += dx) { hp_t e; fp_t a, r; if (++k > n) break; /* guard against dx too small */ e = any(x); a = wk_any(x); if (ISNAN(a) || ISNAN(e)) { ++n_nan; continue; } else if (ISINF(a) || ISINF(e)) { ++n_inf; continue; } else ++n_test; r = (fp_t)((e == ZERO) ? a : ((hp_t)a - e) / e); r = FABS(r); rsum += r; if (r / oneulp() <= HALF) ++n_exact; if (r > rmax) { amax = a; rmax = r; xmax = x; save(x, r, a, e); } else save(x, r, a, e); } report(name, "[uniform]", lo, hi, n_test, xmax, amax, rmax, rsum); if (n_inf > 0) PRINTF("Skipped %d Inf results\n", n_inf); if (n_nan > 0) PRINTF("Skipped %d NaN results\n", n_nan); n_finite += n_test; (void)fflush((FILE *)NULL); } void test_wk_any_rand(const char * name, fp_t wk_any(fp_t x), hp_t any(hp_t x), fp_t lo, fp_t hi, int n) { fp_t amax, rmax, rsum, xmax; int k, n_inf, n_nan, n_test; amax = -MAXNORMAL; xmax = lo; rmax = -MAXNORMAL; rsum = ZERO; n_inf = 0; n_nan = 0; n_test = 0; for (k = 0; k < n; ++k) { hp_t e; fp_t a, r, x; x = randlab(lo, hi); e = any(x); a = wk_any(x); if (isnan(a) || isnan(e)) { ++n_nan; continue; } else if (ISINF(a) || ISINF(e)) { ++n_inf; continue; } else ++n_test; r = (fp_t)((e == ZERO) ? a : ((hp_t)a - e) / e); r = FABS(r); rsum += r; if (r / oneulp() <= HALF) ++n_exact; if (r > rmax) { amax = a; rmax = r; xmax = x; save(x, r, a, e); } else save(x, r, a, e); } report(name, "[log random]", lo, hi, n_test, xmax, amax, rmax, rsum); if (n_inf > 0) PRINTF("Skipped %d Inf results\n", n_inf); if (n_nan > 0) PRINTF("Skipped %d NaN results\n", n_nan); n_finite += n_test; (void)fflush((FILE *)NULL); } void test_wk_acos(void) { begin("wk_acos"); test_wk_any("wk_acos", wk_acos, ACOS, ZERO, ONE, N); test_wk_any("wk_acos", wk_acos, ACOS, -ONE, -ZERO, N); test_wk_any_rand("wk_acos", wk_acos, ACOS, MINSUBNORMAL, MINNORMAL, N); test_wk_any_rand("wk_acos", wk_acos, ACOS, -MINNORMAL, -MINSUBNORMAL, N); test_wk_any_rand("wk_acos", wk_acos, ACOS, MINNORMAL, ONE, N); test_wk_any_rand("wk_acos", wk_acos, ACOS, -ONE, -MINNORMAL, N); end(); } void test_wk_acos_series(void) { fp_t xmax; xmax = FP(31.0) / FP(32.0); /* see wk_acos_asin_series() comments */ begin("wk_acos_series"); test_wk_any("wk_acos_series", wk_acos_series, ACOS, ZERO, xmax, N); test_wk_any("wk_acos_series", wk_acos_series, ACOS, -xmax, -ZERO, N); test_wk_any_rand("wk_acos_series", wk_acos_series, ACOS, MINSUBNORMAL, MINNORMAL, N); test_wk_any_rand("wk_acos_series", wk_acos_series, ACOS, -MINNORMAL, -MINSUBNORMAL, N); test_wk_any_rand("wk_acos_series", wk_acos_series, ACOS, MINNORMAL, xmax, N); test_wk_any_rand("wk_acos_series", wk_acos_series, ACOS, -xmax, -MINNORMAL, N); end(); } void test_wk_acosh(void) { begin("wk_acosh"); test_wk_any("wk_acosh", wk_acosh, ACOSH, ONE, ONE + FP(0x1p-12), N); test_wk_any("wk_acosh", wk_acosh, ACOSH, ONE, ONE + FP(0x1p-6), N); test_wk_any("wk_acosh", wk_acosh, ACOSH, ONE, ONE + FP(0x1p-4), N); test_wk_any("wk_acosh", wk_acosh, ACOSH, ONE, ONE + FP(0x1p-2), N); test_wk_any("wk_acosh", wk_acosh, ACOSH, ONE, TWO, N); test_wk_any("wk_acosh", wk_acosh, ACOSH, TWO, FP(40.0), N); test_wk_any("wk_acosh", wk_acosh, ACOSH, FP(40.0), FP(74.0), N); test_wk_any("wk_acosh", wk_acosh, ACOSH, FP(74.0), FP(1098.0), N); test_wk_any_rand("wk_acosh", wk_acosh, ACOSH, ONE, SQRT(MAXNORMAL), N); test_wk_any_rand("wk_acosh", wk_acosh, ACOSH, SQRT(MAXNORMAL), HALF * MAXNORMAL, N); test_wk_any_rand("wk_acosh", wk_acosh, ACOSH, HALF * MAXNORMAL, MAXNORMAL, N); end(); } void test_wk_acosh_alt(void) { begin("wk_acosh_alt"); test_wk_any("wk_acosh_alt", wk_acosh_alt, ACOSH, ONE, ONE + FP(0x1p-12), N); test_wk_any("wk_acosh_alt", wk_acosh_alt, ACOSH, ONE, ONE + FP(0x1p-6), N); test_wk_any("wk_acosh_alt", wk_acosh_alt, ACOSH, ONE, ONE + FP(0x1p-4), N); test_wk_any("wk_acosh_alt", wk_acosh_alt, ACOSH, ONE, ONE + FP(0x1p-2), N); test_wk_any("wk_acosh_alt", wk_acosh_alt, ACOSH, ONE, TWO, N); test_wk_any("wk_acosh_alt", wk_acosh_alt, ACOSH, TWO, FP(40.0), N); test_wk_any("wk_acosh_alt", wk_acosh_alt, ACOSH, FP(40.0), FP(74.0), N); test_wk_any("wk_acosh_alt", wk_acosh_alt, ACOSH, FP(74.0), FP(1098.0), N); test_wk_any_rand("wk_acosh_alt", wk_acosh_alt, ACOSH, ONE, SQRT(MAXNORMAL), N); test_wk_any_rand("wk_acosh_alt", wk_acosh_alt, ACOSH, SQRT(MAXNORMAL), HALF * MAXNORMAL, N); test_wk_any_rand("wk_acosh_alt", wk_acosh_alt, ACOSH, HALF * MAXNORMAL, MAXNORMAL, N); end(); } void test_wk_asin(void) { begin("wk_asin"); test_wk_any("wk_asin", wk_asin, ASIN, ZERO, ONE, N); test_wk_any("wk_asin", wk_asin, ASIN, -ONE, -ZERO, N); test_wk_any_rand("wk_asin", wk_asin, ASIN, MINSUBNORMAL, MINNORMAL, N); test_wk_any_rand("wk_asin", wk_asin, ASIN, -MINNORMAL, -MINSUBNORMAL, N); test_wk_any_rand("wk_asin", wk_asin, ASIN, MINNORMAL, ONE, N); test_wk_any_rand("wk_asin", wk_asin, ASIN, -ONE, -MINNORMAL, N); end(); } void test_wk_asin_series(void) { fp_t xmax; xmax = FP(31.0) / FP(32.0); /* see wk_acos_asin_series() comments */ begin("wk_asin_series"); test_wk_any("wk_asin_series", wk_asin_series, ASIN, ZERO, xmax, N); test_wk_any("wk_asin_series", wk_asin_series, ASIN, -xmax, -ZERO, N); test_wk_any_rand("wk_asin_series", wk_asin_series, ASIN, MINSUBNORMAL, MINNORMAL, N); test_wk_any_rand("wk_asin_series", wk_asin_series, ASIN, -MINNORMAL, -MINSUBNORMAL, N); test_wk_any_rand("wk_asin_series", wk_asin_series, ASIN, MINNORMAL, xmax, N); test_wk_any_rand("wk_asin_series", wk_asin_series, ASIN, -xmax, -MINNORMAL, N); end(); } void test_wk_asinh(void) { begin("wk_asinh"); test_wk_any("wk_asinh", wk_asinh, ASINH, ZERO, ONE, N); test_wk_any("wk_asinh", wk_asinh, ASINH, -ONE, -ZERO, N); test_wk_any_rand("wk_asinh", wk_asinh, ASINH, MINSUBNORMAL, MINNORMAL, N); test_wk_any_rand("wk_asinh", wk_asinh, ASINH, -MINNORMAL, -MINSUBNORMAL, N); test_wk_any_rand("wk_asinh", wk_asinh, ASINH, MINNORMAL, ONE, N); test_wk_any_rand("wk_asinh", wk_asinh, ASINH, -ONE, -MINNORMAL, N); test_wk_any("wk_asinh", wk_asinh, ASINH, ONE, TWO, N); test_wk_any("wk_asinh", wk_asinh, ASINH, TWO, FP(10.0), N); test_wk_any("wk_asinh", wk_asinh, ASINH, FP(10.0), FP(74.0), N); test_wk_any("wk_asinh", wk_asinh, ASINH, FP(74.0), FP(1098.0), N); test_wk_any("wk_asinh", wk_asinh, ASINH, -TWO, -ONE, N); test_wk_any("wk_asinh", wk_asinh, ASINH, -FP(10.0), -TWO, N); test_wk_any("wk_asinh", wk_asinh, ASINH, -FP(74.0), -FP(10.0), N); test_wk_any("wk_asinh", wk_asinh, ASINH, -FP(1098.0), -FP(74.0), N); test_wk_any_rand("wk_asinh", wk_asinh, ASINH, ONE, FP(100.0), N); test_wk_any_rand("wk_asinh", wk_asinh, ASINH, -FP(100.0), -ONE, N); test_wk_any_rand("wk_asinh", wk_asinh, ASINH, ONE, MAXNORMAL, N); test_wk_any_rand("wk_asinh", wk_asinh, ASINH, -MAXNORMAL, -ONE, N); end(); } void test_wk_asinh_alt(void) { begin("wk_asinh_alt"); test_wk_any("wk_asinh_alt", wk_asinh_alt, ASINH, ZERO, ONE, N); test_wk_any("wk_asinh_alt", wk_asinh_alt, ASINH, -ONE, -ZERO, N); test_wk_any_rand("wk_asinh_alt", wk_asinh_alt, ASINH, MINSUBNORMAL, MINNORMAL, N); test_wk_any_rand("wk_asinh_alt", wk_asinh_alt, ASINH, -MINNORMAL, -MINSUBNORMAL, N); test_wk_any_rand("wk_asinh_alt", wk_asinh_alt, ASINH, MINNORMAL, ONE, N); test_wk_any_rand("wk_asinh_alt", wk_asinh_alt, ASINH, -ONE, -MINNORMAL, N); test_wk_any("wk_asinh_alt", wk_asinh_alt, ASINH, ONE, TWO, N); test_wk_any("wk_asinh_alt", wk_asinh_alt, ASINH, TWO, FP(10.0), N); test_wk_any("wk_asinh_alt", wk_asinh_alt, ASINH, FP(10.0), FP(74.0), N); test_wk_any("wk_asinh_alt", wk_asinh_alt, ASINH, FP(74.0), FP(1098.0), N); test_wk_any("wk_asinh_alt", wk_asinh_alt, ASINH, -TWO, -ONE, N); test_wk_any("wk_asinh_alt", wk_asinh_alt, ASINH, -FP(10.0), -TWO, N); test_wk_any("wk_asinh_alt", wk_asinh_alt, ASINH, -FP(74.0), -FP(10.0), N); test_wk_any("wk_asinh_alt", wk_asinh_alt, ASINH, -FP(1098.0), -FP(74.0), N); test_wk_any_rand("wk_asinh_alt", wk_asinh_alt, ASINH, ONE, FP(100.0), N); test_wk_any_rand("wk_asinh_alt", wk_asinh_alt, ASINH, -FP(100.0), -ONE, N); test_wk_any_rand("wk_asinh_alt", wk_asinh_alt, ASINH, ONE, MAXNORMAL, N); test_wk_any_rand("wk_asinh_alt", wk_asinh_alt, ASINH, -MAXNORMAL, -ONE, N); end(); } void test_wk_atan(void) { begin("wk_atan"); begin("wk_atan"); test_wk_any("wk_atan", wk_atan, ATAN, ZERO, ONE, N); test_wk_any("wk_atan", wk_atan, ATAN, -ONE, -ZERO, N); test_wk_any_rand("wk_atan", wk_atan, ATAN, MINSUBNORMAL, MINNORMAL, N); test_wk_any_rand("wk_atan", wk_atan, ATAN, -MINNORMAL, -MINSUBNORMAL, N); test_wk_any_rand("wk_atan", wk_atan, ATAN, MINNORMAL, ONE, N); test_wk_any_rand("wk_atan", wk_atan, ATAN, -ONE, -MINNORMAL, N); test_wk_any_rand("wk_atan", wk_atan, ATAN, ONE, FP(10.0), N); test_wk_any_rand("wk_atan", wk_atan, ATAN, -FP(10.0), -ONE, N); test_wk_any_rand("wk_atan", wk_atan, ATAN, FP(10.0), MAXNORMAL, N); test_wk_any_rand("wk_atan", wk_atan, ATAN, -MAXNORMAL, -FP(10.0), N); end(); } void test_wk_atan_series(void) { fp_t xmax; xmax = FOURTH + FP(0x1p-6); /* see wk_atan_series() comments */ begin("wk_atan_series"); test_wk_any("wk_atan_series", wk_atan_series, ATAN, FOURTH, xmax, N); test_wk_any_rand("wk_atan_series", wk_atan_series, ATAN, FOURTH, xmax, N); test_wk_any("wk_atan_series", wk_atan_series, ATAN, -xmax, -FOURTH, N); test_wk_any_rand("wk_atan_series", wk_atan_series, ATAN, -xmax, -FOURTH, N); xmax = SQRT(TWO) - ONE; test_wk_any("wk_atan_series", wk_atan_series, ATAN, ZERO, xmax, N); test_wk_any("wk_atan_series", wk_atan_series, ATAN, -xmax, -ZERO, N); test_wk_any_rand("wk_atan_series", wk_atan_series, ATAN, MINSUBNORMAL, MINNORMAL, N); test_wk_any_rand("wk_atan_series", wk_atan_series, ATAN, -MINNORMAL, -MINSUBNORMAL, N); test_wk_any_rand("wk_atan_series", wk_atan_series, ATAN, MINNORMAL, xmax, N); test_wk_any_rand("wk_atan_series", wk_atan_series, ATAN, -xmax, -MINNORMAL, N); end(); } void test_wk_atanh(void) { begin("wk_atanh"); test_wk_any("wk_atanh", wk_atanh, ATANH, ZERO, ONE, N); test_wk_any("wk_atanh", wk_atanh, ATANH, -ONE, -ZERO, N); test_wk_any_rand("wk_atanh", wk_atanh, ATANH, MINSUBNORMAL, MINNORMAL, N); test_wk_any_rand("wk_atanh", wk_atanh, ATANH, -MINNORMAL, -MINSUBNORMAL, N); test_wk_any_rand("wk_atanh", wk_atanh, ATANH, MINNORMAL, ONE, N); test_wk_any_rand("wk_atanh", wk_atanh, ATANH, -ONE, -MINNORMAL, N); end(); } void test_wk_cos(void) { begin("wk_cos"); test_wk_any("wk_cos", wk_cos, COS, ZERO, ONE, N); test_wk_any("wk_cos", wk_cos, COS, -ONE, -ZERO, N); test_wk_any("wk_cos", wk_cos, COS, ONE, PI_HALF, N); test_wk_any("wk_cos", wk_cos, COS, -PI_HALF, -ONE, N); test_wk_any("wk_cos", wk_cos, COS, MINSUBNORMAL, MINNORMAL, N); test_wk_any("wk_cos", wk_cos, COS, -MINNORMAL, -MINSUBNORMAL, N); test_wk_any_rand("wk_cos", wk_cos, COS, MINSUBNORMAL, MINNORMAL, N); test_wk_any_rand("wk_cos", wk_cos, COS, -MINNORMAL, -MINSUBNORMAL, N); test_wk_any_rand("wk_cos", wk_cos, COS, MINNORMAL, ONE, N); test_wk_any_rand("wk_cos", wk_cos, COS, -ONE, -MINNORMAL, N); test_wk_any_rand("wk_cos", wk_cos, COS, ONE, PI_HALF, N); test_wk_any_rand("wk_cos", wk_cos, COS, -PI_HALF, -ONE, N); end(); } void test_wk_cos_series(void) { begin("wk_cos_series"); test_wk_any("wk_cos_series", wk_cos_series, COS, ZERO, PI_FOURTH, N); test_wk_any("wk_cos_series", wk_cos_series, COS, -PI_FOURTH, -ZERO, N); test_wk_any("wk_cos_series", wk_cos_series, COS, MINSUBNORMAL, MINNORMAL, N); test_wk_any("wk_cos_series", wk_cos_series, COS, -MINNORMAL, -MINSUBNORMAL, N); test_wk_any_rand("wk_cos_series", wk_cos_series, COS, MINSUBNORMAL, MINNORMAL, N); test_wk_any_rand("wk_cos_series", wk_cos_series, COS, -MINNORMAL, -MINSUBNORMAL, N); test_wk_any_rand("wk_cos_series", wk_cos_series, COS, MINNORMAL, PI_FOURTH, N); test_wk_any_rand("wk_cos_series", wk_cos_series, COS, -PI_FOURTH, -MINNORMAL, N); end(); } void test_wk_cosh(void) { fp_t xmax; xmax = LOG(MAXNORMAL) - FP(5.0); begin("wk_cosh"); test_wk_any("wk_cosh", wk_cosh, COSH, ZERO, ONE, N); test_wk_any("wk_cosh", wk_cosh, COSH, -ONE, -ZERO, N); test_wk_any("wk_cosh", wk_cosh, COSH, ONE, xmax, N); test_wk_any("wk_cosh", wk_cosh, COSH, -xmax, -ONE, N); test_wk_any_rand("wk_cosh", wk_cosh, COSH, MINSUBNORMAL, MINNORMAL, N); test_wk_any_rand("wk_cosh", wk_cosh, COSH, -MINNORMAL, -MINSUBNORMAL, N); test_wk_any_rand("wk_cosh", wk_cosh, COSH, MINNORMAL, ONE, N); test_wk_any_rand("wk_cosh", wk_cosh, COSH, -ONE, -MINNORMAL, N); test_wk_any_rand("wk_cosh", wk_cosh, COSH, ONE, xmax, N); test_wk_any_rand("wk_cosh", wk_cosh, COSH, -xmax, -ONE, N); test_wk_any("wk_cosh", wk_cosh, COSH, -MINNORMAL, -MINSUBNORMAL, N); test_wk_any_rand("wk_cosh", wk_cosh, COSH, -MINNORMAL, -MINSUBNORMAL, N); test_wk_any("wk_cosh", wk_cosh, COSH, MINSUBNORMAL, MINNORMAL, N); test_wk_any_rand("wk_cosh", wk_cosh, COSH, MINSUBNORMAL, MINNORMAL, N); end(); } void test_wk_sin(void) { begin("wk_sin"); test_wk_any("wk_sin", wk_sin, SIN, ZERO, ONE, N); test_wk_any("wk_sin", wk_sin, SIN, -ONE, -ZERO, N); test_wk_any("wk_sin", wk_sin, SIN, ONE, PI_HALF, N); test_wk_any("wk_sin", wk_sin, SIN, -PI_HALF, -ONE, N); test_wk_any("wk_sin", wk_sin, SIN, MINSUBNORMAL, MINNORMAL, N); test_wk_any("wk_sin", wk_sin, SIN, -MINNORMAL, -MINSUBNORMAL, N); test_wk_any_rand("wk_sin", wk_sin, SIN, MINSUBNORMAL, MINNORMAL, N); test_wk_any_rand("wk_sin", wk_sin, SIN, -MINNORMAL, -MINSUBNORMAL, N); test_wk_any_rand("wk_sin", wk_sin, SIN, MINNORMAL, ONE, N); test_wk_any_rand("wk_sin", wk_sin, SIN, -ONE, -MINNORMAL, N); test_wk_any_rand("wk_sin", wk_sin, SIN, ONE, PI_HALF, N); test_wk_any_rand("wk_sin", wk_sin, SIN, -PI_HALF, -ONE, N); end(); } void test_wk_sin_series(void) { /* NB: series limited to x in [-PI/4, +PI/4] */ begin("wk_sin_series"); test_wk_any("wk_sin_series", wk_sin_series, SIN, ZERO, PI_EIGHTH, N); test_wk_any("wk_sin_series", wk_sin_series, SIN, -PI_EIGHTH, -ZERO, N); test_wk_any("wk_sin_series", wk_sin_series, SIN, PI_EIGHTH, PI_FOURTH, N); test_wk_any("wk_sin_series", wk_sin_series, SIN, -PI_FOURTH, -PI_EIGHTH, N); test_wk_any_rand("wk_sin_series", wk_sin_series, SIN, MINSUBNORMAL, MINNORMAL, N); test_wk_any_rand("wk_sin_series", wk_sin_series, SIN, -MINNORMAL, -MINSUBNORMAL, N); test_wk_any_rand("wk_sin_series", wk_sin_series, SIN, MINNORMAL, PI_FOURTH, N); test_wk_any_rand("wk_sin_series", wk_sin_series, SIN, -PI_FOURTH, -MINNORMAL, N); end(); } void test_wk_sinh(void) { fp_t xmax; xmax = LOG(MAXNORMAL) - FP(5.0); begin("wk_sinh"); test_wk_any("wk_sinh", wk_sinh, SINH, ZERO, ONE, N); test_wk_any("wk_sinh", wk_sinh, SINH, -ONE, -ZERO, N); test_wk_any("wk_sinh", wk_sinh, SINH, ONE, xmax, N); test_wk_any("wk_sinh", wk_sinh, SINH, -xmax, -ONE, N); test_wk_any_rand("wk_sinh", wk_sinh, SINH, MINSUBNORMAL, MINNORMAL, N); test_wk_any_rand("wk_sinh", wk_sinh, SINH, -MINNORMAL, -MINSUBNORMAL, N); test_wk_any_rand("wk_sinh", wk_sinh, SINH, MINNORMAL, ONE, N); test_wk_any_rand("wk_sinh", wk_sinh, SINH, -ONE, -MINNORMAL, N); test_wk_any_rand("wk_sinh", wk_sinh, SINH, ONE, xmax, N); test_wk_any_rand("wk_sinh", wk_sinh, SINH, -xmax, -ONE, N); test_wk_any("wk_sinh", wk_sinh, SINH, -MINNORMAL, -MINSUBNORMAL, N); test_wk_any_rand("wk_sinh", wk_sinh, SINH, -MINNORMAL, -MINSUBNORMAL, N); test_wk_any("wk_sinh", wk_sinh, SINH, MINSUBNORMAL, MINNORMAL, N); test_wk_any_rand("wk_sinh", wk_sinh, SINH, MINSUBNORMAL, MINNORMAL, N); end(); } void test_wk_tan(void) { begin("wk_tan"); test_wk_any("wk_tan", wk_tan, TAN, MINSUBNORMAL, MINNORMAL, N); test_wk_any("wk_tan", wk_tan, TAN, -MINNORMAL, -MINSUBNORMAL, N); test_wk_any("wk_tan", wk_tan, TAN, ZERO, ONE, N); test_wk_any("wk_tan", wk_tan, TAN, -ONE, -ZERO, N); test_wk_any("wk_tan", wk_tan, TAN, ONE, PI_HALF, N); test_wk_any("wk_tan", wk_tan, TAN, -PI_HALF, -ONE, N); test_wk_any_rand("wk_tan", wk_tan, TAN, MINSUBNORMAL, MINNORMAL, N); test_wk_any_rand("wk_tan", wk_tan, TAN, -MINNORMAL, -MINSUBNORMAL, N); test_wk_any_rand("wk_tan", wk_tan, TAN, MINNORMAL, ONE, N); test_wk_any_rand("wk_tan", wk_tan, TAN, -ONE, -MINNORMAL, N); test_wk_any_rand("wk_tan", wk_tan, TAN, ONE, PI_HALF, N); test_wk_any_rand("wk_tan", wk_tan, TAN, -PI_HALF, -ONE, N); end(); } void test_wk_tan_series(void) { begin("wk_tan_series"); test_wk_any("wk_tan_series", wk_tan_series, TAN, MINSUBNORMAL, MINNORMAL, N); test_wk_any("wk_tan_series", wk_tan_series, TAN, -MINNORMAL, -MINSUBNORMAL, N); test_wk_any("wk_tan_series", wk_tan_series, TAN, ZERO, PI_EIGHTH, N); test_wk_any("wk_tan_series", wk_tan_series, TAN, -PI_EIGHTH, -ZERO, N); test_wk_any_rand("wk_tan_series", wk_tan_series, TAN, MINSUBNORMAL, MINNORMAL, N); test_wk_any_rand("wk_tan_series", wk_tan_series, TAN, -MINNORMAL, -MINSUBNORMAL, N); test_wk_any_rand("wk_tan_series", wk_tan_series, TAN, MINNORMAL, PI_EIGHTH, N); test_wk_any_rand("wk_tan_series", wk_tan_series, TAN, -PI_EIGHTH, -MINNORMAL, N); end(); } void test_wk_tan_series_32(void) { #if P == 24 begin("wk_tan_series_32"); test_wk_any("wk_tan_series_32", wk_tan_series_32, TAN, MINSUBNORMAL, MINNORMAL, N); test_wk_any("wk_tan_series_32", wk_tan_series_32, TAN, -MINNORMAL, -MINSUBNORMAL, N); test_wk_any("wk_tan_series_32", wk_tan_series_32, TAN, ZERO, PI_EIGHTH, N); test_wk_any("wk_tan_series_32", wk_tan_series_32, TAN, -PI_EIGHTH, -ZERO, N); test_wk_any_rand("wk_tan_series_32", wk_tan_series_32, TAN, MINSUBNORMAL, MINNORMAL, N); test_wk_any_rand("wk_tan_series_32", wk_tan_series_32, TAN, -MINNORMAL, -MINSUBNORMAL, N); test_wk_any_rand("wk_tan_series_32", wk_tan_series_32, TAN, MINNORMAL, PI_EIGHTH, N); test_wk_any_rand("wk_tan_series_32", wk_tan_series_32, TAN, -PI_EIGHTH, -MINNORMAL, N); end(); #endif /* P == 24 */ } void test_wk_tan_series_64(void) { #if P == 53 begin("wk_tan_series_64"); test_wk_any("wk_tan_series_64", wk_tan_series_64, TAN, MINSUBNORMAL, MINNORMAL, N); test_wk_any("wk_tan_series_64", wk_tan_series_64, TAN, -MINNORMAL, -MINSUBNORMAL, N); test_wk_any("wk_tan_series_64", wk_tan_series_64, TAN, ZERO, PI_EIGHTH, N); test_wk_any("wk_tan_series_64", wk_tan_series_64, TAN, -PI_EIGHTH, -ZERO, N); test_wk_any_rand("wk_tan_series_64", wk_tan_series_64, TAN, MINSUBNORMAL, MINNORMAL, N); test_wk_any_rand("wk_tan_series_64", wk_tan_series_64, TAN, -MINNORMAL, -MINSUBNORMAL, N); test_wk_any_rand("wk_tan_series_64", wk_tan_series_64, TAN, MINNORMAL, PI_EIGHTH, N); test_wk_any_rand("wk_tan_series_64", wk_tan_series_64, TAN, -PI_EIGHTH, -MINNORMAL, N); end(); #endif /* P == 53 */ } void test_wk_tan_series_80(void) { #if P == 64 begin("wk_tan_series_80"); test_wk_any("wk_tan_series_80", wk_tan_series_80, TAN, MINSUBNORMAL, MINNORMAL, N); test_wk_any("wk_tan_series_80", wk_tan_series_80, TAN, -MINNORMAL, -MINSUBNORMAL, N); test_wk_any("wk_tan_series_80", wk_tan_series_80, TAN, ZERO, PI_EIGHTH, N); test_wk_any("wk_tan_series_80", wk_tan_series_80, TAN, -PI_EIGHTH, -ZERO, N); test_wk_any_rand("wk_tan_series_80", wk_tan_series_80, TAN, MINSUBNORMAL, MINNORMAL, N); test_wk_any_rand("wk_tan_series_80", wk_tan_series_80, TAN, -MINNORMAL, -MINSUBNORMAL, N); test_wk_any_rand("wk_tan_series_80", wk_tan_series_80, TAN, MINNORMAL, PI_EIGHTH, N); test_wk_any_rand("wk_tan_series_80", wk_tan_series_80, TAN, -PI_EIGHTH, -MINNORMAL, N); end(); #endif /* P == 64 */ } void test_wk_tan_series_128(void) { #if P == 113 begin("tan_series_128"); test_wk_any("wk_tan_series_128", wk_tan_series_128, TAN, MINSUBNORMAL, MINNORMAL, N); test_wk_any("wk_tan_series_128", wk_tan_series_128, TAN, -MINNORMAL, -MINSUBNORMAL, N); test_wk_any("wk_tan_series_128", wk_tan_series_128, TAN, ZERO, PI_EIGHTH, N); test_wk_any("wk_tan_series_128", wk_tan_series_128, TAN, -PI_EIGHTH, -ZERO, N); test_wk_any_rand("wk_tan_series_128", wk_tan_series_128, TAN, MINSUBNORMAL, MINNORMAL, N); test_wk_any_rand("wk_tan_series_128", wk_tan_series_128, TAN, -MINNORMAL, -MINSUBNORMAL, N); test_wk_any_rand("wk_tan_series_128", wk_tan_series_128, TAN, MINNORMAL, PI_EIGHTH, N); test_wk_any_rand("wk_tan_series_128", wk_tan_series_128, TAN, -PI_EIGHTH, -MINNORMAL, N); end(); #endif /* P == 113 */ } void test_wk_tanh(void) { fp_t xmax; xmax = LOG(MAXNORMAL) - FP(5.0); begin("wk_tanh"); test_wk_any("wk_tanh", wk_tanh, TANH, ZERO, ONE, N); test_wk_any("wk_tanh", wk_tanh, TANH, -ONE, -ZERO, N); test_wk_any("wk_tanh", wk_tanh, TANH, ONE, FP(30.0), N); test_wk_any("wk_tanh", wk_tanh, TANH, -FP(30.0), -ONE, N); test_wk_any("wk_tanh", wk_tanh, TANH, ONE, xmax, N); test_wk_any("wk_tanh", wk_tanh, TANH, -xmax, -ONE, N); test_wk_any_rand("wk_tanh", wk_tanh, TANH, MINSUBNORMAL, MINNORMAL, N); test_wk_any_rand("wk_tanh", wk_tanh, TANH, -MINNORMAL, -MINSUBNORMAL, N); test_wk_any_rand("wk_tanh", wk_tanh, TANH, MINNORMAL, ONE, N); test_wk_any_rand("wk_tanh", wk_tanh, TANH, -ONE, -MINNORMAL, N); test_wk_any_rand("wk_tanh", wk_tanh, TANH, ONE, xmax, N); test_wk_any_rand("wk_tanh", wk_tanh, TANH, -xmax, -ONE, N); test_wk_any("wk_tanh", wk_tanh, TANH, -MINNORMAL, -MINSUBNORMAL, N); test_wk_any_rand("wk_tanh", wk_tanh, TANH, -MINNORMAL, -MINSUBNORMAL, N); test_wk_any("wk_tanh", wk_tanh, TANH, MINSUBNORMAL, MINNORMAL, N); test_wk_any_rand("wk_tanh", wk_tanh, TANH, MINSUBNORMAL, MINNORMAL, N); end(); } int main(int argc, char* argv[]) { /* usage: a.out [alternate N value] [plotdataflag] */ int status; static const long int INITSEED = 0x7eadbeefcafefeedL; save_data = 0; /* normally 0 */ if (argc > 1) N = atoi(argv[1]); if (argc > 2) save_data = atoi(argv[2]); if (argc > 3) suffix = argv[3]; rounding_mode = (argc > 4) ? argv[4] : "FE_TONEAREST"; if (strcmp(rounding_mode, "FE_DOWNWARD") == 0) status = fesetround(FE_DOWNWARD); else if (strcmp(rounding_mode, "FE_TONEAREST") == 0) status = fesetround(FE_TONEAREST); else if (strcmp(rounding_mode, "FE_TOWARDZERO") == 0) status = fesetround(FE_TOWARDZERO); else if (strcmp(rounding_mode, "FE_UPWARD") == 0) status = fesetround(FE_UPWARD); else status = fesetround(FE_TONEAREST); srand48(INITSEED); PRINTF("Rounding mode = %s%s\n", (status == 0) ? rounding_mode : "FE_TONEAREST", (status == 0) ? " [selected]" : " [default]"); PRINTF("N = %d\n", (int)N); PRINTF("P = %d\n", (int)P); PRINTF("sizeof(fp_t) = %d\n", (int)sizeof(fp_t)); PRINTF("EPS = %.36Lg = %La\n", (long double)EPS, (long double)EPS); PRINTF("__TINY__ = %.36Lg = %La\n", (long double)__TINY__, (long double)__TINY__); PRINTF("save_data = %d\n", save_data); PRINTF("suffix = \"%s\"\n", suffix == (const char *)NULL ? "" : suffix); PRINTF("\n"); show_options(); test_signed_zero(); #if 0 test_wk_acos_series(); test_wk_asin_series(); test_wk_atan_series(); test_wk_cos_series(); test_wk_sin_series(); test_wk_tan_series(); return (EXIT_SUCCESS); #endif test_PI_HALF_minus_w(); test_PI_HALF_minus_w_alt(); test_wk_acos(); test_wk_acos_series(); test_wk_asin(); test_wk_asin_series(); test_wk_atan(); test_wk_atan_series(); test_wk_acosh(); test_wk_acosh_alt(); test_wk_asinh(); test_wk_asinh_alt(); test_wk_atanh(); test_wk_cos(); test_wk_cos_series(); test_wk_sin(); test_wk_sin_series(); test_wk_tan(); test_wk_tan_series(); test_wk_tan_series_32(); test_wk_tan_series_64(); test_wk_tan_series_80(); test_wk_tan_series_128(); test_wk_cosh(); test_wk_sinh(); test_wk_tanh(); end(); return (EXIT_SUCCESS); }