update from dtk-template - clangd :) (#66)

* update from dtk-template and start work towards using clangd

* include <a> -> "a"

* Update build.yml

* remove/add non-trivial class in union warning
This commit is contained in:
Elijah Thomas
2024-10-16 15:36:02 -04:00
committed by GitHub
parent 4543ec129e
commit 26af4db82d
2230 changed files with 18282 additions and 8154 deletions
+5 -5
View File
@@ -1,4 +1,4 @@
#include <DynamicLink.h>
#include "DynamicLink.h"
DynamicModuleControlBase *DynamicModuleControlBase::mFirst;
DynamicModuleControlBase *DynamicModuleControlBase::mLast;
@@ -311,10 +311,10 @@ int DynamicModuleControl::getModuleSize() const {
const char *DynamicModuleControl::getModuleTypeString() const {
static const char *REL_LOAD_TYPES[4] = {
"????",
"MEM",
"ARAM",
"DVD",
"????",
"MEM",
"ARAM",
"DVD",
};
return REL_LOAD_TYPES[mResourceType & 3];
}
@@ -0,0 +1,65 @@
#ifndef MSL_ALGORITHM_H_
#define MSL_ALGORITHM_H_
namespace std {
template <class ForwardIterator, class T>
ForwardIterator lower_bound(ForwardIterator first, ForwardIterator last, const T &val);
template <class ForwardIterator, class T>
ForwardIterator upper_bound(ForwardIterator first, ForwardIterator last, const T &val);
template <class InputIt, class UnaryPredicate>
InputIt find_if(InputIt first, InputIt last, UnaryPredicate p);
/*
template<class OutputIt, class Size, int A2>
struct __fill_n {
OutputIt fill_n(OutputIt first, Size count, const unsigned long& value);
};
template<>
unsigned long* __fill_n<unsigned long, long, 0>::fill_n(unsigned long* first, long count, const unsigned long& value) {
for (; count > 0; count--) {
*first++ = value;
}
return first;
}
template<class OutputIt, class Size, class T>
OutputIt fill_n(OutputIt first, Size count, const T& value) {
return __fill_n::fill_n(first, count, value);
}
template<class ForwardIt, class T>
void __fill(ForwardIt first, ForwardIt last, const T& value, std::random_access_iterator_tag param_3) {
fill_n(first, last - first, value);
}
*/
template <class ForwardIt, class T>
void fill(ForwardIt first, ForwardIt last, const T &value) {
for (; first != last; ++first) {
*first = value;
}
}
template <class InputIt, class OutputIt>
inline OutputIt copy(InputIt first, InputIt last, OutputIt d_first) {
for (; first < last; ++first, ++d_first) {
*d_first = *first;
}
return d_first;
}
template <class BidirIt1, class BidirIt2>
inline BidirIt2 copy_backward(BidirIt1 first, BidirIt1 last, BidirIt2 d_last) {
while (first != last) {
*(--d_last) = *(--last);
}
return d_last;
}
} // namespace std
#endif
@@ -0,0 +1,89 @@
#ifndef MSL_BITSET_H_
#define MSL_BITSET_H_
#include "algorithm.h"
#include "stdio.h"
#include "stdlib.h"
namespace std {
// TODO: where does this go?
inline void __msl_error(const char *param_0) {
fprintf(stderr, param_0);
abort();
}
template <size_t N>
class __bitset_base {
public:
__bitset_base();
bool test(size_t pos) const;
bool any() const;
void set(size_t pos, bool val);
void reset(size_t pos);
private:
size_t data[N];
};
template <size_t N>
__bitset_base<N>::__bitset_base() {
std::fill(data, data + N, 0);
}
template <size_t N>
bool __bitset_base<N>::test(size_t pos) const {
size_t i = pos / (sizeof(size_t) * 8);
size_t mask = 1 << (pos % (sizeof(size_t) * 8));
return data[i] & mask;
}
template <size_t N>
void __bitset_base<N>::set(size_t pos, bool val) {
size_t i = pos / (sizeof(size_t) * 8);
size_t mask = 1 << (pos % (sizeof(size_t) * 8));
if (val) {
data[i] |= mask;
} else {
data[i] &= ~mask;
}
}
template <size_t N>
void __bitset_base<N>::reset(size_t pos) {
size_t i = pos / (sizeof(size_t) * 8);
size_t mask = 1 << (pos % (sizeof(size_t) * 8));
data[i] &= ~mask;
}
template <size_t N>
class bitset : private __bitset_base<(N - 1) / (sizeof(size_t) * 8) + 1> {
public:
typedef __bitset_base<(N - 1) / (sizeof(size_t) * 8) + 1> base;
bitset(){};
void set(size_t pos, bool val) {
if (pos >= N) {
__msl_error("index out of range of bitset::set");
}
base::set(pos, val);
}
void reset(size_t pos) {
if (pos >= N) {
__msl_error("index out of range of bitset::reset");
}
base::reset(pos);
}
bool test(size_t pos) const {
if (pos >= N) {
__msl_error("index out of range of bitset::test");
}
return base::test(pos);
}
bool any() const;
};
} // namespace std
#endif
@@ -0,0 +1,8 @@
#ifndef MSL_FUNCTIONAL_H_
#define MSL_FUNCTIONAL_H_
namespace std {
template <class T> struct less {};
} // namespace std
#endif
@@ -0,0 +1,21 @@
#ifndef MSL_ITERATOR_H_
#define MSL_ITERATOR_H_
namespace std {
template< class InputIt, class Distance >
inline void advance( InputIt& it, Distance n) {
while (n > 0) {
--n;
++it;
}
}
// This needs to be defined with gcc concepts or something similar. Workaround.
template< class InputIt, class Distance >
inline void advance_pointer( InputIt& it, Distance n) {
it += n;
}
}
#endif
@@ -0,0 +1,24 @@
#ifndef MSL_MEMORY_H_
#define MSL_MEMORY_H_
namespace std {
template<class ForwardIt, class Size, class T>
inline ForwardIt uninitialized_fill_n(ForwardIt first, Size count, const T& value) {
for (; count > 0; ++first, (void) --count) {
*first = value;
}
return first;
}
template<class InputIt, class NoThrowForwardIt>
inline NoThrowForwardIt uninitialized_copy(InputIt first, InputIt last, NoThrowForwardIt d_first) {
for (; first != last; ++first, ++d_first) {
*d_first = *first;
}
return d_first;
}
}
#endif
@@ -0,0 +1,17 @@
#ifndef MSL_UTILITY_H_
#define MSL_UTILITY_H_
namespace std {
template <class T1, class T2>
struct pair {
T1 first;
T2 second;
pair() {
first = T1();
second = T2();
}
};
} // namespace std
#endif
@@ -0,0 +1,18 @@
#ifndef _MSL_COMMON_FILE_POS_H
#define _MSL_COMMON_FILE_POS_H
#include "ansi_files.h"
#ifdef __cplusplus
extern "C" {
#endif
int fseek(FILE* file, unsigned long offset, int mode);
int _fseek(FILE* file, fpos_t offset, int mode);
long ftell(FILE* file);
#ifdef __cplusplus
}
#endif
#endif /* _MSL_COMMON_FILE_POS_H */
@@ -0,0 +1,17 @@
#ifndef _MSL_COMMON_ABORT_EXIT_H
#define _MSL_COMMON_ABORT_EXIT_H
#ifdef __cplusplus
extern "C" {
#endif
void exit(int status);
void abort(void);
extern void (*__stdio_exit)(void);
#ifdef __cplusplus
};
#endif
#endif /* _MSL_COMMON_ABORT_EXIT_H */
@@ -0,0 +1,15 @@
#ifndef _MSL_COMMON_ALLOC_H
#define _MSL_COMMON_ALLOC_H
#ifdef __cplusplus
extern "C" {
#endif
void free(void* ptr);
#ifdef __cplusplus
}
#endif
#endif /* _MSL_COMMON_ALLOC_H */
@@ -0,0 +1,129 @@
#ifndef _MSL_COMMON_ANSI_FILES_H
#define _MSL_COMMON_ANSI_FILES_H
#include "stddef.h"
#ifdef __cplusplus
extern "C" {
#endif
#define SEEK_SET 0
#define SEEK_CUR 1
#define SEEK_END 2
typedef unsigned long __file_handle;
typedef unsigned long fpos_t;
#ifndef __cplusplus
typedef unsigned short wchar_t;
#endif
#define set_error(file) \
do { \
(file)->file_state.error = 1; \
(file)->buffer_length = 0; \
} while (0)
enum __file_kinds {
__closed_file,
__disk_file,
__console_file,
__unavailable_file,
};
enum __file_orientation {
/* 0x0 */ UNORIENTED,
/* 0x1 */ CHAR_ORIENTED,
/* 0x2 */ WIDE_ORIENTED,
};
typedef struct _file_modes {
unsigned int open_mode : 2;
unsigned int io_mode : 3;
unsigned int buffer_mode : 2;
unsigned int file_kind : 3;
unsigned int file_orientation : 2;
unsigned int binary_io : 1;
} file_modes;
enum __io_modes {
__read = 1,
__write = 2,
__read_write = 3,
__append = 4,
};
enum __io_states {
__neutral,
__writing,
__reading,
__rereading,
};
enum __io_results {
__no_io_error,
__io_error,
__io_EOF,
};
typedef struct _file_states {
unsigned int io_state : 3;
unsigned int free_buffer : 1;
unsigned char eof;
unsigned char error;
} file_states;
typedef void (*__idle_proc)(void);
typedef int (*__pos_proc)(__file_handle file, fpos_t* position, int mode, __idle_proc idle_proc);
typedef int (*__io_proc)(__file_handle file, unsigned char* buff, size_t* count,
__idle_proc idle_proc);
typedef int (*__close_proc)(__file_handle file);
typedef struct _FILE {
/* 0x00 */ __file_handle handle;
/* 0x04 */ file_modes file_mode;
/* 0x08 */ file_states file_state;
/* 0x0C */ unsigned char is_dynamically_allocated;
/* 0x0D */ char char_buffer;
/* 0x0E */ char char_buffer_overflow;
/* 0x0F */ char ungetc_buffer[2];
/* 0x12 */ wchar_t ungetc_wide_buffer[2];
/* 0x18 */ unsigned long position;
/* 0x1C */ unsigned char* buffer;
/* 0x20 */ unsigned long buffer_size;
/* 0x24 */ unsigned char* buffer_ptr;
/* 0x28 */ unsigned long buffer_length;
/* 0x2C */ unsigned long buffer_alignment;
/* 0x30 */ unsigned long save_buffer_length;
/* 0x34 */ unsigned long buffer_position;
/* 0x38 */ __pos_proc position_fn;
/* 0x3C */ __io_proc read_fn;
/* 0x40 */ __io_proc write_fn;
/* 0x44 */ __close_proc close_fn;
/* 0x48 */ __idle_proc idle_fn;
/* 0x4C */ struct _FILE* next_file;
} FILE;
typedef struct _files {
FILE _stdin;
FILE _stdout;
FILE _stderr;
FILE empty;
} files;
#define _IONBF 0
#define _IOLBF 1
#define _IOFBF 2
extern files __files;
extern int __close_console(__file_handle file);
extern int __write_console(__file_handle file, unsigned char* buf, size_t* count, __idle_proc idle_fn);
extern int __read_console(__file_handle file, unsigned char* buf, size_t* count, __idle_proc idle_fn);
unsigned int __flush_all(void);
void __close_all(void);
#ifdef __cplusplus
};
#endif
#endif /* _MSL_COMMON_ANSI_FILES_H */
@@ -0,0 +1,35 @@
#ifndef _MSL_COMMON_ANSI_FP_H
#define _MSL_COMMON_ANSI_FP_H
#define SIGDIGLEN 36
typedef struct decimal {
char sign;
char unk1;
short exp;
struct {
unsigned char length;
unsigned char text[36];
unsigned char unk41;
} sig;
} decimal;
typedef struct decform {
char style;
char unk1;
short digits;
} decform;
/* void __ull2dec(decimal*, u64);
void __timesdec(decimal*, const decimal*, const decimal*);
void __str2dec(decimal*, const char*, short);
void __two_exp(decimal*, s32);
BOOL __equals_dec(const decimal*, const decimal*);
BOOL __less_dec(const decimal*, const decimal*);
void __minus_dec(decimal*, const decimal*, const decimal*);
void __num2dec_internal(decimal*, f64);
void __num2dec(const decform*, f64, decimal*);
f64 __dec2num(const decimal*); */
#endif
@@ -0,0 +1,19 @@
#ifndef _MSL_COMMON_ARITH_H
#define _MSL_COMMON_ARITH_H
#ifdef __cplusplus
extern "C" {
#endif
typedef struct {
int quot; /* quotient */
int rem; /* remainder */
} div_t;
div_t div(int numerator, int denominator);
#ifdef __cplusplus
}
#endif
#endif /* _MSL_COMMON_ARITH_H */
@@ -0,0 +1,11 @@
#ifndef _MSL_COMMON_BUFFER_IO_H
#define _MSL_COMMON_BUFFER_IO_H
#include "ansi_files.h"
enum { __align_buffer, __dont_align_buffer };
void __prep_buffer(FILE* file);
int __flush_buffer(FILE* file, size_t* bytes_flushed);
#endif /* _MSL_COMMON_BUFFER_IO_H */
@@ -0,0 +1,17 @@
#ifndef _MSL_COMMON_CHAR_IO_H
#define _MSL_COMMON_CHAR_IO_H
#include "ansi_files.h"
#ifdef __cplusplus
extern "C" {
#endif
int fputs(const char *str, FILE *stream);
int __put_char(int c, FILE *stream);
#ifdef __cplusplus
}
#endif
#endif /* _MSL_COMMON_CHAR_IO_H */
@@ -0,0 +1,23 @@
#ifndef STD_CMATH_H_
#define STD_CMATH_H_
namespace std {
/* float fabs(float num) {
return ::fabsf(num);
} */
inline float fabsf(float num) {
return ::fabsf(num);
}
inline float sqrt(float x) {
return ::sqrtf(x);
}
inline float abs(float x) {
return ::fabsf(x);
}
} // namespace std
#endif
@@ -0,0 +1,30 @@
#ifndef _MSL_COMMON_CRITICAL_REGIONS_H
#define _MSL_COMMON_CRITICAL_REGIONS_H
#ifdef __cplusplus
extern "C" {
#endif
enum critical_regions {
atexit_funcs_access,
malloc_pool_access,
stdin_access,
stdout_access,
stderr_access,
files_access,
console_status_access,
signal_funcs_access,
thread_access,
num_critical_regions
};
void __init_critical_regions(void);
void __kill_critical_regions(void);
void __begin_critical_region(int region);
void __end_critical_region(int region);
#ifdef __cplusplus
}
#endif
#endif
@@ -0,0 +1,16 @@
#ifndef STD_CMATH_H_
#define STD_CMATH_H_
#include "string.h"
namespace std {
inline size_t strlen(const char* str) {
return ::strlen(str);
}
inline char* strcpy(char* dest, const char* src) {
return ::strcpy(dest, src);
}
} // namespace std
#endif
@@ -0,0 +1,46 @@
#ifndef _MSL_COMMON_CTYPE_H
#define _MSL_COMMON_CTYPE_H
#ifdef __cplusplus
extern "C" {
#endif
#define EOF -1L
extern unsigned char __ctype_map[];
extern unsigned char __lower_map[];
extern unsigned char __upper_map[];
#define __control_char 0x01
#define __motion_char 0x02
#define __space_char 0x04
#define __punctuation 0x08
#define __digit 0x10
#define __hex_digit 0x20
#define __lower_case 0x40
#define __upper_case 0x80
#define __letter (__lower_case | __upper_case)
#define __alphanumeric (__letter | __digit)
#define __graphic (__alphanumeric | __punctuation)
#define __printable (__graphic | __space_char)
#define __whitespace (__motion_char | __space_char)
#define __control (__motion_char | __control_char)
#define __zero_fill(c) ((int)(unsigned char)(c))
int tolower(int);
inline int isalpha(int c) { return (int)(__ctype_map[(unsigned char)c] & __letter); }
inline int isdigit(int c) { return (int)(__ctype_map[(unsigned char)c] & __digit); }
inline int isspace(int c) { return (int)(__ctype_map[(unsigned char)c] & __whitespace); }
inline int isupper(int c) { return (int)(__ctype_map[(unsigned char)c] & __upper_case); }
inline int isxdigit(int c) { return (int)(__ctype_map[(unsigned char)c] & __hex_digit); }
// added underscore to avoid naming conflicts
inline int _tolower(int c) { return (c == -1 ? -1 : (int)__lower_map[(unsigned char)c]); }
inline int toupper(int c) { return (c == -1 ? -1 : (int)__upper_map[(unsigned char)c]); }
#ifdef __cplusplus
}
#endif
#endif /* _MSL_COMMON_CTYPE_H */
@@ -0,0 +1,17 @@
#ifndef _MSL_COMMON_DIRECT_IO_H
#define _MSL_COMMON_DIRECT_IO_H
#include "ansi_files.h"
#ifdef __cplusplus
extern "C" {
#endif
size_t __fwrite(const void* buffer, size_t size, size_t count, FILE* stream);
size_t fwrite(const void* buffer, size_t size, size_t count, FILE* stream);
#ifdef __cplusplus
}
#endif
#endif /* _MSL_COMMON_DIRECT_IO_H */
@@ -0,0 +1,20 @@
#ifndef MSL_COMMON_SRC_ERRNO_H
#define MSL_COMMON_SRC_ERRNO_H
#ifdef __cplusplus
extern "C" {
#endif
#define ENOERR 0
#define EDOM 33
#define ERANGE 34
#define EFPOS 40
#define ESIGPARM 36
extern int errno;
#ifdef __cplusplus
}
#endif
#endif /* MSL_COMMON_SRC_ERRNO_H */
@@ -0,0 +1,15 @@
#ifndef _MSL_COMMON_EXTRAS_H
#define _MSL_COMMON_EXTRAS_H
#ifdef __cplusplus
extern "C" {
#endif
int strnicmp(const char* str1, const char* str2, int n);
int stricmp(const char* str1, const char* str2);
#ifdef __cplusplus
}
#endif
#endif /* _MSL_COMMON_EXTRAS_H */
@@ -0,0 +1,18 @@
#ifndef _MSL_COMMON_FILE_IO_H
#define _MSL_COMMON_FILE_IO_H
#include "ansi_files.h"
#ifdef __cplusplus
extern "C" {
#endif
int __msl_strnicmp(const char* str1, const char* str2, int n);
int fflush(FILE* file);
int fclose(FILE* file);
#ifdef __cplusplus
}
#endif
#endif /* _MSL_COMMON_FILE_IO_H */
@@ -0,0 +1,91 @@
#ifndef _MSL_COMMON_FLOAT_H
#define _MSL_COMMON_FLOAT_H
#include "fdlibm.h"
#define FP_SNAN 0
#define FP_QNAN 1
#define FP_INFINITE 2
#define FP_ZERO 3
#define FP_NORMAL 4
#define FP_SUBNORMAL 5
#define FP_NAN FP_QNAN
#define fpclassify(x) ((sizeof(x) == sizeof(float)) ? __fpclassifyf(x) : __fpclassifyd(x))
#define signbit(x) ((sizeof(x) == sizeof(float)) ? __signbitf(x) : __signbitd(x))
#define isfinite(x) ((fpclassify(x) > 2))
#define __signbitf(x) ((int)(__HI(x) & 0x80000000))
// TODO: OK?
#define __signbitd(x) ((int)(__HI(x) & 0x80000000))
extern unsigned long __float_nan[];
extern unsigned long __float_huge[];
extern unsigned long __float_max[];
extern unsigned long __float_epsilon[];
inline int __fpclassifyf(float __value) {
unsigned long integer = *(unsigned long *)&__value;
switch (integer & 0x7f800000) {
case 0x7f800000:
if ((integer & 0x7fffff) != 0) {
return FP_QNAN;
}
return FP_INFINITE;
case 0:
if ((integer & 0x7fffff) != 0) {
return FP_SUBNORMAL;
}
return FP_ZERO;
}
return FP_NORMAL;
}
inline int __fpclassifyd(double __value) {
switch (__HI(__value) & 0x7ff00000) {
case 0x7ff00000: {
if ((__HI(__value) & 0x000fffff) || (__LO(__value) & 0xffffffff)) {
return FP_QNAN;
} else {
return FP_INFINITE;
}
break;
}
case 0: {
if ((__HI(__value) & 0x000fffff) || (__LO(__value) & 0xffffffff)) {
return FP_SUBNORMAL;
} else {
return FP_ZERO;
}
break;
}
}
return FP_NORMAL;
}
#define FLT_MANT_DIG 24
#define FLT_DIG 6
#define FLT_MIN_EXP (-125)
#define FLT_MIN_10_EXP (-37)
#define FLT_MAX_EXP 128
#define FLT_MAX_10_EXP 38
// #define FLT_MAX (*(float*) __float_max)
// #define FLT_EPSILON (*(float*) __float_epsilon)
#define FLT_MAX (3.402823466e+38f)
#define FLT_EPSILON (1.192092896e-07f)
#define DBL_MANT_DIG 53
#define DBL_DIG 15
#define DBL_MIN_EXP (-1021)
#define DBL_MIN_10_EXP (-308)
#define DBL_MAX_EXP 1024
#define DBL_MAX_10_EXP 308
#endif /* _MSL_COMMON_FLOAT_H */
@@ -0,0 +1,102 @@
#ifndef _STD_LIMITS_H
#define _STD_LIMITS_H
#ifdef __cplusplus
extern "C" {
#endif
#define CHAR_BIT 8
#define SCHAR_MIN (-0x7F - 1)
#define SCHAR_MAX 0x7F
#define UCHAR_MAX 0xFF
#define CHAR_MIN 0
#define CHAR_MAX SCHAR_MAX
#define SHRT_MIN (-0x7FFF - 1)
#define SHRT_MAX 0x7FFF
#define USHRT_MAX 0xFFFF
#define INT_MIN (-0x7FFFFFFF - 1)
#define INT_MAX 0x7FFFFFFF
#define UINT_MAX 0xFFFFFFFF
#define LONG_MIN (-0x7FFFFFFFL - 1)
#define LONG_MAX 0x7FFFFFFFL
#define ULONG_MAX 0xFFFFFFFFUL
#define LLONG_MIN (-0x7FFFFFFFFFFFFFFFLL - 1)
#define LLONG_MAX 0x7FFFFFFFFFFFFFFFLL
#define ULLONG_MAX 0xFFFFFFFFFFFFFFFFULL
#ifdef __cplusplus
}
namespace std {
template <typename T>
class numeric_limits {
public:
inline static T min();
inline static T max();
};
template <>
class numeric_limits<char> {
public:
inline static char min() { return -0x80; }
inline static char max() { return 0x7F; }
};
template <>
class numeric_limits<short> {
public:
inline static short min() { return -0x8000; }
inline static short max() { return 0x7FFF; }
};
template <>
class numeric_limits<int> {
public:
inline static int min() { return -0x80000000; }
inline static int max() { return 0x7FFFFFFF; }
};
template <>
class numeric_limits<long> {
public:
inline static long min() { return -0x80000000; }
inline static long max() { return 0x7FFFFFFF; }
};
template <>
class numeric_limits<unsigned char> {
public:
inline static unsigned char min() { return 0x0; }
inline static unsigned char max() { return 0xFF; }
};
template <>
class numeric_limits<unsigned short> {
public:
inline static unsigned short min() { return 0x0; }
inline static unsigned short max() { return 0xFFFF; }
};
template <>
class numeric_limits<unsigned int> {
public:
inline static unsigned int min() { return 0x0; }
inline static unsigned int max() { return 0xFFFFFFFF; }
};
template <>
class numeric_limits<unsigned long> {
public:
inline static unsigned long min() { return 0x0; }
inline static unsigned long max() { return 0xFFFFFFFF; }
};
} // namespace std
#endif
#endif
@@ -0,0 +1,114 @@
#ifndef MSL_MATH_H_
#define MSL_MATH_H_
#include "float.h"
#define NAN (*(float *)__float_nan)
#define HUGE_VALF (*(float *)__float_huge)
#define M_PI 3.14159265358979323846f
#define M_SQRT3 1.73205f
#define DEG_TO_RAD(degrees) (degrees * (M_PI / 180.0f))
#define RAD_TO_DEG(radians) (radians / (180.0f / M_PI))
#ifdef __cplusplus
extern "C" {
#endif
int abs(int);
int labs(int);
double acos(double);
float acosf(float);
double asin(double);
double atan(double);
double atan2(double, double);
double ceil(double);
double copysign(double, double);
double cos(double);
float cosf(float);
double exp(double);
double ceil(double);
float ceilf(float);
double frexp(double, int *);
double ldexp(double, int);
double modf(double, double *);
double pow(double, double);
double sin(double);
float sinf(float);
double tan(double);
float tanf(float);
double floor(double);
float floorf(float);
double fmod(double, double);
float fmodf(float, float);
float fmodff(float, float *);
extern float __fabsf(float);
inline double fabs(double f) { return __fabs(f); }
inline double fabsf2(float f) { return __fabsf(f); }
inline float fabsf(float f) { return fabsf2(f); }
inline float fmodf(float f1, float f2) { return fmod(f1, f2); }
inline float modff(float x, float *iptr) {
float frac;
double intg;
frac = modf((double)x, &intg);
*iptr = intg;
return frac;
}
inline double sqrt_step(double tmpd, float mag) {
return tmpd * 0.5 * (3.0 - mag * (tmpd * tmpd));
}
extern inline float sqrtf(float x) {
const double _half = .5;
const double _three = 3.0;
volatile float y;
if (x > 0.0f) {
double guess = __frsqrte((double)x); // returns an approximation to
guess =
_half * guess * (_three - guess * guess * x); // now have 12 sig bits
guess =
_half * guess * (_three - guess * guess * x); // now have 24 sig bits
guess =
_half * guess * (_three - guess * guess * x); // now have 32 sig bits
y = (float)(x * guess);
return y;
}
return x;
}
extern inline double sqrt(double x) {
if (x > 0.0) {
double guess = __frsqrte(x); /* returns an approximation to */
guess = .5 * guess * (3.0 - guess * guess * x); /* now have 8 sig bits */
guess = .5 * guess * (3.0 - guess * guess * x); /* now have 16 sig bits */
guess = .5 * guess * (3.0 - guess * guess * x); /* now have 32 sig bits */
guess = .5 * guess * (3.0 - guess * guess * x); /* now have > 53 sig bits */
return x * guess;
} else if (x == 0) {
return 0;
} else if (x) {
return NAN;
}
return HUGE_VALF;
}
inline float atan2f(float y, float x) { return (float)atan2(y, x); }
// these are duplicated due to sinf/cosf having a symbol, but
// still being used as inlines elsewhere
inline float i_sinf(float x) { return sin(x); }
inline float i_cosf(float x) { return cos(x); }
inline float i_tanf(float x) { return tanf(x); }
#ifdef __cplusplus
};
#endif
#endif
@@ -0,0 +1,16 @@
#ifndef _MSL_COMMON_MBSTRING_H
#define _MSL_COMMON_MBSTRING_H
#include "wchar_io.h"
#ifdef __cplusplus
extern "C" {
#endif
size_t wcstombs(char* dst, const wchar_t* src, size_t n);
#ifdef __cplusplus
}
#endif
#endif /* _MSL_COMMON_MBSTRING_H */
@@ -0,0 +1,19 @@
#ifndef _MSL_COMMON_MEM_FUNCS_H
#define _MSL_COMMON_MEM_FUNCS_H
#include "stddef.h"
#ifdef __cplusplus
extern "C" {
#endif
void __copy_longs_rev_unaligned(void* dst, const void* src, size_t n);
void __copy_longs_unaligned(void* dst, const void* src, size_t n);
void __copy_longs_rev_aligned(void* dst, const void* src, size_t n);
void __copy_longs_aligned(void* dst, const void* src, size_t n);
#ifdef __cplusplus
}
#endif
#endif /* _MSL_COMMON_MEM_FUNCS_H */
@@ -0,0 +1,14 @@
#ifndef _MSL_COMMON_MISC_IO_H
#define _MSL_COMMON_MISC_IO_H
#ifdef __cplusplus
extern "C" {
#endif
void __stdio_atexit(void);
#ifdef __cplusplus
}
#endif
#endif /* _MSL_COMMON_MISC_IO_H */
@@ -0,0 +1,10 @@
#ifndef MSL_NEW_H_
#define MSL_NEW_H_
#include "stddef.h"
inline void *operator new(size_t size, void *ptr) {
return ptr;
}
#endif
@@ -0,0 +1,22 @@
#ifndef _MSL_COMMON_PRINTF_H
#define _MSL_COMMON_PRINTF_H
#include "__va_arg.h"
#include "ansi_files.h"
#ifdef __cplusplus
extern "C" {
#endif
int fprintf(FILE *stream, const char *format, ...);
int printf(const char *format, ...);
int sprintf(char *s, const char *format, ...);
int snprintf(char *s, size_t n, const char *format, ...);
int vsnprintf(char *s, size_t n, const char *format, va_list arg);
int vprintf(const char *format, va_list arg);
#ifdef __cplusplus
}
#endif
#endif /* _MSL_COMMON_PRINTF_H */
@@ -0,0 +1,42 @@
#ifndef _MSL_COMMON_SCANF_H
#define _MSL_COMMON_SCANF_H
#include "ansi_files.h"
#ifdef __cplusplus
extern "C" {
#endif
enum __ReadProcActions { __GetAChar, __UngetAChar, __TestForError };
enum __WReadProcActions { __GetAwChar, __UngetAwChar, __TestForwcsError };
typedef struct {
char* CharStr;
size_t MaxCharCount;
size_t CharsWritten;
} __OutStrCtrl;
typedef struct {
char* NextChar;
int NullCharDetected;
} __InStrCtrl;
typedef struct {
wchar_t* wCharStr;
size_t MaxCharCount;
size_t CharsWritten;
} __wOutStrCtrl;
typedef struct {
wchar_t* wNextChar;
int wNullCharDetected;
} __wInStrCtrl;
int __StringRead(void* str, int ch, int behavior);
#ifdef __cplusplus
}
#endif
#endif /* _MSL_COMMON_SCANF_H */
@@ -0,0 +1,16 @@
#ifndef _MSL_COMMON_SIGNAL_H
#define _MSL_COMMON_SIGNAL_H
#ifdef __cplusplus
extern "C" {
#endif
typedef void (*__signal_func_ptr)(int);
int raise(int sig);
#ifdef __cplusplus
}
#endif
#endif /* _MSL_COMMON_SIGNAL_H */
@@ -0,0 +1,26 @@
#ifndef _STDDEF_H_
#define _STDDEF_H_
#ifdef __cplusplus
extern "C" {
#endif
#if defined __INTELLISENSE__
typedef unsigned int size_t;
typedef int ptrdiff_t;
#else
typedef unsigned long size_t;
typedef long ptrdiff_t;
#endif
#define offsetof(type, member) ((size_t) & (((type *)0)->member))
#ifndef NULL
#define NULL (0)
#endif
#ifdef __cplusplus
};
#endif
#endif
@@ -0,0 +1,14 @@
#ifndef __STDC_STDINT_H__
#define __STDC_STDINT_H__
typedef unsigned char uint8_t;
typedef unsigned short uint16_t;
typedef unsigned long uint32_t;
typedef unsigned long long uint64_t;
typedef signed char int8_t;
typedef signed short int16_t;
typedef signed long int32_t;
typedef signed long long int64_t;
#endif // __STDC_STDINT_H__
@@ -0,0 +1,13 @@
#ifndef MSL_STDIO_H_
#define MSL_STDIO_H_
#include "char_io.h"
#include "extras.h"
#include "file_io.h"
#include "printf.h"
#define stdin (&__files._stdin)
#define stdout (&__files._stdout)
#define stderr (&__files._stderr)
#endif
@@ -0,0 +1,8 @@
#ifndef MSL_STDLIB_H_
#define MSL_STDLIB_H_
#include "abort_exit.h"
#include "arith.h"
#include "strtoul.h"
#endif
@@ -0,0 +1,30 @@
#ifndef _MSL_COMMON_STRING_H
#define _MSL_COMMON_STRING_H
#include "stddef.h"
#ifdef __cplusplus
extern "C" {
#endif
void *memcpy(void *dst, const void *src, size_t n);
void *memset(void *dst, int val, size_t n);
int memcmp(const void *lhs, const void *rhs, size_t count);
void *__memrchr(const void *ptr, int ch, size_t count);
void *memchr(const void *ptr, int ch, size_t count);
void *memmove(void *dst, const void *src, size_t n);
char *strrchr(const char *str, int c);
char *strchr(const char *str, int c);
int strncmp(const char *str1, const char *str2, size_t n);
int strcmp(const char *str1, const char *str2);
char *strcat(char *dst, const char *src);
char *strncpy(char *dst, const char *src, size_t n);
char *strcpy(char *dst, const char *src);
int sscanf(const char *buffer, const char *format, ...);
size_t strlen(const char *str);
#ifdef __cplusplus
}
#endif
#endif /* _MSL_COMMON_STRING_H */
@@ -0,0 +1,17 @@
#ifndef _MSL_COMMON_STRTOUL_H
#define _MSL_COMMON_STRTOUL_H
#ifdef __cplusplus
extern "C" {
#endif
long strtol(const char* str, char** end, int base);
unsigned long strtoul(const char* str, char** end, int base);
unsigned long __strtoul(int base, int max_width, int (*ReadProc)(void*, int, int), void* ReadProcArg, int* chars_scanned, int* negative,
int* overflow);
#ifdef __cplusplus
}
#endif
#endif /* _MSL_COMMON_STRTOUL_H */
@@ -0,0 +1,19 @@
#ifndef MSL_WCHAR_H
#define MSL_WCHAR_H
#include "common.h"
#include "limits.h"
#include "mbstring.h"
#include "wprintf.h"
#include "wstring.h"
#ifdef __cplusplus
extern "C" {
#endif
#define WCHAR_MIN SHRT_MIN
#define WCHAR_MAX USHRT_MAX
#ifdef __cplusplus
}
#endif
#endif
@@ -0,0 +1,12 @@
#ifndef _MSL_COMMON_WCHAR_IO_H
#define _MSL_COMMON_WCHAR_IO_H
#include "ansi_files.h"
#ifndef __cplusplus
typedef unsigned short wchar_t;
#endif
int fwide(FILE* file, int mode);
#endif /* _MSL_COMMON_WCHAR_IO_H */
@@ -0,0 +1,17 @@
#ifndef MSL_WPRINTF_H
#define MSL_WPRINTF_H
#include "__va_arg.h"
#include "common.h"
#ifdef __cplusplus
extern "C" {
#endif
int vswprintf(wchar_t *, size_t, const wchar_t *, va_list);
int swprintf(wchar_t *, size_t, const wchar_t *, ...);
#ifdef __cplusplus
}
#endif
#endif
@@ -0,0 +1,18 @@
#ifndef MSL_WSTRING_H
#define MSL_WSTRING_H
#include "common.h"
#ifdef __cplusplus
extern "C" {
#endif
size_t wcslen(const wchar_t *);
wchar_t *wcscpy(wchar_t *, const wchar_t *);
wchar_t *wcsncpy(wchar_t *, const wchar_t *, size_t);
wchar_t *wcscat(wchar_t *, const wchar_t *);
int wcscmp(const wchar_t *, const wchar_t *);
wchar_t *wcschr(const wchar_t *, wchar_t);
#ifdef __cplusplus
}
#endif
#endif
@@ -0,0 +1,147 @@
#include "FILE_POS.h"
#include "critical_regions.h"
#include "errno.h"
/* 80365E90-80365F74 3607D0 00E4+00 0/0 2/2 0/0 .text ftell */
long ftell(FILE* stream) {
int retval;
__begin_critical_region(stdin_access);
retval = (long)_ftell(stream);
__end_critical_region(stdin_access);
return retval;
}
int _ftell(FILE* file) {
int charsInUndoBuffer = 0;
int position;
unsigned char tmp_kind = file->file_mode.file_kind;
if (!(tmp_kind == __disk_file || tmp_kind == __console_file) || file->file_state.error) {
errno = 0x28;
return -1;
}
if (file->file_state.io_state == __neutral)
return (file->position);
position = file->buffer_position + (file->buffer_ptr - file->buffer);
if (file->file_state.io_state >= __rereading) {
charsInUndoBuffer = file->file_state.io_state - __rereading + 1;
position -= charsInUndoBuffer;
}
if (!file->file_mode.binary_io) {
int n = file->buffer_ptr - file->buffer - charsInUndoBuffer;
unsigned char* p = (unsigned char*)file->buffer;
while (n--)
if (*p++ == '\n')
position++;
}
return (position);
}
/* 80365C20-80365E90 360560 0270+00 1/1 0/0 0/0 .text _fseek */
int _fseek(FILE* file, unsigned long offset, int whence) {
int bufferCode;
int pos;
int adjust;
unsigned long state;
int buffLen;
unsigned char* ptr;
if (file->file_mode.file_kind != 1 || file->file_state.error != 0) {
errno = 0x28;
return -1;
}
if (file->file_state.io_state == 1) {
if (__flush_buffer(file, NULL) != 0) {
file->file_state.error = 1;
file->buffer_length = 0;
errno = 0x28;
return -1;
}
}
if (whence == SEEK_CUR) {
whence = SEEK_SET;
adjust = 0;
if ((file->file_mode.file_kind != 1 && file->file_mode.file_kind != 2) ||
file->file_state.error != 0)
{
errno = 0x28;
pos = -1;
} else {
state = file->file_state.io_state;
if (state == 0) {
pos = file->position;
} else {
pos = file->buffer_position;
ptr = file->buffer;
buffLen = (file->buffer_ptr - ptr);
pos += buffLen;
if ((state >= 3)) {
adjust = (state - 2);
pos -= adjust;
}
if (file->file_mode.binary_io == 0) {
int i;
for (i = (buffLen - adjust); i != 0; i--) {
unsigned char c = *ptr;
ptr++;
if (c == 10) {
pos++;
}
}
}
}
}
offset += pos;
}
if ((whence != SEEK_END) && (file->file_mode.io_mode != 3) &&
(file->file_state.io_state == 2 || file->file_state.io_state == 3))
{
if ((offset >= file->position) || !(offset >= file->buffer_position)) {
file->file_state.io_state = 0;
} else {
file->buffer_ptr = file->buffer + (offset - file->buffer_position);
file->buffer_length = file->position - offset;
file->file_state.io_state = 2;
}
} else {
file->file_state.io_state = 0;
}
if (file->file_state.io_state == 0) {
if (file->position_fn != NULL &&
(int)file->position_fn(file->handle, &offset, whence, file->idle_fn))
{
file->file_state.error = 1;
file->buffer_length = 0;
errno = 0x28;
return -1;
} else {
file->file_state.eof = 0;
file->position = offset;
file->buffer_length = 0;
}
}
return 0;
}
/* 80365BB4-80365C20 3604F4 006C+00 0/0 2/2 0/0 .text fseek */
int fseek(FILE* stream, unsigned long offset, int whence) {
int code;
__begin_critical_region(stdin_access);
code = _fseek(stream, offset, whence); // 0 if successful, -1 if error
__end_critical_region(stdin_access);
return code;
}
@@ -0,0 +1,77 @@
#include "abort_exit.h"
#include "critical_regions.h"
#include "stddef.h"
extern void (*_dtors[])(void);
/* 8044D440-8044D540 07A160 0100+00 2/2 0/0 0/0 .bss __atexit_funcs */
static void (*__atexit_funcs[64])(void);
/* 8045199C-804519A0 000E9C 0004+00 2/2 0/0 0/0 .sbss __console_exit */
static void (*__console_exit)(void);
/* 80451998-8045199C 000E98 0004+00 1/1 1/1 0/0 .sbss __stdio_exit */
void (*__stdio_exit)(void);
/* 80451994-80451998 000E94 0004+00 2/2 0/0 0/0 .sbss __atexit_curr_func */
static int __atexit_curr_func;
/* 80451990-80451994 000E90 0004+00 2/2 0/0 0/0 .sbss __aborting */
static int __aborting;
/* 80362ABC-80362B58 35D3FC 009C+00 0/0 9/9 0/0 .text abort */
void abort(void) {
raise(1);
__aborting = 1;
__begin_critical_region(atexit_funcs_access);
while (__atexit_curr_func > 0)
__atexit_funcs[--__atexit_curr_func]();
__end_critical_region(atexit_funcs_access);
__kill_critical_regions();
if (__console_exit != NULL) {
__console_exit();
__console_exit = NULL;
}
_ExitProcess();
}
/* 803629CC-80362ABC 35D30C 00F0+00 0/0 2/2 0/0 .text exit */
void exit(int status) {
int i;
void (**dtor)(void);
if (!__aborting) {
__begin_critical_region(atexit_funcs_access);
__end_critical_region(atexit_funcs_access);
__destroy_global_chain();
dtor = _dtors;
while (*dtor != NULL) {
(*dtor)();
dtor++;
}
if (__stdio_exit != NULL) {
__stdio_exit();
__stdio_exit = NULL;
}
}
__begin_critical_region(atexit_funcs_access);
while (__atexit_curr_func > 0)
__atexit_funcs[--__atexit_curr_func]();
__end_critical_region(atexit_funcs_access);
__kill_critical_regions();
if (__console_exit != NULL) {
__console_exit();
__console_exit = NULL;
}
_ExitProcess();
}
@@ -0,0 +1,334 @@
#include "alloc.h"
#include "critical_regions.h"
typedef struct Block {
struct Block* prev;
struct Block* next;
unsigned long max_size;
unsigned long size;
} Block;
typedef struct SubBlock {
unsigned long size;
Block* block;
struct SubBlock* prev;
struct SubBlock* next;
} SubBlock;
struct FixSubBlock;
typedef struct FixBlock {
struct FixBlock* prev_;
struct FixBlock* next_;
unsigned long client_size_;
struct FixSubBlock* start_;
unsigned long n_allocated_;
} FixBlock;
typedef struct FixSubBlock {
FixBlock* block_;
struct FixSubBlock* next_;
} FixSubBlock;
typedef struct FixStart {
FixBlock* tail_;
FixBlock* head_;
} FixStart;
typedef struct __mem_pool_obj {
Block* start_;
FixStart fix_start[6];
} __mem_pool_obj;
typedef struct __mem_pool {
void* reserved[14];
} __mem_pool;
typedef long tag_word;
typedef struct block_header {
tag_word tag;
struct block_header* prev;
struct block_header* next;
} block_header;
typedef struct list_header {
block_header* rover;
block_header header;
} list_header;
typedef struct heap_header {
struct heap_header* prev;
struct heap_header* next;
} heap_header;
struct mem_pool_obj;
typedef void* (*sys_alloc_ptr)(unsigned long, struct mem_pool_obj*);
typedef void (*sys_free_ptr)(void*, struct mem_pool_obj*);
typedef struct pool_options {
sys_alloc_ptr sys_alloc_func;
sys_free_ptr sys_free_func;
unsigned long min_heap_size;
int always_search_first;
} pool_options;
typedef struct mem_pool_obj {
list_header free_list;
pool_options options;
heap_header* heap_list;
void* userData;
} mem_pool_obj;
mem_pool_obj __malloc_pool;
static int initialized = 0;
static SubBlock* SubBlock_merge_prev(SubBlock*, SubBlock**);
static void SubBlock_merge_next(SubBlock*, SubBlock**);
static const unsigned long fix_pool_sizes[] = {4, 12, 20, 36, 52, 68};
#define SubBlock_size(ths) ((ths)->size & 0xFFFFFFF8)
#define SubBlock_block(ths) ((Block*)((unsigned long)((ths)->block) & ~0x1))
#define Block_size(ths) ((ths)->size & 0xFFFFFFF8)
#define Block_start(ths) (*(SubBlock**)((char*)(ths) + Block_size((ths)) - sizeof(unsigned long)))
#define SubBlock_set_free(ths) \
unsigned long this_size = SubBlock_size((ths)); \
(ths)->size &= ~0x2; \
*(unsigned long*)((char*)(ths) + this_size) &= ~0x4; \
*(unsigned long*)((char*)(ths) + this_size - sizeof(unsigned long)) = this_size
#define SubBlock_is_free(ths) !((ths)->size & 2)
#define SubBlock_set_size(ths, sz) \
(ths)->size &= ~0xFFFFFFF8; \
(ths)->size |= (sz) & 0xFFFFFFF8; \
if (SubBlock_is_free((ths))) \
*(unsigned long*)((char*)(ths) + (sz) - sizeof(unsigned long)) = (sz)
#define SubBlock_from_pointer(ptr) ((SubBlock*)((char*)(ptr)-8))
#define FixSubBlock_from_pointer(ptr) ((FixSubBlock*)((char*)(ptr)-4))
#define FixBlock_client_size(ths) ((ths)->client_size_)
#define FixSubBlock_size(ths) (FixBlock_client_size((ths)->block_))
#define classify(ptr) (*(unsigned long*)((char*)(ptr) - sizeof(unsigned long)) & 1)
#define __msize_inline(ptr) \
(!classify(ptr) ? FixSubBlock_size(FixSubBlock_from_pointer(ptr)) : \
SubBlock_size(SubBlock_from_pointer(ptr)) - 8)
#define Block_empty(ths) \
(_sb = (SubBlock*)((char*)(ths) + 16)), \
SubBlock_is_free(_sb) && SubBlock_size(_sb) == Block_size((ths)) - 24
void __sys_free();
static inline SubBlock* SubBlock_merge_prev(SubBlock* ths, SubBlock** start) {
unsigned long prevsz;
SubBlock* p;
if (!(ths->size & 0x04)) {
prevsz = *(unsigned long*)((char*)ths - sizeof(unsigned long));
if (prevsz & 0x2)
return ths;
p = (SubBlock*)((char*)ths - prevsz);
SubBlock_set_size(p, prevsz + SubBlock_size(ths));
if (*start == ths)
*start = (*start)->next;
ths->next->prev = ths->prev;
ths->next->prev->next = ths->next;
return p;
}
return ths;
}
static inline void SubBlock_merge_next(SubBlock* pBlock, SubBlock** pStart) {
SubBlock* next_sub_block;
unsigned long this_cur_size;
next_sub_block = (SubBlock*)((char*)pBlock + (pBlock->size & 0xFFFFFFF8));
if (!(next_sub_block->size & 2)) {
this_cur_size = (pBlock->size & 0xFFFFFFF8) + (next_sub_block->size & 0xFFFFFFF8);
pBlock->size &= ~0xFFFFFFF8;
pBlock->size |= this_cur_size & 0xFFFFFFF8;
if (!(pBlock->size & 2)) {
*(unsigned long*)((char*)(pBlock) + (this_cur_size)-4) = (this_cur_size);
}
if (!(pBlock->size & 2)) {
*(unsigned long*)((char*)pBlock + this_cur_size) &= ~4;
} else {
*(unsigned long*)((char*)pBlock + this_cur_size) |= 4;
}
if (*pStart == next_sub_block) {
*pStart = (*pStart)->next;
}
if (*pStart == next_sub_block) {
*pStart = 0;
}
next_sub_block->next->prev = next_sub_block->prev;
next_sub_block->prev->next = next_sub_block->next;
}
}
inline void Block_link(Block* ths, SubBlock* sb) {
SubBlock** st;
SubBlock_set_free(sb);
st = &Block_start(ths);
if (*st != 0) {
sb->prev = (*st)->prev;
sb->prev->next = sb;
sb->next = *st;
(*st)->prev = sb;
*st = sb;
*st = SubBlock_merge_prev(*st, st);
SubBlock_merge_next(*st, st);
} else {
*st = sb;
sb->prev = sb;
sb->next = sb;
}
if (ths->max_size < SubBlock_size(*st))
ths->max_size = SubBlock_size(*st);
}
static inline Block* __unlink(__mem_pool_obj* pool_obj, Block* bp) {
Block* result = bp->next;
if (result == bp) {
result = 0;
}
if (pool_obj->start_ == bp) {
pool_obj->start_ = result;
}
if (result != 0) {
result->prev = bp->prev;
result->prev->next = result;
}
bp->next = 0;
bp->prev = 0;
return result;
}
/* 80362D78-8036300C 35D6B8 0294+00 2/2 0/0 0/0 .text deallocate_from_var_pools */
static void deallocate_from_var_pools(__mem_pool_obj* pool_obj, void* ptr) {
SubBlock* sb = SubBlock_from_pointer(ptr);
SubBlock* _sb;
Block* bp = SubBlock_block(sb);
Block_link(bp, sb);
if (Block_empty(bp)) {
__unlink(pool_obj, bp);
__sys_free(bp);
}
}
inline void __init_pool_obj(__mem_pool* pool_obj) {
memset(pool_obj, 0, sizeof(__mem_pool_obj));
}
static inline __mem_pool* get_malloc_pool(void) {
static __mem_pool protopool;
static unsigned char init = 0;
if (!init) {
__init_pool_obj(&protopool);
init = 1;
}
return &protopool;
}
/* 80362C20-80362D78 35D560 0158+00 1/1 0/0 0/0 .text deallocate_from_fixed_pools */
void deallocate_from_fixed_pools(__mem_pool_obj* pool_obj, void* ptr, unsigned long size) {
unsigned long i = 0;
FixSubBlock* p;
FixBlock* b;
FixStart* fs;
while (size > fix_pool_sizes[i]) {
++i;
}
fs = &pool_obj->fix_start[i];
p = FixSubBlock_from_pointer(ptr);
b = p->block_;
if (b->start_ == 0 && fs->head_ != b) {
if (fs->tail_ == b) {
fs->head_ = fs->head_->prev_;
fs->tail_ = fs->tail_->prev_;
} else {
b->prev_->next_ = b->next_;
b->next_->prev_ = b->prev_;
b->next_ = fs->head_;
b->prev_ = b->next_->prev_;
b->prev_->next_ = b;
b->next_->prev_ = b;
fs->head_ = b;
}
}
p->next_ = b->start_;
b->start_ = p;
if (--b->n_allocated_ == 0) {
if (fs->head_ == b) {
fs->head_ = b->next_;
}
if (fs->tail_ == b) {
fs->tail_ = b->prev_;
}
b->prev_->next_ = b->next_;
b->next_->prev_ = b->prev_;
if (fs->head_ == b) {
fs->head_ = 0;
}
if (fs->tail_ == b) {
fs->tail_ = 0;
}
deallocate_from_var_pools(pool_obj, b);
}
}
/* 80362BC8-80362C20 35D508 0058+00 1/1 0/0 0/0 .text __pool_free */
void __pool_free(__mem_pool* pool, void* ptr) {
__mem_pool_obj* pool_obj;
unsigned long size;
if (ptr == 0) {
return;
}
pool_obj = (__mem_pool_obj*)pool;
size = __msize_inline(ptr);
if (size <= 68) {
deallocate_from_fixed_pools(pool_obj, ptr, size);
} else {
deallocate_from_var_pools(pool_obj, ptr);
}
}
/* 80362B58-80362BC8 35D498 0070+00 0/0 2/2 0/0 .text free */
void free(void* ptr) {
__begin_critical_region(malloc_pool_access);
__pool_free(get_malloc_pool(), ptr);
__end_critical_region(malloc_pool_access);
}
@@ -0,0 +1,158 @@
#include "ansi_files.h"
#include "critical_regions.h"
/* 8044D778-8044D878 07A498 0100+00 1/0 0/0 0/0 .bss stdin_buff */
static unsigned char stdin_buff[0x100];
/* 8044D678-8044D778 07A398 0100+00 1/0 0/0 0/0 .bss stdout_buff */
static unsigned char stdout_buff[0x100];
/* 8044D578-8044D678 07A298 0100+00 1/0 0/0 0/0 .bss stderr_buff */
static unsigned char stderr_buff[0x100];
/* 803D29B0-803D2AF0 -00001 0140+00 3/2 15/15 0/0 .data __files */
extern files __files = {
{
0,
0,
1,
1,
2,
0,
0,
0,
0,
0,
0,
0,
0,
0,
{0, 0},
{0, 0},
0,
stdin_buff,
0x00000100,
stdin_buff,
0,
0,
0,
0,
NULL,
__read_console,
__write_console,
__close_console,
NULL,
&__files._stdout,
},
{
1,
0,
2,
1,
2,
0,
0,
0,
0,
0,
0,
0,
0,
0,
{0, 0},
{0, 0},
0,
stdout_buff,
0x00000100,
stdout_buff,
0,
0,
0,
0,
NULL,
__read_console,
__write_console,
__close_console,
NULL,
&__files._stderr,
},
{
2,
0,
2,
0,
2,
0,
0,
0,
0,
0,
0,
0,
0,
0,
{0, 0},
{0, 0},
0,
stderr_buff,
0x00000100,
stderr_buff,
0,
0,
0,
0,
NULL,
__read_console,
__write_console,
__close_console,
NULL,
&__files.empty,
},
{
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, {0, 0},
{0, 0}, 0, NULL, 0x00000000, NULL, 0, 0, 0, 0, NULL, NULL, NULL, NULL, NULL, NULL,
},
};
/* 8036307C-80363124 35D9BC 00A8+00 0/0 1/1 0/0 .text __close_all */
void __close_all(void) {
FILE* file = &__files._stdin;
FILE* last_file;
__begin_critical_region(stdin_access);
while (file != NULL) {
if (file->file_mode.file_kind != __closed_file) {
fclose(file);
}
last_file = file;
file = file->next_file;
if (last_file->is_dynamically_allocated) {
free(last_file);
} else {
last_file->file_mode.file_kind = __unavailable_file;
if (file != NULL && file->is_dynamically_allocated) {
last_file->next_file = NULL;
}
}
}
__end_critical_region(stdin_access);
}
/* 8036300C-8036307C 35D94C 0070+00 0/0 2/2 0/0 .text __flush_all */
unsigned int __flush_all(void) {
unsigned int ret = 0;
FILE* file = &__files._stdin;
while (file) {
if (file->file_mode.file_kind != __closed_file && fflush(file)) {
ret = -1;
}
file = file->next_file;
}
return ret;
}
@@ -0,0 +1,30 @@
#include "arith.h"
/* 803650D0-803650E0 35FA10 0010+00 0/0 66/66 225/225 .text abs */
int abs(int n) {
if (n < 0)
return (-n);
else
return (n);
}
/* 80365078-803650D0 35F9B8 0058+00 0/0 1/1 0/0 .text div */
div_t div(int numerator, int denominator) {
div_t ret;
int i = 1;
int j = 1;
if (numerator < 0 ) {
numerator = -numerator;
i = -1;
}
if (denominator < 0) {
denominator = -denominator;
j = -1;
}
ret.quot = (numerator / denominator) * (i * j);
ret.rem = numerator * i - j * (ret.quot * denominator);
return ret;
}
@@ -0,0 +1,43 @@
#include "buffer_io.h"
void __convert_from_newlines(unsigned char* buf, size_t* n) {}
void __convert_to_newlines(unsigned char* buf, size_t* n) {}
/* 803651A4-803651D8 35FAE4 0034+00 0/0 2/2 0/0 .text __prep_buffer */
void __prep_buffer(FILE* file) {
file->buffer_ptr = file->buffer;
file->buffer_length = file->buffer_size;
file->buffer_length -= file->position & file->buffer_alignment;
file->buffer_position = file->position;
}
/* 803650E0-803651A4 35FA20 00C4+00 0/0 5/5 0/0 .text __flush_buffer */
int __flush_buffer(FILE* file, size_t* bytes_flushed) {
size_t buffer_len;
int ioresult;
buffer_len = file->buffer_ptr - file->buffer;
if (buffer_len) {
file->buffer_length = buffer_len;
if (!file->file_mode.binary_io)
__convert_from_newlines(file->buffer, (size_t*)&file->buffer_length);
ioresult = (*file->write_fn)(file->handle, file->buffer, (size_t*)&file->buffer_length,
file->idle_fn);
if (bytes_flushed)
*bytes_flushed = file->buffer_length;
if (ioresult)
return (ioresult);
file->position += file->buffer_length;
}
__prep_buffer(file);
return 0;
}
@@ -0,0 +1,98 @@
#include "char_io.h"
#include "critical_regions.h"
#include "misc_io.h"
#include "wchar_io.h"
/* 803652AC-80365464 35FBEC 01B8+00 1/1 0/0 0/0 .text __put_char */
int __put_char(int c, FILE* stream) {
int ret;
int file_kind = stream->file_mode.file_kind;
stream->buffer_length = 0;
if (stream->file_state.error != 0 || file_kind == __closed_file) {
return -1;
}
if (file_kind == __console_file) {
__stdio_atexit();
}
if (stream->file_state.io_state == __neutral && (stream->file_mode.io_mode & __write)) {
if ((stream->file_mode.io_mode & __append) && fseek(stream, 0, 2) != 0) {
return 0;
}
stream->file_state.io_state = __writing;
__prep_buffer(stream);
}
if (stream->file_state.io_state != __writing) {
stream->file_state.error = 1;
ret = -1;
stream->buffer_length = 0;
} else if ((stream->file_mode.buffer_mode == 2 ||
stream->buffer_size ==
(unsigned int)stream->buffer_ptr - (unsigned int)stream->buffer) &&
__flush_buffer(stream, NULL) != 0)
{
stream->file_state.error = 1;
ret = -1;
stream->buffer_length = 0;
} else {
stream->buffer_length--;
*stream->buffer_ptr++ = c;
if (stream->file_mode.buffer_mode != 2) {
if ((stream->file_mode.buffer_mode == 0 || c == 10) &&
__flush_buffer(stream, NULL) != 0)
{
stream->file_state.error = 1;
ret = -1;
stream->buffer_length = 0;
goto exit;
}
stream->buffer_length = 0;
}
ret = c & 0xFF;
}
exit:
return ret;
}
/* 803651D8-803652AC 35FB18 00D4+00 0/0 1/1 0/0 .text fputs */
int fputs(const char* s, FILE* stream) {
char c;
int var_r3;
unsigned long len;
int ret = 0;
__begin_critical_region(stdin_access);
while (c = *s++, c != 0) {
if (fwide(stream, -1) >= 0) {
var_r3 = -1;
} else {
len = stream->buffer_length;
stream->buffer_length = len - 1;
if (len != 0) {
char* buf = (char*)stream->buffer_ptr;
stream->buffer_ptr++;
*buf = var_r3 = c & 0xFF;
} else {
var_r3 = __put_char(c, stream);
}
}
if (var_r3 == -1) {
ret = -1;
break;
}
}
__end_critical_region(stdin_access);
return ret;
}
@@ -0,0 +1,82 @@
#include "ctype.h"
#define ctrl __control_char
#define motn __motion_char
#define spac __space_char
#define punc __punctuation
#define digi __digit
#define hexd __hex_digit
#define lowc __lower_case
#define uppc __upper_case
#define dhex (hexd | digi)
#define uhex (hexd | uppc)
#define lhex (hexd | lowc)
/* 803D2C18-803D2D18 02FD38 0100+00 0/0 3/3 0/0 .data __ctype_map */
unsigned char __ctype_map[256] = {
// clang-format off
ctrl, ctrl, ctrl, ctrl, ctrl, ctrl, ctrl, ctrl, ctrl, motn, motn, motn, motn, motn, ctrl, ctrl,
ctrl, ctrl, ctrl, ctrl, ctrl, ctrl, ctrl, ctrl, ctrl, ctrl, ctrl, ctrl, ctrl, ctrl, ctrl, ctrl,
spac, punc, punc, punc, punc, punc, punc, punc, punc, punc, punc, punc, punc, punc, punc, punc,
dhex, dhex, dhex, dhex, dhex, dhex, dhex, dhex, dhex, dhex, punc, punc, punc, punc, punc, punc,
punc, uhex, uhex, uhex, uhex, uhex, uhex, uppc, uppc, uppc, uppc, uppc, uppc, uppc, uppc, uppc,
uppc, uppc, uppc, uppc, uppc, uppc, uppc, uppc, uppc, uppc, uppc, punc, punc, punc, punc, punc,
punc, lhex, lhex, lhex, lhex, lhex, lhex, lowc, lowc, lowc, lowc, lowc, lowc, lowc, lowc, lowc,
lowc, lowc, lowc, lowc, lowc, lowc, lowc, lowc, lowc, lowc, lowc, punc, punc, punc, punc, ctrl,
// clang-format on
};
/* 803D2D18-803D2E18 02FE38 0100+00 1/1 3/3 0/0 .data __lower_map */
unsigned char __lower_map[256] = {
// clang-format off
0x00, 0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08, 0x09, 0x0A, 0x0B, 0x0C, 0x0D, 0x0E, 0x0F,
0x10, 0x11, 0x12, 0x13, 0x14, 0x15, 0x16, 0x17, 0x18, 0x19, 0x1A, 0x1B, 0x1C, 0x1D, 0x1E, 0x1F,
' ', '!', '"', '#', '$', '%', '&', '\'', '(', ')', '*', '+', ',', '-', '.', '/',
'0', '1', '2', '3', '4', '5', '6', '7', '8', '9', ':', ';', '<', '=', '>', '?',
'@', 'a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'i', 'j', 'k', 'l', 'm', 'n', 'o',
'p', 'q', 'r', 's', 't', 'u', 'v', 'w', 'x', 'y', 'z', '[', '\\', ']', '^', '_',
'`', 'a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'i', 'j', 'k', 'l', 'm', 'n', 'o',
'p', 'q', 'r', 's', 't', 'u', 'v', 'w', 'x', 'y', 'z', '{', '|', '}', '~', 0x7F,
0x80, 0x81, 0x82, 0x83, 0x84, 0x85, 0x86, 0x87, 0x88, 0x89, 0x8A, 0x8B, 0x8C, 0x8D, 0x8E, 0x8F,
0x90, 0x91, 0x92, 0x93, 0x94, 0x95, 0x96, 0x97, 0x98, 0x99, 0x9A, 0x9B, 0x9C, 0x9D, 0x9E, 0x9F,
0xA0, 0xA1, 0xA2, 0xA3, 0xA4, 0xA5, 0xA6, 0xA7, 0xA8, 0xA9, 0xAA, 0xAB, 0xAC, 0xAD, 0xAE, 0xAF,
0xB0, 0xB1, 0xB2, 0xB3, 0xB4, 0xB5, 0xB6, 0xB7, 0xB8, 0xB9, 0xBA, 0xBB, 0xBC, 0xBD, 0xBE, 0xBF,
0xC0, 0xC1, 0xC2, 0xC3, 0xC4, 0xC5, 0xC6, 0xC7, 0xC8, 0xC9, 0xCA, 0xCB, 0xCC, 0xCD, 0xCE, 0xCF,
0xD0, 0xD1, 0xD2, 0xD3, 0xD4, 0xD5, 0xD6, 0xD7, 0xD8, 0xD9, 0xDA, 0xDB, 0xDC, 0xDD, 0xDE, 0xDF,
0xE0, 0xE1, 0xE2, 0xE3, 0xE4, 0xE5, 0xE6, 0xE7, 0xE8, 0xE9, 0xEA, 0xEB, 0xEC, 0xED, 0xEE, 0xEF,
0xF0, 0xF1, 0xF2, 0xF3, 0xF4, 0xF5, 0xF6, 0xF7, 0xF8, 0xF9, 0xFA, 0xFB, 0xFC, 0xFD, 0xFE, 0xFF,
// clang-format on
};
/* 80365470-80365494 35FDB0 0024+00 0/0 4/4 0/0 .text tolower */
int tolower(int __c) {
{
if (__c == -1) {
return 0xffffffff;
}
return (unsigned int)__lower_map[__c & 0xff];
}
}
/* ############################################################################################## */
/* 803D2E18-803D2F18 02FF38 0100+00 0/0 1/1 0/0 .data __upper_map */
unsigned char __upper_map[256] = {
0x00, 0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08, 0x09, 0x0A, 0x0B, 0x0C, 0x0D, 0x0E, 0x0F,
0x10, 0x11, 0x12, 0x13, 0x14, 0x15, 0x16, 0x17, 0x18, 0x19, 0x1A, 0x1B, 0x1C, 0x1D, 0x1E, 0x1F,
' ', '!', '"', '#', '$', '%', '&', '\'', '(', ')', '*', '+', ',', '-', '.', '/',
'0', '1', '2', '3', '4', '5', '6', '7', '8', '9', ':', ';', '<', '=', '>', '?',
'@', 'A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J', 'K', 'L', 'M', 'N', 'O',
'P', 'Q', 'R', 'S', 'T', 'U', 'V', 'W', 'X', 'Y', 'Z', '[', '\\', ']', '^', '_',
'`', 'A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J', 'K', 'L', 'M', 'N', 'O',
'P', 'Q', 'R', 'S', 'T', 'U', 'V', 'W', 'X', 'Y', 'Z', '{', '|', '}', '~', 0x7F,
0x80, 0x81, 0x82, 0x83, 0x84, 0x85, 0x86, 0x87, 0x88, 0x89, 0x8A, 0x8B, 0x8C, 0x8D, 0x8E, 0x8F,
0x90, 0x91, 0x92, 0x93, 0x94, 0x95, 0x96, 0x97, 0x98, 0x99, 0x9A, 0x9B, 0x9C, 0x9D, 0x9E, 0x9F,
0xA0, 0xA1, 0xA2, 0xA3, 0xA4, 0xA5, 0xA6, 0xA7, 0xA8, 0xA9, 0xAA, 0xAB, 0xAC, 0xAD, 0xAE, 0xAF,
0xB0, 0xB1, 0xB2, 0xB3, 0xB4, 0xB5, 0xB6, 0xB7, 0xB8, 0xB9, 0xBA, 0xBB, 0xBC, 0xBD, 0xBE, 0xBF,
0xC0, 0xC1, 0xC2, 0xC3, 0xC4, 0xC5, 0xC6, 0xC7, 0xC8, 0xC9, 0xCA, 0xCB, 0xCC, 0xCD, 0xCE, 0xCF,
0xD0, 0xD1, 0xD2, 0xD3, 0xD4, 0xD5, 0xD6, 0xD7, 0xD8, 0xD9, 0xDA, 0xDB, 0xDC, 0xDD, 0xDE, 0xDF,
0xE0, 0xE1, 0xE2, 0xE3, 0xE4, 0xE5, 0xE6, 0xE7, 0xE8, 0xE9, 0xEA, 0xEB, 0xEC, 0xED, 0xEE, 0xEF,
0xF0, 0xF1, 0xF2, 0xF3, 0xF4, 0xF5, 0xF6, 0xF7, 0xF8, 0xF9, 0xFA, 0xFB, 0xFC, 0xFD, 0xFE, 0xFF,
// clang-format on
};
@@ -0,0 +1,121 @@
#include "direct_io.h"
#include "critical_regions.h"
#include "wchar_io.h"
/* 80365494-803657A0 35FDD4 030C+00 1/1 0/0 0/0 .text __fwrite */
size_t fwrite(const void* buffer, size_t size, size_t count, FILE* stream) {
size_t retval;
__begin_critical_region(stdin_access);
retval = __fwrite(buffer, size, count, stream);
__end_critical_region(stdin_access);
return (retval);
}
/* 803657A0-8036581C 3600E0 007C+00 0/0 1/1 0/0 .text fwrite */
size_t __fwrite(const void* buffer, size_t size, size_t count, FILE* stream) {
unsigned char* write_ptr;
size_t num_bytes, bytes_to_go, bytes_written;
int ioresult, always_buffer;
#ifndef __NO_WIDE_CHAR
if (fwide(stream, 0) == 0)
fwide(stream, -1);
#endif
bytes_to_go = size * count;
if (!bytes_to_go || stream->file_state.error || stream->file_mode.file_kind == __closed_file)
return 0;
if (stream->file_mode.file_kind == __console_file)
__stdio_atexit();
always_buffer = !stream->file_mode.binary_io || stream->file_mode.buffer_mode == _IOFBF ||
stream->file_mode.buffer_mode == _IOLBF;
if (stream->file_state.io_state == __neutral) {
if (stream->file_mode.io_mode & __write) {
if (stream->file_mode.io_mode & __append) {
if (fseek(stream, 0, SEEK_END))
return 0;
}
stream->file_state.io_state = __writing;
__prep_buffer(stream);
}
}
if (stream->file_state.io_state != __writing) {
set_error(stream);
return 0;
}
write_ptr = (unsigned char*)buffer;
bytes_written = 0;
if (bytes_to_go && (stream->buffer_ptr != stream->buffer || always_buffer)) {
stream->buffer_length = stream->buffer_size - (stream->buffer_ptr - stream->buffer);
do {
unsigned char* newline = NULL;
num_bytes = stream->buffer_length;
if (num_bytes > bytes_to_go)
num_bytes = bytes_to_go;
if (stream->file_mode.buffer_mode == _IOLBF && num_bytes)
if ((newline = (unsigned char*)__memrchr(write_ptr, '\n', num_bytes)) != NULL)
num_bytes = newline + 1 - write_ptr;
if (num_bytes) {
memcpy(stream->buffer_ptr, write_ptr, num_bytes);
write_ptr += num_bytes;
bytes_written += num_bytes;
bytes_to_go -= num_bytes;
stream->buffer_ptr += num_bytes;
stream->buffer_length -= num_bytes;
}
if (!stream->buffer_length || newline != NULL ||
(stream->file_mode.buffer_mode == _IONBF))
{
ioresult = __flush_buffer(stream, NULL);
if (ioresult) {
set_error(stream);
bytes_to_go = 0;
break;
}
}
} while (bytes_to_go && always_buffer);
}
if (bytes_to_go && !always_buffer) {
unsigned char* save_buffer = stream->buffer;
size_t save_size = stream->buffer_size;
stream->buffer = write_ptr;
stream->buffer_size = bytes_to_go;
stream->buffer_ptr = write_ptr + bytes_to_go;
if (__flush_buffer(stream, &num_bytes) != __no_io_error)
set_error(stream);
bytes_written += num_bytes;
stream->buffer = save_buffer;
stream->buffer_size = save_size;
__prep_buffer(stream);
stream->buffer_length = 0;
}
if (stream->file_mode.buffer_mode != _IOFBF)
stream->buffer_length = 0;
return ((bytes_written + size - 1) / size);
}
@@ -0,0 +1,4 @@
#include "errno.h"
/* 804519A8-804519B0 000EA8 0004+04 0/0 6/6 0/0 .sbss errno */
int errno;
@@ -0,0 +1,27 @@
#include "extras.h"
#include "ctype.h"
/* 8036CA94-8036CB20 3673D4 008C+00 0/0 1/1 0/0 .text stricmp */
int stricmp(const char* str1, const char* str2) {
char a_var;
char b_var;
do {
b_var = _tolower(*str1++);
a_var = _tolower(*str2++);
if (b_var < a_var) {
return -1;
}
if (b_var > a_var) {
return 1;
}
} while (b_var != 0);
return 0;
}
/* 8036CA74-8036CA94 3673B4 0020+00 0/0 1/1 0/0 .text strnicmp */
int strnicmp(const char* str1, const char* str2, int n) {
return __msl_strnicmp(str1, str2, n);
}
@@ -0,0 +1,95 @@
#include "file_io.h"
#include "ctype.h"
/* 803659F8-80365BB4 360338 01BC+00 0/0 1/1 0/0 .text fclose */
int fclose(FILE* file) {
int flush_result, close_result;
if (file == NULL)
return (-1);
if (file->file_mode.file_kind == __closed_file)
return (0);
flush_result = fflush(file);
close_result = (*file->close_fn)(file->handle);
file->file_mode.file_kind = __closed_file;
file->handle = NULL;
if (file->file_state.free_buffer)
free((FILE*)file->buffer);
return ((flush_result || close_result) ? -1 : 0);
}
/* 803658C0-803659F8 360200 0138+00 0/0 4/4 0/0 .text fflush */
int fflush(FILE* file) {
int pos;
if (file == NULL) {
return __flush_all();
}
if (file->file_state.error != 0 || file->file_mode.file_kind == __closed_file) {
return -1;
}
if (file->file_mode.io_mode == __read) {
return 0;
}
if (file->file_state.io_state >= __rereading) {
file->file_state.io_state = __reading;
}
if (file->file_state.io_state == __reading) {
file->buffer_length = 0;
}
if (file->file_state.io_state != __writing) {
file->file_state.io_state = __neutral;
return 0;
}
if (file->file_mode.file_kind != __disk_file) {
pos = 0;
} else {
pos = ftell(file);
}
if (__flush_buffer(file, 0) != 0) {
file->file_state.error = 1;
file->buffer_length = 0;
return -1;
}
file->file_state.io_state = __neutral;
file->position = pos;
file->buffer_length = 0;
return 0;
}
/* 8036581C-803658C0 36015C 00A4+00 0/0 1/1 0/0 .text __msl_strnicmp */
int __msl_strnicmp(const char* str1, const char* str2, int n) {
int i;
char c1, c2;
for (i = 0; i < n; i++) {
c1 = _tolower(*str1++);
c2 = _tolower(*str2++);
if (c1 < c2) {
return -1;
}
if (c1 > c2) {
return 1;
}
if (c1 == '\0') {
return 0;
}
}
return 0;
}
@@ -0,0 +1,13 @@
#include "float.h"
/* 80450AE0-80450AE4 000560 0004+00 0/0 204/204 1060/1060 .sdata __float_nan */
unsigned long __float_nan[] = {0x7FFFFFFF};
/* 80450AE4-80450AE8 000564 0004+00 0/0 1/1 0/0 .sdata __float_huge */
unsigned long __float_huge[] = {0x7F800000};
/* 80450AE8-80450AEC 000568 0004+00 0/0 18/18 14/14 .sdata __float_max */
unsigned long __float_max[] = {0x7F7FFFFF};
/* 80450AEC-80450AF0 00056C 0004+00 0/0 28/28 0/0 .sdata __float_epsilon */
unsigned long __float_epsilon[] = {0x34000000};
@@ -0,0 +1,65 @@
#include "mbstring.h"
#include "string.h"
inline static int unicode_to_UTF8(char* s, wchar_t wchar) {
int number_of_bytes;
wchar_t wide_char;
char* target_ptr;
char first_byte_mark[4] = { 0x00, 0x00, 0xc0, 0xe0 };
if (!s)
return (0);
wide_char = wchar;
if (wide_char < 0x0080)
number_of_bytes = 1;
else if (wide_char < 0x0800)
number_of_bytes = 2;
else
number_of_bytes = 3;
target_ptr = s + number_of_bytes;
switch (number_of_bytes) {
case 3:
*--target_ptr = (wide_char & 0x003f) | 0x80;
wide_char >>= 6;
case 2:
*--target_ptr = (wide_char & 0x003f) | 0x80;
wide_char >>= 6;
case 1:
*--target_ptr = wide_char | first_byte_mark[number_of_bytes];
}
return number_of_bytes;
}
inline int wctomb(char* s, wchar_t wchar) { return (unicode_to_UTF8(s, wchar)); }
/* 80365F74-8036608C 3608B4 0118+00 0/0 1/1 0/0 .text wcstombs */
size_t wcstombs(char* s, const wchar_t* pwcs, size_t n) {
int chars_written = 0;
int result;
char temp[3];
wchar_t* source;
if (!s || !pwcs)
return (0);
source = (wchar_t*)pwcs;
while (chars_written <= n) {
if (!*source) {
*(s + chars_written) = '\0';
break;
} else {
result = wctomb(temp, *source++);
if ((chars_written + result) <= n) {
strncpy(s + chars_written, temp, result);
chars_written += result;
} else
break;
}
}
return chars_written;
}
@@ -0,0 +1,85 @@
#include "string.h"
/* 80366130-803661FC 360A70 00CC+00 0/0 2/2 0/0 .text memmove */
void* memmove(void* dst, const void* src, size_t n) {
unsigned char* csrc;
unsigned char* cdst;
int reverse = (unsigned int)src < (unsigned int)dst;
if (n >= 32) {
if (((unsigned int)dst ^ (unsigned int)src) & 3) {
if (!reverse) {
__copy_longs_unaligned(dst, src, n);
} else {
__copy_longs_rev_unaligned(dst, src, n);
}
} else {
if (!reverse) {
__copy_longs_aligned(dst, src, n);
} else {
__copy_longs_rev_aligned(dst, src, n);
}
}
return dst;
} else {
if (!reverse) {
csrc = ((unsigned char*)src) - 1;
cdst = ((unsigned char*)dst) - 1;
n++;
while (--n > 0) {
*++cdst = *++csrc;
}
} else {
csrc = (unsigned char*)src + n;
cdst = (unsigned char*)dst + n;
n++;
while (--n > 0) {
*--cdst = *--csrc;
}
}
}
return dst;
}
/* 80366104-80366130 360A44 002C+00 0/0 1/1 0/0 .text memchr */
void* memchr(const void* ptr, int ch, size_t count) {
const unsigned char* p;
unsigned long v = (ch & 0xff);
for (p = (unsigned char*)ptr - 1, count++; --count;)
if ((*++p & 0xff) == v)
return (void*)p;
return NULL;
}
/* 803660D8-80366104 360A18 002C+00 0/0 1/1 0/0 .text __memrchr */
void* __memrchr(const void* ptr, int ch, size_t count) {
const unsigned char* p;
unsigned long v = (ch & 0xff);
for (p = (unsigned char*)ptr + count, count++; --count;)
if (*--p == v)
return (void*)p;
return NULL;
}
/* 8036608C-803660D8 3609CC 004C+00 0/0 19/19 5/5 .text memcmp */
int memcmp(const void* lhs, const void* rhs, size_t count) {
const unsigned char* p1;
const unsigned char* p2;
for (p1 = (const unsigned char*)lhs - 1, p2 = (const unsigned char*)rhs - 1, count++; --count;)
if (*++p1 != *++p2)
return ((*p1 < *p2) ? -1 : +1);
return 0;
}
@@ -0,0 +1,221 @@
#include "mem_funcs.h"
#define cps ((unsigned char*)src)
#define cpd ((unsigned char*)dst)
#define lps ((unsigned long*)src)
#define lpd ((unsigned long*)dst)
#define deref_auto_inc(p) *++(p)
/* 80366410-803664CC 360D50 00BC+00 0/0 1/1 0/0 .text __copy_longs_aligned */
void __copy_longs_aligned(void* dst, const void* src, size_t n) {
unsigned long i;
i = (-(unsigned long)dst) & 3;
cps = ((unsigned char*)src) - 1;
cpd = ((unsigned char*)dst) - 1;
if (i) {
n -= i;
do
deref_auto_inc(cpd) = deref_auto_inc(cps);
while (--i);
}
lps = ((unsigned long*)(cps + 1)) - 1;
lpd = ((unsigned long*)(cpd + 1)) - 1;
i = n >> 5;
if (i)
do {
deref_auto_inc(lpd) = deref_auto_inc(lps);
deref_auto_inc(lpd) = deref_auto_inc(lps);
deref_auto_inc(lpd) = deref_auto_inc(lps);
deref_auto_inc(lpd) = deref_auto_inc(lps);
deref_auto_inc(lpd) = deref_auto_inc(lps);
deref_auto_inc(lpd) = deref_auto_inc(lps);
deref_auto_inc(lpd) = deref_auto_inc(lps);
deref_auto_inc(lpd) = deref_auto_inc(lps);
} while (--i);
i = (n & 31) >> 2;
if (i)
do
deref_auto_inc(lpd) = deref_auto_inc(lps);
while (--i);
cps = ((unsigned char*)(lps + 1)) - 1;
cpd = ((unsigned char*)(lpd + 1)) - 1;
n &= 3;
if (n)
do
deref_auto_inc(cpd) = deref_auto_inc(cps);
while (--n);
return;
}
/* 80366368-80366410 360CA8 00A8+00 0/0 1/1 0/0 .text __copy_longs_rev_aligned */
void __copy_longs_rev_aligned(void* dst, const void* src, size_t n) {
unsigned long i;
cps = ((unsigned char*)src) + n;
cpd = ((unsigned char*)dst) + n;
i = ((unsigned long)cpd) & 3;
if (i) {
n -= i;
do
*--cpd = *--cps;
while (--i);
}
i = n >> 5;
if (i)
do {
*--lpd = *--lps;
*--lpd = *--lps;
*--lpd = *--lps;
*--lpd = *--lps;
*--lpd = *--lps;
*--lpd = *--lps;
*--lpd = *--lps;
*--lpd = *--lps;
} while (--i);
i = (n & 31) >> 2;
if (i)
do
*--lpd = *--lps;
while (--i);
n &= 3;
if (n)
do
*--cpd = *--cps;
while (--n);
return;
}
/* 803662A8-80366368 360BE8 00C0+00 0/0 1/1 0/0 .text __copy_longs_unaligned */
void __copy_longs_unaligned(void* dst, const void* src, size_t n) {
unsigned long i, v1, v2;
unsigned int src_offset, left_shift, right_shift;
i = (-(unsigned long)dst) & 3;
cps = ((unsigned char*)src) - 1;
cpd = ((unsigned char*)dst) - 1;
if (i) {
n -= i;
do
deref_auto_inc(cpd) = deref_auto_inc(cps);
while (--i);
}
src_offset = ((unsigned int)(cps + 1)) & 3;
left_shift = src_offset << 3;
right_shift = 32 - left_shift;
cps -= src_offset;
lps = ((unsigned long*)(cps + 1)) - 1;
lpd = ((unsigned long*)(cpd + 1)) - 1;
i = n >> 3;
v1 = deref_auto_inc(lps);
do {
v2 = deref_auto_inc(lps);
deref_auto_inc(lpd) = (v1 << left_shift) | (v2 >> right_shift);
v1 = deref_auto_inc(lps);
deref_auto_inc(lpd) = (v2 << left_shift) | (v1 >> right_shift);
} while (--i);
if (n & 4) {
v2 = deref_auto_inc(lps);
deref_auto_inc(lpd) = (v1 << left_shift) | (v2 >> right_shift);
}
cps = ((unsigned char*)(lps + 1)) - 1;
cpd = ((unsigned char*)(lpd + 1)) - 1;
n &= 3;
if (n) {
cps -= 4 - src_offset;
do
deref_auto_inc(cpd) = deref_auto_inc(cps);
while (--n);
}
return;
}
/* 803661FC-803662A8 360B3C 00AC+00 0/0 1/1 0/0 .text __copy_longs_rev_unaligned */
void __copy_longs_rev_unaligned(void* dst, const void* src, size_t n) {
unsigned long i, v1, v2;
unsigned int src_offset, left_shift, right_shift;
cps = ((unsigned char*)src) + n;
cpd = ((unsigned char*)dst) + n;
i = ((unsigned long)cpd) & 3;
if (i) {
n -= i;
do
*--cpd = *--cps;
while (--i);
}
src_offset = ((unsigned int)cps) & 3;
left_shift = src_offset << 3;
right_shift = 32 - left_shift;
cps += 4 - src_offset;
i = n >> 3;
v1 = *--lps;
do {
v2 = *--lps;
*--lpd = (v2 << left_shift) | (v1 >> right_shift);
v1 = *--lps;
*--lpd = (v1 << left_shift) | (v2 >> right_shift);
} while (--i);
if (n & 4) {
v2 = *--lps;
*--lpd = (v2 << left_shift) | (v1 >> right_shift);
}
n &= 3;
if (n) {
cps += src_offset;
do
*--cpd = *--cps;
while (--n);
}
return;
}
@@ -0,0 +1,8 @@
#include "misc_io.h"
#include "abort_exit.h"
#include "ansi_files.h"
/* 803664CC-803664DC 360E0C 0010+00 0/0 2/2 0/0 .text __stdio_atexit */
void __stdio_atexit(void) {
__stdio_exit = __close_all;
}
File diff suppressed because it is too large Load Diff
@@ -0,0 +1,34 @@
#include "scanf.h"
/* 8036878C-8036881C 3630CC 0090+00 0/0 2/2 0/0 .text __StringRead */
int __StringRead(void* pPtr, int ch, int act) {
char ret;
__InStrCtrl* Iscp = (__InStrCtrl*)pPtr;
switch (act) {
case __GetAChar:
ret = *(Iscp->NextChar);
if (ret == '\0') {
Iscp->NullCharDetected = 1;
return -1;
} else {
Iscp->NextChar++;
return (unsigned char)ret;
}
case __UngetAChar:
if (Iscp->NullCharDetected == 0) {
Iscp->NextChar--;
} else {
Iscp->NullCharDetected = 0;
}
return ch;
case __TestForError:
return Iscp->NullCharDetected;
}
return 0;
}
@@ -0,0 +1,38 @@
#include "signal.h"
#include "critical_regions.h"
#define __std(ref) ref
/* 8044D878-8044D890 07A598 0018+00 1/1 0/0 0/0 .bss signal_funcs */
__signal_func_ptr signal_funcs[6];
/* 8036881C-803688DC 36315C 00C0+00 0/0 1/1 0/0 .text raise */
int raise(int sig) {
__signal_func_ptr signal_func;
if (sig < 1 || sig > 6) {
return -1;
}
__begin_critical_region(stderr_access);
signal_func = signal_funcs[sig - 1];
if (signal_func != ((__std(__signal_func_ptr))1)) {
signal_funcs[sig - 1] = ((__std(__signal_func_ptr))0);
}
__end_critical_region(stderr_access);
if (signal_func == ((__std(__signal_func_ptr))1) ||
(signal_func == ((__std(__signal_func_ptr))0) && sig == 1))
{
return 0;
}
if (signal_func == ((__std(__signal_func_ptr))0)) {
exit(0);
}
(*signal_func)(sig);
return 0;
}
@@ -0,0 +1,246 @@
#include "string.h"
#define K1 0x80808080
#define K2 0xFEFEFEFF
/* 80368BE4-80368C00 363524 001C+00 0/0 28/28 8/8 .text strlen */
size_t strlen(const char* str) {
size_t len = -1;
unsigned char* p = (unsigned char*)str - 1;
do {
len++;
} while (*++p);
return len;
}
/* 80368B2C-80368BE4 36346C 00B8+00 0/0 131/131 13/13 .text strcpy */
char* strcpy(char* dst, const char* src) {
register unsigned char *destb, *fromb;
register unsigned long w, t, align;
fromb = (unsigned char*)src;
destb = (unsigned char*)dst;
if ((align = ((int)fromb & 3)) != ((int)destb & 3)) {
goto bytecopy;
}
if (align) {
if ((*destb = *fromb) == 0) {
return dst;
}
for (align = 3 - align; align; align--) {
if ((*(++destb) = *(++fromb)) == 0) {
return dst;
}
}
++destb;
++fromb;
}
w = *((int*)(fromb));
t = w + K2;
t &= K1;
if (t) {
goto bytecopy;
}
--((int*)(destb));
do {
*(++((int*)(destb))) = w;
w = *(++((int*)(fromb)));
t = w + K2;
t &= K1;
if (t) {
goto adjust;
}
} while (1);
adjust:
++((int*)(destb));
bytecopy:
if ((*destb = *fromb) == 0) {
return dst;
}
do {
if ((*(++destb) = *(++fromb)) == 0) {
return dst;
}
} while (1);
return dst;
}
/* 80368AE8-80368B2C 363428 0044+00 0/0 9/9 1/1 .text strncpy */
char* strncpy(char* dst, const char* src, size_t n) {
const unsigned char* p = (const unsigned char*)src - 1;
unsigned char* q = (unsigned char*)dst - 1;
n++;
while (--n) {
if (!(*++q = *++p)) {
while (--n) {
*++q = 0;
}
break;
}
}
return dst;
}
/* 80368ABC-80368AE8 3633FC 002C+00 0/0 20/20 8/8 .text strcat */
char* strcat(char* dst, const char* src) {
const unsigned char* p = (unsigned char*)src - 1;
unsigned char* q = (unsigned char*)dst - 1;
while (*++q) {
}
q--;
while (*++q = *++p) {
}
return dst;
}
/* 80368994-80368ABC 3632D4 0128+00 0/0 155/155 279/279 .text strcmp */
int strcmp(const char* str1, const char* str2) {
register unsigned char* left = (unsigned char*)str1;
register unsigned char* right = (unsigned char*)str2;
unsigned long align, l1, r1, x;
l1 = *left;
r1 = *right;
if (l1 - r1) {
return l1 - r1;
}
if ((align = ((int)left & 3)) != ((int)right & 3)) {
goto bytecopy;
}
if (align) {
if (l1 == 0) {
return 0;
}
for (align = 3 - align; align; align--) {
l1 = *(++left);
r1 = *(++right);
if (l1 - r1) {
return l1 - r1;
}
if (l1 == 0) {
return 0;
}
}
left++;
right++;
}
l1 = *(int*)left;
r1 = *(int*)right;
x = l1 + K2;
if (x & K1) {
goto adjust;
}
while (l1 == r1) {
l1 = *(++((int*)(left)));
r1 = *(++((int*)(right)));
x = l1 + K2;
if (x & K1) {
goto adjust;
}
}
if (l1 > r1) {
return 1;
}
return -1;
adjust:
l1 = *left;
r1 = *right;
if (l1 - r1) {
return l1 - r1;
}
bytecopy:
if (l1 == 0) {
return 0;
}
do {
l1 = *(++left);
r1 = *(++right);
if (l1 - r1) {
return l1 - r1;
}
if (l1 == 0) {
return 0;
}
} while (1);
}
/* 80368954-80368994 363294 0040+00 0/0 6/6 0/0 .text strncmp */
int strncmp(const char* str1, const char* str2, size_t n) {
const unsigned char* p1 = (unsigned char*)str1 - 1;
const unsigned char* p2 = (unsigned char*)str2 - 1;
unsigned long c1, c2;
n++;
while (--n) {
if ((c1 = *++p1) != (c2 = *++p2)) {
return c1 - c2;
} else if (c1 == 0) {
break;
}
}
return 0;
}
/* 80368924-80368954 363264 0030+00 0/0 3/3 0/0 .text strchr */
char* strchr(const char* str, int c) {
const unsigned char* p = (unsigned char*)str - 1;
unsigned long chr = (c & 0xFF);
unsigned long ch;
while (ch = *++p) {
if (ch == chr) {
return (char*)p;
}
}
return chr ? NULL : (char*)p;
}
/* 803688DC-80368924 36321C 0048+00 0/0 1/1 0/0 .text strrchr */
char* strrchr(const char* str, int c) {
const unsigned char* p = (unsigned char*)str - 1;
const unsigned char* q = NULL;
unsigned long chr = (c & 0xFF);
unsigned long ch;
while (ch = *++p) {
if (ch == chr) {
q = p;
}
}
if (q != NULL) {
return (char*)q;
}
return chr ? NULL : (char*)p;
}
@@ -0,0 +1,203 @@
#include "strtoul.h"
#include "ctype.h"
#include "errno.h"
#include "limits.h"
#include "scanf.h"
enum scan_states {
start = 0x01,
check_for_zero = 0x02,
leading_zero = 0x04,
need_digit = 0x08,
digit_loop = 0x10,
finished = 0x20,
failure = 0x40
};
#define final_state(scan_state) (scan_state & (finished | failure))
#define success(scan_state) (scan_state & (leading_zero | digit_loop | finished))
#define fetch() (count++, (*ReadProc)(ReadProcArg, 0, __GetAChar))
#define unfetch(c) (*ReadProc)(ReadProcArg, c, __UngetAChar)
/* 80368D9C-80369114 3636DC 0378+00 3/2 0/0 0/0 .text __strtoul */
unsigned long __strtoul(int base, int max_width, int (*ReadProc)(void*, int, int),
void* ReadProcArg, int* chars_scanned, int* negative, int* overflow) {
int scan_state = start;
int count = 0;
int spaces = 0;
unsigned long value = 0;
unsigned long value_max = 0;
int c;
*negative = *overflow = 0;
if (base < 0 || base == 1 || base > 36 || max_width < 1) {
scan_state = failure;
} else {
c = fetch();
}
if (base != 0) {
value_max = ULONG_MAX / base;
}
while (count <= max_width && c != -1 && !final_state(scan_state)) {
switch (scan_state) {
case start:
if (isspace(c)) {
c = fetch();
count--;
spaces++;
break;
}
if (c == '+') {
c = fetch();
} else if (c == '-') {
c = fetch();
*negative = 1;
}
scan_state = check_for_zero;
break;
case check_for_zero:
if (base == 0 || base == 16) {
if (c == '0') {
scan_state = leading_zero;
c = fetch();
break;
}
}
scan_state = need_digit;
break;
case 4:
if (c == 'X' || c == 'x') {
base = 16;
scan_state = need_digit;
c = fetch();
break;
}
if (base == 0) {
base = 8;
}
scan_state = digit_loop;
break;
case need_digit:
case digit_loop:
if (base == 0) {
base = 10;
}
if (!value_max) {
value_max = ULONG_MAX / base;
}
if (isdigit(c)) {
if ((c -= '0') >= base) {
if (scan_state == digit_loop) {
scan_state = finished;
} else {
scan_state = failure;
}
c += '0';
break;
}
} else if (!isalpha(c) || (toupper(c) - 'A' + 10) >= base) {
if (scan_state == digit_loop) {
scan_state = finished;
} else {
scan_state = failure;
}
break;
} else {
c = toupper(c) - 'A' + 10;
}
if (value > value_max) {
*overflow = 1;
}
value *= base;
if (c > (ULONG_MAX - value)) {
*overflow = 1;
}
value += c;
scan_state = digit_loop;
c = fetch();
break;
}
}
if (!success(scan_state)) {
count = 0;
value = 0;
*chars_scanned = 0;
} else {
count--;
*chars_scanned = count + spaces;
}
unfetch(c);
return value;
}
/* 80368CF0-80368D9C 363630 00AC+00 0/0 2/2 0/0 .text strtoul */
unsigned long strtoul(const char* str, char** end, int base) {
unsigned long value;
int count, negative, overflow;
__InStrCtrl isc;
isc.NextChar = (char*)str;
isc.NullCharDetected = 0;
value = __strtoul(base, 0x7FFFFFFF, &__StringRead, (void*)&isc, &count, &negative, &overflow);
if (end) {
*end = (char*)str + count;
}
if (overflow) {
value = ULONG_MAX;
errno = 0x22;
} else if (negative) {
value = -value;
}
return value;
}
/* 80368C00-80368CF0 363540 00F0+00 0/0 3/3 0/0 .text strtol */
long strtol(const char* str, char** end, int base) {
unsigned long uvalue;
long svalue;
int count, negative, overflow;
__InStrCtrl isc;
isc.NextChar = (char*)str;
isc.NullCharDetected = 0;
uvalue = __strtoul(base, 0x7FFFFFFF, &__StringRead, (void*)&isc, &count, &negative, &overflow);
if (end) {
*end = (char*)str + count;
}
if (overflow || (!negative && uvalue > LONG_MAX) || (negative && uvalue > -LONG_MIN)) {
svalue = (negative ? -LONG_MIN : LONG_MAX);
errno = ERANGE;
} else {
svalue = (negative ? (long)-uvalue : (long)uvalue);
}
return svalue;
}
@@ -0,0 +1,22 @@
#include "wchar_io.h"
/* 80369114-8036919C 363A54 0088+00 0/0 5/5 0/0 .text fwide */
int fwide(FILE* file, int mode) {
if (file == NULL || file->file_mode.file_kind == __closed_file) {
return 0;
}
switch (file->file_mode.file_orientation) {
case UNORIENTED:
if (mode > 0) {
file->file_mode.file_orientation = WIDE_ORIENTED;
} else if (mode < 0) {
file->file_mode.file_orientation = CHAR_ORIENTED;
}
return mode;
case WIDE_ORIENTED:
return 1;
case CHAR_ORIENTED:
return -1;
}
}
@@ -0,0 +1,104 @@
/* @(#)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 "math.h"
#ifdef __STDC__
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 */
#ifdef __STDC__
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 */
}
return NAN; /* 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);
}
}
@@ -0,0 +1,115 @@
/* @(#)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 "math.h"
#ifdef __STDC__
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 */
#ifdef __STDC__
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;
return NAN; /* 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;
}
@@ -0,0 +1,142 @@
/* @(#)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"
#ifdef __STDC__
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 */
#ifdef __STDC__
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(-,-) */
}
}
@@ -0,0 +1,161 @@
/* @(#)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"
#ifdef __STDC__
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 */
#ifdef __STDC__
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;
}
}
@@ -0,0 +1,166 @@
/* @(#)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"
#ifdef __STDC__
static const double one = 1.0, Zero[] = {
0.0,
-0.0,
};
#else
static double one = 1.0, Zero[] = {
0.0,
-0.0,
};
#endif
#ifdef __STDC__
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;
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;
while (n--) {
hz = hx - hy;
lz = lx - ly;
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 (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 */
}
@@ -0,0 +1,407 @@
//========================================================================
//
// e_pow.c
//
// Part of the standard mathematical function library
//
//========================================================================
//####ECOSGPLCOPYRIGHTBEGIN####
// -------------------------------------------
// This file is part of eCos, the Embedded Configurable Operating System.
// Copyright (C) 1998, 1999, 2000, 2001, 2002 Red Hat, Inc.
//
// eCos is free software; you can redistribute it and/or modify it under
// the terms of the GNU General Public License as published by the Free
// Software Foundation; either version 2 or (at your option) any later version.
//
// eCos is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or
// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
// for more details.
//
// You should have received a copy of the GNU General Public License along
// with eCos; if not, write to the Free Software Foundation, Inc.,
// 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA.
//
// As a special exception, if other files instantiate templates or use macros
// or inline functions from this file, or you compile this file and link it
// with other works to produce a work based on this file, this file does not
// by itself cause the resulting work to be covered by the GNU General Public
// License. However the source code for this file must still be made available
// in accordance with section (3) of the GNU General Public License.
//
// This exception does not invalidate any other reasons why a work based on
// this file might be covered by the GNU General Public License.
//
// Alternative licenses for eCos may be arranged by contacting Red Hat, Inc.
// at http://sources.redhat.com/ecos/ecos-license/
// -------------------------------------------
//####ECOSGPLCOPYRIGHTEND####
//========================================================================
//#####DESCRIPTIONBEGIN####
//
// Author(s): jlarmour
// Contributors:
// Date: 2001-07-20
// Purpose:
// Description:
// Usage:
//
//####DESCRIPTIONEND####
//
//========================================================================
// CONFIGURATION
/* @(#)e_pow.c 5.1 93/09/24 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/* __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 "types.h"
#include "errno.h"
#include "math.h"
#ifndef _DOUBLE_IS_32BITS
#ifdef __STDC__
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*/
#ifdef __STDC__
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;
double qqq; // necessary temp
int i0, i1, i, j, k, yisint, n;
int hx, hy, ix, iy;
u32 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);
qqq = ax; /*x is +-0,+-inf,+-1*/
/* special value of x */
if (lx == 0) {
if (ix == 0x7ff00000 || ix == 0 || ix == 0x3ff00000) {
z = qqq; /*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;
}
}
/* (x<0)**(non-int) is NaN */
/* CYGNUS LOCAL: This used to be
if((((hx>>31)+1)|yisint)==0) return (x-x)/(x-x);
but ANSI C says a right shift of a signed negative quantity is
implementation defined. */
if (((((int)hx >> 31) + 1) | yisint) == 0) {
errno = 33;
return (double)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|<sqrt(3/2) */
else if (j < 0xBB67A)
k = 1; /* |x|<sqrt(3) */
else {
k = 0;
n += 1;
ix -= 0x00100000;
}
__HI(ax) = ix;
/* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
u = ax - bp[k]; /* bp[0]=1.0, bp[1]=1.5 */
v = one / (ax + bp[k]);
s = u * v;
s_h = s;
__LO(s_h) = 0;
/* t_h=ax+bp[k] High */
t_h = zero;
__HI(t_h) = ((ix >> 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 = ldexp(z, n); /* subnormal output */
else
__HI(z) += (n << 20);
return s * z;
}
#endif /* defined(_DOUBLE_IS_32BITS) */
// EOF e_pow.c
@@ -0,0 +1,179 @@
/* @(#)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 "math.h"
/*
* Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
*/
#ifdef __STDC__
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,
};
#ifdef __STDC__
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)
*/
#ifdef __STDC__
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 */
#ifdef __STDC__
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;
r = t - fn * pio2_2;
w = fn * pio2_2t - ((t - r) - fn * pio2_2);
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;
}
@@ -0,0 +1,461 @@
/* @(#)e_sqrt.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_sqrt(x)
* Return correctly rounded sqrt.
* ------------------------------------------
* | Use the hardware sqrt if you have one |
* ------------------------------------------
* Method:
* Bit by bit method using integer arithmetic. (Slow, but portable)
* 1. Normalization
* Scale x to y in [1,4) with even powers of 2:
* find an integer k such that 1 <= (y=x*2^(2k)) < 4, then
* sqrt(x) = 2^k * sqrt(y)
* 2. Bit by bit computation
* Let q = sqrt(y) truncated to i bit after binary point (q = 1),
* i 0
* i+1 2
* s = 2*q , and y = 2 * ( y - q ). (1)
* i i i i
*
* To compute q from q , one checks whether
* i+1 i
*
* -(i+1) 2
* (q + 2 ) <= y. (2)
* i
* -(i+1)
* If (2) is false, then q = q ; otherwise q = q + 2 .
* i+1 i i+1 i
*
* With some algebric manipulation, it is not difficult to see
* that (2) is equivalent to
* -(i+1)
* s + 2 <= y (3)
* i i
*
* The advantage of (3) is that s and y can be computed by
* i i
* the following recurrence formula:
* if (3) is false
*
* s = s , y = y ; (4)
* i+1 i i+1 i
*
* otherwise,
* -i -(i+1)
* s = s + 2 , y = y - s - 2 (5)
* i+1 i i+1 i i
*
* One may easily use induction to prove (4) and (5).
* Note. Since the left hand side of (3) contain only i+2 bits,
* it does not necessary to do a full (53-bit) comparison
* in (3).
* 3. Final rounding
* After generating the 53 bits result, we compute one more bit.
* Together with the remainder, we can decide whether the
* result is exact, bigger than 1/2ulp, or less than 1/2ulp
* (it will never equal to 1/2ulp).
* The rounding mode can be detected by checking whether
* huge + tiny is equal to huge, and whether huge - tiny is
* equal to huge for some floating point number "huge" and "tiny".
*
* Special cases:
* sqrt(+-0) = +-0 ... exact
* sqrt(inf) = inf
* sqrt(-ve) = NaN ... with invalid signal
* sqrt(NaN) = NaN ... with invalid signal for signaling NaN
*
* Other methods : see the appended file at the end of the program below.
*---------------
*/
#include "errno.h"
#include "math.h"
#ifdef __STDC__
static const double one = 1.0, tiny = 1.0e-300;
#else
static double one = 1.0, tiny = 1.0e-300;
#endif
#ifdef __STDC__
double __ieee754_sqrt(double x)
#else
double __ieee754_sqrt(x) double x;
#endif
{
double z;
int sign = (int)0x80000000;
unsigned r, t1, s1, ix1, q1;
int ix0, s0, q, m, t, i;
ix0 = __HI(x); /* high word of x */
ix1 = __LO(x); /* low word of x */
/* take care of Inf and NaN */
if ((ix0 & 0x7ff00000) == 0x7ff00000) {
errno = 33;
return x * x + x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf
sqrt(-inf)=sNaN */
}
/* take care of zero */
if (ix0 <= 0) {
if (((ix0 & (~sign)) | ix1) == 0)
return x; /* sqrt(+-0) = +-0 */
else if (ix0 < 0) {
errno = 33;
return NAN;
} /* sqrt(-ve) = sNaN */
}
/* normalize x */
m = (ix0 >> 20);
if (m == 0) { /* subnormal x */
while (ix0 == 0) {
m -= 21;
ix0 |= (ix1 >> 11);
ix1 <<= 21;
}
for (i = 0; (ix0 & 0x00100000) == 0; i++)
ix0 <<= 1;
m -= i - 1;
ix0 |= (ix1 >> (32 - i));
ix1 <<= i;
}
m -= 1023; /* unbias exponent */
ix0 = (ix0 & 0x000fffff) | 0x00100000;
if (m & 1) { /* odd m, double x to make it even */
ix0 += ix0 + ((ix1 & sign) >> 31);
ix1 += ix1;
}
m >>= 1; /* m = [m/2] */
/* generate sqrt(x) bit by bit */
ix0 += ix0 + ((ix1 & sign) >> 31);
ix1 += ix1;
q = q1 = s0 = s1 = 0; /* [q,q1] = sqrt(x) */
r = 0x00200000; /* r = moving bit from right to left */
while (r != 0) {
t = s0 + r;
if (t <= ix0) {
s0 = t + r;
ix0 -= t;
q += r;
}
ix0 += ix0 + ((ix1 & sign) >> 31);
ix1 += ix1;
r >>= 1;
}
r = sign;
while (r != 0) {
t1 = s1 + r;
t = s0;
if ((t < ix0) || ((t == ix0) && (t1 <= ix1))) {
s1 = t1 + r;
if (((t1 & sign) == sign) && (s1 & sign) == 0)
s0 += 1;
ix0 -= t;
if (ix1 < t1)
ix0 -= 1;
ix1 -= t1;
q1 += r;
}
ix0 += ix0 + ((ix1 & sign) >> 31);
ix1 += ix1;
r >>= 1;
}
/* use floating add to find out rounding direction */
if ((ix0 | ix1) != 0) {
z = one - tiny; /* trigger inexact flag */
if (z >= one) {
z = one + tiny;
if (q1 == (unsigned)0xffffffff) {
q1 = 0;
q += 1;
} else if (z > one) {
if (q1 == (unsigned)0xfffffffe)
q += 1;
q1 += 2;
} else
q1 += (q1 & 1);
}
}
ix0 = (q >> 1) + 0x3fe00000;
ix1 = q1 >> 1;
if ((q & 1) == 1)
ix1 |= sign;
ix0 += (m << 20);
__HI(z) = ix0;
__LO(z) = ix1;
return z;
}
/*
Other methods (use floating-point arithmetic)
-------------
(This is a copy of a drafted paper by Prof W. Kahan
and K.C. Ng, written in May, 1986)
Two algorithms are given here to implement sqrt(x)
(IEEE double precision arithmetic) in software.
Both supply sqrt(x) correctly rounded. The first algorithm (in
Section A) uses newton iterations and involves four divisions.
The second one uses reciproot iterations to avoid division, but
requires more multiplications. Both algorithms need the ability
to chop results of arithmetic operations instead of round them,
and the INEXACT flag to indicate when an arithmetic operation
is executed exactly with no roundoff error, all part of the
standard (IEEE 754-1985). The ability to perform shift, add,
subtract and logical AND operations upon 32-bit words is needed
too, though not part of the standard.
A. sqrt(x) by Newton Iteration
(1) Initial approximation
Let x0 and x1 be the leading and the trailing 32-bit words of
a floating point number x (in IEEE double format) respectively
1 11 52 ...widths
------------------------------------------------------
x: |s| e | f |
------------------------------------------------------
msb lsb msb lsb ...order
------------------------ ------------------------
x0: |s| e | f1 | x1: | f2 |
------------------------ ------------------------
By performing shifts and subtracts on x0 and x1 (both regarded
as integers), we obtain an 8-bit approximation of sqrt(x) as
follows.
k := (x0>>1) + 0x1ff80000;
y0 := k - T1[31&(k>>15)]. ... y ~ sqrt(x) to 8 bits
Here k is a 32-bit integer and T1[] is an integer array containing
correction terms. Now magically the floating value of y (y's
leading 32-bit word is y0, the value of its trailing word is 0)
approximates sqrt(x) to almost 8-bit.
Value of T1:
static int T1[32]= {
0, 1024, 3062, 5746, 9193, 13348, 18162, 23592,
29598, 36145, 43202, 50740, 58733, 67158, 75992, 85215,
83599, 71378, 60428, 50647, 41945, 34246, 27478, 21581,
16499, 12183, 8588, 5674, 3403, 1742, 661, 130,};
(2) Iterative refinement
Apply Heron's rule three times to y, we have y approximates
sqrt(x) to within 1 ulp (Unit in the Last Place):
y := (y+x/y)/2 ... almost 17 sig. bits
y := (y+x/y)/2 ... almost 35 sig. bits
y := y-(y-x/y)/2 ... within 1 ulp
Remark 1.
Another way to improve y to within 1 ulp is:
y := (y+x/y) ... almost 17 sig. bits to 2*sqrt(x)
y := y - 0x00100006 ... almost 18 sig. bits to sqrt(x)
2
(x-y )*y
y := y + 2* ---------- ...within 1 ulp
2
3y + x
This formula has one division fewer than the one above; however,
it requires more multiplications and additions. Also x must be
scaled in advance to avoid spurious overflow in evaluating the
expression 3y*y+x. Hence it is not recommended uless division
is slow. If division is very slow, then one should use the
reciproot algorithm given in section B.
(3) Final adjustment
By twiddling y's last bit it is possible to force y to be
correctly rounded according to the prevailing rounding mode
as follows. Let r and i be copies of the rounding mode and
inexact flag before entering the square root program. Also we
use the expression y+-ulp for the next representable floating
numbers (up and down) of y. Note that y+-ulp = either fixed
point y+-1, or multiply y by nextafter(1,+-inf) in chopped
mode.
I := FALSE; ... reset INEXACT flag I
R := RZ; ... set rounding mode to round-toward-zero
z := x/y; ... chopped quotient, possibly inexact
If(not I) then { ... if the quotient is exact
if(z=y) {
I := i; ... restore inexact flag
R := r; ... restore rounded mode
return sqrt(x):=y.
} else {
z := z - ulp; ... special rounding
}
}
i := TRUE; ... sqrt(x) is inexact
If (r=RN) then z=z+ulp ... rounded-to-nearest
If (r=RP) then { ... round-toward-+inf
y = y+ulp; z=z+ulp;
}
y := y+z; ... chopped sum
y0:=y0-0x00100000; ... y := y/2 is correctly rounded.
I := i; ... restore inexact flag
R := r; ... restore rounded mode
return sqrt(x):=y.
(4) Special cases
Square root of +inf, +-0, or NaN is itself;
Square root of a negative number is NaN with invalid signal.
B. sqrt(x) by Reciproot Iteration
(1) Initial approximation
Let x0 and x1 be the leading and the trailing 32-bit words of
a floating point number x (in IEEE double format) respectively
(see section A). By performing shifs and subtracts on x0 and y0,
we obtain a 7.8-bit approximation of 1/sqrt(x) as follows.
k := 0x5fe80000 - (x0>>1);
y0:= k - T2[63&(k>>14)]. ... y ~ 1/sqrt(x) to 7.8 bits
Here k is a 32-bit integer and T2[] is an integer array
containing correction terms. Now magically the floating
value of y (y's leading 32-bit word is y0, the value of
its trailing word y1 is set to zero) approximates 1/sqrt(x)
to almost 7.8-bit.
Value of T2:
static int T2[64]= {
0x1500, 0x2ef8, 0x4d67, 0x6b02, 0x87be, 0xa395, 0xbe7a, 0xd866,
0xf14a, 0x1091b,0x11fcd,0x13552,0x14999,0x15c98,0x16e34,0x17e5f,
0x18d03,0x19a01,0x1a545,0x1ae8a,0x1b5c4,0x1bb01,0x1bfde,0x1c28d,
0x1c2de,0x1c0db,0x1ba73,0x1b11c,0x1a4b5,0x1953d,0x18266,0x16be0,
0x1683e,0x179d8,0x18a4d,0x19992,0x1a789,0x1b445,0x1bf61,0x1c989,
0x1d16d,0x1d77b,0x1dddf,0x1e2ad,0x1e5bf,0x1e6e8,0x1e654,0x1e3cd,
0x1df2a,0x1d635,0x1cb16,0x1be2c,0x1ae4e,0x19bde,0x1868e,0x16e2e,
0x1527f,0x1334a,0x11051,0xe951, 0xbe01, 0x8e0d, 0x5924, 0x1edd,};
(2) Iterative refinement
Apply Reciproot iteration three times to y and multiply the
result by x to get an approximation z that matches sqrt(x)
to about 1 ulp. To be exact, we will have
-1ulp < sqrt(x)-z<1.0625ulp.
... set rounding mode to Round-to-nearest
y := y*(1.5-0.5*x*y*y) ... almost 15 sig. bits to 1/sqrt(x)
y := y*((1.5-2^-30)+0.5*x*y*y)... about 29 sig. bits to 1/sqrt(x)
... special arrangement for better accuracy
z := x*y ... 29 bits to sqrt(x), with z*y<1
z := z + 0.5*z*(1-z*y) ... about 1 ulp to sqrt(x)
Remark 2. The constant 1.5-2^-30 is chosen to bias the error so that
(a) the term z*y in the final iteration is always less than 1;
(b) the error in the final result is biased upward so that
-1 ulp < sqrt(x) - z < 1.0625 ulp
instead of |sqrt(x)-z|<1.03125ulp.
(3) Final adjustment
By twiddling y's last bit it is possible to force y to be
correctly rounded according to the prevailing rounding mode
as follows. Let r and i be copies of the rounding mode and
inexact flag before entering the square root program. Also we
use the expression y+-ulp for the next representable floating
numbers (up and down) of y. Note that y+-ulp = either fixed
point y+-1, or multiply y by nextafter(1,+-inf) in chopped
mode.
R := RZ; ... set rounding mode to round-toward-zero
switch(r) {
case RN: ... round-to-nearest
if(x<= z*(z-ulp)...chopped) z = z - ulp; else
if(x<= z*(z+ulp)...chopped) z = z; else z = z+ulp;
break;
case RZ:case RM: ... round-to-zero or round-to--inf
R:=RP; ... reset rounding mod to round-to-+inf
if(x<z*z ... rounded up) z = z - ulp; else
if(x>=(z+ulp)*(z+ulp) ...rounded up) z = z+ulp;
break;
case RP: ... round-to-+inf
if(x>(z+ulp)*(z+ulp)...chopped) z = z+2*ulp; else
if(x>z*z ...chopped) z = z+ulp;
break;
}
Remark 3. The above comparisons can be done in fixed point. For
example, to compare x and w=z*z chopped, it suffices to compare
x1 and w1 (the trailing parts of x and w), regarding them as
two's complement integers.
...Is z an exact square root?
To determine whether z is an exact square root of x, let z1 be the
trailing part of z, and also let x0 and x1 be the leading and
trailing parts of x.
If ((z1&0x03ffffff)!=0) ... not exact if trailing 26 bits of z!=0
I := 1; ... Raise Inexact flag: z is not exact
else {
j := 1 - [(x0>>20)&1] ... j = logb(x) mod 2
k := z1 >> 26; ... get z's 25-th and 26-th
fraction bits
I := i or (k&j) or ((k&(j+j+1))!=(x1&3));
}
R:= r ... restore rounded mode
return sqrt(x):=z.
If multiplication is cheaper then the foregoing red tape, the
Inexact flag can be evaluated by
I := i;
I := (z*z!=x) or I.
Note that z*z can overwrite I; this value must be sensed if it is
True.
Remark 4. If z*z = x exactly, then bit 25 to bit 0 of z1 must be
zero.
--------------------
z1: | f2 |
--------------------
bit 31 bit 0
Further more, bit 27 and 26 of z1, bit 0 and 1 of x1, and the odd
or even of logb(x) have the following relations:
-------------------------------------------------
bit 27,26 of z1 bit 1,0 of x1 logb(x)
-------------------------------------------------
00 00 odd and even
01 01 even
10 10 odd
10 00 even
11 01 even
-------------------------------------------------
(4) Special cases (see (4) of Section A).
*/
@@ -0,0 +1,92 @@
/* @(#)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"
#ifdef __STDC__
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 */
#ifdef __STDC__
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));
}
}
@@ -0,0 +1,352 @@
/* @(#)k_rem_pio2.c 1.2 95/01/04 */
/* $Id: k_rem_pio2.c,v 1.2.14.1 2002/01/31 15:24:13 ceciliar Exp $ */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/*
* __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 = ldexp(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 ldexp(), 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"
#ifdef __STDC__
static const int init_jk[] = { 2, 3, 4, 6 };
/* initial value for jk */ /*- cc 020130 -*/
#else
static int init_jk[] = { 2, 3, 4, 6 }; /*- cc 020130 -*/
#endif
#ifdef __STDC__
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 */
};
#ifdef __STDC__
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 */
#ifdef __STDC__
int __kernel_rem_pio2(double* x, double* y, int e0, int nx, int prec, const int* ipio2) /*- cc 020130 -*/
#else
int __kernel_rem_pio2(x, y, e0, nx, prec, ipio2) /*- cc 020130 -*/
double x[],
y[];
int e0, nx, prec;
int ipio2[]; /*- cc 020130 -*/
#endif
{
int jz, jx, jv, jp, jk, carry, n, iq[20], i, j, k, m, q0, ih; /*- cc 020130 -*/
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)); /*- cc 020130 -*/
iq[i] = (int)(z - two24 * fw); /*- cc 020130 -*/
z = q[j - 1] + fw;
}
/* compute n */
z = ldexp(z, q0); /* actual value of z */
z -= 8.0 * floor(z * 0.125); /* trim off integer >= 8 */
n = (int)z; /*- cc 020130 -*/
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 -= ldexp(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 = ldexp(z, -q0);
if (z >= two24) {
fw = (double)((int)(twon24 * z)); /*- cc 020130 -*/
iq[jz] = (int)(z - two24 * fw); /*- cc 020130 -*/
jz += 1;
q0 += 24;
iq[jz] = (int)fw; /*- cc 020130 -*/
} else
iq[jz] = (int)z; /*- cc 020130 -*/
}
/* convert integer "bit" chunk to floating-point value */
fw = ldexp(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;
}
@@ -0,0 +1,79 @@
/* @(#)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"
#ifdef __STDC__
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 */
#ifdef __STDC__
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);
}
@@ -0,0 +1,180 @@
//===========================================================================
//
// k_tan.c
//
// Part of the standard mathematical function library
//
//===========================================================================
//####ECOSGPLCOPYRIGHTBEGIN####
// -------------------------------------------
// This file is part of eCos, the Embedded Configurable Operating System.
// Copyright (C) 1998, 1999, 2000, 2001, 2002 Red Hat, Inc.
//
// eCos is free software; you can redistribute it and/or modify it under
// the terms of the GNU General Public License as published by the Free
// Software Foundation; either version 2 or (at your option) any later version.
//
// eCos is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or
// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
// for more details.
//
// You should have received a copy of the GNU General Public License along
// with eCos; if not, write to the Free Software Foundation, Inc.,
// 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA.
//
// As a special exception, if other files instantiate templates or use macros
// or inline functions from this file, or you compile this file and link it
// with other works to produce a work based on this file, this file does not
// by itself cause the resulting work to be covered by the GNU General Public
// License. However the source code for this file must still be made available
// in accordance with section (3) of the GNU General Public License.
//
// This exception does not invalidate any other reasons why a work based on
// this file might be covered by the GNU General Public License.
//
// Alternative licenses for eCos may be arranged by contacting Red Hat, Inc.
// at http://sources.redhat.com/ecos/ecos-license/
// -------------------------------------------
//####ECOSGPLCOPYRIGHTEND####
//===========================================================================
//#####DESCRIPTIONBEGIN####
//
// Author(s): jlarmour
// Contributors: jlarmour
// Date: 1998-02-13
// Purpose:
// Description:
// Usage:
//
//####DESCRIPTIONEND####
//
//===========================================================================
// Derived from code with the following copyright
/* @(#)k_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.
* ====================================================
*/
/* __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 "math.h"
static const double one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
pio4 = 7.85398163397448278999e-01, /* 0x3FE921FB, 0x54442D18 */
pio4lo = 3.06161699786838301793e-17, /* 0x3C81A626, 0x33145C07 */
T[] = {
3.33333333333334091986e-01, /* 0x3FD55555, 0x55555563 */
1.33333333333201242699e-01, /* 0x3FC11111, 0x1110FE7A */
5.39682539762260521377e-02, /* 0x3FABA1BA, 0x1BB341FE */
2.18694882948595424599e-02, /* 0x3F9664F4, 0x8406D637 */
8.86323982359930005737e-03, /* 0x3F8226E3, 0xE96E8493 */
3.59207910759131235356e-03, /* 0x3F6D6D22, 0xC9560328 */
1.45620945432529025516e-03, /* 0x3F57DBC8, 0xFEE08315 */
5.88041240820264096874e-04, /* 0x3F4344D8, 0xF2F26501 */
2.46463134818469906812e-04, /* 0x3F3026F7, 0x1A8D1068 */
7.81794442939557092300e-05, /* 0x3F147E88, 0xA03792A6 */
7.14072491382608190305e-05, /* 0x3F12B80F, 0x32F0A7E9 */
-1.85586374855275456654e-05, /* 0xBEF375CB, 0xDB605373 */
2.59073051863633712884e-05, /* 0x3EFB2A70, 0x74BF7AD4 */
};
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) {
double ret = fabs(x);
return one / ret;
} else
return (iy == 1) ? x : -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);
}
}
// EOF k_tan.c
@@ -0,0 +1,142 @@
/* @(#)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"
#ifdef __STDC__
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 */
};
#ifdef __STDC__
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 */
};
#ifdef __STDC__
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 */
};
#ifdef __STDC__
static const double
#else
static double
#endif
one
= 1.0,
huge = 1.0e300;
#ifdef __STDC__
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;
}
}
@@ -0,0 +1,89 @@
/* @(#)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"
#ifdef __STDC__
static const double huge = 1.0e300;
#else
static double huge = 1.0e300;
#endif
#ifdef __STDC__
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;
}
@@ -0,0 +1,29 @@
/* @(#)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"
#ifdef __STDC__
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;
}
@@ -0,0 +1,81 @@
/* @(#)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"
#ifdef __STDC__
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);
}
}
}
@@ -0,0 +1,88 @@
/* @(#)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"
#ifdef __STDC__
static const double huge = 1.0e300;
#else
static double huge = 1.0e300;
#endif
#ifdef __STDC__
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;
}
@@ -0,0 +1,56 @@
/* @(#)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"
#ifdef __STDC__
static const double
#else
static double
#endif
two54 = 1.80143985094819840000e+16; /* 0x43500000, 0x00000000 */
#ifdef __STDC__
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;
}
@@ -0,0 +1,58 @@
/* @(#)s_ldexp.c 1.2 95/01/04 */
/* $Id: s_ldexp.c,v 1.3.14.1 2002/01/31 15:24:14 ceciliar Exp $ */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#include "math.h" /* for isfinite macro */
static const double
two54
= 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
big = 1.0e+300, tiny = 1.0e-300;
double ldexp(double x, int n)
{
long k, hx, lx; /*- cc 020130 -*/
if (!isfinite(x) || x == 0.0)
return x;
hx = __HI(x);
lx = __LO(x);
k = (hx & 0x7ff00000) >> 20; /* extract exponent */
if (k == 0) { /* 0 or subnormal x */
if ((lx | (hx & 0x7fffffff)) == 0)
return x; /* +-0 */
x *= two54;
hx = __HI(x);
k = ((hx & 0x7ff00000) >> 20) - 54;
if (n < -50000)
return tiny * x; /*underflow*/
}
if (k == 0x7ff)
return x + x; /* NaN or Inf */
k = k + n;
if (k > 0x7fe)
return big * copysign(big, x); /* overflow */
if (k > 0) /* normal result */
{
__HI(x) = (hx & 0x800fffff) | (k << 20);
return x;
}
if (k <= -54)
if (n > 50000) /* in case integer overflow in n+k */
return big * copysign(big, x); /*overflow*/
else
return tiny * copysign(tiny, x); /*underflow*/
k += 54; /* subnormal result */
__HI(x) = (hx & 0x800fffff) | (k << 20);
return x * twom54;
}
@@ -0,0 +1,78 @@
/* @(#)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"
#ifdef __STDC__
static const double one = 1.0;
#else
static double one = 1.0;
#endif
#ifdef __STDC__
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;
}
}
}
@@ -0,0 +1,81 @@
/* @(#)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"
#ifdef __STDC__
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]);
}
}
}
@@ -0,0 +1,72 @@
/* @(#)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"
#ifdef __STDC__
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 */
}
}
@@ -0,0 +1,6 @@
#include "fdlibm.h"
/* 8036C6E0-8036C700 367020 0020+00 0/0 2/2 3/3 .text acos */
double acos(double x) {
return __ieee754_acos(x);
}
@@ -0,0 +1,6 @@
#include "fdlibm.h"
/* 8036C700-8036C720 367040 0020+00 0/0 2/2 0/0 .text asin */
double asin(double x) {
return __ieee754_asin(x);
}
@@ -0,0 +1,6 @@
#include "fdlibm.h"
/* 8036C720-8036C740 367060 0020+00 0/0 6/6 0/0 .text atan2 */
double atan2(double x, double y) {
return __ieee754_atan2(x, y);
}
@@ -0,0 +1,6 @@
#include "fdlibm.h"
/* 8036C740-8036C760 367080 0020+00 0/0 1/1 0/0 .text exp */
double exp(double x) {
return __ieee754_exp(x);
}
@@ -0,0 +1,6 @@
#include "fdlibm.h"
/* 8036C760-8036C780 3670A0 0020+00 0/0 8/8 0/0 .text fmod */
double fmod(double x, double y) {
return __ieee754_fmod(x, y);
}
@@ -0,0 +1,6 @@
#include "fdlibm.h"
/* 8036C780-8036C7A0 3670C0 0020+00 0/0 3/3 24/24 .text pow */
double pow(double x, double y) {
return __ieee754_pow(x, y);
}
@@ -0,0 +1,6 @@
#include "fdlibm.h"
/* 8036CA54-8036CA74 367394 0020+00 0/0 8/8 1/1 .text sqrt */
double sqrt(double x) {
return __ieee754_sqrt(x);
}
@@ -0,0 +1,228 @@
#ifndef FDLIBM_H
#define FDLIBM_H
/* @(#)fdlibm.h 1.5 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.
* ====================================================
*/
#ifdef __cplusplus
extern "C" {
#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
#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
// TODO: should __STDC__ actually be defined?
//#ifdef __STDC__
#define __P(p) p
//#else
//#define __P(p) ()
//#endif
/*
* ANSI/POSIX
*/
extern int signgam;
#define MAXFLOAT ((float)3.40282346638528860e+38)
enum fdversion { fdlibm_ieee = -1, fdlibm_svid, fdlibm_xopen, fdlibm_posix };
#define _LIB_VERSION_TYPE enum fdversion
#define _LIB_VERSION _fdlib_version
/* if global variable _LIB_VERSION is not desirable, one may
* change the following to be a constant by:
* #define _LIB_VERSION_TYPE const enum version
* In that case, after one initializes the value _LIB_VERSION (see
* s_lib_version.c) during compile time, it cannot be modified
* in the middle of a program
*/
extern _LIB_VERSION_TYPE _LIB_VERSION;
#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;
};
#define HUGE MAXFLOAT
/*
* set X_TLOSS = pi*2**52, which is possibly defined in <values.h"
* (one may replace the following line by "#include "values.h"")
*/
#define X_TLOSS 1.41484755040568800000e+16
#define DOMAIN 1
#define SING 2
#define OVERFLOW 3
#define UNDERFLOW 4
#define TLOSS 5
#define PLOSS 6
/*
* ANSI/POSIX
*/
extern double acos __P((double));
extern double asin __P((double));
extern double atan __P((double));
extern double atan2 __P((double, double));
extern double cos __P((double));
extern double sin __P((double));
extern double tan __P((double));
extern double cosh __P((double));
extern double sinh __P((double));
extern double tanh __P((double));
extern double exp __P((double));
extern double frexp __P((double, int*));
extern double ldexp __P((double, int));
extern double log __P((double));
extern double log10 __P((double));
extern double modf __P((double, double*));
extern double pow __P((double, double));
extern double sqrt __P((double));
extern double ceil __P((double));
extern double fabs __P((double));
extern double floor __P((double));
extern double fmod __P((double, double));
extern double erf __P((double));
extern double erfc __P((double));
extern double gamma __P((double));
extern double hypot __P((double, double));
extern int isnan __P((double));
extern int finite __P((double));
extern double j0 __P((double));
extern double j1 __P((double));
extern double jn __P((int, double));
extern double lgamma __P((double));
extern double y0 __P((double));
extern double y1 __P((double));
extern double yn __P((int, double));
extern double acosh __P((double));
extern double asinh __P((double));
extern double atanh __P((double));
extern double cbrt __P((double));
extern double logb __P((double));
extern double nextafter __P((double, double));
extern double remainder __P((double, double));
#ifdef _SCALB_INT
extern double scalb __P((double, int));
#else
extern double scalb __P((double, double));
#endif
extern int matherr __P((struct exception*));
/*
* IEEE Test Vector
*/
extern double significand __P((double));
/*
* Functions callable from C, intended to support IEEE arithmetic.
*/
extern double copysign __P((double, double));
extern int ilogb __P((double));
extern double rint __P((double));
extern double scalbn __P((double, int));
/*
* BSD math library entry points
*/
extern double expm1 __P((double));
extern double log1p __P((double));
/*
* Reentrant version of gamma & lgamma; passes signgam back by reference
* 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*));
#endif /* _REENTRANT */
/* ieee style elementary functions */
extern double __ieee754_sqrt __P((double));
extern double __ieee754_acos __P((double));
extern double __ieee754_acosh __P((double));
extern double __ieee754_log __P((double));
extern double __ieee754_atanh __P((double));
extern double __ieee754_asin __P((double));
extern double __ieee754_atan2 __P((double, double));
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 __P((double));
extern double __ieee754_gamma __P((double));
extern double __ieee754_log10 __P((double));
extern double __ieee754_sinh __P((double));
extern double __ieee754_hypot __P((double, double));
extern double __ieee754_j0 __P((double));
extern double __ieee754_j1 __P((double));
extern double __ieee754_y0 __P((double));
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*));
#ifdef _SCALB_INT
extern double __ieee754_scalb __P((double, int));
#else
extern double __ieee754_scalb __P((double, double));
#endif
/* fdlibm kernel function */
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*));
#ifdef __cplusplus
};
#endif
#endif /* FDLIBM_H */
@@ -0,0 +1,573 @@
#include "ansi_fp.h"
#include "ansi_fp.h"
#include "limits.h"
#include "float.h"
static int __count_trailing_zerol(unsigned long x) {
int result = 0;
int bits_not_checked = sizeof(unsigned long) * CHAR_BIT;
int n = bits_not_checked / 2;
int mask_size = n;
unsigned long mask = (~0UL) >> (bits_not_checked - n);
while (bits_not_checked) {
if (!(x & mask)) {
result += mask_size;
x >>= mask_size;
bits_not_checked -= mask_size;
} else if (mask == 1) {
break;
}
if (n > 1) {
n /= 2;
}
if (mask > 1) {
mask >>= n;
mask_size -= n;
}
}
return result;
}
static int __count_trailing_zero(double x) {
unsigned long* l = (unsigned long*)&x;
if (l[1] != 0) {
return __count_trailing_zerol(l[1]);
}
return (int)(sizeof(unsigned long) * CHAR_BIT + __count_trailing_zerol(l[0] | 0x00100000));
}
static int __must_round(const decimal* d, int digits) {
unsigned char const* i = d->sig.text + digits;
if (*i > 5) {
return 1;
}
if (*i < 5) {
return -1;
}
{
unsigned char const* e = d->sig.text + d->sig.length;
for (i++; i < e; i++) {
if (*i != 0) {
return 1;
}
}
}
if (d->sig.text[digits - 1] & 1) {
return 1;
}
return -1;
}
static void __dorounddecup(decimal* d, int digits) {
unsigned char* b = d->sig.text;
unsigned char* i = b + digits - 1;
while (1) {
if (*i < 9) {
*i += 1;
break;
}
if (i == b) {
*i = 1;
d->exp++;
break;
}
*i-- = 0;
}
}
static void __rounddec(decimal* d, int digits) {
if (digits > 0 && digits < d->sig.length) {
int unkBool = __must_round(d, digits);
d->sig.length = digits;
if (unkBool >= 0) {
__dorounddecup(d, digits);
}
}
}
void __ull2dec(decimal* result, unsigned long long val) {
result->sign = 0;
if (val == 0) {
result->exp = 0;
result->sig.length = 1;
result->sig.text[0] = 0;
return;
}
if (val < 0) {
val = -val;
result->sign = 1;
}
result->sig.length = 0;
for (; val != 0; val /= 10) {
result->sig.text[result->sig.length++] = (unsigned char)(val % 10);
}
{
unsigned char* i = result->sig.text;
unsigned char* j = result->sig.text + result->sig.length;
for (; i < --j; ++i) {
unsigned char t = *i;
*i = *j;
*j = t;
}
}
result->exp = result->sig.length - 1;
}
/* 80364E00-80365078 35F740 0278+00 2/2 0/0 0/0 .text __timesdec */
void __timesdec(decimal* result, const decimal* x, const decimal* y) {
unsigned long accumulator = 0;
unsigned char mantissa[SIGDIGLEN * 2];
int i = x->sig.length + y->sig.length - 1;
unsigned char* pDigit;
unsigned char* ip = mantissa + i + 1;
unsigned char* ep = ip;
result->sign = 0;
for (; i > 0; i--) {
int k = y->sig.length - 1;
int j = i - k - 1;
int l;
int t;
const unsigned char* jp;
const unsigned char* kp;
if (j < 0) {
j = 0;
k = i - 1;
}
jp = x->sig.text + j;
kp = y->sig.text + k;
l = k + 1;
t = x->sig.length - j;
if (l > t)
l = t;
for (; l > 0; l--, jp++, kp--) {
accumulator += *jp * *kp;
}
*--ip = (unsigned char)(accumulator % 10);
accumulator /= 10;
}
result->exp = (short)(x->exp + y->exp);
if (accumulator) {
*--ip = (unsigned char)(accumulator);
result->exp++;
}
for (i = 0; i < SIGDIGLEN && ip < ep; i++, ip++) {
result->sig.text[i] = *ip;
}
result->sig.length = (unsigned char)(i);
if (ip < ep && *ip >= 5) {
if (*ip == 5) {
unsigned char* jp = ip + 1;
for (; jp < ep; jp++) {
if (*jp != 0)
goto round;
}
if ((ip[-1] & 1) == 0)
return;
}
round:
__dorounddecup(result, result->sig.length);
}
}
void __str2dec(decimal* d, const char* s, short exp) {
int i;
d->exp = exp;
d->sign = 0;
for (i = 0; i < SIGDIGLEN && *s;) {
d->sig.text[i++] = *s++ - '0';
}
d->sig.length = i;
if (*s != 0) {
if (*s < 5)
return;
if (*s > 5)
goto round;
{
const char* p = s + 1;
for (; *p != 0; p++) {
if (*p != '0')
goto round;
}
if ((d->sig.text[i - 1] & 1) == 0)
return;
}
round:
__dorounddecup(d, d->sig.length);
}
}
static const char* const unused = "179769313486231580793729011405303420";
/* 8036367C-80364E00 35DFBC 1784+00 2/1 0/0 0/0 .text __two_exp */
void __two_exp(decimal* result, long exp) {
switch (exp) {
case -64:
__str2dec(result, "542101086242752217003726400434970855712890625", -20);
return;
case -53:
__str2dec(result, "11102230246251565404236316680908203125", -16);
return;
case -32:
__str2dec(result, "23283064365386962890625", -10);
return;
case -16:
__str2dec(result, "152587890625", -5);
return;
case -8:
__str2dec(result, "390625", -3);
return;
case -7:
__str2dec(result, "78125", -3);
return;
case -6:
__str2dec(result, "15625", -2);
return;
case -5:
__str2dec(result, "3125", -2);
return;
case -4:
__str2dec(result, "625", -2);
return;
case -3:
__str2dec(result, "125", -1);
return;
case -2:
__str2dec(result, "25", -1);
return;
case -1:
__str2dec(result, "5", -1);
return;
case 0:
__str2dec(result, "1", 0);
return;
case 1:
__str2dec(result, "2", 0);
return;
case 2:
__str2dec(result, "4", 0);
return;
case 3:
__str2dec(result, "8", 0);
return;
case 4:
__str2dec(result, "16", 1);
return;
case 5:
__str2dec(result, "32", 1);
return;
case 6:
__str2dec(result, "64", 1);
return;
case 7:
__str2dec(result, "128", 2);
return;
case 8:
__str2dec(result, "256", 2);
return;
}
{
decimal x2, temp;
__two_exp(&x2, exp / 2);
__timesdec(result, &x2, &x2);
if (exp & 1) {
temp = *result;
if (exp > 0) {
__str2dec(&x2, "2", 0);
} else {
__str2dec(&x2, "5", -1);
}
__timesdec(result, &temp, &x2);
}
}
}
int __equals_dec(const decimal* x, const decimal* y) {
if (x->sig.text[0] == 0) {
if (y->sig.text[0] == 0)
return 1;
return 0;
}
if (y->sig.text[0] == 0) {
if (x->sig.text[0] == 0)
return 1;
return 0;
}
if (x->exp == y->exp) {
int i;
int l = x->sig.length;
if (l > y->sig.length) {
l = y->sig.length;
}
for (i = 0; i < l; i++) {
if (x->sig.text[i] != y->sig.text[i]) {
return 0;
}
}
if (l == x->sig.length) {
for (; i < y->sig.length; ++i) {
if (y->sig.text[i] != 0) {
return 0;
}
}
} else {
for (; i < x->sig.length; ++i) {
if (x->sig.text[i] != 0) {
return 0;
}
}
}
return 1;
}
return 0;
}
int __less_dec(const decimal* x, const decimal* y) {
if (x->sig.text[0] == 0) {
if (y->sig.text[0] != 0)
return 1;
return 0;
}
if (y->sig.text[0] == 0) {
return 0;
}
if (x->exp == y->exp) {
int i;
int l = x->sig.length;
if (l > y->sig.length) {
l = y->sig.length;
}
for (i = 0; i < l; i++) {
if (x->sig.text[i] < y->sig.text[i]) {
return 1;
} else if (y->sig.text[i] < x->sig.text[i]) {
return 0;
}
}
if (l == x->sig.length) {
for (; i < y->sig.length; i++) {
if (y->sig.text[i] != 0) {
return 1;
}
}
}
return 0;
}
return x->exp < y->exp;
}
void __minus_dec(decimal* z, const decimal* x, const decimal* y) {
int zlen, dexp;
unsigned char *ib, *i, *ie;
unsigned char const *jb, *j, *jn;
*z = *x;
if (y->sig.text[0] == 0)
return;
zlen = z->sig.length;
if (zlen < y->sig.length)
zlen = y->sig.length;
dexp = z->exp - y->exp;
zlen += dexp;
if (zlen > SIGDIGLEN)
zlen = SIGDIGLEN;
while (z->sig.length < zlen) {
z->sig.text[z->sig.length++] = 0;
}
ib = z->sig.text;
i = ib + zlen;
if (y->sig.length + dexp < zlen) {
i = ib + (y->sig.length + dexp);
}
jb = y->sig.text;
j = jb + (i - ib - dexp);
jn = j;
while (i > ib && j > jb) {
i--;
j--;
if (*i < *j) {
unsigned char* k = i - 1;
while (*k == 0)
k--;
while (k != i) {
--*k;
*++k += 10;
}
}
*i -= *j;
}
if (jn - jb < y->sig.length) {
int round_down = 0;
if (*jn < 5)
round_down = 1;
else if (*jn == 5) {
unsigned char const* ibPtr = y->sig.text + y->sig.length;
for (j = jn + 1; j < ibPtr; j++) {
if (*j != 0)
goto done;
}
i = ib + (jn - jb) + dexp - 1;
if (*i & 1)
round_down = 1;
}
if (round_down) {
if (*i < 1) {
unsigned char* k = i - 1;
while (*k == 0)
k--;
while (k != i) {
--*k;
*++k += 10;
}
}
*i -= 1;
}
}
done:
for (i = ib; *i == 0; ++i) {
}
if (i > ib) {
unsigned char dl = (unsigned char)(i - ib);
z->exp -= dl;
ie = ib + z->sig.length;
for (; i < ie; ++i, ++ib)
*ib = *i;
z->sig.length -= dl;
}
ib = z->sig.text;
for (i = ib + z->sig.length; i > ib;) {
i--;
if (*i != 0)
break;
}
z->sig.length = (unsigned char)(i - ib + 1);
}
/* 803632C8-8036367C 35DC08 03B4+00 1/1 0/0 0/0 .text __num2dec_internal */
void __num2dec_internal(decimal* d, double x) {
signed char sign = (signed char)(signbit(x) != 0);
if (x == 0) {
d->sign = sign;
d->exp = 0;
d->sig.length = 1;
d->sig.text[0] = 0;
return;
}
if (!isfinite(x)) {
d->sign = sign;
d->exp = 0;
d->sig.length = 1;
d->sig.text[0] = fpclassify(x) == 1 ? 'N' : 'I';
return;
}
if (sign != 0) {
x = -x;
}
{
int exp;
double frac = frexp(x, &exp);
long num_bits_extract = DBL_MANT_DIG - __count_trailing_zero(frac);
double integer;
decimal int_d, pow2_d;
__two_exp(&pow2_d, exp - num_bits_extract);
frac = modf(ldexp(frac, num_bits_extract), &integer);
__ull2dec(&int_d, (unsigned long long)integer);
__timesdec(d, &int_d, &pow2_d);
d->sign = sign;
}
}
/* 80363124-803632C8 35DA64 01A4+00 0/0 2/2 0/0 .text __num2dec */
void __num2dec(const decform* form, double x, decimal* d) {
short digits = form->digits;
int i;
__num2dec_internal(d, x);
if (d->sig.text[0] > 9) {
return;
}
if (digits > SIGDIGLEN) {
digits = SIGDIGLEN;
}
__rounddec(d, digits);
while (d->sig.length < digits) {
d->sig.text[d->sig.length++] = 0;
}
d->exp -= d->sig.length - 1;
for (i = 0; i < d->sig.length; i++) {
d->sig.text[i] += '0';
}
}
@@ -0,0 +1,6 @@
#ifndef PPC_EABI_SRC_CRITICAL_REGIONSGAMECUBE_H
#define PPC_EABI_SRC_CRITICAL_REGIONSGAMECUBE_H
#endif /* PPC_EABI_SRC_CRITICAL_REGIONSGAMECUBE_H */
@@ -0,0 +1,6 @@
#ifndef PPC_EABI_SRC_MATH_PPC_H
#define PPC_EABI_SRC_MATH_PPC_H
#endif /* PPC_EABI_SRC_MATH_PPC_H */
@@ -0,0 +1,6 @@
#ifndef PPC_EABI_SRC_UART_CONSOLE_IO_GCN_H
#define PPC_EABI_SRC_UART_CONSOLE_IO_GCN_H
#include "ansi_files.h"
#endif /* PPC_EABI_SRC_UART_CONSOLE_IO_GCN_H */

Some files were not shown because too many files have changed in this diff Show More