More floating point...

In particular:
- update math.h, time.h:
  - add: expm1(), pow(), cosh(), sinh(), tanh(), difftime()
  - improve (using assembly): frexp(), ldexp(), modf()
- re-initialize x87 FPU in DOS implementation of system()
- add Mandelbrot set test/demo
- use floating point in VESALFB test/demo
- update tests, binaries and libraries

license.txt

@@ -1,4 +1,4 @@
-Copyright (c) 2012-2015, Alexey Frunze
+Copyright (c) 2012-2016, Alexey Frunze
 All rights reserved.
 
 Redistribution and use in source and binary forms, with or without

v0100/include/math.h

@@ -40,6 +40,10 @@ double exp(double);
 float expf(float);
 double exp2(double);
 float exp2f(float);
+double expm1(double);
+float expm1f(float);
+double pow(double, double);
+float powf(float, float);
 double log(double);
 float logf(float);
 double log2(double);
@@ -62,14 +66,12 @@ double asin(double);
 float asinf(float);
 double fmod(double, double);
 float fmodf(float, float);
-
-/*
-// C89 TBD!!!
-double cosh(double x);
-double sinh(double x);
-double tanh(double x);
-double pow(double x, double y);
-*/
+double cosh(double);
+float coshf(float);
+double sinh(double);
+float sinhf(float);
+double tanh(double);
+float tanhf(float);
 
 /*
 // C99 TBD???

v0100/include/time.h

@@ -61,6 +61,7 @@ struct tm* gmtime(time_t*);
 struct tm* localtime(time_t*);
 time_t mktime(struct tm*);
 char* ctime(time_t*);
+double difftime(time_t, time_t);
 #endif
 char* asctime(struct tm*);
 size_t strftime(char*, size_t, char*, struct tm*);

v0100/srclib/cosh.c

@@ -0,0 +1,21 @@
+/*
+  Copyright (c) 2016, Alexey Frunze
+  2-clause BSD license.
+*/
+
+#ifdef __SMALLER_C_32__
+
+#include <math.h>
+
+float coshf(float x)
+{
+  x = x * 1.44269504f/*log2(e)*/;
+  return exp2f(x - 1) + exp2f(-x - 1);
+}
+
+double cosh(double x)
+{
+  return coshf(x);
+}
+
+#endif

v0100/srclib/difftime.c

@@ -0,0 +1,17 @@
+/*
+  Copyright (c) 2016, Alexey Frunze
+  2-clause BSD license.
+*/
+
+#ifdef __SMALLER_C_32__
+
+#include <time.h>
+
+double difftime(time_t time1, time_t time0)
+{
+  if (time1 >= time0)
+    return (double)((unsigned long)time1 - time0);
+  return -(double)((unsigned long)time0 - time1);
+}
+
+#endif

v0100/srclib/frexp.c

@@ -25,6 +25,8 @@ float frexpf(float value, int* e)
 }
 */
 
+#ifdef UNUSED
+
 float frexpf(float value, int* e)
 {
   union
@@ -61,6 +63,49 @@ float frexpf(float value, int* e)
   return u.f;
 }
 
