Match and link ansi_fp.c

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"),

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); }

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

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

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;
 }

src/Runtime/ansi_fp.c

@@ -0,0 +1,601 @@
+#include <ansi_fp.h>
+#include <ctype.h>
+#include <float.h>
+#include <math.h>
+
+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;
+  }
+}