up

.gitmodules

@@ -1,9 +1,6 @@
 [submodule "3rd/lz4/lz4"]
 	path = 3rd/lz4/lz4
 	url = https://github.com/lz4/lz4
-[submodule "3rd/dmath/dmath"]
-	path = 3rd/dmath/dmath
-	url = https://github.com/pocketpy/dmath
 [submodule "3rd/periphery/c-periphery"]
 	path = 3rd/periphery/c-periphery
 	url = https://github.com/vsergeev/c-periphery.git

3rd/dmath/dmath

@@ -1 +0,0 @@
-Subproject commit cc7cdccae657b0f81a1e234eb9d4ed4c458bda0f

CMakeLists.txt

@@ -62,7 +62,6 @@ else()
 endif()
 
 if(PK_ENABLE_DETERMINISM)
-    add_subdirectory(3rd/dmath/dmath)
     add_definitions(-DPK_ENABLE_DETERMINISM=1)
 else()
     add_definitions(-DPK_ENABLE_DETERMINISM=0)
@@ -149,24 +148,12 @@ target_include_directories(
     ${CMAKE_CURRENT_LIST_DIR}/include
 )
 
-if(PK_ENABLE_DETERMINISM)
-    target_link_libraries(${PROJECT_NAME} dmath)
-    if(MSVC)
-        target_link_options(${PROJECT_NAME} PRIVATE /FORCE:MULTIPLE)
-    endif()
-endif()
-
 if(PK_ENABLE_THREADS)
     find_package(Threads REQUIRED)
     target_link_libraries(${PROJECT_NAME} Threads::Threads)
 endif()
 
 if(UNIX AND NOT APPLE)
-    if(NOT PK_ENABLE_DETERMINISM)
-        # use platform libm
-        target_link_libraries(${PROJECT_NAME} m)
-    endif()
-
     if(PK_ENABLE_OS)
         target_link_libraries(${PROJECT_NAME} dl)
     endif()

include/pocketpy/common/algorithm.h

@@ -1,6 +1,7 @@
 #pragma once
 
 #include <stdbool.h>
+#include "pocketpy/common/utils.h"
 
 #define c11__less(a, b) ((a) < (b))
 

include/pocketpy/config.h

@@ -52,12 +52,6 @@
 // This is the maximum character length of a module path
 #define PK_MAX_MODULE_PATH_LEN      63
 
-// This is some math constants
-#define PK_M_PI                     3.1415926535897932384
-#define PK_M_E                      2.7182818284590452354
-#define PK_M_DEG2RAD                0.017453292519943295
-#define PK_M_RAD2DEG                57.29577951308232
-
 // Hash table load factor (smaller ones mean less collision but more memory)
 // For class instance
 #define PK_INST_ATTR_LOAD_FACTOR    0.67f

scripts/gen_sin_table.py

@@ -0,0 +1,25 @@
+import math
+
+def generate_sin_table(num_points = 4096 + 1):
+    y_step = 1.0 / (num_points - 1)
+    sin_table = []
+    for i in range(num_points):
+        y = i * y_step
+        x = math.asin(y)
+        x = max(0, min(x, math.pi / 2))
+        sin_table.append(x)
+    return sin_table
+
+if __name__ == "__main__":
+    sin_table = generate_sin_table()
+    print(len(sin_table))
+    # print(sin_table)
+
+    filepath = 'include/pocketpy/common/_sin_table.h'
+    with open(filepath, 'w') as f:
+        f.write('static const double _sin_table[] = {')
+        for i, val in enumerate(sin_table):
+            if i % 8 == 0:
+                f.write('\n    ')
+            f.write(str(val) + ', ')
+        f.write('\n};\n')

src/bindings/py_number.c

@@ -1,8 +1,8 @@
 #include "pocketpy/interpreter/vm.h"
 #include "pocketpy/common/sstream.h"
+#include "pocketpy/common/dmath.h"
 #include "pocketpy/pocketpy.h"
 