+#else
+
+#ifdef __HUGE__
+#define xbp "bp"
+#define xbx "bx"
+#else
+#define xbp "ebp"
+#define xbx "ebx"
+#endif
+
+float frexpf(float value, int* e)
+{
+  asm
+  (
+  "mov    eax, ["xbp"+8]\n"
+  "mov    ebx, ["xbp"+12]\n"
+#ifdef __HUGE__
+  "ror    ebx, 4\n"
+  "mov    ds, bx\n"
+  "shr    ebx, 28\n"
+#endif
+  "mov    dword ["xbx"], 0\n" // preset *e=0 for +/-0.0, +/-INF, NAN
+  "mov    ecx, eax\n"
+  "shl    ecx, 1\n" // shift the sign out
+  "jz     .done\n" // done if +/-0.0
+  "cmp    ecx, 0xFF000000\n"
+  "jnc    .done\n" // done if +/-INF or NAN
+  "fld    dword ["xbp"+8]\n"
+  "fxtract\n"
+  "push   dword 0x3F000000\n" // push 0.5f
+  "fld    dword ["xbp"-4]\n"
+  "pop    eax\n"
+  "fmulp\n"                  // scale significand down from [1,2) to [0.5,1)
+  "fstp   dword ["xbp"+8]\n"
+  "mov    eax, ["xbp"+8]\n"
+  "fistp  dword ["xbx"]\n"
+  "inc    dword ["xbx"]\n" // correct exponent
+  ".done:"
+  );
+}
+
+#endif
+
 double frexp(double value, int* e)
 {
   return frexpf(value, e);

v0100/srclib/lcdh.txt

@@ -157,6 +157,7 @@ mktime.c
 asctime.c
 ctime.c
 strftime.c
+difftime.c
 
 frexp.c
 ldexp.c
@@ -170,6 +171,7 @@ sqrt.c
 hypot.c
 exp.c
 exp2.c
+expm1.c
 log.c
 log2.c
 log10.c
@@ -181,6 +183,10 @@ atan.c
 asin.c
 acos.c
 fmod.c
+pow.c
+cosh.c
+sinh.c
+tanh.c
 cvtdf.c
 
 lngfu.c

v0100/srclib/lcdp.txt

@@ -158,6 +158,7 @@ mktime.c
 asctime.c
 ctime.c
 strftime.c
+difftime.c
 
 frexp.c
 ldexp.c
@@ -171,6 +172,7 @@ sqrt.c
 hypot.c
 exp.c
 exp2.c
+expm1.c
 log.c
 log2.c
 log10.c
@@ -182,6 +184,10 @@ atan.c
 asin.c
 acos.c
 fmod.c
+pow.c
+cosh.c
+sinh.c
+tanh.c
 cvtdf.c
 
 lngfu.c

v0100/srclib/lcl.txt

@@ -157,6 +157,7 @@ mktime.c
 asctime.c
 ctime.c
 strftime.c
+difftime.c
 
 frexp.c
 ldexp.c
@@ -170,6 +171,7 @@ sqrt.c
 hypot.c
 exp.c
 exp2.c
+expm1.c
 log.c
 log2.c
 log10.c
@@ -181,6 +183,10 @@ atan.c
 asin.c
 acos.c
 fmod.c
+pow.c
+cosh.c
+sinh.c
+tanh.c
 cvtdf.c
 
 lngfu.c

v0100/srclib/lcw.txt

@@ -159,6 +159,7 @@ mktime.c
 asctime.c
 ctime.c
 strftime.c
+difftime.c
 
 frexp.c
 ldexp.c
@@ -172,6 +173,7 @@ sqrt.c
 hypot.c
 exp.c
 exp2.c
+expm1.c
 log.c
 log2.c
 log10.c
@@ -183,6 +185,10 @@ atan.c
 asin.c
 acos.c
 fmod.c
+pow.c
+cosh.c
+sinh.c
+tanh.c
 cvtdf.c
 
 lngfu.c

v0100/srclib/ldexp.c

@@ -5,8 +5,6 @@
 
 #ifdef __SMALLER_C_32__
 
-#include <limits.h>
-
 /*
 float ldexpf(float value, int e)
 {
@@ -25,6 +23,10 @@ float ldexpf(float value, int e)
 }
 */
 
+#ifdef UNUSED
+
+#include <limits.h>
+
 float ldexpf(float value, int e)
 {
   union
@@ -92,6 +94,29 @@ float ldexpf(float value, int e)
   return u.f;
 }
 
+#else
+
+#ifdef __HUGE__
+#define xbp "bp"
+#else
+#define xbp "ebp"
+#endif
+
+float ldexpf(float value, int e)
+{
+  asm
+  (
+  "fild  dword ["xbp"+12]\n"
+  "fld   dword ["xbp"+8]\n"
+  "fscale\n"
+  "fstp  st1\n"
+  "fstp  dword ["xbp"+8]\n"
+  "mov   eax, ["xbp"+8]"
+  );
+}
+
+#endif
+
 double ldexp(double value, int e)
 {
   return ldexpf(value, e);

v0100/srclib/modf.c

@@ -5,6 +5,8 @@
 
 #ifdef __SMALLER_C_32__
 
+#ifdef UNUSED
+
 float modff(float value, float* iptr)
 {
   union
@@ -63,6 +65,66 @@ float modff(float value, float* iptr)
   return u.f;
 }
 
+#else
+
+#ifdef __HUGE__
+#define xbp "bp"
+#define xsp "sp"
+#define xbx "bx"
+#else
+#define xbp "ebp"
+#define xsp "esp"
+#define xbx "ebx"
+#endif
+
+float modff(float value, float* iptr)
+{
+  asm
+  (
+  "mov    eax, ["xbp"+8]\n"
+  "mov    ebx, ["xbp"+12]\n"
+#ifdef __HUGE__
+  "ror    ebx, 4\n"
+  "mov    ds, bx\n"
+  "shr    ebx, 28\n"
+#endif
+  "mov    ["xbx"], eax\n" // preset *iptr=value for +/-0.0, +/-INF, NAN
+  "mov    ecx, eax\n"
+  "shl    ecx, 1\n" // shift the sign out
+  "jz     .done\n" // done if +/-0.0
+  "cmp    ecx, 0xFF000000\n"
+  "jb     .finite\n"
+  "ja     .done\n" // done if NAN
+  "and    eax, 0x80000000\n" // return +0.0 for +INF, -0.0 for -INF
+  "jmp    .done\n"
+
+  ".finite:\n"
+  "sub    "xsp", 4\n"
+  "fnstcw ["xbp"-2]\n" // save rounding
+  "mov    ax, ["xbp"-2]\n"
+  "mov    ah, 0x0c\n" // rounding towards 0 (AKA truncate)
+  "mov    ["xbp"-4], ax\n"
+  "fld    dword ["xbp"+8]\n"
+  "fld    st0\n"
+  "fldcw  ["xbp"-4]\n"
+  "frndint\n" // trunc(value)
+  "fst    dword ["xbx"]\n" // *iptr = trunc(value)
+  "fldcw  ["xbp"-2]\n" // restore rounding
+  "add    "xsp", 4\n"
+  "fsubp\n" // value - trunc(value)
+  "fstp   dword ["xbp"+8]\n"
+  "mov    eax, ["xbp"+8]\n" // eax = value - trunc(value)
+  // almost done, copy the sign of *iptr into the returned value...
+  "mov    ebx, dword ["xbx"]\n" // ebx = trunc(value)
+  "shl    eax, 1\n"
+  "rcl    ebx, 1\n"
+  "rcr    eax, 1\n"
+  ".done:"
+  );
+}
+
+#endif
+
 double modf(double value, double* iptr)
 {
   return modff(value, iptr);