From bfdb9ceb4782d93166f3eb4b9432097474b7dfc9 Mon Sep 17 00:00:00 2001 From: Nathan Date: Mon, 1 Jul 2024 07:00:46 -0600 Subject: [PATCH] Match fdlibm; re-organize flags.json (#276) * MSL_C float.c * MSL_C bool_def.h * Re-organize flag lists * Bring in fdlibm sources Not fully matching, some modifications need to be made * Pick off most fdlibm matches * Finish fdlibm --- config/SZBE69_B8/flags.json | 109 +++--- config/SZBE69_B8/objects.json | 67 +++- .../MSL/MSL_C/MSL_Common/bool_def.h | 13 + .../MSL/MSL_C/MSL_Common/float.c | 7 + .../MSL/MSL_C/MSL_Common/float.h | 20 + .../MSL/MSL_C/MSL_Common/math_api.h | 17 - .../Math/Double_precision/e_acos.c | 112 ++++++ .../Math/Double_precision/e_asin.c | 120 ++++++ .../Math/Double_precision/e_atan2.c | 144 ++++++++ .../Math/Double_precision/e_exp.c | 162 ++++++++ .../Math/Double_precision/e_fmod.c | 178 +++++++++ .../Math/Double_precision/e_log.c | 158 ++++++++ .../Math/Double_precision/e_log10.c | 99 +++++ .../Math/Double_precision/e_pow.c | 343 +++++++++++++++++ .../Math/Double_precision/e_rem_pio2.c | 186 ++++++++++ .../Math/{ => Double_precision}/fdlibm.h | 93 ++--- .../Math/Double_precision/k_cos.c | 93 +++++ .../Math/Double_precision/k_rem_pio2.c | 349 ++++++++++++++++++ .../Math/Double_precision/k_sin.c | 80 ++++ .../Math/Double_precision/k_tan.c | 133 +++++++ .../Math/Double_precision/readme | 262 +++++++++++++ .../Math/Double_precision/s_atan.c | 143 +++++++ .../Math/Double_precision/s_ceil.c | 91 +++++ .../Math/Double_precision/s_copysign.c | 31 ++ .../Math/Double_precision/s_cos.c | 83 +++++ .../Math/Double_precision/s_floor.c | 90 +++++ .../Math/Double_precision/s_frexp.c | 58 +++ .../Math/Double_precision/s_ldexp.c | 70 ++++ .../Math/Double_precision/s_modf.c | 80 ++++ .../Math/Double_precision/s_rint.c | 92 +++++ .../Math/Double_precision/s_sin.c | 83 +++++ .../Math/Double_precision/s_tan.c | 74 ++++ .../Math/Double_precision/w_acos.c | 41 +- .../Math/Double_precision/w_asin.c | 42 ++- .../Math/Double_precision/w_atan2.c | 42 ++- .../Math/Double_precision/w_exp.c | 51 ++- .../Math/Double_precision/w_fmod.c | 41 +- .../Math/Double_precision/w_log.c | 41 +- .../Math/Double_precision/w_log10.c | 44 ++- .../Math/Double_precision/w_pow.c | 62 +++- .../PowerPC_EABI_Support/MSL/MSL_C/errno.h | 2 +- src/sdk/PowerPC_EABI_Support/MSL/MSL_C/math.h | 147 +------- .../PowerPC_EABI_Support/MSL/MSL_C/stdbool.h | 9 +- src/sdk/RVL_SDK/revolution/math.h | 13 +- 44 files changed, 3866 insertions(+), 309 deletions(-) create mode 100644 src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common/bool_def.h create mode 100644 src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common/float.c create mode 100644 src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common/float.h create mode 100644 src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_acos.c create mode 100644 src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_asin.c create mode 100644 src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_atan2.c create mode 100644 src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_exp.c create mode 100644 src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_fmod.c create mode 100644 src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_log.c create mode 100644 src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_log10.c create mode 100644 src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_pow.c create mode 100644 src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_rem_pio2.c rename src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/{ => Double_precision}/fdlibm.h (77%) create mode 100644 src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/k_cos.c create mode 100644 src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/k_rem_pio2.c create mode 100644 src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/k_sin.c create mode 100644 src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/k_tan.c create mode 100644 src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/readme create mode 100644 src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_atan.c create mode 100644 src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_ceil.c create mode 100644 src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_copysign.c create mode 100644 src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_cos.c create mode 100644 src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_floor.c create mode 100644 src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_frexp.c create mode 100644 src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_ldexp.c create mode 100644 src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_modf.c create mode 100644 src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_rint.c create mode 100644 src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_sin.c create mode 100644 src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_tan.c diff --git a/config/SZBE69_B8/flags.json b/config/SZBE69_B8/flags.json index 6d894e94..117fd47c 100644 --- a/config/SZBE69_B8/flags.json +++ b/config/SZBE69_B8/flags.json @@ -28,29 +28,6 @@ "-d MILO_DEBUG" ] }, - "runtime": { - "base": "base", - "flags": [ - "-str reuse,pool,readonly", - "-pragma \"merge_float_consts on\"", - - "-use_lmw_stmw on", - "-common off", - - "-lang=c99", - "-O4,p", - "-inline auto", - "-func_align 4" - ] - }, - "rvl_sdk": { - "base": "base", - "flags": [ - "-func_align 16", - "-O4,p", - "-Isrc/sdk/PowerPC_EABI_Support/MSL" - ] - }, "main": { "base": "base", "flags": [ @@ -82,15 +59,11 @@ "-ipa file" ] }, - "system": { - "base": "band3", - "flags": [] - }, "network": { "base": "band3", "flags": [] }, - "json_c": { + "network/json_c": { "base": "base", "flags": [ "-sdata 2", @@ -106,7 +79,7 @@ "-ipa file" ] }, - "quazal": { + "network/Platform": { "base": "base", "flags": [ "-sdata 2", @@ -120,29 +93,54 @@ "-inline auto" ] }, - "zlib": { + "system": { + "base": "band3", + "flags": [] + }, + "system/flex": { + "base": "system/zlib", + "flags": [] + }, + "system/speex": { + "base": "system/zlib", + "flags": [ + "-fp_contract on" + ] + }, + "system/soundtouch": { "base": "base", "flags": [ "-sdata 2", "-sdata2 2", - - "-pool on", "-str reuse,pool", "-pragma \"merge_float_consts on\"", + "-fp_contract on", - "-lang=c99", "-O4,p", + "-lang=c++", "-inline auto" ] }, - "zlib/zutil": { + "system/tomcrypt": { + "base": "system/zlib", + "flags": [ + "-inline off" + ] + }, + "system/vorbis": { + "base": "system/zlib", + "flags": [ + "-fp_contract on" + ] + }, + "system/zlib": { "base": "base", "flags": [ "-sdata 2", "-sdata2 2", "-pool on", - "-str nopool", + "-str reuse,pool", "-pragma \"merge_float_consts on\"", "-lang=c99", @@ -150,43 +148,34 @@ "-inline auto" ] }, - "soundtouch": { + "sdk/MSL": { "base": "base", "flags": [ - "-sdata 2", - "-sdata2 2", - "-str reuse,pool", - "-pragma \"merge_float_consts on\"", - "-fp_contract on", + "-str reuse,pool,readonly", + "-use_lmw_stmw on", + "-common off", + + "-lang=c99", "-O4,p", - "-lang=c++", - "-inline auto" - ] - }, - "libtomcrypt": { - "base": "zlib", - "flags": [ - "-inline off" + "-inline auto", + "-func_align 4" ] }, - "speex": { - "base": "zlib", + "sdk/MSL/fdlibm": { + "base": "sdk/MSL", "flags": [ - "-fp_contract on" + "-d _IEEE_LIBM" ] }, - "flex": { - "base": "zlib", - "flags": [] - }, - "vorbis": { - "base": "zlib", + "sdk/RVL_SDK": { + "base": "base", "flags": [ - "-fp_contract on" + "-func_align 16", + "-O4,p" ] }, - "bt": { + "sdk/RVL_SDK/bte": { "base": "base", "flags": [ "-sdata 2", diff --git a/config/SZBE69_B8/objects.json b/config/SZBE69_B8/objects.json index c3e662ac..3d9dcb43 100644 --- a/config/SZBE69_B8/objects.json +++ b/config/SZBE69_B8/objects.json @@ -83,7 +83,7 @@ }, "network/json_c": { "mw_version": "Wii/1.3", - "cflags": "json_c", + "cflags": "network/json_c", "objects": { "network/net/json-c/arraylist.c": "Matching", "network/net/json-c/debug.c": "Matching", @@ -95,7 +95,7 @@ }, "network/Platform": { "mw_version": "Wii/1.3", - "cflags": "quazal", + "cflags": "network/Platform", "objects": { "network/Platform/CriticalSection.cpp": "Matching", "network/Platform/DateTime.cpp": "NonMatching", @@ -562,14 +562,14 @@ }, "system/flex": { "mw_version": "Wii/1.3", - "cflags": "flex", + "cflags": "system/flex", "objects": { "system/obj/DataFlex.c": "NonMatching" } }, "system/speex": { "mw_version": "Wii/1.3", - "cflags": "speex", + "cflags": "system/speex", "objects": { "system/speex/libspeex/bits.c": "NonMatching", "system/speex/libspeex/cb_search.c": "Matching", @@ -598,7 +598,7 @@ }, "system/soundtouch": { "mw_version": "Wii/1.3", - "cflags": "soundtouch", + "cflags": "system/soundtouch", "objects": { "system/synthwii/soundtouch/source/SoundTouch/AAFilter.cpp": "NonMatching", "system/synthwii/soundtouch/source/SoundTouch/FIFOSampleBuffer.cpp": "NonMatching", @@ -610,7 +610,7 @@ }, "system/tomcrypt": { "mw_version": "Wii/1.3", - "cflags": "libtomcrypt", + "cflags": "system/tomcrypt", "objects": { "system/synth/tomcrypt/aes.c": "Matching", "system/synth/tomcrypt/crypt.c": "Matching", @@ -619,7 +619,7 @@ }, "system/vorbis": { "mw_version": "Wii/1.3", - "cflags": "vorbis", + "cflags": "system/vorbis", "objects": { "system/oggvorbis/analysis.c": "NonMatching", "system/oggvorbis/bitrate.c": "Matching", @@ -645,7 +645,7 @@ }, "system/zlib": { "mw_version": "Wii/1.3", - "cflags": "zlib", + "cflags": "system/zlib", "objects": { "system/zlib/adler32.c": "Matching", "system/zlib/crc32.c": "Matching", @@ -654,7 +654,7 @@ "system/zlib/inflate.c": "Matching", "system/zlib/inftrees.c": "Matching", "system/zlib/inffast.c": "Matching", - "system/zlib/zutil.c": { "status": "Matching", "cflags": "zlib/zutil" } + "system/zlib/zutil.c": { "status": "Matching", "extra_cflags": ["-str nopool"] } } }, "lib/binkwii": { @@ -666,13 +666,47 @@ }, "sdk/MSL_C": { "mw_version": "Wii/1.3", - "cflags": "runtime", + "cflags": "sdk/MSL", "objects": { + "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common/alloc.c": "NonMatching", "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common/arith.c": "Matching", "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common/buffer_io.c": "NonMatching", "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common/errno.c": "Matching", + "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common/float.c": "Matching", "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common/rand.c": "Matching", + "sdk/PowerPC_EABI_Support/MSL/MSL_C/PPC_EABI/time.dolphin.c": "Matching", + "sdk/PowerPC_EABI_Support/MSL/MSL_C/PPC_EABI/math_ppc.c": "Matching" + } + }, + "sdk/MSL_C/fdlibm": { + "mw_version": "Wii/1.3", + "cflags": "sdk/MSL/fdlibm", + "objects": { + "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_acos.c": "Matching", + "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_asin.c": "Matching", + "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_atan2.c": "Matching", + "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_exp.c": "Matching", + "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_fmod.c": "Matching", + "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_log.c": "Matching", + "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_log10.c": "Matching", + "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_pow.c": "Matching", + "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_rem_pio2.c": "Matching", + "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/k_cos.c": "Matching", + "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/k_rem_pio2.c": "Matching", + "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/k_sin.c": "Matching", + "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/k_tan.c": "Matching", + "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_atan.c": "Matching", + "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_ceil.c": "Matching", + "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_copysign.c": "Matching", + "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_cos.c": "Matching", + "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_floor.c": "Matching", + "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_frexp.c": "Matching", + "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_ldexp.c": "Matching", + "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_modf.c": "Matching", + "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_rint.c": "Matching", + "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_sin.c": "Matching", + "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_tan.c": "Matching", "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/w_acos.c": "Matching", "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/w_asin.c": "Matching", "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/w_atan2.c": "Matching", @@ -680,36 +714,33 @@ "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/w_fmod.c": "Matching", "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/w_log.c": "Matching", "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/w_log10.c": "Matching", - "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/w_pow.c": "Matching", - - "sdk/PowerPC_EABI_Support/MSL/MSL_C/PPC_EABI/time.dolphin.c": "Matching", - "sdk/PowerPC_EABI_Support/MSL/MSL_C/PPC_EABI/math_ppc.c": "Matching" + "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/w_pow.c": "Matching" } }, "sdk/MSL_C++": { "mw_version": "Wii/1.3", - "cflags": "runtime", + "cflags": "sdk/MSL", "objects": { } }, "sdk/Runtime": { "mw_version": "Wii/1.3", - "cflags": "runtime", + "cflags": "sdk/MSL", "objects": { "sdk/PowerPC_EABI_Support/Runtime/Gecko_setjmp.cpp": "Matching" } }, "sdk/MetroTRK": { "mw_version": "Wii/1.3", - "cflags": "runtime", + "cflags": "sdk/MSL", "objects": { } }, "sdk/RVL_SDK": { "mw_version": "Wii/1.3", - "cflags": "rvl_sdk", + "cflags": "sdk/RVL_SDK", "objects": { "sdk/RVL_SDK/revolution/gx/GXBump.c": "NonMatching", "sdk/RVL_SDK/revolution/gx/GXDisplayList.c": "Matching", diff --git a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common/bool_def.h b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common/bool_def.h new file mode 100644 index 00000000..ea84630b --- /dev/null +++ b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common/bool_def.h @@ -0,0 +1,13 @@ +#ifndef _MSL_BOOL_DEF_H_ +#define _MSL_BOOL_DEF_H_ + +#ifndef __cplusplus + +typedef unsigned char bool; + +#define false 0 +#define true 1 + +#endif + +#endif diff --git a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common/float.c b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common/float.c new file mode 100644 index 00000000..d2c4b8c5 --- /dev/null +++ b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common/float.c @@ -0,0 +1,7 @@ +// Seems they forgot the include, these are stored in .sdata here +// but elsewhere they're are read as if they're in .data +// #include "MSL_Common/float.h" + +int __float_nan = 0x7FFFFFFF; +float __float_huge = 0x1.0P128F; +double __double_huge = 0x1.0P1024; diff --git a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common/float.h b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common/float.h new file mode 100644 index 00000000..e5ba266d --- /dev/null +++ b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common/float.h @@ -0,0 +1,20 @@ +#ifndef MSL_COMMON_FLOAT_H +#define MSL_COMMON_FLOAT_H + +#ifdef __cplusplus +extern "C" { +#endif + +__declspec(section ".data") extern int __float_nan; +__declspec(section ".data") extern float __float_huge; +__declspec(section ".data") extern double __double_huge; + +#define INFINITY (__float_huge) +#define NAN (*(float *)&__float_nan) +#define HUGE_VAL (__double_huge) + +#ifdef __cplusplus +} +#endif + +#endif diff --git a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common/math_api.h b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common/math_api.h index 6720c05e..1ddbdeff 100644 --- a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common/math_api.h +++ b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common/math_api.h @@ -5,23 +5,6 @@ extern "C" { #endif -#if defined(i386) || defined(i486) || defined(intel) || defined(x86) || defined(i86pc) \ - || defined(__alpha) || defined(__osf__) -#define __LITTLE_ENDIAN -#endif - -#ifdef __LITTLE_ENDIAN -#define __HI(x) *(1 + (int *)&x) -#define __LO(x) *(int *)&x -#define __HIp(x) *(1 + (int *)x) -#define __LOp(x) *(int *)x -#else -#define __HI(x) *(int *)&x -#define __LO(x) *(1 + (int *)&x) -#define __HIp(x) *(int *)x -#define __LOp(x) *(1 + (int *)x) -#endif - typedef float float_t; typedef double double_t; diff --git a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_acos.c b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_acos.c new file mode 100644 index 00000000..437e84e9 --- /dev/null +++ b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_acos.c @@ -0,0 +1,112 @@ + +/* @(#)e_acos.c 1.3 95/01/18 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* __ieee754_acos(x) + * Method : + * acos(x) = pi/2 - asin(x) + * acos(-x) = pi/2 + asin(x) + * For |x|<=0.5 + * acos(x) = pi/2 - (x + x*x^2*R(x^2)) (see asin.c) + * For x>0.5 + * acos(x) = pi/2 - (pi/2 - 2asin(sqrt((1-x)/2))) + * = 2asin(sqrt((1-x)/2)) + * = 2s + 2s*z*R(z) ...z=(1-x)/2, s=sqrt(z) + * = 2f + (2c + 2s*z*R(z)) + * where f=hi part of s, and c = (z-f*f)/(s+f) is the correction term + * for f so that f+c ~ sqrt(z). + * For x<-0.5 + * acos(x) = pi - 2asin(sqrt((1-|x|)/2)) + * = pi - 0.5*(s+s*z*R(z)), where z=(1-|x|)/2,s=sqrt(z) + * + * Special cases: + * if x is NaN, return x itself; + * if |x|>1, return NaN with invalid signal. + * + * Function needed: sqrt + */ + +#include "fdlibm.h" +#include +#include + +#if defined(__STDC__) || defined(__cplusplus) +static const double +#else +static double +#endif + one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ + pi = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */ + pio2_hi = 1.57079632679489655800e+00, /* 0x3FF921FB, 0x54442D18 */ + pio2_lo = 6.12323399573676603587e-17, /* 0x3C91A626, 0x33145C07 */ + pS0 = 1.66666666666666657415e-01, /* 0x3FC55555, 0x55555555 */ + pS1 = -3.25565818622400915405e-01, /* 0xBFD4D612, 0x03EB6F7D */ + pS2 = 2.01212532134862925881e-01, /* 0x3FC9C155, 0x0E884455 */ + pS3 = -4.00555345006794114027e-02, /* 0xBFA48228, 0xB5688F3B */ + pS4 = 7.91534994289814532176e-04, /* 0x3F49EFE0, 0x7501B288 */ + pS5 = 3.47933107596021167570e-05, /* 0x3F023DE1, 0x0DFDF709 */ + qS1 = -2.40339491173441421878e+00, /* 0xC0033A27, 0x1C8A2D4B */ + qS2 = 2.02094576023350569471e+00, /* 0x40002AE5, 0x9C598AC8 */ + qS3 = -6.88283971605453293030e-01, /* 0xBFE6066C, 0x1B8D0159 */ + qS4 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */ + +#if defined(__STDC__) || defined(__cplusplus) +double __ieee754_acos(double x) +#else +double __ieee754_acos(x) +double x; +#endif +{ + double z, p, q, r, w, s, c, df; + int hx, ix; + hx = __HI(x); + ix = hx & 0x7fffffff; + if (ix >= 0x3ff00000) { /* |x| >= 1 */ + if (((ix - 0x3ff00000) | __LO(x)) == 0) { /* |x|==1 */ + if (hx > 0) + return 0.0; /* acos(1) = 0 */ + else + return pi + 2.0 * pio2_lo; /* acos(-1)= pi */ + } + errno = EDOM; + return NAN; + // return (x - x) / (x - x); /* acos(|x|>1) is NaN */ + } + if (ix < 0x3fe00000) { /* |x| < 0.5 */ + if (ix <= 0x3c600000) + return pio2_hi + pio2_lo; /*if|x|<2**-57*/ + z = x * x; + p = z * (pS0 + z * (pS1 + z * (pS2 + z * (pS3 + z * (pS4 + z * pS5))))); + q = one + z * (qS1 + z * (qS2 + z * (qS3 + z * qS4))); + r = p / q; + return pio2_hi - (x - (pio2_lo - x * r)); + } else if (hx < 0) { /* x < -0.5 */ + z = (one + x) * 0.5; + p = z * (pS0 + z * (pS1 + z * (pS2 + z * (pS3 + z * (pS4 + z * pS5))))); + q = one + z * (qS1 + z * (qS2 + z * (qS3 + z * qS4))); + s = sqrt(z); + r = p / q; + w = r * s - pio2_lo; + return pi - 2.0 * (s + w); + } else { /* x > 0.5 */ + z = (one - x) * 0.5; + s = sqrt(z); + df = s; + __LO(df) = 0; + c = (z - df * df) / (s + df); + p = z * (pS0 + z * (pS1 + z * (pS2 + z * (pS3 + z * (pS4 + z * pS5))))); + q = one + z * (qS1 + z * (qS2 + z * (qS3 + z * qS4))); + r = p / q; + w = r * s + c; + return 2.0 * (df + w); + } +} diff --git a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_asin.c b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_asin.c new file mode 100644 index 00000000..13e1b59e --- /dev/null +++ b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_asin.c @@ -0,0 +1,120 @@ + +/* @(#)e_asin.c 1.3 95/01/18 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* __ieee754_asin(x) + * Method : + * Since asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ... + * we approximate asin(x) on [0,0.5] by + * asin(x) = x + x*x^2*R(x^2) + * where + * R(x^2) is a rational approximation of (asin(x)-x)/x^3 + * and its remez error is bounded by + * |(asin(x)-x)/x^3 - R(x^2)| < 2^(-58.75) + * + * For x in [0.5,1] + * asin(x) = pi/2-2*asin(sqrt((1-x)/2)) + * Let y = (1-x), z = y/2, s := sqrt(z), and pio2_hi+pio2_lo=pi/2; + * then for x>0.98 + * asin(x) = pi/2 - 2*(s+s*z*R(z)) + * = pio2_hi - (2*(s+s*z*R(z)) - pio2_lo) + * For x<=0.98, let pio4_hi = pio2_hi/2, then + * f = hi part of s; + * c = sqrt(z) - f = (z-f*f)/(s+f) ...f+c=sqrt(z) + * and + * asin(x) = pi/2 - 2*(s+s*z*R(z)) + * = pio4_hi+(pio4-2s)-(2s*z*R(z)-pio2_lo) + * = pio4_hi+(pio4-2f)-(2s*z*R(z)-(pio2_lo+2c)) + * + * Special cases: + * if x is NaN, return x itself; + * if |x|>1, return NaN with invalid signal. + * + */ + +#include "fdlibm.h" +#include +#include + +#if defined(__STDC__) || defined(__cplusplus) +static const double +#else +static double +#endif + one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ + huge = 1.000e+300, pio2_hi = 1.57079632679489655800e+00, /* 0x3FF921FB, 0x54442D18 */ + pio2_lo = 6.12323399573676603587e-17, /* 0x3C91A626, 0x33145C07 */ + pio4_hi = 7.85398163397448278999e-01, /* 0x3FE921FB, 0x54442D18 */ + /* coefficient for R(x^2) */ + pS0 = 1.66666666666666657415e-01, /* 0x3FC55555, 0x55555555 */ + pS1 = -3.25565818622400915405e-01, /* 0xBFD4D612, 0x03EB6F7D */ + pS2 = 2.01212532134862925881e-01, /* 0x3FC9C155, 0x0E884455 */ + pS3 = -4.00555345006794114027e-02, /* 0xBFA48228, 0xB5688F3B */ + pS4 = 7.91534994289814532176e-04, /* 0x3F49EFE0, 0x7501B288 */ + pS5 = 3.47933107596021167570e-05, /* 0x3F023DE1, 0x0DFDF709 */ + qS1 = -2.40339491173441421878e+00, /* 0xC0033A27, 0x1C8A2D4B */ + qS2 = 2.02094576023350569471e+00, /* 0x40002AE5, 0x9C598AC8 */ + qS3 = -6.88283971605453293030e-01, /* 0xBFE6066C, 0x1B8D0159 */ + qS4 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */ + +#if defined(__STDC__) || defined(__cplusplus) +double __ieee754_asin(double x) +#else +double __ieee754_asin(x) +double x; +#endif +{ + double t, w, p, q, c, r, s; + int hx, ix; + hx = __HI(x); + ix = hx & 0x7fffffff; + if (ix >= 0x3ff00000) { /* |x|>= 1 */ + if (((ix - 0x3ff00000) | __LO(x)) == 0) + /* asin(1)=+-pi/2 with inexact */ + return x * pio2_hi + x * pio2_lo; + errno = EDOM; + return NAN; + // return (x - x) / (x - x); /* asin(|x|>1) is NaN */ + } else if (ix < 0x3fe00000) { /* |x|<0.5 */ + if (ix < 0x3e400000) { /* if |x| < 2**-27 */ + if (huge + x > one) + return x; /* return x with inexact if x!=0*/ + } else + t = x * x; + p = t * (pS0 + t * (pS1 + t * (pS2 + t * (pS3 + t * (pS4 + t * pS5))))); + q = one + t * (qS1 + t * (qS2 + t * (qS3 + t * qS4))); + w = p / q; + return x + x * w; + } + /* 1> |x|>= 0.5 */ + w = one - __fabs(x); + t = w * 0.5; + p = t * (pS0 + t * (pS1 + t * (pS2 + t * (pS3 + t * (pS4 + t * pS5))))); + q = one + t * (qS1 + t * (qS2 + t * (qS3 + t * qS4))); + s = sqrt(t); + if (ix >= 0x3FEF3333) { /* if |x| > 0.975 */ + w = p / q; + t = pio2_hi - (2.0 * (s + s * w) - pio2_lo); + } else { + w = s; + __LO(w) = 0; + c = (t - w * w) / (s + w); + r = p / q; + p = 2.0 * s * r - (pio2_lo - 2.0 * c); + q = pio4_hi - 2.0 * w; + t = pio4_hi - (p - q); + } + if (hx > 0) + return t; + else + return -t; +} diff --git a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_atan2.c b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_atan2.c new file mode 100644 index 00000000..7501857d --- /dev/null +++ b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_atan2.c @@ -0,0 +1,144 @@ + +/* @(#)e_atan2.c 1.3 95/01/18 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + * + */ + +/* __ieee754_atan2(y,x) + * Method : + * 1. Reduce y to positive by atan2(y,x)=-atan2(-y,x). + * 2. Reduce x to positive by (if x and y are unexceptional): + * ARG (x+iy) = arctan(y/x) ... if x > 0, + * ARG (x+iy) = pi - arctan[y/(-x)] ... if x < 0, + * + * Special cases: + * + * ATAN2((anything), NaN ) is NaN; + * ATAN2(NAN , (anything) ) is NaN; + * ATAN2(+-0, +(anything but NaN)) is +-0 ; + * ATAN2(+-0, -(anything but NaN)) is +-pi ; + * ATAN2(+-(anything but 0 and NaN), 0) is +-pi/2; + * ATAN2(+-(anything but INF and NaN), +INF) is +-0 ; + * ATAN2(+-(anything but INF and NaN), -INF) is +-pi; + * ATAN2(+-INF,+INF ) is +-pi/4 ; + * ATAN2(+-INF,-INF ) is +-3pi/4; + * ATAN2(+-INF, (anything but,0,NaN, and INF)) is +-pi/2; + * + * Constants: + * The hexadecimal values are the intended ones for the following + * constants. The decimal values may be used, provided that the + * compiler will convert from decimal to binary accurately enough + * to produce the hexadecimal values shown. + */ + +#include "fdlibm.h" + +#if defined(__STDC__) || defined(__cplusplus) +static const double +#else +static double +#endif + tiny = 1.0e-300, + zero = 0.0, pi_o_4 = 7.8539816339744827900E-01, /* 0x3FE921FB, 0x54442D18 */ + pi_o_2 = 1.5707963267948965580E+00, /* 0x3FF921FB, 0x54442D18 */ + pi = 3.1415926535897931160E+00, /* 0x400921FB, 0x54442D18 */ + pi_lo = 1.2246467991473531772E-16; /* 0x3CA1A626, 0x33145C07 */ + +#if defined(__STDC__) || defined(__cplusplus) +double __ieee754_atan2(double y, double x) +#else +double __ieee754_atan2(y, x) +double y, x; +#endif +{ + double z; + int k, m, hx, hy, ix, iy; + unsigned lx, ly; + + hx = __HI(x); + ix = hx & 0x7fffffff; + lx = __LO(x); + hy = __HI(y); + iy = hy & 0x7fffffff; + ly = __LO(y); + if (((ix | ((lx | -lx) >> 31)) > 0x7ff00000) + || ((iy | ((ly | -ly) >> 31)) > 0x7ff00000)) /* x or y is NaN */ + return x + y; + if ((hx - 0x3ff00000 | lx) == 0) + return atan(y); /* x=1.0 */ + m = ((hy >> 31) & 1) | ((hx >> 30) & 2); /* 2*sign(x)+sign(y) */ + + /* when y = 0 */ + if ((iy | ly) == 0) { + switch (m) { + case 0: + case 1: + return y; /* atan(+-0,+anything)=+-0 */ + case 2: + return pi + tiny; /* atan(+0,-anything) = pi */ + case 3: + return -pi - tiny; /* atan(-0,-anything) =-pi */ + } + } + /* when x = 0 */ + if ((ix | lx) == 0) + return (hy < 0) ? -pi_o_2 - tiny : pi_o_2 + tiny; + + /* when x is INF */ + if (ix == 0x7ff00000) { + if (iy == 0x7ff00000) { + switch (m) { + case 0: + return pi_o_4 + tiny; /* atan(+INF,+INF) */ + case 1: + return -pi_o_4 - tiny; /* atan(-INF,+INF) */ + case 2: + return 3.0 * pi_o_4 + tiny; /*atan(+INF,-INF)*/ + case 3: + return -3.0 * pi_o_4 - tiny; /*atan(-INF,-INF)*/ + } + } else { + switch (m) { + case 0: + return zero; /* atan(+...,+INF) */ + case 1: + return -zero; /* atan(-...,+INF) */ + case 2: + return pi + tiny; /* atan(+...,-INF) */ + case 3: + return -pi - tiny; /* atan(-...,-INF) */ + } + } + } + /* when y is INF */ + if (iy == 0x7ff00000) + return (hy < 0) ? -pi_o_2 - tiny : pi_o_2 + tiny; + + /* compute y/x */ + k = (iy - ix) >> 20; + if (k > 60) + z = pi_o_2 + 0.5 * pi_lo; /* |y/x| > 2**60 */ + else if (hx < 0 && k < -60) + z = 0.0; /* |y|/x < -2**60 */ + else + z = atan(__fabs(y / x)); /* safe to do y/x */ + switch (m) { + case 0: + return z; /* atan(+,+) */ + case 1: + __HI(z) ^= 0x80000000; + return z; /* atan(-,+) */ + case 2: + return pi - (z - pi_lo); /* atan(+,-) */ + default: /* case 3 */ + return (z - pi_lo) - pi; /* atan(-,-) */ + } +} diff --git a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_exp.c b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_exp.c new file mode 100644 index 00000000..ab2d7fcb --- /dev/null +++ b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_exp.c @@ -0,0 +1,162 @@ + +/* @(#)e_exp.c 1.6 04/04/22 */ +/* + * ==================================================== + * Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved. + * + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* __ieee754_exp(x) + * Returns the exponential of x. + * + * Method + * 1. Argument reduction: + * Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658. + * Given x, find r and integer k such that + * + * x = k*ln2 + r, |r| <= 0.5*ln2. + * + * Here r will be represented as r = hi-lo for better + * accuracy. + * + * 2. Approximation of exp(r) by a special rational function on + * the interval [0,0.34658]: + * Write + * R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ... + * We use a special Remes algorithm on [0,0.34658] to generate + * a polynomial of degree 5 to approximate R. The maximum error + * of this polynomial approximation is bounded by 2**-59. In + * other words, + * R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5 + * (where z=r*r, and the values of P1 to P5 are listed below) + * and + * | 5 | -59 + * | 2.0+P1*z+...+P5*z - R(z) | <= 2 + * | | + * The computation of exp(r) thus becomes + * 2*r + * exp(r) = 1 + ------- + * R - r + * r*R1(r) + * = 1 + r + ----------- (for better accuracy) + * 2 - R1(r) + * where + * 2 4 10 + * R1(r) = r - (P1*r + P2*r + ... + P5*r ). + * + * 3. Scale back to obtain exp(x): + * From step 1, we have + * exp(x) = 2^k * exp(r) + * + * Special cases: + * exp(INF) is INF, exp(NaN) is NaN; + * exp(-INF) is 0, and + * for finite argument, only exp(0)=1 is exact. + * + * Accuracy: + * according to an error analysis, the error is always less than + * 1 ulp (unit in the last place). + * + * Misc. info. + * For IEEE double + * if x > 7.09782712893383973096e+02 then exp(x) overflow + * if x < -7.45133219101941108420e+02 then exp(x) underflow + * + * Constants: + * The hexadecimal values are the intended ones for the following + * constants. The decimal values may be used, provided that the + * compiler will convert from decimal to binary accurately enough + * to produce the hexadecimal values shown. + */ + +#include "fdlibm.h" + +#if defined(__STDC__) || defined(__cplusplus) +static const double +#else +static double +#endif +one = 1.0, +halF[2] = {0.5,-0.5,}, +huge = 1.0e+300, +twom1000= 9.33263618503218878990e-302, /* 2**-1000=0x01700000,0*/ +o_threshold= 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */ +u_threshold= -7.45133219101941108420e+02, /* 0xc0874910, 0xD52D3051 */ +ln2HI[2] ={ 6.93147180369123816490e-01, /* 0x3fe62e42, 0xfee00000 */ + -6.93147180369123816490e-01,},/* 0xbfe62e42, 0xfee00000 */ +ln2LO[2] ={ 1.90821492927058770002e-10, /* 0x3dea39ef, 0x35793c76 */ + -1.90821492927058770002e-10,},/* 0xbdea39ef, 0x35793c76 */ +invln2 = 1.44269504088896338700e+00, /* 0x3ff71547, 0x652b82fe */ +P1 = 1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */ +P2 = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */ +P3 = 6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */ +P4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */ +P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */ + +#if defined(__STDC__) || defined(__cplusplus) +double __ieee754_exp(double x) /* default IEEE double exp */ +#else +double __ieee754_exp(x) /* default IEEE double exp */ +double x; +#endif +{ + double y, hi, lo, c, t; + int k, xsb; + unsigned hx; + + hx = __HI(x); /* high word of x */ + xsb = (hx >> 31) & 1; /* sign bit of x */ + hx &= 0x7fffffff; /* high word of |x| */ + + /* filter out non-finite argument */ + if (hx >= 0x40862E42) { /* if |x|>=709.78... */ + if (hx >= 0x7ff00000) { + if (((hx & 0xfffff) | __LO(x)) != 0) + return x + x; /* NaN */ + else + return (xsb == 0) ? x : 0.0; /* exp(+-inf)={inf,0} */ + } + if (x > o_threshold) + return huge * huge; /* overflow */ + if (x < u_threshold) + return twom1000 * twom1000; /* underflow */ + } + + /* argument reduction */ + if (hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */ + if (hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */ + hi = x - ln2HI[xsb]; + lo = ln2LO[xsb]; + k = 1 - xsb - xsb; + } else { + k = (int)(invln2 * x + halF[xsb]); + t = k; + hi = x - t * ln2HI[0]; /* t*ln2HI is exact here */ + lo = t * ln2LO[0]; + } + x = hi - lo; + } else if (hx < 0x3e300000) { /* when |x|<2**-28 */ + if (huge + x > one) + return one + x; /* trigger inexact */ + } else + k = 0; + + /* x is now in primary range */ + t = x * x; + c = x - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5)))); + if (k == 0) + return one - ((x * c) / (c - 2.0) - x); + else + y = one - ((lo - (x * c) / (2.0 - c)) - hi); + if (k >= -1021) { + __HI(y) += (k << 20); /* add k to y's exponent */ + return y; + } else { + __HI(y) += ((k + 1000) << 20); /* add k to y's exponent */ + return y * twom1000; + } +} diff --git a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_fmod.c b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_fmod.c new file mode 100644 index 00000000..59348156 --- /dev/null +++ b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_fmod.c @@ -0,0 +1,178 @@ + +/* @(#)e_fmod.c 1.3 95/01/18 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* + * __ieee754_fmod(x,y) + * Return x mod y in exact arithmetic + * Method: shift and subtract + */ + +#include "fdlibm.h" + +#if defined(__STDC__) || defined(__cplusplus) +static const double one = 1.0, + Zero[] = { + 0.0, + -0.0, + }; +#else +static double one = 1.0, + Zero[] = { + 0.0, + -0.0, + }; +#endif + +#if defined(__STDC__) || defined(__cplusplus) +double __ieee754_fmod(double x, double y) +#else +double __ieee754_fmod(x, y) +double x, y; +#endif +{ + int n, hx, hy, hz, ix, iy, sx, i, unk1; + unsigned lx, ly, lz; + + hx = __HI(x); /* high word of x */ + lx = __LO(x); /* low word of x */ + hy = __HI(y); /* high word of y */ + ly = __LO(y); /* low word of y */ + sx = hx & 0x80000000; /* sign of x */ + hx ^= sx; /* |x| */ + hy &= 0x7fffffff; /* |y| */ + + /* purge off exception values */ + if ((hy | ly) == 0 || (hx >= 0x7ff00000) || /* y=0,or x not finite */ + ((hy | ((ly | -ly) >> 31)) > 0x7ff00000)) /* or y is NaN */ + return (x * y) / (x * y); + if (hx <= hy) { + if ((hx < hy) || (lx < ly)) + return x; /* |x|<|y| return x */ + if (lx == ly) + return Zero[(unsigned)sx >> 31]; /* |x|=|y| return x*0*/ + } + + /* determine ix = ilogb(x) */ + if (hx < 0x00100000) { /* subnormal x */ + if (hx == 0) { + for (ix = -1043, i = lx; i > 0; i <<= 1) + ix -= 1; + } else { + for (ix = -1022, i = (hx << 11); i > 0; i <<= 1) + ix -= 1; + } + } else + ix = (hx >> 20) - 1023; + + /* determine iy = ilogb(y) */ + if (hy < 0x00100000) { /* subnormal y */ + if (hy == 0) { + for (iy = -1043, i = ly; i > 0; i <<= 1) + iy -= 1; + } else { + for (iy = -1022, i = (hy << 11); i > 0; i <<= 1) + iy -= 1; + } + } else + iy = (hy >> 20) - 1023; + + /* set up {hx,lx}, {hy,ly} and align y to x */ + if (ix >= -1022) + hx = 0x00100000 | (0x000fffff & hx); + else { /* subnormal x, shift x to normal */ + n = -1022 - ix; + if (n <= 31) { + hx = (hx << n) | (lx >> (32 - n)); + lx <<= n; + } else { + hx = lx << (n - 32); + lx = 0; + } + } + if (iy >= -1022) + hy = 0x00100000 | (0x000fffff & hy); + else { /* subnormal y, shift y to normal */ + n = -1022 - iy; + if (n <= 31) { + hy = (hy << n) | (ly >> (32 - n)); + ly <<= n; + } else { + hy = ly << (n - 32); + ly = 0; + } + } + + /* fix point fmod */ + n = ix - iy; + unk1 = n + 2; + while (n--) { + hz = hx - hy; + lz = lx - ly; + + if (hz == 0 && (lx >> unk1 == ly >> unk1)) { + return Zero[(unsigned)sx >> 31]; + } + if (lx < ly) + hz -= 1; + if (hz < 0) { + hx = hx + hx + (lx >> 31); + lx = lx + lx; + } else { + if ((hz | lz) == 0) /* return sign(x)*0 */ + return Zero[(unsigned)sx >> 31]; + hx = hz + hz + (lz >> 31); + lx = lz + lz; + } + } + hz = hx - hy; + lz = lx - ly; + if (hz == 0 && (lx >> unk1 == ly >> unk1)) { + return Zero[(unsigned)sx >> 31]; + } + if (lx < ly) + hz -= 1; + if (hz >= 0) { + hx = hz; + lx = lz; + } + + /* convert back to floating value and restore the sign */ + if ((hx | lx) == 0) /* return sign(x)*0 */ + return Zero[(unsigned)sx >> 31]; + while (hx < 0x00100000) { /* normalize x */ + hx = hx + hx + (lx >> 31); + lx = lx + lx; + iy -= 1; + } + if (iy >= -1022) { /* normalize output */ + hx = ((hx - 0x00100000) | ((iy + 1023) << 20)); + __HI(x) = hx | sx; + __LO(x) = lx; + } else { /* subnormal output */ + n = -1022 - iy; + if (n <= 20) { + lx = (lx >> n) | ((unsigned)hx << (32 - n)); + hx >>= n; + } else if (n <= 31) { + lx = (hx << (32 - n)) | (lx >> n); + hx = sx; + } else { + lx = hx >> (n - 32); + hx = sx; + } + __HI(x) = hx | sx; + __LO(x) = lx; + x *= one; /* create necessary signal */ + } + return x; /* exact output */ +} diff --git a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_log.c b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_log.c new file mode 100644 index 00000000..ca3413d3 --- /dev/null +++ b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_log.c @@ -0,0 +1,158 @@ + +/* @(#)e_log.c 1.3 95/01/18 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* __ieee754_log(x) + * Return the logrithm of x + * + * Method : + * 1. Argument Reduction: find k and f such that + * x = 2^k * (1+f), + * where sqrt(2)/2 < 1+f < sqrt(2) . + * + * 2. Approximation of log(1+f). + * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s) + * = 2s + 2/3 s**3 + 2/5 s**5 + ....., + * = 2s + s*R + * We use a special Reme algorithm on [0,0.1716] to generate + * a polynomial of degree 14 to approximate R The maximum error + * of this polynomial approximation is bounded by 2**-58.45. In + * other words, + * 2 4 6 8 10 12 14 + * R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s +Lg6*s +Lg7*s + * (the values of Lg1 to Lg7 are listed in the program) + * and + * | 2 14 | -58.45 + * | Lg1*s +...+Lg7*s - R(z) | <= 2 + * | | + * Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2. + * In order to guarantee error in log below 1ulp, we compute log + * by + * log(1+f) = f - s*(f - R) (if f is not too large) + * log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy) + * + * 3. Finally, log(x) = k*ln2 + log(1+f). + * = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo))) + * Here ln2 is split into two floating point number: + * ln2_hi + ln2_lo, + * where n*ln2_hi is always exact for |n| < 2000. + * + * Special cases: + * log(x) is NaN with signal if x < 0 (including -INF) ; + * log(+INF) is +INF; log(0) is -INF with signal; + * log(NaN) is that NaN with no signal. + * + * Accuracy: + * according to an error analysis, the error is always less than + * 1 ulp (unit in the last place). + * + * Constants: + * The hexadecimal values are the intended ones for the following + * constants. The decimal values may be used, provided that the + * compiler will convert from decimal to binary accurately enough + * to produce the hexadecimal values shown. + */ + +#include "fdlibm.h" +#include + +#if defined(__STDC__) || defined(__cplusplus) +static const double +#else +static double +#endif + ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */ + ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */ + two54 = 1.80143985094819840000e+16, /* 43500000 00000000 */ + Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */ + Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */ + Lg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */ + Lg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */ + Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */ + Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */ + Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */ + +static double zero = 0.0; + +#if defined(__STDC__) || defined(__cplusplus) +double __ieee754_log(double x) +#else +double __ieee754_log(x) +double x; +#endif +{ + double hfsq, f, s, z, R, w, t1, t2, dk; + int k, hx, i, j; + unsigned lx; + + hx = __HI(x); /* high word of x */ + lx = __LO(x); /* low word of x */ + + k = 0; + if (hx < 0x00100000) { /* x < 2**-1022 */ + if (((hx & 0x7fffffff) | lx) == 0) + return -two54 / zero; /* log(+-0)=-inf */ + if (hx < 0) { + errno = EDOM; + return (x - x) / zero; /* log(-#) = NaN */ + } + k -= 54; + x *= two54; /* subnormal number, scale up x */ + hx = __HI(x); /* high word of x */ + } + if (hx >= 0x7ff00000) + return x + x; + k += (hx >> 20) - 1023; + hx &= 0x000fffff; + i = (hx + 0x95f64) & 0x100000; + __HI(x) = hx | (i ^ 0x3ff00000); /* normalize x or x/2 */ + k += (i >> 20); + f = x - 1.0; + if ((0x000fffff & (2 + hx)) < 3) { /* |f| < 2**-20 */ + if (f == zero) + if (k == 0) + return zero; + else { + dk = (double)k; + return dk * ln2_hi + dk * ln2_lo; + } + R = f * f * (0.5 - 0.33333333333333333 * f); + if (k == 0) + return f - R; + else { + dk = (double)k; + return dk * ln2_hi - ((R - dk * ln2_lo) - f); + } + } + s = f / (2.0 + f); + dk = (double)k; + z = s * s; + i = hx - 0x6147a; + w = z * z; + j = 0x6b851 - hx; + t1 = w * (Lg2 + w * (Lg4 + w * Lg6)); + t2 = z * (Lg1 + w * (Lg3 + w * (Lg5 + w * Lg7))); + i |= j; + R = t2 + t1; + if (i > 0) { + hfsq = 0.5 * f * f; + if (k == 0) + return f - (hfsq - s * (hfsq + R)); + else + return dk * ln2_hi - ((hfsq - (s * (hfsq + R) + dk * ln2_lo)) - f); + } else { + if (k == 0) + return f - s * (f - R); + else + return dk * ln2_hi - ((s * (f - R) - dk * ln2_lo) - f); + } +} diff --git a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_log10.c b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_log10.c new file mode 100644 index 00000000..a72afa45 --- /dev/null +++ b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_log10.c @@ -0,0 +1,99 @@ + +/* @(#)e_log10.c 1.3 95/01/18 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* __ieee754_log10(x) + * Return the base 10 logarithm of x + * + * Method : + * Let log10_2hi = leading 40 bits of log10(2) and + * log10_2lo = log10(2) - log10_2hi, + * ivln10 = 1/log(10) rounded. + * Then + * n = ilogb(x), + * if(n<0) n = n+1; + * x = scalbn(x,-n); + * log10(x) := n*log10_2hi + (n*log10_2lo + ivln10*log(x)) + * + * Note 1: + * To guarantee log10(10**n)=n, where 10**n is normal, the rounding + * mode must set to Round-to-Nearest. + * Note 2: + * [1/log(10)] rounded to 53 bits has error .198 ulps; + * log10 is monotonic at all binary break points. + * + * Special cases: + * log10(x) is NaN with signal if x < 0; + * log10(+INF) is +INF with no signal; log10(0) is -INF with signal; + * log10(NaN) is that NaN with no signal; + * log10(10**N) = N for N=0,1,...,22. + * + * Constants: + * The hexadecimal values are the intended ones for the following constants. + * The decimal values may be used, provided that the compiler will convert + * from decimal to binary accurately enough to produce the hexadecimal values + * shown. + */ + +#include "fdlibm.h" +#include + +#if defined(__STDC__) || defined(__cplusplus) +static const double +#else +static double +#endif + two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */ + ivln10 = 4.34294481903251816668e-01, /* 0x3FDBCB7B, 0x1526E50E */ + log10_2hi = 3.01029995663611771306e-01, /* 0x3FD34413, 0x509F6000 */ + log10_2lo = 3.69423907715893078616e-13; /* 0x3D59FEF3, 0x11F12B36 */ + +static double zero = 0.0; + +#if defined(__STDC__) || defined(__cplusplus) +double __ieee754_log10(double x) +#else +double __ieee754_log10(x) +double x; +#endif +{ + double y, z; + int i, k, hx; + unsigned lx; + + hx = __HI(x); /* high word of x */ + lx = __LO(x); /* low word of x */ + + k = 0; + if (hx < 0x00100000) { /* x < 2**-1022 */ + if (((hx & 0x7fffffff) | lx) == 0) { + errno = EDOM; + return -two54 / zero; /* log(+-0)=-inf */ + } + if (hx < 0) { + errno = EDOM; + return (x - x) / zero; /* log(-#) = NaN */ + } + k -= 54; + x *= two54; /* subnormal number, scale up x */ + hx = __HI(x); /* high word of x */ + } + if (hx >= 0x7ff00000) + return x + x; + k += (hx >> 20) - 1023; + i = ((unsigned)k & 0x80000000) >> 31; + hx = (hx & 0x000fffff) | ((0x3ff - i) << 20); + y = (double)(k + i); + __HI(x) = hx; + z = y * log10_2lo + ivln10 * __ieee754_log(x); + return z + y * log10_2hi; +} diff --git a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_pow.c b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_pow.c new file mode 100644 index 00000000..d72815bf --- /dev/null +++ b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_pow.c @@ -0,0 +1,343 @@ + +#ifndef lint +static char sccsid[] = "@(#)e_pow.c 1.5 04/04/22 SMI"; +#endif + +/* + * ==================================================== + * Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved. + * + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* __ieee754_pow(x,y) return x**y + * + * n + * Method: Let x = 2 * (1+f) + * 1. Compute and return log2(x) in two pieces: + * log2(x) = w1 + w2, + * where w1 has 53-24 = 29 bit trailing zeros. + * 2. Perform y*log2(x) = n+y' by simulating muti-precision + * arithmetic, where |y'|<=0.5. + * 3. Return x**y = 2**n*exp(y'*log2) + * + * Special cases: + * 1. (anything) ** 0 is 1 + * 2. (anything) ** 1 is itself + * 3. (anything) ** NAN is NAN + * 4. NAN ** (anything except 0) is NAN + * 5. +-(|x| > 1) ** +INF is +INF + * 6. +-(|x| > 1) ** -INF is +0 + * 7. +-(|x| < 1) ** +INF is +0 + * 8. +-(|x| < 1) ** -INF is +INF + * 9. +-1 ** +-INF is NAN + * 10. +0 ** (+anything except 0, NAN) is +0 + * 11. -0 ** (+anything except 0, NAN, odd integer) is +0 + * 12. +0 ** (-anything except 0, NAN) is +INF + * 13. -0 ** (-anything except 0, NAN, odd integer) is +INF + * 14. -0 ** (odd integer) = -( +0 ** (odd integer) ) + * 15. +INF ** (+anything except 0,NAN) is +INF + * 16. +INF ** (-anything except 0,NAN) is +0 + * 17. -INF ** (anything) = -0 ** (-anything) + * 18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer) + * 19. (-anything except 0 and inf) ** (non-integer) is NAN + * + * Accuracy: + * pow(x,y) returns x**y nearly rounded. In particular + * pow(integer,integer) + * always returns the correct integer provided it is + * representable. + * + * Constants : + * The hexadecimal values are the intended ones for the following + * constants. The decimal values may be used, provided that the + * compiler will convert from decimal to binary accurately enough + * to produce the hexadecimal values shown. + */ + +#include "fdlibm.h" +#include +#include + +#if defined(__STDC__) || defined(__cplusplus) +static const double +#else +static double +#endif +bp[] = {1.0, 1.5,}, +dp_h[] = { 0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */ +dp_l[] = { 0.0, 1.35003920212974897128e-08,}, /* 0x3E4CFDEB, 0x43CFD006 */ +zero = 0.0, +one = 1.0, +two = 2.0, +two53 = 9007199254740992.0, /* 0x43400000, 0x00000000 */ +huge = 1.0e300, +tiny = 1.0e-300, + /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */ +L1 = 5.99999999999994648725e-01, /* 0x3FE33333, 0x33333303 */ +L2 = 4.28571428578550184252e-01, /* 0x3FDB6DB6, 0xDB6FABFF */ +L3 = 3.33333329818377432918e-01, /* 0x3FD55555, 0x518F264D */ +L4 = 2.72728123808534006489e-01, /* 0x3FD17460, 0xA91D4101 */ +L5 = 2.30660745775561754067e-01, /* 0x3FCD864A, 0x93C9DB65 */ +L6 = 2.06975017800338417784e-01, /* 0x3FCA7E28, 0x4A454EEF */ +P1 = 1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */ +P2 = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */ +P3 = 6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */ +P4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */ +P5 = 4.13813679705723846039e-08, /* 0x3E663769, 0x72BEA4D0 */ +lg2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */ +lg2_h = 6.93147182464599609375e-01, /* 0x3FE62E43, 0x00000000 */ +lg2_l = -1.90465429995776804525e-09, /* 0xBE205C61, 0x0CA86C39 */ +ovt = 8.0085662595372944372e-0017, /* -(1024-log2(ovfl+.5ulp)) */ +cp = 9.61796693925975554329e-01, /* 0x3FEEC709, 0xDC3A03FD =2/(3ln2) */ +cp_h = 9.61796700954437255859e-01, /* 0x3FEEC709, 0xE0000000 =(float)cp */ +cp_l = -7.02846165095275826516e-09, /* 0xBE3E2FE0, 0x145B01F5 =tail of cp_h*/ +ivln2 = 1.44269504088896338700e+00, /* 0x3FF71547, 0x652B82FE =1/ln2 */ +ivln2_h = 1.44269502162933349609e+00, /* 0x3FF71547, 0x60000000 =24b 1/ln2*/ +ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/ + +#if defined(__STDC__) || defined(__cplusplus) +double __ieee754_pow(double x, double y) +#else +double __ieee754_pow(x, y) +double x, y; +#endif +{ + double z, ax, z_h, z_l, p_h, p_l; + double y1, t1, t2, r, s, t, u, v, w; + int i0, i1, i, j, k, yisint, n; + int hx, hy, ix, iy; + unsigned lx, ly; + + i0 = ((*(int *)&one) >> 29) ^ 1; + i1 = 1 - i0; + hx = __HI(x); + lx = __LO(x); + hy = __HI(y); + ly = __LO(y); + ix = hx & 0x7fffffff; + iy = hy & 0x7fffffff; + + /* y==zero: x**0 = 1 */ + if ((iy | ly) == 0) + return one; + + /* +-NaN return x+y */ + if (ix > 0x7ff00000 || ((ix == 0x7ff00000) && (lx != 0)) || iy > 0x7ff00000 + || ((iy == 0x7ff00000) && (ly != 0))) + return x + y; + + /* determine if y is an odd int when x < 0 + * yisint = 0 ... y is not an integer + * yisint = 1 ... y is an odd int + * yisint = 2 ... y is an even int + */ + yisint = 0; + if (hx < 0) { + if (iy >= 0x43400000) + yisint = 2; /* even integer y */ + else if (iy >= 0x3ff00000) { + k = (iy >> 20) - 0x3ff; /* exponent */ + if (k > 20) { + j = ly >> (52 - k); + if ((j << (52 - k)) == ly) + yisint = 2 - (j & 1); + } else if (ly == 0) { + j = iy >> (20 - k); + if ((j << (20 - k)) == iy) + yisint = 2 - (j & 1); + } + } + } + + /* special value of y */ + if (ly == 0) { + if (iy == 0x7ff00000) { /* y is +-inf */ + if (((ix - 0x3ff00000) | lx) == 0) + return y - y; /* inf**+-1 is NaN */ + else if (ix >= 0x3ff00000) /* (|x|>1)**+-inf = inf,0 */ + return (hy >= 0) ? y : zero; + else /* (|x|<1)**-,+inf = inf,0 */ + return (hy < 0) ? -y : zero; + } + if (iy == 0x3ff00000) { /* y is +-1 */ + if (hy < 0) + return one / x; + else + return x; + } + if (hy == 0x40000000) + return x * x; /* y is 2 */ + if (hy == 0x3fe00000) { /* y is 0.5 */ + if (hx >= 0) /* x >= +0 */ + return sqrt(x); + } + } + + ax = __fabs(x); + /* special value of x */ + if (lx == 0) { + if (ix == 0x7ff00000 || ix == 0 || ix == 0x3ff00000) { + z = ax; /*x is +-0,+-inf,+-1*/ + if (hy < 0) + z = one / z; /* z = (1/|x|) */ + if (hx < 0) { + if (((ix - 0x3ff00000) | yisint) == 0) { + z = (z - z) / (z - z); /* (-1)**non-int is NaN */ + } else if (yisint == 1) + z = -z; /* (x<0)**odd = -(|x|**odd) */ + } + return z; + } + } + + if (((((int)hx >> 31) + 1) | yisint) == 0) { + errno = EDOM; + return NAN; + }; + + /* |y| is huge */ + if (iy > 0x41e00000) { /* if |y| > 2**31 */ + if (iy > 0x43f00000) { /* if |y| > 2**64, must o/uflow */ + if (ix <= 0x3fefffff) + return (hy < 0) ? huge * huge : tiny * tiny; + if (ix >= 0x3ff00000) + return (hy > 0) ? huge * huge : tiny * tiny; + } + /* over/underflow if x is not close to one */ + if (ix < 0x3fefffff) + return (hy < 0) ? huge * huge : tiny * tiny; + if (ix > 0x3ff00000) + return (hy > 0) ? huge * huge : tiny * tiny; + /* now |1-x| is tiny <= 2**-20, suffice to compute + log(x) by x-x^2/2+x^3/3-x^4/4 */ + t = x - 1; /* t has 20 trailing zeros */ + w = (t * t) * (0.5 - t * (0.3333333333333333333333 - t * 0.25)); + u = ivln2_h * t; /* ivln2_h has 21 sig. bits */ + v = t * ivln2_l - w * ivln2; + t1 = u + v; + __LO(t1) = 0; + t2 = v - (t1 - u); + } else { + double s2, s_h, s_l, t_h, t_l; + n = 0; + /* take care subnormal number */ + if (ix < 0x00100000) { + ax *= two53; + n -= 53; + ix = __HI(ax); + } + n += ((ix) >> 20) - 0x3ff; + j = ix & 0x000fffff; + /* determine interval */ + ix = j | 0x3ff00000; /* normalize ix */ + if (j <= 0x3988E) + k = 0; /* |x|> 1) | 0x20000000) + 0x00080000 + (k << 18); + t_l = ax - (t_h - bp[k]); + s_l = v * ((u - s_h * t_h) - s_h * t_l); + /* compute log(ax) */ + s2 = s * s; + r = s2 * s2 * (L1 + s2 * (L2 + s2 * (L3 + s2 * (L4 + s2 * (L5 + s2 * L6))))); + r += s_l * (s_h + s); + s2 = s_h * s_h; + t_h = 3.0 + s2 + r; + __LO(t_h) = 0; + t_l = r - ((t_h - 3.0) - s2); + /* u+v = s*(1+...) */ + u = s_h * t_h; + v = s_l * t_h + t_l * s; + /* 2/(3log2)*(s+...) */ + p_h = u + v; + __LO(p_h) = 0; + p_l = v - (p_h - u); + z_h = cp_h * p_h; /* cp_h+cp_l = 2/(3*log2) */ + z_l = cp_l * p_h + p_l * cp + dp_l[k]; + /* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */ + t = (double)n; + t1 = (((z_h + z_l) + dp_h[k]) + t); + __LO(t1) = 0; + t2 = z_l - (((t1 - t) - dp_h[k]) - z_h); + } + + s = one; /* s (sign of result -ve**odd) = -1 else = 1 */ + if (((((int)hx >> 31) + 1) | (yisint - 1)) == 0) + s = -one; /* (-ve)**(odd int) */ + + /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */ + y1 = y; + __LO(y1) = 0; + p_l = (y - y1) * t1 + y * t2; + p_h = y1 * t1; + z = p_l + p_h; + j = __HI(z); + i = __LO(z); + if (j >= 0x40900000) { /* z >= 1024 */ + if (((j - 0x40900000) | i) != 0) /* if z > 1024 */ + return s * huge * huge; /* overflow */ + else { + if (p_l + ovt > z - p_h) + return s * huge * huge; /* overflow */ + } + } else if ((j & 0x7fffffff) >= 0x4090cc00) { /* z <= -1075 */ + if (((j - 0xc090cc00) | i) != 0) /* z < -1075 */ + return s * tiny * tiny; /* underflow */ + else { + if (p_l <= z - p_h) + return s * tiny * tiny; /* underflow */ + } + } + /* + * compute 2**(p_h+p_l) + */ + i = j & 0x7fffffff; + k = (i >> 20) - 0x3ff; + n = 0; + if (i > 0x3fe00000) { /* if |z| > 0.5, set n = [z+0.5] */ + n = j + (0x00100000 >> (k + 1)); + k = ((n & 0x7fffffff) >> 20) - 0x3ff; /* new k for n */ + t = zero; + __HI(t) = (n & ~(0x000fffff >> k)); + n = ((n & 0x000fffff) | 0x00100000) >> (20 - k); + if (j < 0) + n = -n; + p_h -= t; + } + t = p_l + p_h; + __LO(t) = 0; + u = t * lg2_h; + v = (p_l - (t - p_h)) * lg2 + t * lg2_l; + z = u + v; + w = v - (z - u); + t = z * z; + t1 = z - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5)))); + r = (z * t1) / (t1 - two) - (w + z * w); + z = one - (r - z); + j = __HI(z); + j += (n << 20); + if ((j >> 20) <= 0) + z = scalbn(z, n); /* subnormal output */ + else + __HI(z) += (n << 20); + return s * z; +} diff --git a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_rem_pio2.c b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_rem_pio2.c new file mode 100644 index 00000000..a7e48974 --- /dev/null +++ b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/e_rem_pio2.c @@ -0,0 +1,186 @@ + +/* @(#)e_rem_pio2.c 1.4 95/01/18 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + * + */ + +/* __ieee754_rem_pio2(x,y) + * + * return the remainder of x rem pi/2 in y[0]+y[1] + * use __kernel_rem_pio2() + */ + +#include "fdlibm.h" + +/* + * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi + */ +#if defined(__STDC__) || defined(__cplusplus) +static const int two_over_pi[] = { +#else +static int two_over_pi[] = { +#endif + 0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, 0x95993C, 0x439041, + 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A, 0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, + 0xFE1DEB, 0x1CB129, 0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41, + 0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8, 0x97FFDE, 0x05980F, + 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF, 0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, + 0x7527BA, 0xC7EBE5, 0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08, + 0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3, 0x91615E, 0xE61B08, + 0x659985, 0x5F14A0, 0x68408D, 0xFFD880, 0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, + 0x60E27B, 0xC08C6B, +}; + +#if defined(__STDC__) || defined(__cplusplus) +static const int npio2_hw[] = { +#else +static int npio2_hw[] = { +#endif + 0x3FF921FB, 0x400921FB, 0x4012D97C, 0x401921FB, 0x401F6A7A, 0x4022D97C, 0x4025FDBB, + 0x402921FB, 0x402C463A, 0x402F6A7A, 0x4031475C, 0x4032D97C, 0x40346B9C, 0x4035FDBB, + 0x40378FDB, 0x403921FB, 0x403AB41B, 0x403C463A, 0x403DD85A, 0x403F6A7A, 0x40407E4C, + 0x4041475C, 0x4042106C, 0x4042D97C, 0x4043A28C, 0x40446B9C, 0x404534AC, 0x4045FDBB, + 0x4046C6CB, 0x40478FDB, 0x404858EB, 0x404921FB, +}; + +/* + * invpio2: 53 bits of 2/pi + * pio2_1: first 33 bit of pi/2 + * pio2_1t: pi/2 - pio2_1 + * pio2_2: second 33 bit of pi/2 + * pio2_2t: pi/2 - (pio2_1+pio2_2) + * pio2_3: third 33 bit of pi/2 + * pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3) + */ + +#if defined(__STDC__) || defined(__cplusplus) +static const double +#else +static double +#endif + zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */ + half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */ + two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */ + invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */ + pio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */ + pio2_1t = 6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */ + pio2_2 = 6.07710050630396597660e-11, /* 0x3DD0B461, 0x1A600000 */ + pio2_2t = 2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */ + pio2_3 = 2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */ + pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */ + +#if defined(__STDC__) || defined(__cplusplus) +int __ieee754_rem_pio2(double x, double *y) +#else +int __ieee754_rem_pio2(x, y) +double x, y[]; +#endif +{ + double z, w, t, r, fn; + double tx[3]; + int e0, i, j, nx, n, ix, hx; + + hx = __HI(x); /* high word of x */ + ix = hx & 0x7fffffff; + if (ix <= 0x3fe921fb) /* |x| ~<= pi/4 , no need for reduction */ + { + y[0] = x; + y[1] = 0; + return 0; + } + if (ix < 0x4002d97c) { /* |x| < 3pi/4, special case with n=+-1 */ + if (hx > 0) { + z = x - pio2_1; + if (ix != 0x3ff921fb) { /* 33+53 bit pi is good enough */ + y[0] = z - pio2_1t; + y[1] = (z - y[0]) - pio2_1t; + } else { /* near pi/2, use 33+33+53 bit pi */ + z -= pio2_2; + y[0] = z - pio2_2t; + y[1] = (z - y[0]) - pio2_2t; + } + return 1; + } else { /* negative x */ + z = x + pio2_1; + if (ix != 0x3ff921fb) { /* 33+53 bit pi is good enough */ + y[0] = z + pio2_1t; + y[1] = (z - y[0]) + pio2_1t; + } else { /* near pi/2, use 33+33+53 bit pi */ + z += pio2_2; + y[0] = z + pio2_2t; + y[1] = (z - y[0]) + pio2_2t; + } + return -1; + } + } + if (ix <= 0x413921fb) { /* |x| ~<= 2^19*(pi/2), medium size */ + t = __fabs(x); + n = (int)(t * invpio2 + half); + fn = (double)n; + r = t - fn * pio2_1; + w = fn * pio2_1t; /* 1st round good to 85 bit */ + if (n < 32 && ix != npio2_hw[n - 1]) { + y[0] = r - w; /* quick check no cancellation */ + } else { + j = ix >> 20; + y[0] = r - w; + i = j - (((__HI(y[0])) >> 20) & 0x7ff); + if (i > 16) { /* 2nd iteration needed, good to 118 */ + t = r; + w = fn * pio2_2; + r = t - w; + w = fn * pio2_2t - ((t - r) - w); + y[0] = r - w; + i = j - (((__HI(y[0])) >> 20) & 0x7ff); + if (i > 49) { /* 3rd iteration need, 151 bits acc */ + t = r; /* will cover all possible cases */ + w = fn * pio2_3; + r = t - w; + w = fn * pio2_3t - ((t - r) - w); + y[0] = r - w; + } + } + } + y[1] = (r - y[0]) - w; + if (hx < 0) { + y[0] = -y[0]; + y[1] = -y[1]; + return -n; + } else + return n; + } + /* + * all other (large) arguments + */ + if (ix >= 0x7ff00000) { /* x is inf or NaN */ + y[0] = y[1] = x - x; + return 0; + } + /* set z = scalbn(|x|,ilogb(x)-23) */ + __LO(z) = __LO(x); + e0 = (ix >> 20) - 1046; /* e0 = ilogb(z)-23; */ + __HI(z) = ix - (e0 << 20); + for (i = 0; i < 2; i++) { + tx[i] = (double)((int)(z)); + z = (z - tx[i]) * two24; + } + tx[2] = z; + nx = 3; + while (tx[nx - 1] == zero) + nx--; /* skip zero term */ + n = __kernel_rem_pio2(tx, y, e0, nx, 2, two_over_pi); + if (hx < 0) { + y[0] = -y[0]; + y[1] = -y[1]; + return -n; + } + return n; +} diff --git a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/fdlibm.h b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/fdlibm.h similarity index 77% rename from src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/fdlibm.h rename to src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/fdlibm.h index 13b3587f..61304139 100644 --- a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/fdlibm.h +++ b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/fdlibm.h @@ -1,5 +1,3 @@ -#ifndef MSL_MATH_FDLIBM_H -#define MSL_MATH_FDLIBM_H /* @(#)fdlibm.h 1.5 04/04/22 */ /* @@ -12,32 +10,38 @@ * ==================================================== */ +#ifndef _FDLIBM_H +#define _FDLIBM_H + #ifdef __cplusplus extern "C" { #endif -/* #ifdef __STDC__ */ -#include -#include -/* #endif */ +/* Sometimes it's necessary to define __LITTLE_ENDIAN explicitly + but these catch some common cases. */ + +#if defined(i386) || defined(i486) || defined(intel) || defined(x86) || defined(i86pc) \ + || defined(__alpha) || defined(__osf__) +#define __LITTLE_ENDIAN +#endif -/*#if __option(little_endian) -#define __HIp(x) *(1 + (int*)x) -#define __LOp(x) *(int*)x +#ifdef __LITTLE_ENDIAN +#define __HI(x) *(1 + (int *)&x) +#define __LO(x) *(int *)&x +#define __HIp(x) *(1 + (int *)x) +#define __LOp(x) *(int *)x #else -#define __HIp(x) *(int*)x -#define __LOp(x) *(1 + (int*)x) -#endif*/ +#define __HI(x) *(int *)&x +#define __LO(x) *(1 + (int *)&x) +#define __HIp(x) *(int *)x +#define __LOp(x) *(1 + (int *)x) +#endif -/* TODO: should __STDC__ actually be defined? -#ifdef __STDC__ -*/ +#if defined(__STDC__) || defined(__cplusplus) #define __P(p) p -/* #else #define __P(p) () #endif -*/ /* * ANSI/POSIX @@ -46,12 +50,16 @@ extern "C" { extern int signgam; #define MAXFLOAT ((float)3.40282346638528860e+38) -#define M_PI ((float)3.14159274101257324219) -enum fdversion { fdlibm_ieee = -1, fdlibm_svid, fdlibm_xopen, fdlibm_posix }; +enum fdversion { + fdlibm_ieee = -1, + fdlibm_svid, + fdlibm_xopen, + fdlibm_posix +}; #define _LIB_VERSION_TYPE enum fdversion -#define _LIB_VERSION _fdlib_version +#define _LIB_VERSION _fdlib_version /* if global variable _LIB_VERSION is not desirable, one may * change the following to be a constant by: @@ -62,17 +70,17 @@ enum fdversion { fdlibm_ieee = -1, fdlibm_svid, fdlibm_xopen, fdlibm_posix }; */ extern _LIB_VERSION_TYPE _LIB_VERSION; -#define _IEEE_ fdlibm_ieee -#define _SVID_ fdlibm_svid +#define _IEEE_ fdlibm_ieee +#define _SVID_ fdlibm_svid #define _XOPEN_ fdlibm_xopen #define _POSIX_ fdlibm_posix struct exception { - int type; - char* name; - double arg1; - double arg2; - double retval; + int type; + char *name; + double arg1; + double arg2; + double retval; }; #define HUGE MAXFLOAT @@ -84,12 +92,12 @@ struct exception { #define X_TLOSS 1.41484755040568800000e+16 -#define DOMAIN 1 -#define SING 2 -#define OVERFLOW 3 +#define DOMAIN 1 +#define SING 2 +#define OVERFLOW 3 #define UNDERFLOW 4 -#define TLOSS 5 -#define PLOSS 6 +#define TLOSS 5 +#define PLOSS 6 /* * ANSI/POSIX @@ -107,12 +115,11 @@ extern double sinh __P((double)); extern double tanh __P((double)); extern double exp __P((double)); -extern double frexp __P((double, int*)); +extern double frexp __P((double, int *)); extern double ldexp __P((double, int)); -extern double scalbn __P((double, int)); extern double log __P((double)); extern double log10 __P((double)); -extern double modf __P((double, double*)); +extern double modf __P((double, double *)); extern double pow __P((double, double)); extern double sqrt __P((double)); @@ -149,7 +156,7 @@ extern double scalb __P((double, int)); extern double scalb __P((double, double)); #endif -extern int matherr __P((struct exception*)); +extern int matherr __P((struct exception *)); /* * IEEE Test Vector @@ -175,8 +182,8 @@ extern double log1p __P((double)); * as the second argument; user must allocate space for signgam. */ #ifdef _REENTRANT -extern double gamma_r __P((double, int*)); -extern double lgamma_r __P((double, int*)); +extern double gamma_r __P((double, int *)); +extern double lgamma_r __P((double, int *)); #endif /* _REENTRANT */ /* ieee style elementary functions */ @@ -191,8 +198,8 @@ extern double __ieee754_exp __P((double)); extern double __ieee754_cosh __P((double)); extern double __ieee754_fmod __P((double, double)); extern double __ieee754_pow __P((double, double)); -extern double __ieee754_lgamma_r __P((double, int*)); -extern double __ieee754_gamma_r __P((double, int*)); +extern double __ieee754_lgamma_r __P((double, int *)); +extern double __ieee754_gamma_r __P((double, int *)); extern double __ieee754_lgamma __P((double)); extern double __ieee754_gamma __P((double)); extern double __ieee754_log10 __P((double)); @@ -205,7 +212,7 @@ extern double __ieee754_y1 __P((double)); extern double __ieee754_jn __P((int, double)); extern double __ieee754_yn __P((int, double)); extern double __ieee754_remainder __P((double, double)); -extern int __ieee754_rem_pio2 __P((double, double*)); +extern int __ieee754_rem_pio2 __P((double, double *)); #ifdef _SCALB_INT extern double __ieee754_scalb __P((double, int)); #else @@ -217,10 +224,10 @@ extern double __kernel_standard __P((double, double, int)); extern double __kernel_sin __P((double, double, int)); extern double __kernel_cos __P((double, double)); extern double __kernel_tan __P((double, double, int)); -extern int __kernel_rem_pio2 __P((double*, double*, int, int, int, const int*)); +extern int __kernel_rem_pio2 __P((double *, double *, int, int, int, const int *)); #ifdef __cplusplus -}; +} #endif #endif diff --git a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/k_cos.c b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/k_cos.c new file mode 100644 index 00000000..7f66a06c --- /dev/null +++ b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/k_cos.c @@ -0,0 +1,93 @@ + +/* @(#)k_cos.c 1.3 95/01/18 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* + * __kernel_cos( x, y ) + * kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164 + * Input x is assumed to be bounded by ~pi/4 in magnitude. + * Input y is the tail of x. + * + * Algorithm + * 1. Since cos(-x) = cos(x), we need only to consider positive x. + * 2. if x < 2^-27 (hx<0x3e400000 0), return 1 with inexact if x!=0. + * 3. cos(x) is approximated by a polynomial of degree 14 on + * [0,pi/4] + * 4 14 + * cos(x) ~ 1 - x*x/2 + C1*x + ... + C6*x + * where the remez error is + * + * | 2 4 6 8 10 12 14 | -58 + * |cos(x)-(1-.5*x +C1*x +C2*x +C3*x +C4*x +C5*x +C6*x )| <= 2 + * | | + * + * 4 6 8 10 12 14 + * 4. let r = C1*x +C2*x +C3*x +C4*x +C5*x +C6*x , then + * cos(x) = 1 - x*x/2 + r + * since cos(x+y) ~ cos(x) - sin(x)*y + * ~ cos(x) - x*y, + * a correction term is necessary in cos(x) and hence + * cos(x+y) = 1 - (x*x/2 - (r - x*y)) + * For better accuracy when x > 0.3, let qx = |x|/4 with + * the last 32 bits mask off, and if x > 0.78125, let qx = 0.28125. + * Then + * cos(x+y) = (1-qx) - ((x*x/2-qx) - (r-x*y)). + * Note that 1-qx and (x*x/2-qx) is EXACT here, and the + * magnitude of the latter is at least a quarter of x*x/2, + * thus, reducing the rounding error in the subtraction. + */ + +#include "fdlibm.h" + +#if defined(__STDC__) || defined(__cplusplus) +static const double +#else +static double +#endif + one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ + C1 = 4.16666666666666019037e-02, /* 0x3FA55555, 0x5555554C */ + C2 = -1.38888888888741095749e-03, /* 0xBF56C16C, 0x16C15177 */ + C3 = 2.48015872894767294178e-05, /* 0x3EFA01A0, 0x19CB1590 */ + C4 = -2.75573143513906633035e-07, /* 0xBE927E4F, 0x809C52AD */ + C5 = 2.08757232129817482790e-09, /* 0x3E21EE9E, 0xBDB4B1C4 */ + C6 = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */ + +#if defined(__STDC__) || defined(__cplusplus) +double __kernel_cos(double x, double y) +#else +double __kernel_cos(x, y) +double x, y; +#endif +{ + double a, hz, z, r, qx; + int ix; + ix = __HI(x) & 0x7fffffff; /* ix = |x|'s high word*/ + if (ix < 0x3e400000) { /* if x < 2**27 */ + if (((int)x) == 0) + return one; /* generate inexact */ + } + z = x * x; + r = z * (C1 + z * (C2 + z * (C3 + z * (C4 + z * (C5 + z * C6))))); + if (ix < 0x3FD33333) /* if |x| < 0.3 */ + return one - (0.5 * z - (z * r - x * y)); + else { + if (ix > 0x3fe90000) { /* x > 0.78125 */ + qx = 0.28125; + } else { + __HI(qx) = ix - 0x00200000; /* x/4 */ + __LO(qx) = 0; + } + hz = 0.5 * z - qx; + a = one - qx; + return a - (hz - (z * r - x * y)); + } +} diff --git a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/k_rem_pio2.c b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/k_rem_pio2.c new file mode 100644 index 00000000..ba458843 --- /dev/null +++ b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/k_rem_pio2.c @@ -0,0 +1,349 @@ + +/* @(#)k_rem_pio2.c 1.3 95/01/18 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* + * __kernel_rem_pio2(x,y,e0,nx,prec,ipio2) + * double x[],y[]; int e0,nx,prec; int ipio2[]; + * + * __kernel_rem_pio2 return the last three digits of N with + * y = x - N*pi/2 + * so that |y| < pi/2. + * + * The method is to compute the integer (mod 8) and fraction parts of + * (2/pi)*x without doing the full multiplication. In general we + * skip the part of the product that are known to be a huge integer ( + * more accurately, = 0 mod 8 ). Thus the number of operations are + * independent of the exponent of the input. + * + * (2/pi) is represented by an array of 24-bit integers in ipio2[]. + * + * Input parameters: + * x[] The input value (must be positive) is broken into nx + * pieces of 24-bit integers in double precision format. + * x[i] will be the i-th 24 bit of x. The scaled exponent + * of x[0] is given in input parameter e0 (i.e., x[0]*2^e0 + * match x's up to 24 bits. + * + * Example of breaking a double positive z into x[0]+x[1]+x[2]: + * e0 = ilogb(z)-23 + * z = scalbn(z,-e0) + * for i = 0,1,2 + * x[i] = floor(z) + * z = (z-x[i])*2**24 + * + * + * y[] ouput result in an array of double precision numbers. + * The dimension of y[] is: + * 24-bit precision 1 + * 53-bit precision 2 + * 64-bit precision 2 + * 113-bit precision 3 + * The actual value is the sum of them. Thus for 113-bit + * precison, one may have to do something like: + * + * long double t,w,r_head, r_tail; + * t = (long double)y[2] + (long double)y[1]; + * w = (long double)y[0]; + * r_head = t+w; + * r_tail = w - (r_head - t); + * + * e0 The exponent of x[0] + * + * nx dimension of x[] + * + * prec an integer indicating the precision: + * 0 24 bits (single) + * 1 53 bits (double) + * 2 64 bits (extended) + * 3 113 bits (quad) + * + * ipio2[] + * integer array, contains the (24*i)-th to (24*i+23)-th + * bit of 2/pi after binary point. The corresponding + * floating value is + * + * ipio2[i] * 2^(-24(i+1)). + * + * External function: + * double scalbn(), floor(); + * + * + * Here is the description of some local variables: + * + * jk jk+1 is the initial number of terms of ipio2[] needed + * in the computation. The recommended value is 2,3,4, + * 6 for single, double, extended,and quad. + * + * jz local integer variable indicating the number of + * terms of ipio2[] used. + * + * jx nx - 1 + * + * jv index for pointing to the suitable ipio2[] for the + * computation. In general, we want + * ( 2^e0*x[0] * ipio2[jv-1]*2^(-24jv) )/8 + * is an integer. Thus + * e0-3-24*jv >= 0 or (e0-3)/24 >= jv + * Hence jv = max(0,(e0-3)/24). + * + * jp jp+1 is the number of terms in PIo2[] needed, jp = jk. + * + * q[] double array with integral value, representing the + * 24-bits chunk of the product of x and 2/pi. + * + * q0 the corresponding exponent of q[0]. Note that the + * exponent for q[i] would be q0-24*i. + * + * PIo2[] double precision array, obtained by cutting pi/2 + * into 24 bits chunks. + * + * f[] ipio2[] in floating point + * + * iq[] integer array by breaking up q[] in 24-bits chunk. + * + * fq[] final product of x*(2/pi) in fq[0],..,fq[jk] + * + * ih integer. If >0 it indicates q[] is >= 0.5, hence + * it also indicates the *sign* of the result. + * + */ + +/* + * Constants: + * The hexadecimal values are the intended ones for the following + * constants. The decimal values may be used, provided that the + * compiler will convert from decimal to binary accurately enough + * to produce the hexadecimal values shown. + */ + +#include "fdlibm.h" + +#if defined(__STDC__) || defined(__cplusplus) +static const int init_jk[] = { 2, 3, 4, 6 }; /* initial value for jk */ +#else +static int init_jk[] = { 2, 3, 4, 6 }; +#endif + +#if defined(__STDC__) || defined(__cplusplus) +static const double PIo2[] = { +#else +static double PIo2[] = { +#endif + 1.57079625129699707031e+00, /* 0x3FF921FB, 0x40000000 */ + 7.54978941586159635335e-08, /* 0x3E74442D, 0x00000000 */ + 5.39030252995776476554e-15, /* 0x3CF84698, 0x80000000 */ + 3.28200341580791294123e-22, /* 0x3B78CC51, 0x60000000 */ + 1.27065575308067607349e-29, /* 0x39F01B83, 0x80000000 */ + 1.22933308981111328932e-36, /* 0x387A2520, 0x40000000 */ + 2.73370053816464559624e-44, /* 0x36E38222, 0x80000000 */ + 2.16741683877804819444e-51, /* 0x3569F31D, 0x00000000 */ +}; + +#if defined(__STDC__) || defined(__cplusplus) +static const double +#else +static double +#endif + zero = 0.0, + one = 1.0, two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */ + twon24 = 5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */ + +#if defined(__STDC__) || defined(__cplusplus) +int __kernel_rem_pio2(double *x, double *y, int e0, int nx, int prec, const int *ipio2) +#else +int __kernel_rem_pio2(x, y, e0, nx, prec, ipio2) +double x[], y[]; +int e0, nx, prec; +int ipio2[]; +#endif +{ + int jz, jx, jv, jp, jk, carry, n, iq[20], i, j, k, m, q0, ih; + double z, fw, f[20], fq[20], q[20]; + + /* initialize jk*/ + jk = init_jk[prec]; + jp = jk; + + /* determine jx,jv,q0, note that 3>q0 */ + jx = nx - 1; + jv = (e0 - 3) / 24; + if (jv < 0) + jv = 0; + q0 = e0 - 24 * (jv + 1); + + /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */ + j = jv - jx; + m = jx + jk; + for (i = 0; i <= m; i++, j++) + f[i] = (j < 0) ? zero : (double)ipio2[j]; + + /* compute q[0],q[1],...q[jk] */ + for (i = 0; i <= jk; i++) { + for (j = 0, fw = 0.0; j <= jx; j++) + fw += x[j] * f[jx + i - j]; + q[i] = fw; + } + + jz = jk; +recompute: + /* distill q[] into iq[] reversingly */ + for (i = 0, j = jz, z = q[jz]; j > 0; i++, j--) { + fw = (double)((int)(twon24 * z)); + iq[i] = (int)(z - two24 * fw); + z = q[j - 1] + fw; + } + + /* compute n */ + z = scalbn(z, q0); /* actual value of z */ + z -= 8.0 * floor(z * 0.125); /* trim off integer >= 8 */ + n = (int)z; + z -= (double)n; + ih = 0; + if (q0 > 0) { /* need iq[jz-1] to determine n */ + i = (iq[jz - 1] >> (24 - q0)); + n += i; + iq[jz - 1] -= i << (24 - q0); + ih = iq[jz - 1] >> (23 - q0); + } else if (q0 == 0) + ih = iq[jz - 1] >> 23; + else if (z >= 0.5) + ih = 2; + + if (ih > 0) { /* q > 0.5 */ + n += 1; + carry = 0; + for (i = 0; i < jz; i++) { /* compute 1-q */ + j = iq[i]; + if (carry == 0) { + if (j != 0) { + carry = 1; + iq[i] = 0x1000000 - j; + } + } else + iq[i] = 0xffffff - j; + } + if (q0 > 0) { /* rare case: chance is 1 in 12 */ + switch (q0) { + case 1: + iq[jz - 1] &= 0x7fffff; + break; + case 2: + iq[jz - 1] &= 0x3fffff; + break; + } + } + if (ih == 2) { + z = one - z; + if (carry != 0) + z -= scalbn(one, q0); + } + } + + /* check if recomputation is needed */ + if (z == zero) { + j = 0; + for (i = jz - 1; i >= jk; i--) + j |= iq[i]; + if (j == 0) { /* need recomputation */ + for (k = 1; iq[jk - k] == 0; k++) + ; /* k = no. of terms needed */ + + for (i = jz + 1; i <= jz + k; i++) { /* add q[jz+1] to q[jz+k] */ + f[jx + i] = (double)ipio2[jv + i]; + for (j = 0, fw = 0.0; j <= jx; j++) + fw += x[j] * f[jx + i - j]; + q[i] = fw; + } + jz += k; + goto recompute; + } + } + + /* chop off zero terms */ + if (z == 0.0) { + jz -= 1; + q0 -= 24; + while (iq[jz] == 0) { + jz--; + q0 -= 24; + } + } else { /* break z into 24-bit if necessary */ + z = scalbn(z, -q0); + if (z >= two24) { + fw = (double)((int)(twon24 * z)); + iq[jz] = (int)(z - two24 * fw); + jz += 1; + q0 += 24; + iq[jz] = (int)fw; + } else + iq[jz] = (int)z; + } + + /* convert integer "bit" chunk to floating-point value */ + fw = scalbn(one, q0); + for (i = jz; i >= 0; i--) { + q[i] = fw * (double)iq[i]; + fw *= twon24; + } + + /* compute PIo2[0,...,jp]*q[jz,...,0] */ + for (i = jz; i >= 0; i--) { + for (fw = 0.0, k = 0; k <= jp && k <= jz - i; k++) + fw += PIo2[k] * q[i + k]; + fq[jz - i] = fw; + } + + /* compress fq[] into y[] */ + switch (prec) { + case 0: + fw = 0.0; + for (i = jz; i >= 0; i--) + fw += fq[i]; + y[0] = (ih == 0) ? fw : -fw; + break; + case 1: + case 2: + fw = 0.0; + for (i = jz; i >= 0; i--) + fw += fq[i]; + y[0] = (ih == 0) ? fw : -fw; + fw = fq[0] - fw; + for (i = 1; i <= jz; i++) + fw += fq[i]; + y[1] = (ih == 0) ? fw : -fw; + break; + case 3: /* painful */ + for (i = jz; i > 0; i--) { + fw = fq[i - 1] + fq[i]; + fq[i] += fq[i - 1] - fw; + fq[i - 1] = fw; + } + for (i = jz; i > 1; i--) { + fw = fq[i - 1] + fq[i]; + fq[i] += fq[i - 1] - fw; + fq[i - 1] = fw; + } + for (fw = 0.0, i = jz; i >= 2; i--) + fw += fq[i]; + if (ih == 0) { + y[0] = fq[0]; + y[1] = fq[1]; + y[2] = fw; + } else { + y[0] = -fq[0]; + y[1] = -fq[1]; + y[2] = -fw; + } + } + return n & 7; +} diff --git a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/k_sin.c b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/k_sin.c new file mode 100644 index 00000000..0def29f3 --- /dev/null +++ b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/k_sin.c @@ -0,0 +1,80 @@ + +/* @(#)k_sin.c 1.3 95/01/18 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* __kernel_sin( x, y, iy) + * kernel sin function on [-pi/4, pi/4], pi/4 ~ 0.7854 + * Input x is assumed to be bounded by ~pi/4 in magnitude. + * Input y is the tail of x. + * Input iy indicates whether y is 0. (if iy=0, y assume to be 0). + * + * Algorithm + * 1. Since sin(-x) = -sin(x), we need only to consider positive x. + * 2. if x < 2^-27 (hx<0x3e400000 0), return x with inexact if x!=0. + * 3. sin(x) is approximated by a polynomial of degree 13 on + * [0,pi/4] + * 3 13 + * sin(x) ~ x + S1*x + ... + S6*x + * where + * + * |sin(x) 2 4 6 8 10 12 | -58 + * |----- - (1+S1*x +S2*x +S3*x +S4*x +S5*x +S6*x )| <= 2 + * | x | + * + * 4. sin(x+y) = sin(x) + sin'(x')*y + * ~ sin(x) + (1-x*x/2)*y + * For better accuracy, let + * 3 2 2 2 2 + * r = x *(S2+x *(S3+x *(S4+x *(S5+x *S6)))) + * then 3 2 + * sin(x) = x + (S1*x + (x *(r-y/2)+y)) + */ + +#include "fdlibm.h" + +#if defined(__STDC__) || defined(__cplusplus) +static const double +#else +static double +#endif + half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */ + S1 = -1.66666666666666324348e-01, /* 0xBFC55555, 0x55555549 */ + S2 = 8.33333333332248946124e-03, /* 0x3F811111, 0x1110F8A6 */ + S3 = -1.98412698298579493134e-04, /* 0xBF2A01A0, 0x19C161D5 */ + S4 = 2.75573137070700676789e-06, /* 0x3EC71DE3, 0x57B1FE7D */ + S5 = -2.50507602534068634195e-08, /* 0xBE5AE5E6, 0x8A2B9CEB */ + S6 = 1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */ + +#if defined(__STDC__) || defined(__cplusplus) +double __kernel_sin(double x, double y, int iy) +#else +double __kernel_sin(x, y, iy) +double x, y; +int iy; /* iy=0 if y is zero */ +#endif +{ + double z, r, v; + int ix; + ix = __HI(x) & 0x7fffffff; /* high word of x */ + if (ix < 0x3e400000) /* |x| < 2**-27 */ + { + if ((int)x == 0) + return x; + } /* generate inexact */ + z = x * x; + v = z * x; + r = S2 + z * (S3 + z * (S4 + z * (S5 + z * S6))); + if (iy == 0) + return x + v * (S1 + z * r); + else + return x - ((z * (half * y - v * r) - y) - v * S1); +} diff --git a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/k_tan.c b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/k_tan.c new file mode 100644 index 00000000..8e46bca0 --- /dev/null +++ b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/k_tan.c @@ -0,0 +1,133 @@ +#pragma ident "@(#)k_tan.c 1.5 04/04/22 SMI" + +/* + * ==================================================== + * Copyright 2004 Sun Microsystems, Inc. All Rights Reserved. + * + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* INDENT OFF */ +/* __kernel_tan( x, y, k ) + * kernel tan function on [-pi/4, pi/4], pi/4 ~ 0.7854 + * Input x is assumed to be bounded by ~pi/4 in magnitude. + * Input y is the tail of x. + * Input k indicates whether tan (if k = 1) or -1/tan (if k = -1) is returned. + * + * Algorithm + * 1. Since tan(-x) = -tan(x), we need only to consider positive x. + * 2. if x < 2^-28 (hx<0x3e300000 0), return x with inexact if x!=0. + * 3. tan(x) is approximated by a odd polynomial of degree 27 on + * [0,0.67434] + * 3 27 + * tan(x) ~ x + T1*x + ... + T13*x + * where + * + * |tan(x) 2 4 26 | -59.2 + * |----- - (1+T1*x +T2*x +.... +T13*x )| <= 2 + * | x | + * + * Note: tan(x+y) = tan(x) + tan'(x)*y + * ~ tan(x) + (1+x*x)*y + * Therefore, for better accuracy in computing tan(x+y), let + * 3 2 2 2 2 + * r = x *(T2+x *(T3+x *(...+x *(T12+x *T13)))) + * then + * 3 2 + * tan(x+y) = x + (T1*x + (x *(r+y)+y)) + * + * 4. For x in [0.67434,pi/4], let y = pi/4 - x, then + * tan(x) = tan(pi/4-y) = (1-tan(y))/(1+tan(y)) + * = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y))) + */ + +#include "fdlibm.h" + +static const double T[] = { + 3.33333333333334091986e-01, /* 3FD55555, 55555563 */ + 1.33333333333201242699e-01, /* 3FC11111, 1110FE7A */ + 5.39682539762260521377e-02, /* 3FABA1BA, 1BB341FE */ + 2.18694882948595424599e-02, /* 3F9664F4, 8406D637 */ + 8.86323982359930005737e-03, /* 3F8226E3, E96E8493 */ + 3.59207910759131235356e-03, /* 3F6D6D22, C9560328 */ + 1.45620945432529025516e-03, /* 3F57DBC8, FEE08315 */ + 5.88041240820264096874e-04, /* 3F4344D8, F2F26501 */ + 2.46463134818469906812e-04, /* 3F3026F7, 1A8D1068 */ + 7.81794442939557092300e-05, /* 3F147E88, A03792A6 */ + 7.14072491382608190305e-05, /* 3F12B80F, 32F0A7E9 */ + -1.85586374855275456654e-05, /* BEF375CB, DB605373 */ + 2.59073051863633712884e-05, /* 3EFB2A70, 74BF7AD4 */ +}; + +static const double one = 1.00000000000000000000e+00; /* 3FF00000, 00000000 */ +static const double pio4 = 7.85398163397448278999e-01; /* 3FE921FB, 54442D18 */ +static const double pio4lo = 3.06161699786838301793e-17; /* 3C81A626, 33145C07 */ +/* INDENT ON */ + +double __kernel_tan(double x, double y, int iy) { + double z, r, v, w, s; + int ix, hx; + + hx = __HI(x); /* high word of x */ + ix = hx & 0x7fffffff; /* high word of |x| */ + if (ix < 0x3e300000) { /* x < 2**-28 */ + if ((int)x == 0) { /* generate inexact */ + if (((ix | __LO(x)) | (iy + 1)) == 0) + return one / __fabs(x); + else { + if (iy == 1) + return x; + else { + return -one / x; + } + } + } + } + if (ix >= 0x3FE59428) { /* |x| >= 0.6744 */ + if (hx < 0) { + x = -x; + y = -y; + } + z = pio4 - x; + w = pio4lo - y; + x = z + w; + y = 0.0; + } + z = x * x; + w = z * z; + /* + * Break x^5*(T[1]+x^2*T[2]+...) into + * x^5(T[1]+x^4*T[3]+...+x^20*T[11]) + + * x^5(x^2*(T[2]+x^4*T[4]+...+x^22*[T12])) + */ + r = T[1] + w * (T[3] + w * (T[5] + w * (T[7] + w * (T[9] + w * T[11])))); + v = z * (T[2] + w * (T[4] + w * (T[6] + w * (T[8] + w * (T[10] + w * T[12]))))); + s = z * x; + r = y + z * (s * (r + v) + y); + r += T[0] * s; + w = x + r; + if (ix >= 0x3FE59428) { + v = (double)iy; + return (double)(1 - ((hx >> 30) & 2)) * (v - 2.0 * (x - (w * w / (w + v) - r))); + } + if (iy == 1) + return w; + else { + /* + * if allow error up to 2 ulp, simply return + * -1.0 / (x+r) here + */ + /* compute -1.0 / (x+r) accurately */ + double a, t; + z = w; + __LO(z) = 0; + v = r - (z - x); /* z+v = r+x */ + t = a = -1.0 / w; /* a = -1.0/w */ + __LO(t) = 0; + s = 1.0 + t * z; + return t + a * (s + t * v); + } +} diff --git a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/readme b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/readme new file mode 100644 index 00000000..6031d257 --- /dev/null +++ b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/readme @@ -0,0 +1,262 @@ +Sourced from https://www.netlib.org/fdlibm/ + + ********************************* + * Announcing FDLIBM Version 5.3 * + ********************************* +============================================================ + FDLIBM +============================================================ + developed at Sun Microsystems, Inc. + +What's new in FDLIBM 5.3? + +CONFIGURE + To build FDLIBM, edit the supplied Makefile or create + a local Makefile by running "sh configure" + using the supplied configure script contributed by Nelson Beebe + +BUGS FIXED + + 1. e_pow.c incorrect results when + x is very close to -1.0 and y is very large, e.g. + pow(-1.0000000000000002e+00,4.5035996273704970e+15) = 0 + pow(-9.9999999999999978e-01,4.5035996273704970e+15) = 0 + Correct results are close to -e and -1/e. + + 2. k_tan.c error was > 1 ulp target for FDLIBM + 5.2: Worst error at least 1.45 ulp at + tan(1.7765241907548024E+269) = 1.7733884462610958E+16 + 5.3: Worst error 0.96 ulp + +NOT FIXED YET + + 3. Compiler failure on non-standard code + Statements like + *(1+(int*)&t1) = 0; + are not standard C and cause some optimizing compilers (e.g. GCC) + to generate bad code under optimization. These cases + are to be addressed in the next release. + +FDLIBM (Freely Distributable LIBM) is a C math library +for machines that support IEEE 754 floating-point arithmetic. +In this release, only double precision is supported. + +FDLIBM is intended to provide a reasonably portable (see +assumptions below), reference quality (below one ulp for +major functions like sin,cos,exp,log) math library +(libm.a). For a copy of FDLIBM, please see + http://www.netlib.org/fdlibm/ +or + http://www.validlab.com/software/ + +-------------- +1. ASSUMPTIONS +-------------- +FDLIBM (double precision version) assumes: + a. IEEE 754 style (if not precise compliance) arithmetic; + b. 32 bit 2's complement integer arithmetic; + c. Each double precision floating-point number must be in IEEE 754 + double format, and that each number can be retrieved as two 32-bit + integers through the using of pointer bashing as in the example + below: + + Example: let y = 2.0 + double fp number y: 2.0 + IEEE double format: 0x4000000000000000 + + Referencing y as two integers: + *(int*)&y,*(1+(int*)&y) = {0x40000000,0x0} (on sparc) + {0x0,0x40000000} (on 386) + + Note: Four macros are defined in fdlibm.h to handle this kind of + retrieving: + + __HI(x) the high part of a double x + (sign,exponent,the first 21 significant bits) + __LO(x) the least 32 significant bits of x + __HIp(x) same as __HI except that the argument is a pointer + to a double + __LOp(x) same as __LO except that the argument is a pointer + to a double + + To ensure obtaining correct ordering, one must define __LITTLE_ENDIAN + during compilation for little endian machine (like 386,486). The + default is big endian. + + If the behavior of pointer bashing is undefined, one may hack on the + macro in fdlibm.h. + + d. IEEE exceptions may trigger "signals" as is common in Unix + implementations. + +------------------- +2. EXCEPTION CASES +------------------- +All exception cases in the FDLIBM functions will be mapped +to one of the following four exceptions: + + +-huge*huge, +-tiny*tiny, +-1.0/0.0, +-0.0/0.0 + (overflow) (underflow) (divided-by-zero) (invalid) + +For example, log(0) is a singularity and is thus mapped to + -1.0/0.0 = -infinity. +That is, FDLIBM's log will compute -one/zero and return the +computed value. On an IEEE machine, this will trigger the +divided-by-zero exception and a negative infinity is returned by +default. + +Similarly, exp(-huge) will be mapped to tiny*tiny to generate +an underflow signal. + + +-------------------------------- +3. STANDARD CONFORMANCE WRAPPER +-------------------------------- +The default FDLIBM functions (compiled with -D_IEEE_LIBM flag) +are in "IEEE spirit" (i.e., return the most reasonable result in +floating-point arithmetic). If one wants FDLIBM to comply with +standards like SVID, X/OPEN, or POSIX/ANSI, then one can +create a multi-standard compliant FDLIBM. In this case, each +function in FDLIBM is actually a standard compliant wrapper +function. + +File organization: + 1. For FDLIBM's kernel (internal) function, + File name Entry point + --------------------------- + k_sin.c __kernel_sin + k_tan.c __kernel_tan + --------------------------- + 2. For functions that have no standards conflict + File name Entry point + --------------------------- + s_sin.c sin + s_erf.c erf + --------------------------- + 3. Ieee754 core functions + File name Entry point + --------------------------- + e_exp.c __ieee754_exp + e_sinh.c __ieee754_sinh + --------------------------- + 4. Wrapper functions + File name Entry point + --------------------------- + w_exp.c exp + w_sinh.c sinh + --------------------------- + +Wrapper functions will twist the result of the ieee754 +function to comply to the standard specified by the value +of _LIB_VERSION + if _LIB_VERSION = _IEEE_, return the ieee754 result; + if _LIB_VERSION = _SVID_, return SVID result; + if _LIB_VERSION = _XOPEN_, return XOPEN result; + if _LIB_VERSION = _POSIX_, return POSIX/ANSI result. +(These are macros, see fdlibm.h for their definition.) + + +-------------------------------- +4. HOW TO CREATE FDLIBM's libm.a +-------------------------------- +There are two types of libm.a. One is IEEE only, and the other is +multi-standard compliant (supports IEEE,XOPEN,POSIX/ANSI,SVID). + +To create the IEEE only libm.a, use + make "CFLAGS = -D_IEEE_LIBM" +This will create an IEEE libm.a, which is smaller in size, and +somewhat faster. + +To create a multi-standard compliant libm, use + make "CFLAGS = -D_IEEE_MODE" --- multi-standard fdlibm: default + to IEEE + make "CFLAGS = -D_XOPEN_MODE" --- multi-standard fdlibm: default + to X/OPEN + make "CFLAGS = -D_POSIX_MODE" --- multi-standard fdlibm: default + to POSIX/ANSI + make "CFLAGS = -D_SVID3_MODE" --- multi-standard fdlibm: default + to SVID + + +Here is how one makes a SVID compliant libm. + Make the library by + make "CFLAGS = -D_SVID3_MODE". + The libm.a of FDLIBM will be multi-standard compliant and + _LIB_VERSION is initialized to the value _SVID_ . + + example1: + --------- + main() + { + double y0(); + printf("y0(1e300) = %1.20e\n",y0(1e300)); + exit(0); + } + + % cc example1.c libm.a + % a.out + y0: TLOSS error + y0(1e300) = 0.00000000000000000000e+00 + + +It is possible to change the default standard in multi-standard +fdlibm. Here is an example of how to do it: + example2: + --------- + #include "fdlibm.h" /* must include FDLIBM's fdlibm.h */ + main() + { + double y0(); + _LIB_VERSION = _IEEE_; + printf("IEEE: y0(1e300) = %1.20e\n",y0(1e300)); + _LIB_VERSION = _XOPEN_; + printf("XOPEN y0(1e300) = %1.20e\n",y0(1e300)); + _LIB_VERSION = _POSIX_; + printf("POSIX y0(1e300) = %1.20e\n",y0(1e300)); + _LIB_VERSION = _SVID_; + printf("SVID y0(1e300) = %1.20e\n",y0(1e300)); + exit(0); + } + + % cc example2.c libm.a + % a.out + IEEE: y0(1e300) = -1.36813604503424810557e-151 + XOPEN y0(1e300) = 0.00000000000000000000e+00 + POSIX y0(1e300) = 0.00000000000000000000e+00 + y0: TLOSS error + SVID y0(1e300) = 0.00000000000000000000e+00 + +Note: Here _LIB_VERSION is a global variable. If global variables + are forbidden, then one should modify fdlibm.h to change + _LIB_VERSION to be a global constant. In this case, one + may not change the value of _LIB_VERSION as in example2. + +--------------------------- +5. NOTES ON PORTING FDLIBM +--------------------------- + Care must be taken when installing FDLIBM over existing + libm.a. + All co-existing function prototypes must agree, otherwise + users will encounter mysterious failures. + + So far, the only known likely conflict is the declaration + of the IEEE recommended function scalb: + + double scalb(double,double) (1) SVID3 defined + double scalb(double,int) (2) IBM,DEC,... + + FDLIBM follows Sun definition and use (1) as default. + If one's existing libm.a uses (2), then one may raise + the flags _SCALB_INT during the compilation of FDLIBM + to get the correct function prototype. + (E.g., make "CFLAGS = -D_IEEE_LIBM -D_SCALB_INT".) + NOTE that if -D_SCALB_INT is raised, it won't be SVID3 + conformant. + +-------------- +6. PROBLEMS ? +-------------- +Please send comments and bug reports to the electronic mail address +suggested by: + fdlibm-comments AT sun.com + diff --git a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_atan.c b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_atan.c new file mode 100644 index 00000000..8876d1d3 --- /dev/null +++ b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_atan.c @@ -0,0 +1,143 @@ + +/* @(#)s_atan.c 1.3 95/01/18 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + * + */ + +/* atan(x) + * Method + * 1. Reduce x to positive by atan(x) = -atan(-x). + * 2. According to the integer k=4t+0.25 chopped, t=x, the argument + * is further reduced to one of the following intervals and the + * arctangent of t is evaluated by the corresponding formula: + * + * [0,7/16] atan(x) = t-t^3*(a1+t^2*(a2+...(a10+t^2*a11)...) + * [7/16,11/16] atan(x) = atan(1/2) + atan( (t-0.5)/(1+t/2) ) + * [11/16.19/16] atan(x) = atan( 1 ) + atan( (t-1)/(1+t) ) + * [19/16,39/16] atan(x) = atan(3/2) + atan( (t-1.5)/(1+1.5t) ) + * [39/16,INF] atan(x) = atan(INF) + atan( -1/t ) + * + * Constants: + * The hexadecimal values are the intended ones for the following + * constants. The decimal values may be used, provided that the + * compiler will convert from decimal to binary accurately enough + * to produce the hexadecimal values shown. + */ + +#include "fdlibm.h" + +#if defined(__STDC__) || defined(__cplusplus) +static const double atanhi[] = { +#else +static double atanhi[] = { +#endif + 4.63647609000806093515e-01, /* atan(0.5)hi 0x3FDDAC67, 0x0561BB4F */ + 7.85398163397448278999e-01, /* atan(1.0)hi 0x3FE921FB, 0x54442D18 */ + 9.82793723247329054082e-01, /* atan(1.5)hi 0x3FEF730B, 0xD281F69B */ + 1.57079632679489655800e+00, /* atan(inf)hi 0x3FF921FB, 0x54442D18 */ +}; + +#if defined(__STDC__) || defined(__cplusplus) +static const double atanlo[] = { +#else +static double atanlo[] = { +#endif + 2.26987774529616870924e-17, /* atan(0.5)lo 0x3C7A2B7F, 0x222F65E2 */ + 3.06161699786838301793e-17, /* atan(1.0)lo 0x3C81A626, 0x33145C07 */ + 1.39033110312309984516e-17, /* atan(1.5)lo 0x3C700788, 0x7AF0CBBD */ + 6.12323399573676603587e-17, /* atan(inf)lo 0x3C91A626, 0x33145C07 */ +}; + +#if defined(__STDC__) || defined(__cplusplus) +static const double aT[] = { +#else +static double aT[] = { +#endif + 3.33333333333329318027e-01, /* 0x3FD55555, 0x5555550D */ + -1.99999999998764832476e-01, /* 0xBFC99999, 0x9998EBC4 */ + 1.42857142725034663711e-01, /* 0x3FC24924, 0x920083FF */ + -1.11111104054623557880e-01, /* 0xBFBC71C6, 0xFE231671 */ + 9.09088713343650656196e-02, /* 0x3FB745CD, 0xC54C206E */ + -7.69187620504482999495e-02, /* 0xBFB3B0F2, 0xAF749A6D */ + 6.66107313738753120669e-02, /* 0x3FB10D66, 0xA0D03D51 */ + -5.83357013379057348645e-02, /* 0xBFADDE2D, 0x52DEFD9A */ + 4.97687799461593236017e-02, /* 0x3FA97B4B, 0x24760DEB */ + -3.65315727442169155270e-02, /* 0xBFA2B444, 0x2C6A6C2F */ + 1.62858201153657823623e-02, /* 0x3F90AD3A, 0xE322DA11 */ +}; + +#if defined(__STDC__) || defined(__cplusplus) +static const double +#else +static double +#endif + one = 1.0, + huge = 1.0e300; + +#if defined(__STDC__) || defined(__cplusplus) +double atan(double x) +#else +double atan(x) +double x; +#endif +{ + double w, s1, s2, z; + int ix, hx, id; + + hx = __HI(x); + ix = hx & 0x7fffffff; + if (ix >= 0x44100000) { /* if |x| >= 2^66 */ + if (ix > 0x7ff00000 || (ix == 0x7ff00000 && (__LO(x) != 0))) + return x + x; /* NaN */ + if (hx > 0) + return atanhi[3] + atanlo[3]; + else + return -atanhi[3] - atanlo[3]; + } + if (ix < 0x3fdc0000) { /* |x| < 0.4375 */ + if (ix < 0x3e200000) { /* |x| < 2^-29 */ + if (huge + x > one) + return x; /* raise inexact */ + } + id = -1; + } else { + x = __fabs(x); + if (ix < 0x3ff30000) { /* |x| < 1.1875 */ + if (ix < 0x3fe60000) { /* 7/16 <=|x|<11/16 */ + id = 0; + x = (2.0 * x - one) / (2.0 + x); + } else { /* 11/16<=|x|< 19/16 */ + id = 1; + x = (x - one) / (x + one); + } + } else { + if (ix < 0x40038000) { /* |x| < 2.4375 */ + id = 2; + x = (x - 1.5) / (one + 1.5 * x); + } else { /* 2.4375 <= |x| < 2^66 */ + id = 3; + x = -1.0 / x; + } + } + } + /* end of argument reduction */ + z = x * x; + w = z * z; + /* break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly */ + s1 = z * (aT[0] + w * (aT[2] + w * (aT[4] + w * (aT[6] + w * (aT[8] + w * aT[10]))))); + s2 = w * (aT[1] + w * (aT[3] + w * (aT[5] + w * (aT[7] + w * aT[9])))); + if (id < 0) + return x - x * (s1 + s2); + else { + z = atanhi[id] - ((x * (s1 + s2) - atanlo[id]) - x); + return (hx < 0) ? -z : z; + } +} diff --git a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_ceil.c b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_ceil.c new file mode 100644 index 00000000..027ff4a0 --- /dev/null +++ b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_ceil.c @@ -0,0 +1,91 @@ + +/* @(#)s_ceil.c 1.3 95/01/18 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* + * ceil(x) + * Return x rounded toward -inf to integral value + * Method: + * Bit twiddling. + * Exception: + * Inexact flag raised if x not equal to ceil(x). + */ + +#include "fdlibm.h" + +#if defined(__STDC__) || defined(__cplusplus) +static const double huge = 1.0e300; +#else +static double huge = 1.0e300; +#endif + +#if defined(__STDC__) || defined(__cplusplus) +double ceil(double x) +#else +double ceil(x) +double x; +#endif +{ + int i0, i1, j0; + unsigned i, j; + i0 = __HI(x); + i1 = __LO(x); + j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; + if (j0 < 20) { + if (j0 < 0) { /* raise inexact if x != 0 */ + if (huge + x > 0.0) { /* return 0*sign(x) if |x|<1 */ + if (i0 < 0) { + i0 = 0x80000000; + i1 = 0; + } else if ((i0 | i1) != 0) { + i0 = 0x3ff00000; + i1 = 0; + } + } + } else { + i = (0x000fffff) >> j0; + if (((i0 & i) | i1) == 0) + return x; /* x is integral */ + if (huge + x > 0.0) { /* raise inexact flag */ + if (i0 > 0) + i0 += (0x00100000) >> j0; + i0 &= (~i); + i1 = 0; + } + } + } else if (j0 > 51) { + if (j0 == 0x400) + return x + x; /* inf or NaN */ + else + return x; /* x is integral */ + } else { + i = ((unsigned)(0xffffffff)) >> (j0 - 20); + if ((i1 & i) == 0) + return x; /* x is integral */ + if (huge + x > 0.0) { /* raise inexact flag */ + if (i0 > 0) { + if (j0 == 20) + i0 += 1; + else { + j = i1 + (1 << (52 - j0)); + if (j < i1) + i0 += 1; /* got a carry */ + i1 = j; + } + } + i1 &= (~i); + } + } + __HI(x) = i0; + __LO(x) = i1; + return x; +} diff --git a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_copysign.c b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_copysign.c new file mode 100644 index 00000000..765440ae --- /dev/null +++ b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_copysign.c @@ -0,0 +1,31 @@ + +/* @(#)s_copysign.c 1.3 95/01/18 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* + * copysign(double x, double y) + * copysign(x,y) returns a value with the magnitude of x and + * with the sign bit of y. + */ + +#include "fdlibm.h" + +#if defined(__STDC__) || defined(__cplusplus) +double copysign(double x, double y) +#else +double copysign(x, y) +double x, y; +#endif +{ + __HI(x) = (__HI(x) & 0x7fffffff) | (__HI(y) & 0x80000000); + return x; +} diff --git a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_cos.c b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_cos.c new file mode 100644 index 00000000..053169be --- /dev/null +++ b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_cos.c @@ -0,0 +1,83 @@ + +/* @(#)s_cos.c 1.3 95/01/18 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* cos(x) + * Return cosine function of x. + * + * kernel function: + * __kernel_sin ... sine function on [-pi/4,pi/4] + * __kernel_cos ... cosine function on [-pi/4,pi/4] + * __ieee754_rem_pio2 ... argument reduction routine + * + * Method. + * Let S,C and T denote the sin, cos and tan respectively on + * [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2 + * in [-pi/4 , +pi/4], and let n = k mod 4. + * We have + * + * n sin(x) cos(x) tan(x) + * ---------------------------------------------------------- + * 0 S C T + * 1 C -S -1/T + * 2 -S -C T + * 3 -C S -1/T + * ---------------------------------------------------------- + * + * Special cases: + * Let trig be any of sin, cos, or tan. + * trig(+-INF) is NaN, with signals; + * trig(NaN) is that NaN; + * + * Accuracy: + * TRIG(x) returns trig(x) nearly rounded + */ + +#include "fdlibm.h" + +#if defined(__STDC__) || defined(__cplusplus) +double cos(double x) +#else +double cos(x) +double x; +#endif +{ + double y[2], z = 0.0; + int n, ix; + + /* High word of x. */ + ix = __HI(x); + + /* |x| ~< pi/4 */ + ix &= 0x7fffffff; + if (ix <= 0x3fe921fb) + return __kernel_cos(x, z); + + /* cos(Inf or NaN) is NaN */ + else if (ix >= 0x7ff00000) + return x - x; + + /* argument reduction needed */ + else { + n = __ieee754_rem_pio2(x, y); + switch (n & 3) { + case 0: + return __kernel_cos(y[0], y[1]); + case 1: + return -__kernel_sin(y[0], y[1], 1); + case 2: + return -__kernel_cos(y[0], y[1]); + default: + return __kernel_sin(y[0], y[1], 1); + } + } +} diff --git a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_floor.c b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_floor.c new file mode 100644 index 00000000..a3693613 --- /dev/null +++ b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_floor.c @@ -0,0 +1,90 @@ + +/* @(#)s_floor.c 1.3 95/01/18 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* + * floor(x) + * Return x rounded toward -inf to integral value + * Method: + * Bit twiddling. + * Exception: + * Inexact flag raised if x not equal to floor(x). + */ + +#include "fdlibm.h" + +#if defined(__STDC__) || defined(__cplusplus) +static const double huge = 1.0e300; +#else +static double huge = 1.0e300; +#endif + +#if defined(__STDC__) || defined(__cplusplus) +double floor(double x) +#else +double floor(x) +double x; +#endif +{ + int i0, i1, j0; + unsigned i, j; + i0 = __HI(x); + i1 = __LO(x); + j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; + if (j0 < 20) { + if (j0 < 0) { /* raise inexact if x != 0 */ + if (huge + x > 0.0) { /* return 0*sign(x) if |x|<1 */ + if (i0 >= 0) { + i0 = i1 = 0; + } else if (((i0 & 0x7fffffff) | i1) != 0) { + i0 = 0xbff00000; + i1 = 0; + } + } + } else { + i = (0x000fffff) >> j0; + if (((i0 & i) | i1) == 0) + return x; /* x is integral */ + if (huge + x > 0.0) { /* raise inexact flag */ + if (i0 < 0) + i0 += (0x00100000) >> j0; + i0 &= (~i); + i1 = 0; + } + } + } else if (j0 > 51) { + if (j0 == 0x400) + return x + x; /* inf or NaN */ + else + return x; /* x is integral */ + } else { + i = ((unsigned)(0xffffffff)) >> (j0 - 20); + if ((i1 & i) == 0) + return x; /* x is integral */ + if (huge + x > 0.0) { /* raise inexact flag */ + if (i0 < 0) { + if (j0 == 20) + i0 += 1; + else { + j = i1 + (1 << (52 - j0)); + if (j < i1) + i0 += 1; /* got a carry */ + i1 = j; + } + } + i1 &= (~i); + } + } + __HI(x) = i0; + __LO(x) = i1; + return x; +} diff --git a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_frexp.c b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_frexp.c new file mode 100644 index 00000000..f66992ed --- /dev/null +++ b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_frexp.c @@ -0,0 +1,58 @@ + +/* @(#)s_frexp.c 1.4 95/01/18 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* + * for non-zero x + * x = frexp(arg,&exp); + * return a double fp quantity x such that 0.5 <= |x| <1.0 + * and the corresponding binary exponent "exp". That is + * arg = x*2^exp. + * If arg is inf, 0.0, or NaN, then frexp(arg,&exp) returns arg + * with *exp=0. + */ + +#include "fdlibm.h" + +#if defined(__STDC__) || defined(__cplusplus) +static const double +#else +static double +#endif + two54 = 1.80143985094819840000e+16; /* 0x43500000, 0x00000000 */ + +#if defined(__STDC__) || defined(__cplusplus) +double frexp(double x, int *eptr) +#else +double frexp(x, eptr) +double x; +int *eptr; +#endif +{ + int hx, ix, lx; + hx = __HI(x); + ix = 0x7fffffff & hx; + lx = __LO(x); + *eptr = 0; + if (ix >= 0x7ff00000 || ((ix | lx) == 0)) + return x; /* 0,inf,nan */ + if (ix < 0x00100000) { /* subnormal */ + x *= two54; + hx = __HI(x); + ix = hx & 0x7fffffff; + *eptr = -54; + } + *eptr += (ix >> 20) - 1022; + hx = (hx & 0x800fffff) | 0x3fe00000; + __HI(x) = hx; + return x; +} diff --git a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_ldexp.c b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_ldexp.c new file mode 100644 index 00000000..ce2e92d8 --- /dev/null +++ b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_ldexp.c @@ -0,0 +1,70 @@ + +/* @(#)s_ldexp.c 1.3 95/01/18 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +#include "fdlibm.h" +#include +#include + +#ifdef __STDC__ +static const double +#else +static double +#endif + two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */ + twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */ + huge = 1.0e+300, tiny = 1.0e-300; + +#if defined(__STDC__) || defined(__cplusplus) +double ldexp(double value, int exp) +#else +double ldexp(value, exp) +double value; +int exp; +#endif +{ + int k, hx, lx; + + if (!isfinite(value) || value == 0.0) + return value; + + hx = __HI(value); + lx = __LO(value); + k = (hx & 0x7ff00000) >> 20; /* extract exponent */ + if (k == 0) { /* 0 or subnormal x */ + if ((lx | (hx & 0x7fffffff)) == 0) + return value; /* +-0 */ + value *= two54; + hx = __HI(value); + k = ((hx & 0x7ff00000) >> 20) - 54; + if (exp < -50000) + return tiny * value; /*underflow*/ + } + if (k == 0x7ff) + return value + value; /* NaN or Inf */ + k = k + exp; + if (k > 0x7fe) + return huge * copysign(huge, value); /* overflow */ + if (k > 0) /* normal result */ + { + __HI(value) = (hx & 0x800fffff) | (k << 20); + return value; + } + if (k <= -54) + if (exp > 50000) /* in case integer overflow in n+k */ + return huge * copysign(huge, value); /*overflow*/ + else + return tiny * copysign(tiny, value); /*underflow*/ + k += 54; /* subnormal result */ + __HI(value) = (hx & 0x800fffff) | (k << 20); + return value * twom54; +} diff --git a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_modf.c b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_modf.c new file mode 100644 index 00000000..b3e25f81 --- /dev/null +++ b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_modf.c @@ -0,0 +1,80 @@ + +/* @(#)s_modf.c 1.3 95/01/18 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* + * modf(double x, double *iptr) + * return fraction part of x, and return x's integral part in *iptr. + * Method: + * Bit twiddling. + * + * Exception: + * No exception. + */ + +#include "fdlibm.h" + +#if defined(__STDC__) || defined(__cplusplus) +static const double one = 1.0; +#else +static double one = 1.0; +#endif + +#if defined(__STDC__) || defined(__cplusplus) +double modf(double x, double *iptr) +#else +double modf(x, iptr) +double x, *iptr; +#endif +{ + int i0, i1, j0; + unsigned i; + i0 = __HI(x); /* high x */ + i1 = __LO(x); /* low x */ + j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; /* exponent of x */ + if (j0 < 20) { /* integer part in high x */ + if (j0 < 0) { /* |x|<1 */ + __HIp(iptr) = i0 & 0x80000000; + __LOp(iptr) = 0; /* *iptr = +-0 */ + return x; + } else { + i = (0x000fffff) >> j0; + if (((i0 & i) | i1) == 0) { /* x is integral */ + *iptr = x; + __HI(x) &= 0x80000000; + __LO(x) = 0; /* return +-0 */ + return x; + } else { + __HIp(iptr) = i0 & (~i); + __LOp(iptr) = 0; + return x - *iptr; + } + } + } else if (j0 > 51) { /* no fraction part */ + *iptr = x * one; + __HI(x) &= 0x80000000; + __LO(x) = 0; /* return +-0 */ + return x; + } else { /* fraction part in low x */ + i = ((unsigned)(0xffffffff)) >> (j0 - 20); + if ((i1 & i) == 0) { /* x is integral */ + *iptr = x; + __HI(x) &= 0x80000000; + __LO(x) = 0; /* return +-0 */ + return x; + } else { + __HIp(iptr) = i0; + __LOp(iptr) = i1 & (~i); + return x - *iptr; + } + } +} diff --git a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_rint.c b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_rint.c new file mode 100644 index 00000000..98b76cf7 --- /dev/null +++ b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_rint.c @@ -0,0 +1,92 @@ + +/* @(#)s_rint.c 1.3 95/01/18 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* + * rint(x) + * Return x rounded to integral value according to the prevailing + * rounding mode. + * Method: + * Using floating addition. + * Exception: + * Inexact flag raised if x not equal to rint(x). + */ + +#include "fdlibm.h" + +#if defined(__STDC__) || defined(__cplusplus) +static const double +#else +static double +#endif + TWO52[2] = { + 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */ + -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */ + }; + +#if defined(__STDC__) || defined(__cplusplus) +double rint(double x) +#else +double rint(x) +double x; +#endif +{ + int i0, j0, sx; + unsigned i, i1; + double w, t; + i0 = __HI(x); + sx = (i0 >> 31) & 1; + i1 = __LO(x); + j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; + if (j0 < 20) { + if (j0 < 0) { + if (((i0 & 0x7fffffff) | i1) == 0) + return x; + i1 |= (i0 & 0x0fffff); + i0 &= 0xfffe0000; + i0 |= ((i1 | -i1) >> 12) & 0x80000; + __HI(x) = i0; + w = TWO52[sx] + x; + t = w - TWO52[sx]; + i0 = __HI(t); + __HI(t) = (i0 & 0x7fffffff) | (sx << 31); + return t; + } else { + i = (0x000fffff) >> j0; + if (((i0 & i) | i1) == 0) + return x; /* x is integral */ + i >>= 1; + if (((i0 & i) | i1) != 0) { + if (j0 == 19) + i1 = 0x40000000; + else + i0 = (i0 & (~i)) | ((0x20000) >> j0); + } + } + } else if (j0 > 51) { + if (j0 == 0x400) + return x + x; /* inf or NaN */ + else + return x; /* x is integral */ + } else { + i = ((unsigned)(0xffffffff)) >> (j0 - 20); + if ((i1 & i) == 0) + return x; /* x is integral */ + i >>= 1; + if ((i1 & i) != 0) + i1 = (i1 & (~i)) | ((0x40000000) >> (j0 - 20)); + } + __HI(x) = i0; + __LO(x) = i1; + w = TWO52[sx] + x; + return w - TWO52[sx]; +} diff --git a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_sin.c b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_sin.c new file mode 100644 index 00000000..26b72393 --- /dev/null +++ b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_sin.c @@ -0,0 +1,83 @@ + +/* @(#)s_sin.c 1.3 95/01/18 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* sin(x) + * Return sine function of x. + * + * kernel function: + * __kernel_sin ... sine function on [-pi/4,pi/4] + * __kernel_cos ... cose function on [-pi/4,pi/4] + * __ieee754_rem_pio2 ... argument reduction routine + * + * Method. + * Let S,C and T denote the sin, cos and tan respectively on + * [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2 + * in [-pi/4 , +pi/4], and let n = k mod 4. + * We have + * + * n sin(x) cos(x) tan(x) + * ---------------------------------------------------------- + * 0 S C T + * 1 C -S -1/T + * 2 -S -C T + * 3 -C S -1/T + * ---------------------------------------------------------- + * + * Special cases: + * Let trig be any of sin, cos, or tan. + * trig(+-INF) is NaN, with signals; + * trig(NaN) is that NaN; + * + * Accuracy: + * TRIG(x) returns trig(x) nearly rounded + */ + +#include "fdlibm.h" + +#if defined(__STDC__) || defined(__cplusplus) +double sin(double x) +#else +double sin(x) +double x; +#endif +{ + double y[2], z = 0.0; + int n, ix; + + /* High word of x. */ + ix = __HI(x); + + /* |x| ~< pi/4 */ + ix &= 0x7fffffff; + if (ix <= 0x3fe921fb) + return __kernel_sin(x, z, 0); + + /* sin(Inf or NaN) is NaN */ + else if (ix >= 0x7ff00000) + return x - x; + + /* argument reduction needed */ + else { + n = __ieee754_rem_pio2(x, y); + switch (n & 3) { + case 0: + return __kernel_sin(y[0], y[1], 1); + case 1: + return __kernel_cos(y[0], y[1]); + case 2: + return -__kernel_sin(y[0], y[1], 1); + default: + return -__kernel_cos(y[0], y[1]); + } + } +} diff --git a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_tan.c b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_tan.c new file mode 100644 index 00000000..fa827645 --- /dev/null +++ b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/s_tan.c @@ -0,0 +1,74 @@ + +/* @(#)s_tan.c 1.3 95/01/18 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* tan(x) + * Return tangent function of x. + * + * kernel function: + * __kernel_tan ... tangent function on [-pi/4,pi/4] + * __ieee754_rem_pio2 ... argument reduction routine + * + * Method. + * Let S,C and T denote the sin, cos and tan respectively on + * [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2 + * in [-pi/4 , +pi/4], and let n = k mod 4. + * We have + * + * n sin(x) cos(x) tan(x) + * ---------------------------------------------------------- + * 0 S C T + * 1 C -S -1/T + * 2 -S -C T + * 3 -C S -1/T + * ---------------------------------------------------------- + * + * Special cases: + * Let trig be any of sin, cos, or tan. + * trig(+-INF) is NaN, with signals; + * trig(NaN) is that NaN; + * + * Accuracy: + * TRIG(x) returns trig(x) nearly rounded + */ + +#include "fdlibm.h" + +#if defined(__STDC__) || defined(__cplusplus) +double tan(double x) +#else +double tan(x) +double x; +#endif +{ + double y[2], z = 0.0; + int n, ix; + + /* High word of x. */ + ix = __HI(x); + + /* |x| ~< pi/4 */ + ix &= 0x7fffffff; + if (ix <= 0x3fe921fb) + return __kernel_tan(x, z, 1); + + /* tan(Inf or NaN) is NaN */ + else if (ix >= 0x7ff00000) + return x - x; /* NaN */ + + /* argument reduction needed */ + else { + n = __ieee754_rem_pio2(x, y); + return __kernel_tan(y[0], y[1], 1 - ((n & 1) << 1)); /* 1 -- n even + -1 -- n odd */ + } +} diff --git a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/w_acos.c b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/w_acos.c index 36c54503..0e9813d7 100644 --- a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/w_acos.c +++ b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/w_acos.c @@ -1,2 +1,39 @@ -#include "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/fdlibm.h" -double acos(double x) { return __ieee754_acos(x); } + +/* @(#)w_acos.c 1.3 95/01/18 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* + * wrap_acos(x) + */ + +#include "fdlibm.h" + +#if defined(__STDC__) || defined(__cplusplus) +double acos(double x) /* wrapper acos */ +#else +double acos(x) /* wrapper acos */ +double x; +#endif +{ +#ifdef _IEEE_LIBM + return __ieee754_acos(x); +#else + double z; + z = __ieee754_acos(x); + if (_LIB_VERSION == _IEEE_ || isnan(x)) + return z; + if (fabs(x) > 1.0) { + return __kernel_standard(x, x, 1); /* acos(|x|>1) */ + } else + return z; +#endif +} diff --git a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/w_asin.c b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/w_asin.c index bfb53342..02ad2b43 100644 --- a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/w_asin.c +++ b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/w_asin.c @@ -1,2 +1,40 @@ -#include "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/fdlibm.h" -double asin(double x) { return __ieee754_asin(x); } + +/* @(#)w_asin.c 1.3 95/01/18 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + * + */ + +/* + * wrapper asin(x) + */ + +#include "fdlibm.h" + +#if defined(__STDC__) || defined(__cplusplus) +double asin(double x) /* wrapper asin */ +#else +double asin(x) /* wrapper asin */ +double x; +#endif +{ +#ifdef _IEEE_LIBM + return __ieee754_asin(x); +#else + double z; + z = __ieee754_asin(x); + if (_LIB_VERSION == _IEEE_ || isnan(x)) + return z; + if (fabs(x) > 1.0) { + return __kernel_standard(x, x, 2); /* asin(|x|>1) */ + } else + return z; +#endif +} diff --git a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/w_atan2.c b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/w_atan2.c index d2bf4063..9ba667d1 100644 --- a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/w_atan2.c +++ b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/w_atan2.c @@ -1,2 +1,40 @@ -#include "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/fdlibm.h" -double atan2(double x, double y) { return __ieee754_atan2(x, y); } + +/* @(#)w_atan2.c 1.3 95/01/18 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + * + */ + +/* + * wrapper atan2(y,x) + */ + +#include "fdlibm.h" + +#if defined(__STDC__) || defined(__cplusplus) +double atan2(double y, double x) /* wrapper atan2 */ +#else +double atan2(y, x) /* wrapper atan2 */ +double y, x; +#endif +{ +#ifdef _IEEE_LIBM + return __ieee754_atan2(y, x); +#else + double z; + z = __ieee754_atan2(y, x); + if (_LIB_VERSION == _IEEE_ || isnan(x) || isnan(y)) + return z; + if (x == 0.0 && y == 0.0) { + return __kernel_standard(y, x, 3); /* atan2(+-0,+-0) */ + } else + return z; +#endif +} diff --git a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/w_exp.c b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/w_exp.c index 5a6d74b7..5c3139bd 100644 --- a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/w_exp.c +++ b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/w_exp.c @@ -1,2 +1,49 @@ -#include "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/fdlibm.h" -double exp(double x) { return __ieee754_exp(x); } + +/* @(#)w_exp.c 1.4 04/04/22 */ +/* + * ==================================================== + * Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved. + * + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* + * wrapper exp(x) + */ + +#include "fdlibm.h" + +#if defined(__STDC__) || defined(__cplusplus) +static const double +#else +static double +#endif + o_threshold = 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */ + u_threshold = -7.45133219101941108420e+02; /* 0xc0874910, 0xD52D3051 */ + +#if defined(__STDC__) || defined(__cplusplus) +double exp(double x) /* wrapper exp */ +#else +double exp(x) /* wrapper exp */ +double x; +#endif +{ +#ifdef _IEEE_LIBM + return __ieee754_exp(x); +#else + double z; + z = __ieee754_exp(x); + if (_LIB_VERSION == _IEEE_) + return z; + if (finite(x)) { + if (x > o_threshold) + return __kernel_standard(x, x, 6); /* exp overflow */ + else if (x < u_threshold) + return __kernel_standard(x, x, 7); /* exp underflow */ + } + return z; +#endif +} diff --git a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/w_fmod.c b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/w_fmod.c index 67732982..db05617f 100644 --- a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/w_fmod.c +++ b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/w_fmod.c @@ -1,2 +1,39 @@ -#include "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/fdlibm.h" -double fmod(double x, double y) { return __ieee754_fmod(x,y); } + +/* @(#)w_fmod.c 1.3 95/01/18 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* + * wrapper fmod(x,y) + */ + +#include "fdlibm.h" + +#if defined(__STDC__) || defined(__cplusplus) +double fmod(double x, double y) /* wrapper fmod */ +#else +double fmod(x, y) /* wrapper fmod */ +double x, y; +#endif +{ +#ifdef _IEEE_LIBM + return __ieee754_fmod(x, y); +#else + double z; + z = __ieee754_fmod(x, y); + if (_LIB_VERSION == _IEEE_ || isnan(y) || isnan(x)) + return z; + if (y == 0.0) { + return __kernel_standard(x, y, 27); /* fmod(x,0) */ + } else + return z; +#endif +} diff --git a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/w_log.c b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/w_log.c index 4a2a432f..b1542a9f 100644 --- a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/w_log.c +++ b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/w_log.c @@ -1,2 +1,39 @@ -#include "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/fdlibm.h" -double log(double x) { return __ieee754_log(x); } + +/* @(#)w_log.c 1.3 95/01/18 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* + * wrapper log(x) + */ + +#include "fdlibm.h" + +#if defined(__STDC__) || defined(__cplusplus) +double log(double x) /* wrapper log */ +#else +double log(x) /* wrapper log */ +double x; +#endif +{ +#ifdef _IEEE_LIBM + return __ieee754_log(x); +#else + double z; + z = __ieee754_log(x); + if (_LIB_VERSION == _IEEE_ || isnan(x) || x > 0.0) + return z; + if (x == 0.0) + return __kernel_standard(x, x, 16); /* log(0) */ + else + return __kernel_standard(x, x, 17); /* log(x<0) */ +#endif +} diff --git a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/w_log10.c b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/w_log10.c index 1285a232..54dee48d 100644 --- a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/w_log10.c +++ b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/w_log10.c @@ -1,2 +1,42 @@ -#include "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/fdlibm.h" -double log10(double x) { return __ieee754_log10(x); } + +/* @(#)w_log10.c 1.3 95/01/18 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* + * wrapper log10(X) + */ + +#include "fdlibm.h" + +#if defined(__STDC__) || defined(__cplusplus) +double log10(double x) /* wrapper log10 */ +#else +double log10(x) /* wrapper log10 */ +double x; +#endif +{ +#ifdef _IEEE_LIBM + return __ieee754_log10(x); +#else + double z; + z = __ieee754_log10(x); + if (_LIB_VERSION == _IEEE_ || isnan(x)) + return z; + if (x <= 0.0) { + if (x == 0.0) + return __kernel_standard(x, x, 18); /* log10(0) */ + else + return __kernel_standard(x, x, 19); /* log10(x<0) */ + } else + return z; +#endif +} diff --git a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/w_pow.c b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/w_pow.c index a14efdde..7a2b1331 100644 --- a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/w_pow.c +++ b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/Double_precision/w_pow.c @@ -1,2 +1,60 @@ -#include "sdk/PowerPC_EABI_Support/MSL/MSL_C/MSL_Common_Embedded/Math/fdlibm.h" -double pow(double m, double e) { return __ieee754_pow(m, e); } + + +/* @(#)w_pow.c 1.3 95/01/18 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* + * wrapper pow(x,y) return x**y + */ + +#include "fdlibm.h" + +#if defined(__STDC__) || defined(__cplusplus) +double pow(double x, double y) /* wrapper pow */ +#else +double pow(x, y) /* wrapper pow */ +double x, y; +#endif +{ +#ifdef _IEEE_LIBM + return __ieee754_pow(x, y); +#else + double z; + z = __ieee754_pow(x, y); + if (_LIB_VERSION == _IEEE_ || isnan(y)) + return z; + if (isnan(x)) { + if (y == 0.0) + return __kernel_standard(x, y, 42); /* pow(NaN,0.0) */ + else + return z; + } + if (x == 0.0) { + if (y == 0.0) + return __kernel_standard(x, y, 20); /* pow(0.0,0.0) */ + if (finite(y) && y < 0.0) + return __kernel_standard(x, y, 23); /* pow(0.0,negative) */ + return z; + } + if (!finite(z)) { + if (finite(x) && finite(y)) { + if (isnan(z)) + return __kernel_standard(x, y, 24); /* pow neg**non-int */ + else + return __kernel_standard(x, y, 21); /* pow overflow */ + } + } + if (z == 0.0 && finite(x) && finite(y)) + return __kernel_standard(x, y, 22); /* pow underflow */ + return z; +#endif +} diff --git a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/errno.h b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/errno.h index 62433792..ead68414 100644 --- a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/errno.h +++ b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/errno.h @@ -53,7 +53,7 @@ extern "C" { #define EUNKNOWN 99 /* clang-format on */ -extern int errno; +__declspec(section ".sdata") extern int errno; #ifdef __cplusplus } diff --git a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/math.h b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/math.h index 44380fdb..7dbb80e8 100644 --- a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/math.h +++ b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/math.h @@ -2,18 +2,19 @@ #define _MATH_H #include "MSL_Common/math_api.h" +#include "MSL_Common/float.h" +#include "MSL_Common_Embedded/Math/Double_precision/fdlibm.h" #ifdef __cplusplus extern "C" { #endif -extern int __float_nan[]; -extern int __float_huge[]; -extern int __double_huge[]; +#ifndef __MWERKS__ +/* Get clangd to shut up about __fabs being undefined. */ +#define __fabs(x) (x) +#define __frsqrte(x) (x) +#endif -#define INFINITY (*(float *)__float_huge) -#define NAN (*(float *)__float_nan) -#define HUGE_VAL (*(double *)__double_huge) /* #define HUGE_VALF #define HUGE_VALL @@ -51,10 +52,9 @@ int __fpclassifyd(double); #define MATH_ERRNO 1 #define MATH_ERREXCEPT 2 +#define math_errhandling (MATH_ERRNO | MATH_ERREXCEPT) /* -#define math_errhandling - #define FP_FAST_FMA #define FP_FAST_FMAF #define FP_FAST_FMAL @@ -62,141 +62,12 @@ int __fpclassifyd(double); #define FP_ILOGBNAN */ -#define _DECL_MATH(name) \ - float name##f(float x); \ - double name(double x); \ - long double name##l(long double x) - -#define _DECL_MATHC(name, arg) \ - float name##f(arg); \ - double name(arg); \ - long double name##l(arg) - -#define _DECL_MATHR(name, ret) \ - ret name##f(float x); \ - ret name(double x); \ - ret name##l(long double x) - -#define _DECL_MATH2(name) \ - float name##f(float x, float y); \ - double name(double x, double y); \ - long double name##l(long double x, long double y) - -#define _DECL_MATH2C(name, arg2) \ - float name##f(float x, arg2); \ - double name(double x, arg2); \ - long double name##l(long double x, arg2) - -#define _DECL_MATH3(name) \ - float name##f(float x, float y, float z); \ - double name(double x, double y, double z); \ - long double name##l(long double x, long double y, long double z) - -#define _DECL_MATH3C(name, arg3) \ - float name##f(float x, float y, arg3); \ - double name(double x, double y, arg3); \ - long double name##l(long double x, long double y, arg3) - -_DECL_MATH(fabs); -_DECL_MATH2(fmod); -_DECL_MATH2(remainder); -_DECL_MATH3C(remquo, int *quo); -_DECL_MATH3(fma); -_DECL_MATH2(fmax); -_DECL_MATH2(fmin); -_DECL_MATH2(fdim); - -_DECL_MATHC(nan, const char *str); - -_DECL_MATH(exp); -_DECL_MATH(exp2); -_DECL_MATH(expm1); -_DECL_MATH(log); -_DECL_MATH(log10); -_DECL_MATH(log2); -_DECL_MATH(log1p); - -_DECL_MATH2(pow); -_DECL_MATH(sqrt); -_DECL_MATH(cbrt); -_DECL_MATH2(hypot); - -_DECL_MATH(sin); -_DECL_MATH(cos); -_DECL_MATH(tan); -_DECL_MATH(asin); -_DECL_MATH(acos); -_DECL_MATH(atan); -_DECL_MATH2(atan2); - -_DECL_MATH(sinh); -_DECL_MATH(cosh); -_DECL_MATH(tanh); -_DECL_MATH(asinh); -_DECL_MATH(acosh); -_DECL_MATH(atanh); - -_DECL_MATH(erf); -_DECL_MATH(erfc); -_DECL_MATH(tgamma); -_DECL_MATH(lgamma); - -_DECL_MATH(ceil); -_DECL_MATH(floor); -_DECL_MATH(trunc); -_DECL_MATH(round); -_DECL_MATHR(lround, long); -_DECL_MATHR(llround, long long); - -_DECL_MATH(nearbyint); -_DECL_MATH(rint); -_DECL_MATHR(lrint, long); -_DECL_MATHR(llrint, long long); - -_DECL_MATH2C(frexp, int *exp); -_DECL_MATH2C(ldexp, int exp); - -float modff(float arg, float *iptr); -double modf(double arg, double *iptr); -long double modfl(long double arg, long double *iptr); - -_DECL_MATH2C(scalbn, int exp); -_DECL_MATH2C(scalbln, long exp); -_DECL_MATHR(ilogb, int); -_DECL_MATH(logb); - -_DECL_MATH2(nextafter); -_DECL_MATH2C(nexttoward, long double to); - -_DECL_MATH2(copysign); - -#undef _DECL_MATH -#undef _DECL_MATHC -#undef _DECL_MATHR -#undef _DECL_MATH2 -#undef _DECL_MATH2C -#undef _DECL_MATH3 -#undef _DECL_MATH3C +double nan(const char *str); inline long double fabsl(long double x) { return __fabs((double)x); } -double __ieee754_acos(double); -double __ieee754_fmod(double, double); -double __ieee754_log(double); -double __ieee754_log10(double); -double __ieee754_pow(double, double); -double __ieee754_sqrt(double); -double __ieee754_atan2(double, double); -double __ieee754_asin(double); - -int __ieee754_rem_pio2(double, double *); - -double __kernel_sin(double, double, int); -double __kernel_cos(double, double); -double __kernel_tan(double, double, int); - #ifdef __cplusplus } #endif diff --git a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/stdbool.h b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/stdbool.h index a7fbf7e3..7447f332 100644 --- a/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/stdbool.h +++ b/src/sdk/PowerPC_EABI_Support/MSL/MSL_C/stdbool.h @@ -1,13 +1,6 @@ #ifndef _STDBOOL_H_ #define _STDBOOL_H_ -#ifndef __cplusplus - -typedef unsigned char bool; - -#define false 0 -#define true 1 - -#endif +#include "MSL_Common/bool_def.h" #endif diff --git a/src/sdk/RVL_SDK/revolution/math.h b/src/sdk/RVL_SDK/revolution/math.h index 9b3a679c..d326f918 100644 --- a/src/sdk/RVL_SDK/revolution/math.h +++ b/src/sdk/RVL_SDK/revolution/math.h @@ -2,19 +2,12 @@ #define _REVOSDK_MATH_H #include "types.h" -#include "MSL_C/MSL_Common_Embedded/Math/fdlibm.h" -#include "MSL_C/MSL_Common/math_api.h" +#include "PowerPC_EABI_Support/MSL/MSL_C/math.h" #ifdef __cplusplus extern "C" { #endif // ifdef __cplusplus -#ifndef __MWERKS__ -// Get clangd to shut up about __fabs being undefined. -#define __fabs(x) (x) -#define __frsqrte(x) (x) -#endif - #define FABS(x) (float)__fabs(x) #define fabs(x) __fabs(x) // #define __frsqrtes opword @@ -27,6 +20,10 @@ extern "C" { #define LONG_TAU 6.2831854820251465 +#ifndef M_PI +#define M_PI ((float)3.14159274101257324219) +#endif + inline float cosf(float x){ return (float)cos(x); }