-#include <math.h>
 
 static bool try_castfloat(py_Ref self, double* out) {
     switch(self->type) {
@@ -100,7 +100,7 @@ static bool number__pow__(int argc, py_Ref argv) {
             if(lhs == 0) {
                 return ZeroDivisionError("0.0 cannot be raised to a negative power");
             } else {
-                py_newfloat(py_retval(), pow(lhs, rhs));
+                py_newfloat(py_retval(), dmath_pow(lhs, rhs));
             }
         } else {
             // rhs >= 0
@@ -117,7 +117,7 @@ static bool number__pow__(int argc, py_Ref argv) {
         py_f64 lhs, rhs;
         if(!py_castfloat(&argv[0], &lhs)) return false;
         if(try_castfloat(&argv[1], &rhs)) {
-            py_newfloat(py_retval(), pow(lhs, rhs));
+            py_newfloat(py_retval(), dmath_pow(lhs, rhs));
         } else {
             py_newnotimplemented(py_retval());
         }
@@ -153,7 +153,7 @@ static py_i64 cpy11__fast_mod(py_i64 a, py_i64 b) {
 static void cpy11__float_div_mod(double vx, double wx, double *floordiv, double *mod)
 {
     double div;
-    *mod = fmod(vx, wx);
+    *mod = dmath_fmod(vx, wx);
     /* fmod is typically exact, so vx-mod is *mathematically* an
        exact multiple of wx.  But this is fp arithmetic, and fp
        vx - mod is an approximation; the result is that div may
@@ -172,18 +172,18 @@ static void cpy11__float_div_mod(double vx, double wx, double *floordiv, double
         /* the remainder is zero, and in the presence of signed zeroes
            fmod returns different results across platforms; ensure
            it has the same sign as the denominator. */
-        *mod = copysign(0.0, wx);
+        *mod = dmath_copysign(0.0, wx);
     }
     /* snap quotient to nearest integral value */
     if (div) {
-        *floordiv = floor(div);
+        *floordiv = dmath_floor(div);
         if (div - *floordiv > 0.5) {
             *floordiv += 1.0;
         }
     }
     else {
         /* div is zero - get the same sign as the true quotient */
-        *floordiv = copysign(0.0, vx / wx); /* zero w/ sign of vx/wx */
+        *floordiv = dmath_copysign(0.0, vx / wx); /* zero w/ sign of vx/wx */
     }
 }
 
@@ -471,11 +471,11 @@ static bool float__new__(int argc, py_Ref argv) {
             c11_sv sv = py_tosv(py_arg(1));
 
             if(c11__sveq2(sv, "inf")) {
-                py_newfloat(py_retval(), INFINITY);
+                py_newfloat(py_retval(), DMATH_INFINITY);
                 return true;
             }
             if(c11__sveq2(sv, "-inf")) {
-                py_newfloat(py_retval(), -INFINITY);
+                py_newfloat(py_retval(), -DMATH_INFINITY);
                 return true;
             }
 

src/common/sstream.c

@@ -1,13 +1,13 @@
 #include "pocketpy/common/sstream.h"
 #include "pocketpy/common/str.h"
 #include "pocketpy/common/utils.h"
+#include "pocketpy/common/dmath.h"
 #include "pocketpy/pocketpy.h"
 
 #include <stdarg.h>
 #include <assert.h>
 #include <stdio.h>
 #include <ctype.h>
-#include <math.h>
 
 void c11_sbuf__ctor(c11_sbuf* self) {
     c11_vector__ctor(&self->data, sizeof(char));
@@ -42,11 +42,11 @@ void c11_sbuf__write_i64(c11_sbuf* self, int64_t val) {
 }
 
 void c11_sbuf__write_f64(c11_sbuf* self, double val, int precision) {
-    if(isinf(val)) {
+    if(dmath_isinf(val)) {
         c11_sbuf__write_cstr(self, val > 0 ? "inf" : "-inf");
         return;
     }
-    if(isnan(val)) {
+    if(dmath_isnan(val)) {
         c11_sbuf__write_cstr(self, "nan");
         return;
     }

src/modules/builtins.c

@@ -3,13 +3,12 @@
 #include "pocketpy/objects/codeobject.h"
 #include "pocketpy/pocketpy.h"
 #include "pocketpy/common/utils.h"
+#include "pocketpy/common/dmath.h"
 #include "pocketpy/objects/object.h"
 #include "pocketpy/common/sstream.h"
 #include "pocketpy/interpreter/vm.h"
 #include "pocketpy/common/_generated.h"
 
-#include <math.h>
-
 
 static bool builtins_exit(int argc, py_Ref argv) {
     int code = 0;
@@ -173,7 +172,9 @@ static bool builtins_round(int argc, py_Ref argv) {
             py_newint(py_retval(), (py_i64)(x + offset));
             return true;
         }
-        py_f64 factor = pow(10, ndigits);
+        // py_f64 factor = dmath_exp10(ndigits);
+        py_f64 factor = 1.0;
+        for(int i = 0; i < ndigits; i++) factor *= 10.0;
         py_newfloat(py_retval(), (py_i64)(x * factor + offset) / factor);
         return true;
     }

src/modules/colorcvt.c

@@ -1,40 +1,40 @@
 #include "pocketpy/pocketpy.h"
 
 #include "pocketpy/common/utils.h"
+#include "pocketpy/common/dmath.h"
 #include "pocketpy/objects/object.h"
 #include "pocketpy/common/sstream.h"
 #include "pocketpy/interpreter/vm.h"
-#include <math.h>
 
 // https://bottosson.github.io/posts/gamutclipping/#oklab-to-linear-srgb-conversion
 
 // clang-format off
 static c11_vec3 linear_srgb_to_oklab(c11_vec3 c)
 {
-	float l = 0.4122214708f * c.x + 0.5363325363f * c.y + 0.0514459929f * c.z;
+	double l = 0.4122214708 * c.x + 0.5363325363 * c.y + 0.0514459929 * c.z;
-	float m = 0.2119034982f * c.x + 0.6806995451f * c.y + 0.1073969566f * c.z;
+	double m = 0.2119034982 * c.x + 0.6806995451 * c.y + 0.1073969566 * c.z;
-	float s = 0.0883024619f * c.x + 0.2817188376f * c.y + 0.6299787005f * c.z;
+	double s = 0.0883024619 * c.x + 0.2817188376 * c.y + 0.6299787005 * c.z;
 
-	float l_ = cbrtf(l);
+	double l_ = dmath_cbrt(l);
-	float m_ = cbrtf(m);
+	double m_ = dmath_cbrt(m);
-	float s_ = cbrtf(s);
+	double s_ = dmath_cbrt(s);
 
 	return (c11_vec3){{
-		0.2104542553f * l_ + 0.7936177850f * m_ - 0.0040720468f * s_,
+		0.2104542553 * l_ + 0.7936177850 * m_ - 0.0040720468 * s_,
-		1.9779984951f * l_ - 2.4285922050f * m_ + 0.4505937099f * s_,
+		1.9779984951 * l_ - 2.4285922050 * m_ + 0.4505937099 * s_,
-		0.0259040371f * l_ + 0.7827717662f * m_ - 0.8086757660f * s_,
+		0.0259040371 * l_ + 0.7827717662 * m_ - 0.8086757660 * s_,
 	}};
 }
 
 static c11_vec3 oklab_to_linear_srgb(c11_vec3 c)
 {
-    float l_ = c.x + 0.3963377774f * c.y + 0.2158037573f * c.z;
+    double l_ = c.x + 0.3963377774 * c.y + 0.2158037573 * c.z;
-    float m_ = c.x - 0.1055613458f * c.y - 0.0638541728f * c.z;
+    double m_ = c.x - 0.1055613458 * c.y - 0.0638541728 * c.z;
-    float s_ = c.x - 0.0894841775f * c.y - 1.2914855480f * c.z;
+    double s_ = c.x - 0.0894841775 * c.y - 1.2914855480 * c.z;
 
-    float l = l_ * l_ * l_;
+    double l = l_ * l_ * l_;
-    float m = m_ * m_ * m_;
+    double m = m_ * m_ * m_;
-    float s = s_ * s_ * s_;
+    double s = s_ * s_ * s_;
 
     return (c11_vec3){{
         +4.0767416621f * l - 3.3077115913f * m + 0.2309699292f * s,
@@ -45,12 +45,12 @@ static c11_vec3 oklab_to_linear_srgb(c11_vec3 c)
 
 // clang-format on
 
-static float _gamma_correct_inv(float x) {
+static double _gamma_correct_inv(double x) {
-    return (x <= 0.04045f) ? (x / 12.92f) : powf((x + 0.055f) / 1.055f, 2.4f);
+    return (x <= 0.04045) ? (x / 12.92) : dmath_pow((x + 0.055) / 1.055, 2.4);
 }
 
-static float _gamma_correct(float x) {
+static double _gamma_correct(double x) {
-    return (x <= 0.0031308f) ? (12.92f * x) : (1.055f * powf(x, 1.0f / 2.4f) - 0.055f);
+    return (x <= 0.0031308) ? (12.92 * x) : (1.055 * dmath_pow(x, 1.0 / 2.4) - 0.055);
 }
 
 static c11_vec3 srgb_to_linear_srgb(c11_vec3 c) {
@@ -70,17 +70,17 @@ static c11_vec3 linear_srgb_to_srgb(c11_vec3 c) {
 static c11_vec3 _oklab_to_oklch(c11_vec3 c) {
     c11_vec3 res;
     res.x = c.x;
-    res.y = sqrtf(c.y * c.y + c.z * c.z);
+    res.y = dmath_sqrt(c.y * c.y + c.z * c.z);
-    res.z = fmodf(atan2f(c.z, c.y), 2 * (float)PK_M_PI);
+    res.z = dmath_fmod(dmath_atan2(c.z, c.y), 2 * DMATH_PI);
-    res.z = res.z * PK_M_RAD2DEG;
+    res.z = res.z * DMATH_RAD2DEG;
     return res;
 }
 
 static c11_vec3 _oklch_to_oklab(c11_vec3 c) {
     c11_vec3 res;
     res.x = c.x;
-    res.y = c.y * cosf(c.z * PK_M_DEG2RAD);
+    res.y = c.y * dmath_cos(c.z * DMATH_DEG2RAD);
-    res.z = c.y * sinf(c.z * PK_M_DEG2RAD);
+    res.z = c.y * dmath_sin(c.z * DMATH_DEG2RAD);
     return res;
 }
 
@@ -105,9 +105,9 @@ static c11_vec3 oklch_to_linear_srgb(c11_vec3 c) {
     // fall back to RGB clamping
     candidate = oklab_to_linear_srgb(_oklch_to_oklab(clamped));
     if(!_is_valid_srgb(candidate)) {
-        candidate.x = fmaxf(0.0f, fminf(1.0f, candidate.x));
+        candidate.x = dmath_fmax(0.0, dmath_fmin(1.0, candidate.x));
-        candidate.y = fmaxf(0.0f, fminf(1.0f, candidate.y));
+        candidate.y = dmath_fmax(0.0, dmath_fmin(1.0, candidate.y));
-        candidate.z = fmaxf(0.0f, fminf(1.0f, candidate.z));
+        candidate.z = dmath_fmax(0.0, dmath_fmin(1.0, candidate.z));
         return candidate;
     }
 
@@ -116,7 +116,7 @@ static c11_vec3 oklch_to_linear_srgb(c11_vec3 c) {
     float start = 0.0f;
     float end = c.y;
     float range[2] = {0.0f, 0.4f};
-    float resolution = (range[1] - range[0]) / powf(2, 13);
+    float resolution = (range[1] - range[0]) / dmath_pow(2, 13);
     float _last_good_c = clamped.y;
 
     while(end - start > resolution) {
@@ -142,8 +142,8 @@ static c11_vec3 srgb_to_hsv(c11_vec3 c) {
     float g = c.y;
     float b = c.z;
 
-    float maxc = fmaxf(r, fmaxf(g, b));
+    float maxc = dmath_fmax(r, dmath_fmax(g, b));
-    float minc = fminf(r, fminf(g, b));
+    float minc = dmath_fmin(r, dmath_fmin(g, b));
     float v = maxc;
     if(minc == maxc) {
         return (c11_vec3){
@@ -163,7 +163,7 @@ static c11_vec3 srgb_to_hsv(c11_vec3 c) {
     } else {
         h = 4.0f + gc - rc;
     }
-    h = fmodf(h / 6.0f, 1.0f);
+    h = dmath_fmod(h / 6.0, 1.0);
     return (c11_vec3){
         {h, s, v}
     };

src/modules/easing.c

@@ -1,71 +1,69 @@
+#include "pocketpy/common/dmath.h"
 #include "pocketpy/pocketpy.h"
-#include "pocketpy/interpreter/vm.h"
-
-#include <math.h>
 
 // https://easings.net/
 
-const double kPi = 3.1415926545;
+const double kPi = DMATH_PI;
 
 static double easeLinear(double x) { return x; }
 
-static double easeInSine(double x) { return 1.0 - cos(x * kPi / 2); }
+static double easeInSine(double x) { return 1.0 - dmath_cos(x * kPi / 2); }
 
-static double easeOutSine(double x) { return sin(x * kPi / 2); }
+static double easeOutSine(double x) { return dmath_sin(x * kPi / 2); }
 
-static double easeInOutSine(double x) { return -(cos(kPi * x) - 1) / 2; }
+static double easeInOutSine(double x) { return -(dmath_cos(kPi * x) - 1) / 2; }
 
 static double easeInQuad(double x) { return x * x; }
 
-static double easeOutQuad(double x) { return 1 - pow(1 - x, 2); }
+static double easeOutQuad(double x) { return 1 - dmath_pow(1 - x, 2); }
 
 static double easeInOutQuad(double x) {
     if(x < 0.5) {
         return 2 * x * x;
     } else {
-        return 1 - pow(-2 * x + 2, 2) / 2;
+        return 1 - dmath_pow(-2 * x + 2, 2) / 2;
     }
 }
 
 static double easeInCubic(double x) { return x * x * x; }
 
-static double easeOutCubic(double x) { return 1 - pow(1 - x, 3); }
+static double easeOutCubic(double x) { return 1 - dmath_pow(1 - x, 3); }
 
 static double easeInOutCubic(double x) {
     if(x < 0.5) {
         return 4 * x * x * x;
     } else {
-        return 1 - pow(-2 * x + 2, 3) / 2;
+        return 1 - dmath_pow(-2 * x + 2, 3) / 2;
     }
 }
 
-static double easeInQuart(double x) { return pow(x, 4); }
+static double easeInQuart(double x) { return dmath_pow(x, 4); }
 
-static double easeOutQuart(double x) { return 1 - pow(1 - x, 4); }
+static double easeOutQuart(double x) { return 1 - dmath_pow(1 - x, 4); }
 
 static double easeInOutQuart(double x) {
     if(x < 0.5) {
-        return 8 * pow(x, 4);
+        return 8 * dmath_pow(x, 4);
     } else {
-        return 1 - pow(-2 * x + 2, 4) / 2;
+        return 1 - dmath_pow(-2 * x + 2, 4) / 2;
     }
 }
 
-static double easeInQuint(double x) { return pow(x, 5); }
+static double easeInQuint(double x) { return dmath_pow(x, 5); }
 
-static double easeOutQuint(double x) { return 1 - pow(1 - x, 5); }
+static double easeOutQuint(double x) { return 1 - dmath_pow(1 - x, 5); }
 
 static double easeInOutQuint(double x) {
     if(x < 0.5) {
-        return 16 * pow(x, 5);
+        return 16 * dmath_pow(x, 5);
     } else {
-        return 1 - pow(-2 * x + 2, 5) / 2;
+        return 1 - dmath_pow(-2 * x + 2, 5) / 2;
     }
 }
 
-static double easeInExpo(double x) { return x == 0 ? 0 : pow(2, 10 * x - 10); }
+static double easeInExpo(double x) { return x == 0 ? 0 : dmath_pow(2, 10 * x - 10); }
 
-static double easeOutExpo(double x) { return x == 1 ? 1 : 1 - pow(2, -10 * x); }
+static double easeOutExpo(double x) { return x == 1 ? 1 : 1 - dmath_pow(2, -10 * x); }
 
 static double easeInOutExpo(double x) {
     if(x == 0) {
@@ -73,21 +71,21 @@ static double easeInOutExpo(double x) {
     } else if(x == 1) {
         return 1;
     } else if(x < 0.5) {
-        return pow(2, 20 * x - 10) / 2;
+        return dmath_pow(2, 20 * x - 10) / 2;
     } else {
-        return (2 - pow(2, -20 * x + 10)) / 2;
+        return (2 - dmath_pow(2, -20 * x + 10)) / 2;
     }
 }
 
-static double easeInCirc(double x) { return 1 - sqrt(1 - pow(x, 2)); }
+static double easeInCirc(double x) { return 1 - dmath_sqrt(1 - dmath_pow(x, 2)); }
 
-static double easeOutCirc(double x) { return sqrt(1 - pow(x - 1, 2)); }
+static double easeOutCirc(double x) { return dmath_sqrt(1 - dmath_pow(x - 1, 2)); }
 
 static double easeInOutCirc(double x) {
     if(x < 0.5) {
-        return (1 - sqrt(1 - pow(2 * x, 2))) / 2;
+        return (1 - dmath_sqrt(1 - dmath_pow(2 * x, 2))) / 2;
     } else {
-        return (sqrt(1 - pow(-2 * x + 2, 2)) + 1) / 2;
+        return (dmath_sqrt(1 - dmath_pow(-2 * x + 2, 2)) + 1) / 2;
     }
 }
 
@@ -100,16 +98,16 @@ static double easeInBack(double x) {
 static double easeOutBack(double x) {
     const double c1 = 1.70158;
     const double c3 = c1 + 1;
-    return 1 + c3 * pow(x - 1, 3) + c1 * pow(x - 1, 2);
+    return 1 + c3 * dmath_pow(x - 1, 3) + c1 * dmath_pow(x - 1, 2);
 }
 
 static double easeInOutBack(double x) {
     const double c1 = 1.70158;
     const double c2 = c1 * 1.525;
     if(x < 0.5) {
-        return (pow(2 * x, 2) * ((c2 + 1) * 2 * x - c2)) / 2;
+        return (dmath_pow(2 * x, 2) * ((c2 + 1) * 2 * x - c2)) / 2;
     } else {
-        return (pow(2 * x - 2, 2) * ((c2 + 1) * (x * 2 - 2) + c2) + 2) / 2;
+        return (dmath_pow(2 * x - 2, 2) * ((c2 + 1) * (x * 2 - 2) + c2) + 2) / 2;
     }
 }
 
@@ -120,7 +118,7 @@ static double easeInElastic(double x) {
     } else if(x == 1) {
         return 1;
     } else {
-        return -pow(2, 10 * x - 10) * sin((x * 10 - 10.75) * c4);
+        return -dmath_pow(2, 10 * x - 10) * dmath_sin((x * 10 - 10.75) * c4);
     }
 }
 
@@ -131,7 +129,7 @@ static double easeOutElastic(double x) {
     } else if(x == 1) {
         return 1;
     } else {
-        return pow(2, -10 * x) * sin((x * 10 - 0.75) * c4) + 1;
+        return dmath_pow(2, -10 * x) * dmath_sin((x * 10 - 0.75) * c4) + 1;
     }
 }
 
@@ -142,9 +140,9 @@ static double easeInOutElastic(double x) {
     } else if(x == 1) {
         return 1;
     } else if(x < 0.5) {
-        return -(pow(2, 20 * x - 10) * sin((20 * x - 11.125) * c5)) / 2;
+        return -(dmath_pow(2, 20 * x - 10) * dmath_sin((20 * x - 11.125) * c5)) / 2;
     } else {
-        return (pow(2, -20 * x + 10) * sin((20 * x - 11.125) * c5)) / 2 + 1;
+        return (dmath_pow(2, -20 * x + 10) * dmath_sin((20 * x - 11.125) * c5)) / 2 + 1;
     }
 }
 

src/modules/json.c

@@ -1,10 +1,10 @@
 #include "pocketpy/pocketpy.h"
 
 #include "pocketpy/common/utils.h"
+#include "pocketpy/common/dmath.h"
 #include "pocketpy/objects/object.h"
 #include "pocketpy/common/sstream.h"
 #include "pocketpy/interpreter/vm.h"
-#include <math.h>
 
 static bool json_loads(int argc, py_Ref argv) {
     PY_CHECK_ARGC(1);
@@ -27,9 +27,9 @@ void pk__add_module_json() {
     py_setdict(mod, py_name("true"), py_True());
     py_setdict(mod, py_name("false"), py_False());
     py_TValue tmp;
-    py_newfloat(&tmp, NAN);
+    py_newfloat(&tmp, DMATH_NAN);
     py_setdict(mod, py_name("NaN"), &tmp);
-    py_newfloat(&tmp, INFINITY);
+    py_newfloat(&tmp, DMATH_INFINITY);
     py_setdict(mod, py_name("Infinity"), &tmp);
 
     py_bindfunc(mod, "loads", json_loads);
@@ -104,9 +104,9 @@ static bool json__write_object(c11_sbuf* buf, py_TValue* obj, int indent, int de
         case tp_NoneType: c11_sbuf__write_cstr(buf, "null"); return true;
         case tp_int: c11_sbuf__write_int(buf, obj->_i64); return true;
         case tp_float: {
-            if(isnan(obj->_f64)) {
+            if(dmath_isnan(obj->_f64)) {
                 c11_sbuf__write_cstr(buf, "NaN");
-            } else if(isinf(obj->_f64)) {
+            } else if(dmath_isinf(obj->_f64)) {
                 c11_sbuf__write_cstr(buf, obj->_f64 < 0 ? "-Infinity" : "Infinity");
             } else {
                 c11_sbuf__write_f64(buf, obj->_f64, -1);

src/modules/math.c

@@ -1,13 +1,7 @@
 #include "pocketpy/pocketpy.h"
+#include "pocketpy/common/dmath.h"
 #include "pocketpy/interpreter/vm.h"
 
-#include <math.h>
-
-#if PK_ENABLE_DETERMINISM
-    #ifndef _DMATH_H
-        #error "_DMATH_H not defined"
-    #endif
-#endif
 
 #define ONE_ARG_FUNC(name, func)                                                                   \
     static bool math_##name(int argc, py_Ref argv) {                                               \
@@ -37,10 +31,10 @@
         return true;                                                                               \
     }
 
-ONE_ARG_FUNC(ceil, ceil)
+ONE_ARG_FUNC(ceil, dmath_ceil)
-ONE_ARG_FUNC(fabs, fabs)
+ONE_ARG_FUNC(fabs, dmath_fabs)
-ONE_ARG_FUNC(floor, floor)
+ONE_ARG_FUNC(floor, dmath_floor)
-ONE_ARG_FUNC(trunc, trunc)
+ONE_ARG_FUNC(trunc, dmath_trunc)
 
 static bool math_fsum(int argc, py_Ref argv) {
     PY_CHECK_ARGC(1);
@@ -78,58 +72,58 @@ static bool math_gcd(int argc, py_Ref argv) {
     return true;
 }
 
-ONE_ARG_BOOL_FUNC(isfinite, isfinite)
+ONE_ARG_BOOL_FUNC(isfinite, dmath_isfinite)
-ONE_ARG_BOOL_FUNC(isinf, isinf)
+ONE_ARG_BOOL_FUNC(isinf, dmath_isinf)
-ONE_ARG_BOOL_FUNC(isnan, isnan)
+ONE_ARG_BOOL_FUNC(isnan, dmath_isnan)
 
 static bool math_isclose(int argc, py_Ref argv) {
     PY_CHECK_ARGC(2);
     double a, b;
     if(!py_castfloat(py_arg(0), &a)) return false;
     if(!py_castfloat(py_arg(1), &b)) return false;
-    py_newbool(py_retval(), fabs(a - b) < 1e-9);
+    py_newbool(py_retval(), dmath_fabs(a - b) < 1e-9);
     return true;
 }
 
-ONE_ARG_FUNC(exp, exp)
+ONE_ARG_FUNC(exp, dmath_exp)
 
 static bool math_log(int argc, py_Ref argv) {
     double x;
     if(!py_castfloat(py_arg(0), &x)) return false;
     if(argc == 1) {
-        py_newfloat(py_retval(), log(x));
+        py_newfloat(py_retval(), dmath_log(x));
     } else if(argc == 2) {
         double base;
         if(!py_castfloat(py_arg(1), &base)) return false;
-        py_newfloat(py_retval(), log(x) / log(base));
+        py_newfloat(py_retval(), dmath_log2(x) / dmath_log2(base));
     } else {
         return TypeError("log() takes 1 or 2 arguments");
     }
     return true;
 }
 
-ONE_ARG_FUNC(log2, log2)
+ONE_ARG_FUNC(log2, dmath_log2)
-ONE_ARG_FUNC(log10, log10)
+ONE_ARG_FUNC(log10, dmath_log10)
 
-TWO_ARG_FUNC(pow, pow)
+TWO_ARG_FUNC(pow, dmath_pow)
 
-ONE_ARG_FUNC(sqrt, sqrt)
+ONE_ARG_FUNC(sqrt, dmath_sqrt)
 
-ONE_ARG_FUNC(acos, acos)
+ONE_ARG_FUNC(acos, dmath_acos)
-ONE_ARG_FUNC(asin, asin)
+ONE_ARG_FUNC(asin, dmath_asin)
-ONE_ARG_FUNC(atan, atan)
+ONE_ARG_FUNC(atan, dmath_atan)
 
-ONE_ARG_FUNC(cos, cos)
+ONE_ARG_FUNC(cos, dmath_cos)
-ONE_ARG_FUNC(sin, sin)
+ONE_ARG_FUNC(sin, dmath_sin)
-ONE_ARG_FUNC(tan, tan)
+ONE_ARG_FUNC(tan, dmath_tan)
 
-TWO_ARG_FUNC(atan2, atan2)
+TWO_ARG_FUNC(atan2, dmath_atan2)
 
 static bool math_degrees(int argc, py_Ref argv) {
     PY_CHECK_ARGC(1);
     double x;
     if(!py_castfloat(py_arg(0), &x)) return false;
-    py_newfloat(py_retval(), x * PK_M_RAD2DEG);
+    py_newfloat(py_retval(), x * DMATH_RAD2DEG);
     return true;
 }
 
@@ -137,17 +131,17 @@ static bool math_radians(int argc, py_Ref argv) {
     PY_CHECK_ARGC(1);
     double x;
     if(!py_castfloat(py_arg(0), &x)) return false;
-    py_newfloat(py_retval(), x * PK_M_DEG2RAD);
+    py_newfloat(py_retval(), x * DMATH_DEG2RAD);
     return true;
 }
 
-TWO_ARG_FUNC(fmod, fmod)
+TWO_ARG_FUNC(fmod, dmath_fmod)
-TWO_ARG_FUNC(copysign, copysign)
+TWO_ARG_FUNC(copysign, dmath_copysign)
 
 static bool math_modf(int argc, py_Ref argv) {
     PY_CHECK_ARGC(1);
     double i;
-    double f = modf(py_tofloat(py_arg(0)), &i);
+    double f = dmath_modf(py_tofloat(py_arg(0)), &i);
     py_Ref p = py_newtuple(py_retval(), 2);
     py_newfloat(&p[0], f);
     py_newfloat(&p[1], i);
@@ -169,10 +163,10 @@ static bool math_factorial(int argc, py_Ref argv) {
 void pk__add_module_math() {
     py_Ref mod = py_newmodule("math");
 
-    py_newfloat(py_emplacedict(mod, py_name("pi")), PK_M_PI);
+    py_newfloat(py_emplacedict(mod, py_name("pi")), DMATH_PI);
-    py_newfloat(py_emplacedict(mod, py_name("e")), PK_M_E);
+    py_newfloat(py_emplacedict(mod, py_name("e")), DMATH_E);
-    py_newfloat(py_emplacedict(mod, py_name("inf")), INFINITY);
+    py_newfloat(py_emplacedict(mod, py_name("inf")), DMATH_INFINITY);
-    py_newfloat(py_emplacedict(mod, py_name("nan")), NAN);
+    py_newfloat(py_emplacedict(mod, py_name("nan")), DMATH_NAN);
 
     py_bindfunc(mod, "ceil", math_ceil);
     py_bindfunc(mod, "fabs", math_fabs);

src/modules/vmath.c

@@ -3,10 +3,10 @@
 
 #include "pocketpy/common/sstream.h"
 #include "pocketpy/common/utils.h"
+#include "pocketpy/common/dmath.h"
 #include "pocketpy/interpreter/vm.h"
-#include <math.h>
 
-static bool isclose(float a, float b) { return fabs(a - b) < 1e-4; }
+static bool isclose(float a, float b) { return dmath_fabs(a - b) < 1e-4; }
 
 #define DEFINE_VEC_FIELD(name, T, Tc, field)                                                       \
     static bool name##__##field(int argc, py_Ref argv) {                                           \
@@ -193,7 +193,7 @@ static py_Ref _const(py_Type type, const char* name) {
         float sum = 0;                                                                             \
         for(int i = 0; i < D; i++)                                                                 \
             sum += v.data[i] * v.data[i];                                                          \
-        py_newfloat(py_retval(), sqrtf(sum));                                                      \
+        py_newfloat(py_retval(), dmath_sqrt(sum));                                                \
         return true;                                                                               \
     }                                                                                              \
     static bool vec##D##_length_squared(int argc, py_Ref argv) {                                   \
@@ -223,7 +223,7 @@ static py_Ref _const(py_Type type, const char* name) {
         for(int i = 0; i < D; i++)                                                                 \
             len += self.data[i] * self.data[i];                                                    \
         if(isclose(len, 0)) return ZeroDivisionError("cannot normalize zero vector");              \
-        len = sqrtf(len);                                                                          \
+        len = dmath_sqrt(len);                                                                          \
         c11_vec##D res;                                                                            \
         for(int i = 0; i < D; i++)                                                                 \
             res.data[i] = self.data[i] / len;                                                      \
@@ -344,8 +344,8 @@ static bool vec2_rotate(int argc, py_Ref argv) {
     PY_CHECK_ARGC(2);
     py_f64 radians;
     if(!py_castfloat(&argv[1], &radians)) return false;
-    float cr = cosf(radians);
+    double sr, cr;
-    float sr = sinf(radians);
+    dmath_sincos(radians, &sr, &cr);
     c11_vec2 res;
     res.x = argv[0]._vec2.x * cr - argv[0]._vec2.y * sr;
     res.y = argv[0]._vec2.x * sr + argv[0]._vec2.y * cr;
@@ -357,9 +357,9 @@ static bool vec2_angle_STATIC(int argc, py_Ref argv) {
     PY_CHECK_ARGC(2);
     PY_CHECK_ARG_TYPE(0, tp_vec2);
     PY_CHECK_ARG_TYPE(1, tp_vec2);
-    float val = atan2f(argv[1]._vec2.y, argv[1]._vec2.x) - atan2f(argv[0]._vec2.y, argv[0]._vec2.x);
+    float val = dmath_atan2(argv[1]._vec2.y, argv[1]._vec2.x) - dmath_atan2(argv[0]._vec2.y, argv[0]._vec2.x);
-    if(val > PK_M_PI) val -= 2 * (float)PK_M_PI;
+    if(val > DMATH_PI) val -= 2 * (float)DMATH_PI;
-    if(val < -PK_M_PI) val += 2 * (float)PK_M_PI;
+    if(val < -DMATH_PI) val += 2 * (float)DMATH_PI;
     py_newfloat(py_retval(), val);
     return true;
 }
@@ -398,7 +398,7 @@ static bool vec2_smoothdamp_STATIC(int argc, py_Ref argv) {
     float maxChangeSq = maxChange * maxChange;
     float sqDist = change_x * change_x + change_y * change_y;
     if(sqDist > maxChangeSq) {
-        float mag = sqrtf(sqDist);
+        float mag = dmath_sqrt(sqDist);
         change_x = change_x / mag * maxChange;
         change_y = change_y / mag * maxChange;
     }
@@ -576,8 +576,8 @@ static bool inverse(const c11_mat3x3* m, c11_mat3x3* restrict out) {
 }
 
 static void trs(c11_vec2 t, float r, c11_vec2 s, c11_mat3x3* restrict out) {
-    float cr = cosf(r);
+    double sr, cr;
-    float sr = sinf(r);
+    dmath_sincos(r, &sr, &cr);
     // clang-format off
     *out = (c11_mat3x3){
         ._11 = s.x * cr, ._12 = -s.y * sr, ._13 = t.x,
@@ -733,7 +733,7 @@ static bool mat3x3_t(int argc, py_Ref argv) {
 static bool mat3x3_r(int argc, py_Ref argv) {
     PY_CHECK_ARGC(1);
     c11_mat3x3* ud = py_tomat3x3(argv);
-    float r = atan2f(ud->_21, ud->_11);
+    float r = dmath_atan2(ud->_21, ud->_11);
     py_newfloat(py_retval(), r);
     return true;
 }
@@ -742,8 +742,8 @@ static bool mat3x3_s(int argc, py_Ref argv) {
     PY_CHECK_ARGC(1);
     c11_mat3x3* ud = py_tomat3x3(argv);
     c11_vec2 res;
-    res.x = sqrtf(ud->_11 * ud->_11 + ud->_21 * ud->_21);
+    res.x = dmath_sqrt(ud->_11 * ud->_11 + ud->_21 * ud->_21);
-    res.y = sqrtf(ud->_12 * ud->_12 + ud->_22 * ud->_22);
+    res.y = dmath_sqrt(ud->_12 * ud->_12 + ud->_22 * ud->_22);
     py_newvec2(py_retval(), res);
     return true;
 }

tests/704_math.py

@@ -124,6 +124,7 @@ assert isclose(math.sin(0), 0.0)
 assert isclose(math.sin(math.pi / 2), 1.0)
 assert isclose(math.sin(math.pi), 0.0)
 assert isclose(math.sin(-math.pi / 2), -1.0)
+assert isclose(math.sin(-math.pi / 4), -0.7071067811865476)
 
 # cos - cosine
 assert isclose(math.cos(0), 1.0)

tests/710_cmath.py

@@ -10,7 +10,7 @@ assert isclose(c1 - c2, complex(1, -0.5))
 assert isclose(c1*4, complex(12, 16))
 assert isclose(c1*c2, complex(-12, 21.5))
 assert isclose(c2/c1, complex(0.96, 0.22))
-assert isclose(c2**2, complex(-16.25, 17.99999999999999))
+assert isclose(c2**2, complex(-16.25, 18))
 
 assert 1+2j == complex(1, 2) == 2j+1
 

tests/791_easing.py

@@ -1,12 +1,13 @@
 import easing
 
-def validate(val):
+def validate(k, f, x):
-    assert -2 <= val <= 2.0
+    val = f(x)
-    assert isinstance(val, float)
+    assert (-2 <= val <= 2.0), (k, x, val)
+    assert isinstance(val, float), (k, x, val)
 
 for k, f in easing.__dict__.items():
     if callable(f):
-        validate(f(0.2))
+        validate(k, f, 0.2)
-        validate(f(0.5))
+        validate(k, f, 0.5)
-        validate(f(0.8))
+        validate(k, f, 0.8)
-        validate(f(1.0))
+        validate(k, f, 1.0)