From 9ecf223fc1f94bb51a296fa0e6a6c632742ceb52 Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Thu, 4 Dec 2014 23:32:39 +0530 Subject: [PATCH] Get the ld80 routines from OpenBSD to build on mac and linux. Bump version number and SO major version, since we have introduced new long double APIs. --- Make.inc | 6 +- ld80/Make.files | 12 +- ld80/e_acoshl.c | 2 +- ld80/e_atanhl.c | 2 +- ld80/e_coshl.c | 2 +- ld80/e_expl.c | 2 +- ld80/e_hypotl.c | 2 +- ld80/e_lgammal.c | 3 +- ld80/e_log10l.c | 2 +- ld80/e_log2l.c | 2 +- ld80/e_logl.c | 2 +- ld80/e_powl.c | 2 +- ld80/e_rem_pio2l.h | 20 ++-- ld80/e_sinhl.c | 2 +- ld80/e_tgammal.c | 3 +- ld80/invtrig.h | 2 +- ld80/s_asinhl.c | 2 +- ld80/s_ceill.c | 2 +- ld80/s_erfl.c | 2 +- ld80/s_exp2l.c | 2 +- ld80/s_expm1l.c | 2 +- ld80/s_floorl.c | 2 +- ld80/s_log1pl.c | 2 +- ld80/s_modfl.c | 2 +- ld80/s_nanl.c | 2 +- ld80/s_nextafterl.c | 4 +- ld80/s_nexttoward.c | 2 +- ld80/s_nexttowardf.c | 2 +- ld80/s_remquol.c | 2 +- ld80/s_tanhl.c | 2 +- ld80/s_truncl.c | 4 +- src/Make.files | 2 +- src/math_private.h | 3 +- src/math_private_openbsd.h | 218 +++++++++++++++++++++++++++++++++++++ src/polevll.c | 104 ++++++++++++++++++ 35 files changed, 377 insertions(+), 50 deletions(-) create mode 100644 src/math_private_openbsd.h create mode 100644 src/polevll.c diff --git a/Make.inc b/Make.inc index 83d9612..57e4ef2 100644 --- a/Make.inc +++ b/Make.inc @@ -3,8 +3,8 @@ OS := $(shell uname) # Do not forget to bump SOMINOR when changing VERSION, # and SOMAJOR when breaking ABI in a backward-incompatible way -VERSION = 0.4 -SOMAJOR = 1 +VERSION = 0.5 +SOMAJOR = 2 SOMINOR = 0 DESTDIR = prefix = /usr/local @@ -64,7 +64,7 @@ endif ifeq ($(ARCH),x86_64) override ARCH := amd64 -endif +endif ifneq (,$(findstring MINGW,$(OS))) override OS=WINNT diff --git a/ld80/Make.files b/ld80/Make.files index 506883d..6e31da3 100644 --- a/ld80/Make.files +++ b/ld80/Make.files @@ -1,10 +1,12 @@ $(CUR_SRCS) += invtrig.c \ - e_acoshl.c e_hypotl.c e_powl.c k_tanl.c s_exp2l.c s_nanl.c \ + e_acoshl.c e_powl.c k_tanl.c s_exp2l.c \ e_atanhl.c e_lgammal.c e_sinhl.c s_asinhl.c s_expm1l.c \ - e_coshl.c e_log10l.c e_tgammal.c s_floorl.c s_nextafterl.c \ - e_expl.c e_log2l.c k_cosl.c s_ceill.c s_log1pl.c s_tanhl.c \ - e_fmodl.c e_logl.c k_sinl.c s_erfl.c s_modfl.c s_truncl.c -# s_remquol.c + e_coshl.c e_log10l.c e_tgammal.c \ + e_expl.c e_log2l.c k_cosl.c s_log1pl.c s_tanhl.c \ + e_logl.c k_sinl.c s_erfl.c + +# s_remquol.c e_fmodl.c s_truncl.c +# e_hypotl.c s_floorl.c s_nextafterl.c s_ceill.c s_modfl.c ifneq ($(OS), WINNT) $(CUR_SRCS) += s_nanl.c diff --git a/ld80/e_acoshl.c b/ld80/e_acoshl.c index 0756fad..78fe296 100644 --- a/ld80/e_acoshl.c +++ b/ld80/e_acoshl.c @@ -24,7 +24,7 @@ * acoshl(NaN) is NaN without signal. */ -#include +#include #include "math_private.h" diff --git a/ld80/e_atanhl.c b/ld80/e_atanhl.c index d75a757..86e6df3 100644 --- a/ld80/e_atanhl.c +++ b/ld80/e_atanhl.c @@ -28,7 +28,7 @@ * */ -#include +#include #include "math_private.h" diff --git a/ld80/e_coshl.c b/ld80/e_coshl.c index 25349db..2d0a114 100644 --- a/ld80/e_coshl.c +++ b/ld80/e_coshl.c @@ -31,7 +31,7 @@ * only coshl(0)=1 is exact for finite x. */ -#include "math.h" +#include "openlibm.h" #include "math_private.h" static const long double one = 1.0, half=0.5, huge = 1.0e4900L; diff --git a/ld80/e_expl.c b/ld80/e_expl.c index 2e1008a..8e4fea5 100644 --- a/ld80/e_expl.c +++ b/ld80/e_expl.c @@ -72,7 +72,7 @@ /* Exponential function */ -#include +#include #include "math_private.h" diff --git a/ld80/e_hypotl.c b/ld80/e_hypotl.c index e70d2b1..f69c35c 100644 --- a/ld80/e_hypotl.c +++ b/ld80/e_hypotl.c @@ -42,7 +42,7 @@ * than 1 ulps (units in the last place) */ -#include +#include #include "math_private.h" diff --git a/ld80/e_lgammal.c b/ld80/e_lgammal.c index 04c0aef..0f0970d 100644 --- a/ld80/e_lgammal.c +++ b/ld80/e_lgammal.c @@ -86,9 +86,10 @@ * */ -#include +#include #include "math_private.h" +extern int signgam; static const long double half = 0.5L, diff --git a/ld80/e_log10l.c b/ld80/e_log10l.c index f314fd1..4b0d84f 100644 --- a/ld80/e_log10l.c +++ b/ld80/e_log10l.c @@ -63,7 +63,7 @@ * log domain: x < 0; returns MINLOG */ -#include +#include #include "math_private.h" diff --git a/ld80/e_log2l.c b/ld80/e_log2l.c index e02a34e..ea337c0 100644 --- a/ld80/e_log2l.c +++ b/ld80/e_log2l.c @@ -63,7 +63,7 @@ * log domain: x < 0; returns NAN */ -#include +#include #include "math_private.h" diff --git a/ld80/e_logl.c b/ld80/e_logl.c index 7c1b854..29e8490 100644 --- a/ld80/e_logl.c +++ b/ld80/e_logl.c @@ -63,7 +63,7 @@ * log domain: x < 0; returns NAN */ -#include +#include #include "math_private.h" diff --git a/ld80/e_powl.c b/ld80/e_powl.c index 9581151..128ebab 100644 --- a/ld80/e_powl.c +++ b/ld80/e_powl.c @@ -73,7 +73,7 @@ */ #include -#include +#include #include "math_private.h" diff --git a/ld80/e_rem_pio2l.h b/ld80/e_rem_pio2l.h index 2e7d0e7..65ccaaa 100644 --- a/ld80/e_rem_pio2l.h +++ b/ld80/e_rem_pio2l.h @@ -6,7 +6,7 @@ * * Developed at SunSoft, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== * @@ -17,8 +17,8 @@ //__FBSDID("$FreeBSD: src/lib/msun/ld80/e_rem_pio2l.h,v 1.3 2011/06/18 13:56:33 benl Exp $"); /* ld80 version of __ieee754_rem_pio2l(x,y) - * - * return the remainder of x rem pi/2 in y[0]+y[1] + * + * return the remainder of x rem pi/2 in y[0]+y[1] * use __kernel_rem_pio2() */ @@ -26,7 +26,7 @@ #include "openlibm.h" #include "math_private.h" -#include "fpmath.h" +#include "openlibm.h" #define BIAS (LDBL_MAX_EXP - 1) @@ -102,24 +102,24 @@ __ieee754_rem_pio2l(long double x, long double *y) union IEEEl2bits u2; int ex1; j = ex; - y[0] = r-w; + y[0] = r-w; u2.e = y[0]; ex1 = u2.xbits.expsign & 0x7fff; i = j-ex1; if(i>22) { /* 2nd iteration needed, good to 141 */ t = r; - w = fn*pio2_2; + w = fn*pio2_2; r = t-w; - w = fn*pio2_2t-((t-r)-w); + w = fn*pio2_2t-((t-r)-w); y[0] = r-w; u2.e = y[0]; ex1 = u2.xbits.expsign & 0x7fff; i = j-ex1; if(i>61) { /* 3rd iteration need, 180 bits acc */ t = r; /* will cover all possible cases */ - w = fn*pio2_3; + w = fn*pio2_3; r = t-w; - w = fn*pio2_3t-((t-r)-w); + w = fn*pio2_3t-((t-r)-w); y[0] = r-w; } } @@ -127,7 +127,7 @@ __ieee754_rem_pio2l(long double x, long double *y) y[1] = (r-y[0])-w; return n; } - /* + /* * all other (large) arguments */ if(ex==0x7fff) { /* x is inf or NaN */ diff --git a/ld80/e_sinhl.c b/ld80/e_sinhl.c index b5b0be7..fa1b102 100644 --- a/ld80/e_sinhl.c +++ b/ld80/e_sinhl.c @@ -28,7 +28,7 @@ * only sinhl(0)=0 is exact for finite x. */ -#include +#include #include "math_private.h" diff --git a/ld80/e_tgammal.c b/ld80/e_tgammal.c index e7b9505..571a0a9 100644 --- a/ld80/e_tgammal.c +++ b/ld80/e_tgammal.c @@ -58,9 +58,10 @@ */ #include -#include +#include #include "math_private.h" +extern int signgam; /* tgamma(x+2) = tgamma(x+2) P(x)/Q(x) diff --git a/ld80/invtrig.h b/ld80/invtrig.h index d4dfd5c..489fc81 100644 --- a/ld80/invtrig.h +++ b/ld80/invtrig.h @@ -28,7 +28,7 @@ #include -#include "fpmath.h" +#include #define BIAS (LDBL_MAX_EXP - 1) #define MANH_SIZE LDBL_MANH_SIZE diff --git a/ld80/s_asinhl.c b/ld80/s_asinhl.c index 1dc804b..35c1221 100644 --- a/ld80/s_asinhl.c +++ b/ld80/s_asinhl.c @@ -21,7 +21,7 @@ * := signl(x)*log1pl(|x| + x^2/(1 + sqrtl(1+x^2))) */ -#include +#include #include "math_private.h" diff --git a/ld80/s_ceill.c b/ld80/s_ceill.c index bff5277..663f295 100644 --- a/ld80/s_ceill.c +++ b/ld80/s_ceill.c @@ -19,7 +19,7 @@ * Inexact flag raised if x not equal to ceil(x). */ -#include +#include #include "math_private.h" diff --git a/ld80/s_erfl.c b/ld80/s_erfl.c index 0bb8931..0c8eca5 100644 --- a/ld80/s_erfl.c +++ b/ld80/s_erfl.c @@ -99,7 +99,7 @@ */ -#include +#include #include "math_private.h" diff --git a/ld80/s_exp2l.c b/ld80/s_exp2l.c index d76f298..0c89159 100644 --- a/ld80/s_exp2l.c +++ b/ld80/s_exp2l.c @@ -32,7 +32,7 @@ #include "amd64/bsd_ieeefp.h" -#include "fpmath.h" +#include "openlibm.h" #include "openlibm.h" #include "math_private.h" diff --git a/ld80/s_expm1l.c b/ld80/s_expm1l.c index 06b0539..bd965da 100644 --- a/ld80/s_expm1l.c +++ b/ld80/s_expm1l.c @@ -57,7 +57,7 @@ * */ -#include +#include static const long double MAXLOGL = 1.1356523406294143949492E4L; diff --git a/ld80/s_floorl.c b/ld80/s_floorl.c index a5b7877..77ffb7a 100644 --- a/ld80/s_floorl.c +++ b/ld80/s_floorl.c @@ -19,7 +19,7 @@ * Inexact flag raised if x not equal to floor(x). */ -#include +#include #include "math_private.h" diff --git a/ld80/s_log1pl.c b/ld80/s_log1pl.c index 2fbac5a..decb0c2 100644 --- a/ld80/s_log1pl.c +++ b/ld80/s_log1pl.c @@ -59,7 +59,7 @@ * log domain: x-1 < 0; returns NAN */ -#include +#include #include "math_private.h" diff --git a/ld80/s_modfl.c b/ld80/s_modfl.c index ebfca4b..e06f868 100644 --- a/ld80/s_modfl.c +++ b/ld80/s_modfl.c @@ -20,7 +20,7 @@ * No exception. */ -#include +#include #include "math_private.h" diff --git a/ld80/s_nanl.c b/ld80/s_nanl.c index d56fdd4..7b0fe01 100644 --- a/ld80/s_nanl.c +++ b/ld80/s_nanl.c @@ -28,7 +28,7 @@ #include "openlibm.h" -#include "fpmath.h" +#include "openlibm.h" #include "math_private.h" DLLEXPORT long double diff --git a/ld80/s_nextafterl.c b/ld80/s_nextafterl.c index 0c28b02..2e431c9 100644 --- a/ld80/s_nextafterl.c +++ b/ld80/s_nextafterl.c @@ -17,7 +17,7 @@ * Special cases: */ -#include +#include #include "math_private.h" @@ -87,4 +87,4 @@ nextafterl(long double x, long double y) return x; } -__strong_alias(nexttowardl, nextafterl); +//__strong_alias(nexttowardl, nextafterl); diff --git a/ld80/s_nexttoward.c b/ld80/s_nexttoward.c index a9d6773..d5327ee 100644 --- a/ld80/s_nexttoward.c +++ b/ld80/s_nexttoward.c @@ -17,7 +17,7 @@ * Special cases: */ -#include +#include #include #include "math_private.h" diff --git a/ld80/s_nexttowardf.c b/ld80/s_nexttowardf.c index 3e0a2af..1941bf2 100644 --- a/ld80/s_nexttowardf.c +++ b/ld80/s_nexttowardf.c @@ -10,7 +10,7 @@ * ==================================================== */ -#include +#include #include #include "math_private.h" diff --git a/ld80/s_remquol.c b/ld80/s_remquol.c index 244c105..330f5c6 100644 --- a/ld80/s_remquol.c +++ b/ld80/s_remquol.c @@ -14,7 +14,7 @@ #include #include -#include +#include #include #include "math_private.h" diff --git a/ld80/s_tanhl.c b/ld80/s_tanhl.c index a0b7bd8..bca28bb 100644 --- a/ld80/s_tanhl.c +++ b/ld80/s_tanhl.c @@ -34,7 +34,7 @@ * only tanhl(0)=0 is exact for finite argument. */ -#include +#include #include "math_private.h" diff --git a/ld80/s_truncl.c b/ld80/s_truncl.c index 3d8a1a5..e4aa68c 100644 --- a/ld80/s_truncl.c +++ b/ld80/s_truncl.c @@ -21,10 +21,10 @@ */ #include -#include +//#include #include -#include +#include "openlibm.h" #include #include "math_private.h" diff --git a/src/Make.files b/src/Make.files index 2775cb9..e8a8204 100644 --- a/src/Make.files +++ b/src/Make.files @@ -50,7 +50,7 @@ $(CUR_SRCS) += e_acosl.c e_asinl.c e_atan2l.c e_fmodl.c \ s_frexpl.c s_logbl.c s_nexttoward.c \ s_remquol.c \ s_sinl.c s_sincosl.c s_tanl.c s_truncl.c w_cabsl.c \ - s_nextafterl.c s_rintl.c s_scalbnl.c + s_nextafterl.c s_rintl.c s_scalbnl.c polevll.c # s_cbrtl.c endif diff --git a/src/math_private.h b/src/math_private.h index e728ba8..a8c0595 100644 --- a/src/math_private.h +++ b/src/math_private.h @@ -22,6 +22,7 @@ #include "fpmath.h" #include #include +#include "math_private_openbsd.h" /* * The original fdlibm code used statements like: @@ -228,7 +229,7 @@ do { \ * Common routine to process the arguments to nan(), nanf(), and nanl(). */ void _scan_nan(u_int32_t *__words, int __num_words, const char *__s); - + #ifdef __GNUCLIKE_ASM /* Asm versions of some functions. */ diff --git a/src/math_private_openbsd.h b/src/math_private_openbsd.h new file mode 100644 index 0000000..bd0e5b2 --- /dev/null +++ b/src/math_private_openbsd.h @@ -0,0 +1,218 @@ +/* $OpenBSD: math_private.h,v 1.17 2014/06/02 19:31:17 kettenis Exp $ */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* + * from: @(#)fdlibm.h 5.1 93/09/24 + */ + +#ifndef _MATH_PRIVATE_OPENBSD_H_ +#define _MATH_PRIVATE_OPENBSD_H_ + +#if _IEEE_WORD_ORDER == _BIG_ENDIAN + +typedef union +{ + long double value; + struct { + u_int32_t mswhi; + u_int32_t mswlo; + u_int32_t lswhi; + u_int32_t lswlo; + } parts32; + struct { + u_int64_t msw; + u_int64_t lsw; + } parts64; +} ieee_quad_shape_type; + +#endif + +#if _IEEE_WORD_ORDER == _LITTLE_ENDIAN + +typedef union +{ + long double value; + struct { + u_int32_t lswlo; + u_int32_t lswhi; + u_int32_t mswlo; + u_int32_t mswhi; + } parts32; + struct { + u_int64_t lsw; + u_int64_t msw; + } parts64; +} ieee_quad_shape_type; + +#endif + +/* Get two 64 bit ints from a long double. */ + +#define GET_LDOUBLE_WORDS64(ix0,ix1,d) \ +do { \ + ieee_quad_shape_type qw_u; \ + qw_u.value = (d); \ + (ix0) = qw_u.parts64.msw; \ + (ix1) = qw_u.parts64.lsw; \ +} while (0) + +/* Set a long double from two 64 bit ints. */ + +#define SET_LDOUBLE_WORDS64(d,ix0,ix1) \ +do { \ + ieee_quad_shape_type qw_u; \ + qw_u.parts64.msw = (ix0); \ + qw_u.parts64.lsw = (ix1); \ + (d) = qw_u.value; \ +} while (0) + +/* Get the more significant 64 bits of a long double mantissa. */ + +#define GET_LDOUBLE_MSW64(v,d) \ +do { \ + ieee_quad_shape_type sh_u; \ + sh_u.value = (d); \ + (v) = sh_u.parts64.msw; \ +} while (0) + +/* Set the more significant 64 bits of a long double mantissa from an int. */ + +#define SET_LDOUBLE_MSW64(d,v) \ +do { \ + ieee_quad_shape_type sh_u; \ + sh_u.value = (d); \ + sh_u.parts64.msw = (v); \ + (d) = sh_u.value; \ +} while (0) + +/* Get the least significant 64 bits of a long double mantissa. */ + +#define GET_LDOUBLE_LSW64(v,d) \ +do { \ + ieee_quad_shape_type sh_u; \ + sh_u.value = (d); \ + (v) = sh_u.parts64.lsw; \ +} while (0) + +/* A union which permits us to convert between a long double and + three 32 bit ints. */ + +#if _IEEE_WORD_ORDER == _BIG_ENDIAN + +typedef union +{ + long double value; + struct { +#ifdef __LP64__ + int padh:32; +#endif + int exp:16; + int padl:16; + u_int32_t msw; + u_int32_t lsw; + } parts; +} ieee_extended_shape_type; + +#endif + +#if _IEEE_WORD_ORDER == _LITTLE_ENDIAN + +typedef union +{ + long double value; + struct { + u_int32_t lsw; + u_int32_t msw; + int exp:16; + int padl:16; +#ifdef __LP64__ + int padh:32; +#endif + } parts; +} ieee_extended_shape_type; + +#endif + +/* Get three 32 bit ints from a double. */ + +#define GET_LDOUBLE_WORDS(se,ix0,ix1,d) \ +do { \ + ieee_extended_shape_type ew_u; \ + ew_u.value = (d); \ + (se) = ew_u.parts.exp; \ + (ix0) = ew_u.parts.msw; \ + (ix1) = ew_u.parts.lsw; \ +} while (0) + +/* Set a double from two 32 bit ints. */ + +#define SET_LDOUBLE_WORDS(d,se,ix0,ix1) \ +do { \ + ieee_extended_shape_type iw_u; \ + iw_u.parts.exp = (se); \ + iw_u.parts.msw = (ix0); \ + iw_u.parts.lsw = (ix1); \ + (d) = iw_u.value; \ +} while (0) + +/* Get the more significant 32 bits of a long double mantissa. */ + +#define GET_LDOUBLE_MSW(v,d) \ +do { \ + ieee_extended_shape_type sh_u; \ + sh_u.value = (d); \ + (v) = sh_u.parts.msw; \ +} while (0) + +/* Set the more significant 32 bits of a long double mantissa from an int. */ + +#define SET_LDOUBLE_MSW(d,v) \ +do { \ + ieee_extended_shape_type sh_u; \ + sh_u.value = (d); \ + sh_u.parts.msw = (v); \ + (d) = sh_u.value; \ +} while (0) + +/* Get int from the exponent of a long double. */ + +#define GET_LDOUBLE_EXP(se,d) \ +do { \ + ieee_extended_shape_type ge_u; \ + ge_u.value = (d); \ + (se) = ge_u.parts.exp; \ +} while (0) + +/* Set exponent of a long double from an int. */ + +#define SET_LDOUBLE_EXP(d,se) \ +do { \ + ieee_extended_shape_type se_u; \ + se_u.value = (d); \ + se_u.parts.exp = (se); \ + (d) = se_u.value; \ +} while (0) + +/* + * Common routine to process the arguments to nan(), nanf(), and nanl(). + */ +void _scan_nan(uint32_t *__words, int __num_words, const char *__s); + +/* + * Functions internal to the math package, yet not static. + */ +double __exp__D(double, double); +struct Double __log__D(double); +long double __p1evll(long double, void *, int); +long double __polevll(long double, void *, int); + +#endif /* _MATH_PRIVATE_OPENBSD_H_ */ diff --git a/src/polevll.c b/src/polevll.c new file mode 100644 index 0000000..78cd75c --- /dev/null +++ b/src/polevll.c @@ -0,0 +1,104 @@ +/* $OpenBSD: polevll.c,v 1.2 2013/11/12 20:35:09 martynas Exp $ */ + +/* + * Copyright (c) 2008 Stephen L. Moshier + * + * Permission to use, copy, modify, and distribute this software for any + * purpose with or without fee is hereby granted, provided that the above + * copyright notice and this permission notice appear in all copies. + * + * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES + * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF + * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR + * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES + * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN + * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF + * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. + */ + +/* polevll.c + * p1evll.c + * + * Evaluate polynomial + * + * + * + * SYNOPSIS: + * + * int N; + * long double x, y, coef[N+1], polevl[]; + * + * y = polevll( x, coef, N ); + * + * + * + * DESCRIPTION: + * + * Evaluates polynomial of degree N: + * + * 2 N + * y = C + C x + C x +...+ C x + * 0 1 2 N + * + * Coefficients are stored in reverse order: + * + * coef[0] = C , ..., coef[N] = C . + * N 0 + * + * The function p1evll() assumes that coef[N] = 1.0 and is + * omitted from the array. Its calling arguments are + * otherwise the same as polevll(). + * + * + * SPEED: + * + * In the interest of speed, there are no checks for out + * of bounds arithmetic. This routine is used by most of + * the functions in the library. Depending on available + * equipment features, the user may wish to rewrite the + * program in microcode or assembly language. + * + */ + +#include "openlibm.h" + +#include "math_private.h" + +/* + * Polynomial evaluator: + * P[0] x^n + P[1] x^(n-1) + ... + P[n] + */ +long double +__polevll(long double x, void *PP, int n) +{ + long double y; + long double *P; + + P = (long double *)PP; + y = *P++; + do { + y = y * x + *P++; + } while (--n); + + return (y); +} + +/* + * Polynomial evaluator: + * x^n + P[0] x^(n-1) + P[1] x^(n-2) + ... + P[n] + */ +long double +__p1evll(long double x, void *PP, int n) +{ + long double y; + long double *P; + + P = (long double *)PP; + n -= 1; + y = x + *P++; + do { + y = y * x + *P++; + } while (--n); + + return (y); +}