diff --git a/bake.cmd b/bake.cmd index 58e8bb4..d757f81 100644 --- a/bake.cmd +++ b/bake.cmd @@ -51,4 +51,4 @@ del build\*.obj :skip_crt_compilation echo Compiling test.. -clang test\test4.c ciabatta.lib -std=c11 -lkernel32 -luser32 -lshell32 -nostdlib %CIABATTA_OPTIONS% +clang test\test5.c ciabatta.lib -std=c11 -lkernel32 -luser32 -lshell32 -nostdlib %CIABATTA_OPTIONS% diff --git a/code/math/ieee754.h b/code/math/ieee754.h index a013427..eeef751 100644 --- a/code/math/ieee754.h +++ b/code/math/ieee754.h @@ -9,21 +9,26 @@ // e.g. f for floats #define f_nbits (1+f_ebits+f_mbits) -#define f_emax ((1ULL << (f_ebits-1)) - 1) -#define f_emin (1 - f_emax) +#define f_emax ((1ULL << (f_ebits-1)) - 1) +#define f_emin (1 - f_emax) +#define f_ebias f_emax // Extracting fields from the float -#define f_eoffs (f_mbits) -#define f_soffs (f_mbits+f_ebits) -#define f_emask ((1ULL << f_ebits) - 1) -#define f_mmask ((1ULL << f_mbits) - 1) -#define f_smask 1ULL +#define f_eoffs (f_mbits) +#define f_soffs (f_mbits+f_ebits) +#define f_emask ((1ULL << f_ebits) - 1) +#define f_mmask ((1ULL << f_mbits) - 1) +#define f_smask 1ULL #define f_eval(b) (((b) >> f_eoffs) & f_emask) #define f_sval(b) (((b) >> f_soffs) & f_smask) #define f_mval(b) (((b) >> 0) & f_mmask) #define f_abs(b) ((b) & ~(f_smask << f_soffs)) +#define f_exp(b) (f_eval(b) - f_ebias) -#define b_cons(s,e,m) ((s << f_soffs) | (e << f_eoffs) | (m)) +#define f_qexp(b) (f_eval(b) - f_ebias - f_mbits) +#define f_qman(b) (((b) & f_mmask) | (f_mmask+1)) + +#define b_cons(s,e,m) (((itype)s << f_soffs) | ((itype)e << f_eoffs) | (itype)(m)) // Converting float to integer bits static inline itype suffix(f_bits)(ftype f) { @@ -35,8 +40,8 @@ static inline itype suffix(f_bits)(ftype f) { return u.b; } -static inline itype suffix(f_frombits)(itype b) { - union _u { +static inline ftype suffix(f_frombits)(itype b) { + union _u { ftype f; itype b; } u; @@ -49,7 +54,7 @@ int suffix(_fpclassify)(ftype f) { itype bits = suffix(f_bits)(f); itype exp = f_eval(bits); itype man = f_mval(bits); - if(exp == f_emax) { + if(exp == f_emask) { if(man == 0) return FP_INFINITE; else return FP_NAN; } diff --git a/code/math/math.c b/code/math/math.c index fa1fbe2..5b85c90 100644 --- a/code/math/math.c +++ b/code/math/math.c @@ -16,6 +16,7 @@ enum Ordering { #define f_mbits 23 #define suffix(n) n ## f #include "ieee754.h" +#include "pow.h" #undef suffix #undef f_mbits #undef f_ebits @@ -28,6 +29,7 @@ enum Ordering { #define f_mbits 52 #define suffix(n) n #include "ieee754.h" +#include "pow.h" #undef suffix #undef f_mbits #undef f_ebits @@ -42,7 +44,8 @@ _Static_assert(sizeof(long double) == sizeof(double), #define f_ebits 11 #define f_mbits 52 #define suffix(n) n ## l -//#include "ieee754.h" +#include "ieee754.h" +#include "pow.h" #undef suffix #undef f_mbits #undef f_ebits diff --git a/code/math/pow.h b/code/math/pow.h index e69de29..336ea3d 100644 --- a/code/math/pow.h +++ b/code/math/pow.h @@ -0,0 +1,95 @@ + +#include + +#if !defined(_isqrt_defined) +#define _isqrt_defined +static uint64_t _isqrt(uint64_t num, uint64_t *remp) { + // To find a square root of a number + // We get rid of zero + if(num == 0) { + *remp = 0; + return 0; + } + // Then, starting from the bottom, split num into 2-digit pairs + // and find the top-most non-zero pair + uint64_t i = 0; + while(i != (sizeof(uint64_t)*8) && (num >> i) != 0) { + i += 2; + } + // Then we start taking guesses such that at each step + // sqrt^2 <= number made of consequent pairs of exausted integers + uint64_t sqrt = 0; + uint64_t rem = 0; + // Repeat until remainder is equal to zero: + do { + i -= 2; + // Bring the next two digits of the number to our remainder + rem = (rem << 2) | ((num >> i) & 0x3); + // Find d such that d(2sqrt+d) <= rem + // Since d could be either 0 or 1 we simply check 1, otherwise its 0 + uint64_t d = 1; + uint64_t t = ((sqrt<<2)|1); + if(t <= rem) { + rem -= t; + } + else { + d = 0; + } + // Append the digit to sqrt from the right + sqrt = (sqrt<<1)|d; + } while(i != 0); + + *remp = rem; + return sqrt; +} +#endif + +// For all it's worth this shit is simply equivalent to +// _isqrt((uint64)x) +// I hate porgaming. +ftype suffix(sqrt)(ftype x) { + if(x < 0) { + #if math_errhandling & MATH_ERRNO + errno = EDOM; + #endif + return NAN; + } + if(x == 0 || isinf(x)) { + return x; + } + if(isnan(x)) { + return NAN; + } + itype bits = suffix(f_bits)(x); + itype exp = f_qexp(bits); + itype man = f_qman(bits); + // Get lots of precision by shifting man right by max bits + // and subtracting this from the exponent + itype bit = 0; // highest set-bit of man + while((man >> (bit+1)) != 0) ++bit; + itype prec_shift_n = f_nbits - bit - 3; + man <<= prec_shift_n; + exp -= prec_shift_n; + // Now do the sqrt of 2^exp * man + // If exp is odd then 2^{2k+1}*sqrt(man) = 2^{2k}*sqrt{2*man} + if((2 + (exp % 2)) % 2 != 0) { + man <<= 1; + } + // Take exp sqrt + exp >>= 1; + // Take sqrt of mantissa + uint64_t rem; + man = (itype)_isqrt(man, &rem); + // Now sqrt(x) = 2^exp * man + // we need to normalize this shit + bit = 0; // highest set-bit of man + while((man >> (bit+1)) != 0) ++bit; + exp += bit; + man <<= f_nbits-bit; + exp += f_ebias; + man >>= f_nbits-f_mbits; + man &= f_mmask; + // Cons it back + bits = b_cons(0, exp, man); + return suffix(f_frombits)(bits); +} diff --git a/code/printf.h b/code/printf.h index 2aaf80c..af31298 100644 --- a/code/printf.h +++ b/code/printf.h @@ -1,3 +1,7 @@ + +//TODO: verify printf("%d", 0). From the code it looked like it would print +// an empty string. + // NOTE: this file doesn't exist in a vacuum, it's a template for generating // the formatted print, you should define FMT_CHAR_TYPE before including it inline static int FMT_FUNC_NAME (void *ctx, OutputFunc out, const FMT_CHAR_TYPE *fmt, va_list args) { @@ -83,6 +87,9 @@ inline static int FMT_FUNC_NAME (void *ctx, OutputFunc out, const FMT_CHAR_TYPE } FMT_CHAR_TYPE ch = *fmt++; + const char* characters = "0123456789abcdef"; + if (ch == 'X') characters = "0123456789ABCDEF"; + switch (ch) { case 'c': { const char chr = va_arg(args, int); @@ -90,6 +97,45 @@ inline static int FMT_FUNC_NAME (void *ctx, OutputFunc out, const FMT_CHAR_TYPE full_length ++; break; } + case 'f': + case 'L': { + double d = va_arg(args, double); + + if(isinf(d)) { + out(ctx, sizeof"inf"-1, "inf"); + break; + } + else if(isnan(d)) { + out(ctx, sizeof"nan"-1, "nan"); + break; + } + + if(d < 0) { // TODO: negative zero + out(ctx, 1, "-"); + d = -d; + } + + uint64_t w = (uint64_t)d; + d -= w; + FMT_CHAR_TYPE buffer[20]; + size_t len = sizeof(buffer); + do { + buffer[--len] = characters[w % 10]; + w /= 10; + } while(w != 0); + out(ctx, sizeof(buffer) - (len * sizeof(FMT_CHAR_TYPE)), buffer + len); + + char dot = '.'; + out(ctx, 1, &dot); + + for(int i = 0; i != 6; ++i) { + d *= 10; + int dv = (int)d; + d -= dv; + char digit = characters[dv]; + out(ctx, 1, &digit); + } + } break; case 's': { const FMT_CHAR_TYPE *str = va_arg(args, FMT_CHAR_TYPE*); size_t len = FMT_STRLEN_S(str, precision ? precision : SIZE_MAX); @@ -114,9 +160,6 @@ inline static int FMT_FUNC_NAME (void *ctx, OutputFunc out, const FMT_CHAR_TYPE default: base = 10; break; } - const char* characters = "0123456789abcdef"; - if (ch == 'X') characters = "0123456789ABCDEF"; - uintmax_t i; if (ch == 'd' || ch == 'i') { intmax_t num = 0; diff --git a/code/stdio.c b/code/stdio.c index d97e1e7..b7c782e 100644 --- a/code/stdio.c +++ b/code/stdio.c @@ -2,6 +2,7 @@ #include #include #include +#include #include <_os.h> diff --git a/test/test5.c b/test/test5.c new file mode 100644 index 0000000..856970f --- /dev/null +++ b/test/test5.c @@ -0,0 +1,22 @@ + +#include +#include + +void test_sqrt(float f) { + float s = sqrtf(f); + printf("sqrt of %f is %f\n", f, s); +} + +int main() { + test_sqrt(0.0f); + test_sqrt(1.0f); + test_sqrt(2.0f); + test_sqrt(3.0f); + test_sqrt(4.0f); + test_sqrt(7.0f); + test_sqrt(9.0f); + test_sqrt(16.0f); + test_sqrt(256.0f); + test_sqrt(257.0f); + return 0; +}