From 4b972d2a4e81869086185f18be251fff5be69f54 Mon Sep 17 00:00:00 2001 From: Phillip Stephens Date: Wed, 30 Oct 2024 23:54:08 -0700 Subject: [PATCH] Match and link ansi_fp.c --- configure.py | 2 +- include/Kyoto/Math/CVector3f.hpp | 2 +- libc/float.h | 6 + libc/math.h | 13 +- src/Kyoto/Text/CWordBreakTables.cpp | 9 +- src/Runtime/ansi_fp.c | 601 ++++++++++++++++++++++++++++ 6 files changed, 626 insertions(+), 7 deletions(-) create mode 100644 src/Runtime/ansi_fp.c diff --git a/configure.py b/configure.py index 718ee05a..5b65d8a7 100755 --- a/configure.py +++ b/configure.py @@ -1209,7 +1209,7 @@ config.libs = [ Object(Matching, "Runtime/abort_exit.c"), Object(NonMatching, "Runtime/alloc.c"), Object(Matching, "Runtime/ansi_files.c"), - Object(NonMatching, "Runtime/ansi_fp.c"), + Object(Matching, "Runtime/ansi_fp.c"), Object(Matching, "Runtime/arith.c"), Object(Matching, "Runtime/buffer_io.c"), Object(Matching, "Runtime/ctype.c"), diff --git a/include/Kyoto/Math/CVector3f.hpp b/include/Kyoto/Math/CVector3f.hpp index a02ec7cf..62d87cf8 100644 --- a/include/Kyoto/Math/CVector3f.hpp +++ b/include/Kyoto/Math/CVector3f.hpp @@ -75,7 +75,7 @@ public: const float& operator[](EDimZ) const { return mZ; } float& operator[](int i) { return (&mX)[i]; } - float operator[](int i) const { return (&mX)[i]; } + const float operator[](int i) const { return (&mX)[i]; } bool IsNonZero() const { return mX != 0.f || mY != 0.f || mZ != 0.f; } CVector3f DropZ() const { return CVector3f(mX, mY, 0.f); } diff --git a/libc/float.h b/libc/float.h index ee8b96b6..f3a4f133 100644 --- a/libc/float.h +++ b/libc/float.h @@ -5,12 +5,18 @@ extern "C" { #endif +#define FLT_DIG 6 #define FLT_MAX 3.402823466e+38f #define FLT_EPSILON 1.192092896e-07f #define FLT_MIN 1.175494351e-38f +#define DBL_DIG 6 +#define DBL_MIN 5.8774717e-39 +#define DBL_MAX (* (double *) __double_max) #define DBL_EPSILON 1.1920929e-07 +#define DBL_MANT_DIG 53 + #ifdef __cplusplus } #endif diff --git a/libc/math.h b/libc/math.h index 9d8ce973..8d69b046 100644 --- a/libc/math.h +++ b/libc/math.h @@ -56,6 +56,7 @@ extern _INT32 __float_huge[]; extern _INT32 __float_nan[]; extern _INT32 __double_huge[]; extern _INT32 __extended_huge[]; +extern _INT32 __double_max[]; #define HUGE_VAL (*(double*)__double_huge) #define INFINITY (*(float*)__float_huge) @@ -111,6 +112,8 @@ _MATH_INLINE float powf(float __x, float __y) { return pow(__x, __y); } #define __UHI(x) (*(_UINT32*)&x) #endif +#define signbit(x)((int)(__HI(x)&0x80000000)) + #define FP_NAN 1 #define FP_INFINITE 2 #define FP_ZERO 3 @@ -204,14 +207,16 @@ float sqrtf(float x); double sqrt(double x); #endif -static inline float ldexpf(float x, int exp) { return (float)ldexp((double)x, exp); } -static inline double scalbn(double x, int n) { return ldexp(x, n); } -static inline float scalbnf(float x, int n) { return (float)ldexpf(x, n); } - #ifdef __MWERKS__ #pragma cplusplus reset #endif +static inline float ldexpf(float x, int exp) { return (float)ldexp((double)x, exp); } +double frexp(double, int *exp); +static inline double scalbn(double x, int n) { return ldexp(x, n); } +static inline float scalbnf(float x, int n) { return (float)ldexpf(x, n); } +double nextafter(double, double); + #ifdef __cplusplus } #endif diff --git a/src/Kyoto/Text/CWordBreakTables.cpp b/src/Kyoto/Text/CWordBreakTables.cpp index fc82f5b6..f6b698ea 100644 --- a/src/Kyoto/Text/CWordBreakTables.cpp +++ b/src/Kyoto/Text/CWordBreakTables.cpp @@ -5,6 +5,11 @@ struct CCharacterIdentifier { wchar_t chr; uint rank; + + struct Compare { + bool operator()(const CCharacterIdentifier& ident, wchar_t chr) { return ident.chr == chr; } + bool operator()(wchar_t chr, const CCharacterIdentifier& ident) { return chr == ident.chr; } + }; }; const CCharacterIdentifier gCantBeginChars[63] = { @@ -21,7 +26,7 @@ const CCharacterIdentifier gCantBeginChars[63] = { const CCharacterIdentifier gCantEndChars[89] = { {L'#', 2}, {L'$', 2}, {L'(', 1}, {L'@', 2}, {L'B', 4}, {L'C', 4}, {L'D', 4}, - {L'E', 4}, {L'F', 4}, {L'G', 4}, {L'J', 4}, {L'K', 4}, {L'L', 4}, {L'M', 4}, + {L'F', 4}, {L'G', 4}, {L'H', 4}, {L'J', 4}, {L'K', 4}, {L'L', 4}, {L'M', 4}, {L'N', 4}, {L'P', 4}, {L'Q', 4}, {L'R', 4}, {L'S', 4}, {L'T', 4}, {L'V', 4}, {L'W', 4}, {L'X', 4}, {L'Y', 4}, {L'Z', 4}, {L'b', 4}, {L'c', 4}, {L'd', 4}, {L'f', 4}, {L'g', 4}, {L'h', 4}, {L'j', 4}, {L'k', 4}, {L'l', 4}, {L'm', 4}, @@ -36,4 +41,6 @@ const CCharacterIdentifier gCantEndChars[89] = { }; int CWordBreakTables::GetBeginRank(wchar_t chr) { + //rstl::binary_find(&gCantBeginChars[0], &gCantBeginChars[62], chr, CCharacterIdentifier::Compare()); + return -1; } diff --git a/src/Runtime/ansi_fp.c b/src/Runtime/ansi_fp.c new file mode 100644 index 00000000..59aff744 --- /dev/null +++ b/src/Runtime/ansi_fp.c @@ -0,0 +1,601 @@ +#include +#include +#include +#include + +static int __count_trailing_zerol(unsigned long x) { + int result = 0; + int bits_not_checked = sizeof(unsigned long) * 8; + 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) * 8 + __count_trailing_zerol(l[0] | 0x00100000)); +} + +static int __must_round(const decimal* d, int pos) { + unsigned char const* i = d->sig.text + pos; + 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[pos - 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)++; + 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) + return; + { + int r = __must_round(d, digits); + d->sig.length = (unsigned char)digits; + if (r < 0) + return; + } + __dorounddecup(d, digits); +} + +static void __ull2dec(decimal* result, unsigned long long integer) { + result->sgn = 0; + if (integer == 0) { + result->exp = 0; + result->sig.length = 1; + result->sig.text[0] = 0; + return; + } + if (integer < 0) { + integer = -integer; + result->sgn = 1; + } + result->sig.length = 0; + for (; integer != 0; integer /= 10) + result->sig.text[result->sig.length++] = (unsigned char)(integer % 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 = (short)(result->sig.length - 1); +} + +static void __timesdec(decimal* result, const decimal* x, const decimal* y) { + unsigned long accumulator = 0; + unsigned char mantissa[2 * SIGDIGLEN]; + int i = x->sig.length + y->sig.length - 1; + unsigned char* ip = mantissa + i + 1; + unsigned char* ep = ip; + result->sgn = 0; + for (; i > 0; --i) { + int k = y->sig.length - 1; + int j = i - k - 1; + int l; + int t; + unsigned char const *jp, *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 += (unsigned long)(*jp) * *kp; + *--ip = (unsigned char)(accumulator % 10); + accumulator /= 10; + } + result->exp = (short)(x->exp + y->exp); + if (accumulator > 0) { + *--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); + } +} + +static void __str2dec(decimal* d, const char* s, short exp) { + int i; + d->exp = exp; + d->sgn = 0; + for (i = 0; i < SIGDIGLEN && *s != 0;) + d->sig.text[i++] = (unsigned char)(*s++ - '0'); + d->sig.length = (unsigned char)i; + if (*s != 0) { + if (*s < 5) + return; + { + 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 void __two_exp(decimal* result, short 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; + __two_exp(&x2, (short)(exp / 2)); + __timesdec(result, &x2, &x2); + } + if (exp & 1) { + decimal temp = *result; + if (exp > 0) { + decimal two; + __str2dec(&two, "2", 0); + __timesdec(result, &temp, &two); + } else { + decimal one_half; + __str2dec(&one_half, "5", -1); + __timesdec(result, &temp, &one_half); + } + } +} + +static 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; +} + +static 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; + 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; +} + +static 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* je = y->sig.text + y->sig.length; + for (j = jn + 1; j < je; ++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); +} + +static void __num2dec_internal(decimal* d, double x) { + char sgn = signbit(x) != 0; + if (x == 0.0) { + d->sgn = sgn; + d->exp = 0; + d->sig.length = 1; + d->sig.text[0] = 0; + return; + } + if (!isfinite(x)) { + d->sgn = sgn; + d->exp = 0; + d->sig.length = 1; + d->sig.text[0] = (unsigned char)(isnan(x) ? 'N' : 'I'); + return; + } + + if (sgn) + x = -x; + { + int exp; + double frac = frexp(x, &exp); + short num_bits_extract = (short)(DBL_MANT_DIG - __count_trailing_zero(frac)); + double integer; + decimal int_d, pow2_d; + __two_exp(&pow2_d, (short)(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->sgn = sgn; + } +} + +void __num2dec(const decform* f, double x, decimal* d) { + short digits = f->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'; +} + +double __dec2num(const decimal* d) { + if (d->sig.length <= 0) + return copysign(0.0, d->sgn == 0 ? 1.0 : -1.0); + switch (d->sig.text[0]) { + case '0': + return copysign(0.0, d->sgn == 0 ? 1.0 : -1.0); + case 'I': + return copysign((double)INFINITY, d->sgn == 0 ? 1.0 : -1.0); + case 'N': { + double result; + unsigned long long* ll = (unsigned long long*)&result; + *ll = 0x7FF0000000000000; + if (d->sgn) + *ll |= 0x8000000000000000; + if (d->sig.length == 1) + *ll |= 0x0008000000000000; + else { + unsigned char* p = (unsigned char*)&result + 1; + int placed_non_zero = 0; + int low = 1; + int i; + int e = d->sig.length; + if (e > 14) + e = 14; + for (i = 1; i < e; ++i) { + unsigned char c = d->sig.text[i]; + if (isdigit(c)) + c -= '0'; + else + c = (unsigned char)(tolower(c) - 'a' + 10); + if (c != 0) + placed_non_zero = 1; + if (low) + *p++ |= c; + else + *p = (unsigned char)(c << 4); + low = !low; + } + if (!placed_non_zero) + *ll |= 0x0008000000000000; + } + return result; + } + } + { + static double pow_10[8] = {1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8}; + decimal dec = *d; + unsigned char* i = dec.sig.text; + unsigned char* e = i + dec.sig.length; + double first_guess; + int exponent; + for (; i < e; ++i) + *i -= '0'; + dec.exp += dec.sig.length - 1; + exponent = dec.exp; + i = dec.sig.text; + first_guess = *i++; + while (i < e) { + unsigned long ival = 0; + int j; + double temp1, temp2; + int ndig = (e - i) % 8; + if (ndig == 0) + ndig = 8; + for (j = 0; j < ndig; ++j, ++i) + ival = ival * 10 + *i; + temp1 = first_guess * pow_10[ndig - 1]; + temp2 = temp1 + ival; + if (ival != 0 && temp1 == temp2) + break; + first_guess = temp2; + exponent -= ndig; + } + if (exponent < 0) + first_guess /= pow(5.0, -exponent); + else + first_guess *= pow(5.0, exponent); + first_guess = ldexp(first_guess, exponent); + if (isinf(first_guess)) { + decimal max; + __str2dec(&max, "179769313486231580793729011405303420", 308); + if (__less_dec(&max, &dec)) + goto done; + first_guess = DBL_MAX; + } + { + decimal feedback1; + __num2dec_internal(&feedback1, first_guess); + if (__equals_dec(&feedback1, &dec)) + goto done; + if (__less_dec(&feedback1, &dec)) { + + decimal feedback2, difflow, diffhigh; + double next_guess = nextafter(first_guess, (double)INFINITY); + if (isinf(next_guess)) { + first_guess = next_guess; + goto done; + } + __num2dec_internal(&feedback2, next_guess); + while (__less_dec(&feedback2, &dec)) { + feedback1 = feedback2; + first_guess = next_guess; + next_guess = nextafter(next_guess, (double)INFINITY); + if (isinf(next_guess)) { + first_guess = next_guess; + goto done; + } + __num2dec_internal(&feedback2, next_guess); + } + __minus_dec(&difflow, &dec, &feedback1); + __minus_dec(&diffhigh, &feedback2, &dec); + if (__equals_dec(&difflow, &diffhigh)) { + if (*(unsigned long long*)&first_guess & 1) + first_guess = next_guess; + } else if (!__less_dec(&difflow, &diffhigh)) + first_guess = next_guess; + } else { + decimal feedback2, difflow, diffhigh; + double next_guess = nextafter(first_guess, (double)(-INFINITY)); + __num2dec_internal(&feedback2, next_guess); + while (__less_dec(&dec, &feedback2)) { + feedback1 = feedback2; + first_guess = next_guess; + next_guess = nextafter(next_guess, (double)(-INFINITY)); + __num2dec_internal(&feedback2, next_guess); + } + __minus_dec(&difflow, &dec, &feedback2); + __minus_dec(&diffhigh, &feedback1, &dec); + if (__equals_dec(&difflow, &diffhigh)) { + if (*(unsigned long long*)&first_guess & 1) + first_guess = next_guess; + } else if (__less_dec(&difflow, &diffhigh)) + first_guess = next_guess; + } + } + done: + if (dec.sgn) + first_guess = -first_guess; + return first_guess; + } +}