Merge c048ec9faf4bfe02da004dce69471897933c3617 into 05c0337ba071129dabd167f33f48f1eb6cd75b47

Problem: When an exception occurs in a WITH block, __exit__ was not called,
preventing proper cleanup of context managers.

Solution:
1. Wrap WITH block body in try-except structure
2. On normal exit: call __exit__(None, None, None)
3. On exception: call __exit__ with exception info before re-raising

Changes:
- compiler.c: Wrap WITH body in try-except, ensure __exit__ called in both paths
- ceval.c: Update OP_WITH_EXIT to accept three arguments (exc_type, exc_val, exc_tb)
- tests/520_context.py: Add test to verify __exit__ called on exceptions

.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()

CMakeOptions.txt

@@ -8,7 +8,7 @@ endif()
 # system features
 option(PK_ENABLE_OS "" ON)
 option(PK_ENABLE_THREADS "" ON)
-option(PK_ENABLE_DETERMINISM "" OFF)
+option(PK_ENABLE_DETERMINISM "" ON)
 option(PK_ENABLE_WATCHDOG "" OFF)
 option(PK_ENABLE_CUSTOM_SNAME "" OFF)
 option(PK_ENABLE_MIMALLOC "" OFF)

build_android_libs.sh

@@ -16,7 +16,6 @@ cmake \
     -DPK_BUILD_SHARED_LIB=ON \
     -DCMAKE_BUILD_TYPE=Release \
     -DPK_ENABLE_OS=OFF \
-    -DPK_ENABLE_DETERMINISM=ON \
     -DPK_BUILD_MODULE_LZ4=ON \
     -DPK_BUILD_MODULE_CUTE_PNG=ON \
     -DPK_BUILD_MODULE_MSGPACK=ON

build_darwin_libs.sh

@@ -8,7 +8,6 @@ cd build
 
 FLAGS="-DPK_BUILD_STATIC_LIB=ON \
     -DPK_ENABLE_OS=OFF \
-    -DPK_ENABLE_DETERMINISM=ON \
     -DPK_BUILD_MODULE_LZ4=ON \
     -DPK_BUILD_MODULE_CUTE_PNG=ON \
     -DPK_BUILD_MODULE_MSGPACK=ON \

build_ios_libs.sh

@@ -10,7 +10,6 @@ FLAGS="-DCMAKE_TOOLCHAIN_FILE=3rd/ios.toolchain.cmake \
     -DDEPLOYMENT_TARGET=13.0 \
     -DPK_BUILD_STATIC_LIB=ON \
     -DPK_ENABLE_OS=OFF \
-    -DPK_ENABLE_DETERMINISM=ON \
     -DPK_BUILD_MODULE_LZ4=ON \
     -DPK_BUILD_MODULE_CUTE_PNG=ON \
     -DPK_BUILD_MODULE_MSGPACK=ON \

cmake_build.py

@@ -20,7 +20,7 @@ assert config in ['Debug', 'Release', 'RelWithDebInfo']
 
 os.chdir("build")
 
-code = os.system(f"cmake .. -DPK_ENABLE_MIMALLOC=ON -DPK_ENABLE_DETERMINISM=ON -DCMAKE_BUILD_TYPE={config} {extra_flags}")
+code = os.system(f"cmake .. -DPK_ENABLE_MIMALLOC=ON -DCMAKE_BUILD_TYPE={config} {extra_flags}")
 assert code == 0
 code = os.system(f"cmake --build . --config {config} -j 4")
 assert code == 0

docs/modules/colorcvt.md

@@ -1,10 +0,0 @@
----
-icon: package
-label: colorcvt
----
-
-Provide color conversion functions.
-
-#### Source code
-
-:::code source="../../include/typings/colorcvt.pyi" :::

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/common/dmath.h

@@ -0,0 +1,47 @@
+#pragma once
+
+#define DMATH_INFINITY ((float)((1e+300 * 1e+300)))
+#define DMATH_NAN ((float)(DMATH_INFINITY * 0.0F))
+#define DMATH_PI 3.1415926535897932384
+#define DMATH_E 2.7182818284590452354
+#define DMATH_DEG2RAD 0.017453292519943295
+#define DMATH_RAD2DEG 57.29577951308232
+#define DMATH_EPSILON 1e-10
+#define DMATH_LOG2_E 1.4426950408889634
+
+double dmath_exp2(double x);
+double dmath_log2(double x);
+
+double dmath_exp(double x);
+double dmath_exp10(double x);
+double dmath_log(double x);
+double dmath_log10(double x);
+double dmath_pow(double base, double exp);
+double dmath_sqrt(double x);
+double dmath_cbrt(double x);
+
+void dmath_sincos(double x, double* sin, double* cos);
+double dmath_sin(double x);
+double dmath_cos(double x);
+double dmath_tan(double x);
+double dmath_asin(double x);
+double dmath_acos(double x);
+double dmath_atan(double x);
+double dmath_atan2(double y, double x);
+
+int dmath_isinf(double x);
+int dmath_isnan(double x);
+int dmath_isnormal(double x);
+int dmath_isfinite(double x);
+
+double dmath_fmod(double x, double y);
+double dmath_copysign(double x, double y);
+
+double dmath_fabs(double x);
+double dmath_ceil(double x);
+double dmath_floor(double x);
+double dmath_trunc(double x);
+double dmath_modf(double x, double* intpart);
+
+double dmath_fmin(double x, double y);
+double dmath_fmax(double x, double y);

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

plugins/flutter/pocketpy/src/CMakeLists.txt

@@ -12,7 +12,6 @@ set(CMAKE_POSITION_INDEPENDENT_CODE ON)
 include(FetchContent)
 
 set(PK_ENABLE_OS OFF CACHE BOOL "" FORCE)
-set(PK_ENABLE_DETERMINISM ON CACHE BOOL "" FORCE)
 set(PK_BUILD_MODULE_LZ4 ON CACHE BOOL "" FORCE)
 set(PK_BUILD_MODULE_CUTE_PNG ON CACHE BOOL "" FORCE)
 set(PK_BUILD_MODULE_MSGPACK ON CACHE BOOL "" FORCE)

python/typing.py

@@ -57,3 +57,5 @@ assert_never = lambda x: x
 
 TypedDict = dict
 NotRequired = _PLACEHOLDER
+
+cast = lambda _, val: val

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/_generated.c

@@ -11,7 +11,7 @@ const char kPythonLibs_functools[] = "class cache:\n    def __init__(self, f):\n
 const char kPythonLibs_heapq[] = "# Heap queue algorithm (a.k.a. priority queue)\ndef heappush(heap, item):\n    \"\"\"Push item onto heap, maintaining the heap invariant.\"\"\"\n    heap.append(item)\n    _siftdown(heap, 0, len(heap)-1)\n\ndef heappop(heap):\n    \"\"\"Pop the smallest item off the heap, maintaining the heap invariant.\"\"\"\n    lastelt = heap.pop()    # raises appropriate IndexError if heap is empty\n    if heap:\n        returnitem = heap[0]\n        heap[0] = lastelt\n        _siftup(heap, 0)\n        return returnitem\n    return lastelt\n\ndef heapreplace(heap, item):\n    \"\"\"Pop and return the current smallest value, and add the new item.\n\n    This is more efficient than heappop() followed by heappush(), and can be\n    more appropriate when using a fixed-size heap.  Note that the value\n    returned may be larger than item!  That constrains reasonable uses of\n    this routine unless written as part of a conditional replacement:\n\n        if item > heap[0]:\n            item = heapreplace(heap, item)\n    \"\"\"\n    returnitem = heap[0]    # raises appropriate IndexError if heap is empty\n    heap[0] = item\n    _siftup(heap, 0)\n    return returnitem\n\ndef heappushpop(heap, item):\n    \"\"\"Fast version of a heappush followed by a heappop.\"\"\"\n    if heap and heap[0] < item:\n        item, heap[0] = heap[0], item\n        _siftup(heap, 0)\n    return item\n\ndef heapify(x):\n    \"\"\"Transform list into a heap, in-place, in O(len(x)) time.\"\"\"\n    n = len(x)\n    # Transform bottom-up.  The largest index there's any point to looking at\n    # is the largest with a child index in-range, so must have 2*i + 1 < n,\n    # or i < (n-1)/2.  If n is even = 2*j, this is (2*j-1)/2 = j-1/2 so\n    # j-1 is the largest, which is n//2 - 1.  If n is odd = 2*j+1, this is\n    # (2*j+1-1)/2 = j so j-1 is the largest, and that's again n//2-1.\n    for i in reversed(range(n//2)):\n        _siftup(x, i)\n\n# 'heap' is a heap at all indices >= startpos, except possibly for pos.  pos\n# is the index of a leaf with a possibly out-of-order value.  Restore the\n# heap invariant.\ndef _siftdown(heap, startpos, pos):\n    newitem = heap[pos]\n    # Follow the path to the root, moving parents down until finding a place\n    # newitem fits.\n    while pos > startpos:\n        parentpos = (pos - 1) >> 1\n        parent = heap[parentpos]\n        if newitem < parent:\n            heap[pos] = parent\n            pos = parentpos\n            continue\n        break\n    heap[pos] = newitem\n\ndef _siftup(heap, pos):\n    endpos = len(heap)\n    startpos = pos\n    newitem = heap[pos]\n    # Bubble up the smaller child until hitting a leaf.\n    childpos = 2*pos + 1    # leftmost child position\n    while childpos < endpos:\n        # Set childpos to index of smaller child.\n        rightpos = childpos + 1\n        if rightpos < endpos and not heap[childpos] < heap[rightpos]:\n            childpos = rightpos\n        # Move the smaller child up.\n        heap[pos] = heap[childpos]\n        pos = childpos\n        childpos = 2*pos + 1\n    # The leaf at pos is empty now.  Put newitem there, and bubble it up\n    # to its final resting place (by sifting its parents down).\n    heap[pos] = newitem\n    _siftdown(heap, startpos, pos)";
 const char kPythonLibs_linalg[] = "from vmath import *";
 const char kPythonLibs_operator[] = "# https://docs.python.org/3/library/operator.html#mapping-operators-to-functions\n\ndef le(a, b): return a <= b\ndef lt(a, b): return a < b\ndef ge(a, b): return a >= b\ndef gt(a, b): return a > b\ndef eq(a, b): return a == b\ndef ne(a, b): return a != b\n\ndef and_(a, b): return a & b\ndef or_(a, b): return a | b\ndef xor(a, b): return a ^ b\ndef invert(a): return ~a\ndef lshift(a, b): return a << b\ndef rshift(a, b): return a >> b\n\ndef is_(a, b): return a is b\ndef is_not(a, b): return a is not b\ndef not_(a): return not a\ndef truth(a): return bool(a)\ndef contains(a, b): return b in a\n\ndef add(a, b): return a + b\ndef sub(a, b): return a - b\ndef mul(a, b): return a * b\ndef truediv(a, b): return a / b\ndef floordiv(a, b): return a // b\ndef mod(a, b): return a % b\ndef pow(a, b): return a ** b\ndef neg(a): return -a\ndef matmul(a, b): return a @ b\n\ndef getitem(a, b): return a[b]\ndef setitem(a, b, c): a[b] = c\ndef delitem(a, b): del a[b]\n\ndef iadd(a, b): a += b; return a\ndef isub(a, b): a -= b; return a\ndef imul(a, b): a *= b; return a\ndef itruediv(a, b): a /= b; return a\ndef ifloordiv(a, b): a //= b; return a\ndef imod(a, b): a %= b; return a\n# def ipow(a, b): a **= b; return a\n# def imatmul(a, b): a @= b; return a\ndef iand(a, b): a &= b; return a\ndef ior(a, b): a |= b; return a\ndef ixor(a, b): a ^= b; return a\ndef ilshift(a, b): a <<= b; return a\ndef irshift(a, b): a >>= b; return a\n\nclass attrgetter:\n    def __init__(self, attr):\n        self.attr = attr\n    def __call__(self, obj):\n        return getattr(obj, self.attr)\n";
-const char kPythonLibs_typing[] = "class _Placeholder:\n    def __init__(self, *args, **kwargs):\n        pass\n    def __getitem__(self, *args):\n        return self\n    def __call__(self, *args, **kwargs):\n        return self\n    def __and__(self, other):\n        return self\n    def __or__(self, other):\n        return self\n    def __xor__(self, other):\n        return self\n\n\n_PLACEHOLDER = _Placeholder()\n\nSequence = _PLACEHOLDER\nList = _PLACEHOLDER\nDict = _PLACEHOLDER\nTuple = _PLACEHOLDER\nSet = _PLACEHOLDER\nAny = _PLACEHOLDER\nUnion = _PLACEHOLDER\nOptional = _PLACEHOLDER\nCallable = _PLACEHOLDER\nType = _PLACEHOLDER\nTypeAlias = _PLACEHOLDER\nNewType = _PLACEHOLDER\n\nClassVar = _PLACEHOLDER\n\nLiteral = _PLACEHOLDER\nLiteralString = _PLACEHOLDER\n\nIterable = _PLACEHOLDER\nGenerator = _PLACEHOLDER\nIterator = _PLACEHOLDER\n\nHashable = _PLACEHOLDER\n\nTypeVar = _PLACEHOLDER\nSelf = _PLACEHOLDER\n\nProtocol = object\nGeneric = object\nNever = object\n\nTYPE_CHECKING = False\n\n# decorators\noverload = lambda x: x\nfinal = lambda x: x\n\n# exhaustiveness checking\nassert_never = lambda x: x\n\nTypedDict = dict\nNotRequired = _PLACEHOLDER\n";
+const char kPythonLibs_typing[] = "class _Placeholder:\n    def __init__(self, *args, **kwargs):\n        pass\n    def __getitem__(self, *args):\n        return self\n    def __call__(self, *args, **kwargs):\n        return self\n    def __and__(self, other):\n        return self\n    def __or__(self, other):\n        return self\n    def __xor__(self, other):\n        return self\n\n\n_PLACEHOLDER = _Placeholder()\n\nSequence = _PLACEHOLDER\nList = _PLACEHOLDER\nDict = _PLACEHOLDER\nTuple = _PLACEHOLDER\nSet = _PLACEHOLDER\nAny = _PLACEHOLDER\nUnion = _PLACEHOLDER\nOptional = _PLACEHOLDER\nCallable = _PLACEHOLDER\nType = _PLACEHOLDER\nTypeAlias = _PLACEHOLDER\nNewType = _PLACEHOLDER\n\nClassVar = _PLACEHOLDER\n\nLiteral = _PLACEHOLDER\nLiteralString = _PLACEHOLDER\n\nIterable = _PLACEHOLDER\nGenerator = _PLACEHOLDER\nIterator = _PLACEHOLDER\n\nHashable = _PLACEHOLDER\n\nTypeVar = _PLACEHOLDER\nSelf = _PLACEHOLDER\n\nProtocol = object\nGeneric = object\nNever = object\n\nTYPE_CHECKING = False\n\n# decorators\noverload = lambda x: x\nfinal = lambda x: x\n\n# exhaustiveness checking\nassert_never = lambda x: x\n\nTypedDict = dict\nNotRequired = _PLACEHOLDER\n\ncast = lambda _, val: val\n";
 
 const char* load_kPythonLib(const char* name) {
     if (strchr(name, '.') != NULL) return NULL;

src/common/dmath.c

@@ -0,0 +1,706 @@
+#include "pocketpy/common/dmath.h"
+#include "pocketpy/common/algorithm.h"
+#include "pocketpy/common/_log_spline_tbl.h"
+#include <stdint.h>
+
+union Float64Bits {
+    double f;
+    uint64_t i;
+};
+
+/* IEEE 754 double precision floating point data manipulation */
+typedef union 
+{
+    double   f;
+    uint64_t u;
+    struct {int32_t  i0,i1;} s;
+}  udi_t;
+
+// https://github.com/akohlmey/fastermath/blob/master/src/exp.c#L63
+double dmath_exp2(double x) {
+    if (x > 1000) return DMATH_INFINITY;
+    if (x < -1000) return 0;
+	if (dmath_isnan(x)) return DMATH_NAN;
+    
+    const int FM_DOUBLE_BIAS = 1023;
+
+    static const double fm_exp2_q[] = {
+    /*  1.00000000000000000000e0, */
+        2.33184211722314911771e2,
+        4.36821166879210612817e3
+    };
+    static const double fm_exp2_p[] = {
+        2.30933477057345225087e-2,
+        2.02020656693165307700e1,
+        1.51390680115615096133e3
+    };
+
+    double   ipart, fpart, px, qx;
+    udi_t    epart;
+
+    ipart = dmath_floor(x+0.5);
+    fpart = x - ipart;
+
+    // FM_DOUBLE_INIT_EXP(epart,ipart);
+    epart.s.i0 = 0;
+    epart.s.i1 = (((int) ipart) + FM_DOUBLE_BIAS) << 20;
+
+    x = fpart*fpart;
+
+    px =        fm_exp2_p[0];
+    px = px*x + fm_exp2_p[1];
+    qx =    x + fm_exp2_q[0];
+    px = px*x + fm_exp2_p[2];
+    qx = qx*x + fm_exp2_q[1];
+
+    px = px * fpart;
+
+    x = 1.0 + 2.0*(px/(qx-px));
+    return epart.f*x;
+}
+
+double dmath_log2(double x) {
+	if(x < 0) return DMATH_NAN;
+	if(x == 0) return -DMATH_INFINITY;
+	if(x == DMATH_INFINITY) return DMATH_INFINITY;
+	if(dmath_isnan(x)) return DMATH_NAN;
+
+    const double fm_log_dinv =  4.09600000000000000000e+03;
+    const double fm_log_dsq6 =  9.93410746256510361521e-09;
+
+    const int FM_DOUBLE_BIAS = 1023;
+    const int FM_DOUBLE_EMASK = 2146435072;
+    const int FM_DOUBLE_MBITS = 20;
+    const int FM_DOUBLE_MMASK = 1048575;
+    const int FM_DOUBLE_EZERO = 1072693248;
+
+    const int FM_SPLINE_SHIFT = 8;
+
+    udi_t val;
+    double a,b,y;
+    int32_t hx, ipart;
+
+    val.f = x;
+    hx = val.s.i1;
+    
+    /* extract exponent and subtract bias */
+    ipart = (((hx & FM_DOUBLE_EMASK) >> FM_DOUBLE_MBITS) - FM_DOUBLE_BIAS);
+
+    /* mask out exponent to get the prefactor to 2**ipart */
+    hx &= FM_DOUBLE_MMASK;
+    val.s.i1 = hx | FM_DOUBLE_EZERO;
+    x = val.f;
+
+    /* table index */
+    hx >>= FM_SPLINE_SHIFT;
+
+    /* compute x value matching table index */
+    val.s.i0 = 0;
+    val.s.i1 = FM_DOUBLE_EZERO | (hx << FM_SPLINE_SHIFT);
+    b = (x - val.f) * fm_log_dinv;
+    a = 1.0 - b;
+
+    /* evaluate spline */
+    y = a * fm_log_q1[hx] + b * fm_log_q1[hx+1];
+    a = (a*a*a-a) * fm_log_q2[hx];
+    b = (b*b*b-b) * fm_log_q2[hx+1];
+    y += (a + b) * fm_log_dsq6;
+
+    return ((double) ipart) + (y * DMATH_LOG2_E);
+}
+
+double dmath_exp(double x) {
+    return dmath_exp2(x * DMATH_LOG2_E); // log2(e)
+}
+
+double dmath_exp10(double x) {
+    return dmath_exp2(x * 3.321928094887362); // log2(10)
+}
+
+double dmath_log(double x) {
+    return dmath_log2(x) / DMATH_LOG2_E; // log2(e)
+}
+
+double dmath_log10(double x) {
+    return dmath_log2(x) / 3.321928094887362; // log2(10)
+}
+
+double dmath_pow(double base, double exp) {
+    int exp_int = (int)exp;
+    if(exp_int == exp) {
+        if(exp_int == 0) return 1;
+        if(exp_int < 0) {
+            base = 1 / base;
+            exp_int = -exp_int;
+        }
+        double res = 1;
+        while(exp_int > 0) {
+            if(exp_int & 1) res *= base;
+            base *= base;
+            exp_int >>= 1;
+        }
+        return res;
+    }
+    if (base > 0) {
+        return dmath_exp(exp * dmath_log(base));
+    }
+    if (base == 0) {
+        if (exp > 0) return 0;
+        if (exp == 0) return 1;
+    }
+    return DMATH_NAN;
+}
+
+double dmath_sqrt(double x) {
+    return dmath_pow(x, 0.5);
+}
+
+double dmath_cbrt(double x) {
+    return dmath_pow(x, 1.0 / 3.0);
+}
+
+// https://github.com/kraj/musl/blob/kraj/master/src/math/sincos.c
+static double __sin(double x, double y, int iy)
+{
+static const double
+S1  = -1.66666666666666324348e-01, /* 0xBFC55555, 0x55555549 */
+S2  =  8.33333333332248946124e-03, /* 0x3F811111, 0x1110F8A6 */
+S3  = -1.98412698298579493134e-04, /* 0xBF2A01A0, 0x19C161D5 */
+S4  =  2.75573137070700676789e-06, /* 0x3EC71DE3, 0x57B1FE7D */
+S5  = -2.50507602534068634195e-08, /* 0xBE5AE5E6, 0x8A2B9CEB */
+S6  =  1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */
+
+	double z,r,v,w;
+
+	z = x*x;
+	w = z*z;
+	r = S2 + z*(S3 + z*S4) + z*w*(S5 + z*S6);
+	v = z*x;
+	if (iy == 0)
+		return x + v*(S1 + z*r);
+	else
+		return x - ((z*(0.5*y - v*r) - y) - v*S1);
+}
+
+static double __cos(double x, double y)
+{
+static const double
+C1  =  4.16666666666666019037e-02, /* 0x3FA55555, 0x5555554C */
+C2  = -1.38888888888741095749e-03, /* 0xBF56C16C, 0x16C15177 */
+C3  =  2.48015872894767294178e-05, /* 0x3EFA01A0, 0x19CB1590 */
+C4  = -2.75573143513906633035e-07, /* 0xBE927E4F, 0x809C52AD */
+C5  =  2.08757232129817482790e-09, /* 0x3E21EE9E, 0xBDB4B1C4 */
+C6  = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */
+
+	double hz,z,r,w;
+
+	z  = x*x;
+	w  = z*z;
+	r  = z*(C1+z*(C2+z*C3)) + w*w*(C4+z*(C5+z*C6));
+	hz = 0.5*z;
+	w  = 1.0-hz;
+	return w + (((1.0-w)-hz) + (z*r-x*y));
+}
+
+int __rem_pio2(double x, double *y)
+{
+static const double
+toint   = 1.5/2.22044604925031308085e-16,
+pio4    = 0x1.921fb54442d18p-1,
+invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
+pio2_1  = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */
+pio2_1t = 6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */
+pio2_2  = 6.07710050630396597660e-11, /* 0x3DD0B461, 0x1A600000 */
+pio2_2t = 2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */
+pio2_3  = 2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */
+pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
+
+	union Float64Bits u = { .f = x };
+	double z,w,t,r,fn;
+	double tx[3],ty[2];
+	uint32_t ix;
+	int sign, n, ex, ey, i;
+
+	sign = u.i>>63;
+	ix = u.i>>32 & 0x7fffffff;
+	if (ix <= 0x400f6a7a) {  /* |x| ~<= 5pi/4 */
+		if ((ix & 0xfffff) == 0x921fb)  /* |x| ~= pi/2 or 2pi/2 */
+			goto medium;  /* cancellation -- use medium case */
+		if (ix <= 0x4002d97c) {  /* |x| ~<= 3pi/4 */
+			if (!sign) {
+				z = x - pio2_1;  /* one round good to 85 bits */
+				y[0] = z - pio2_1t;
+				y[1] = (z-y[0]) - pio2_1t;
+				return 1;
+			} else {
+				z = x + pio2_1;
+				y[0] = z + pio2_1t;
+				y[1] = (z-y[0]) + pio2_1t;
+				return -1;
+			}
+		} else {
+			if (!sign) {
+				z = x - 2*pio2_1;
+				y[0] = z - 2*pio2_1t;
+				y[1] = (z-y[0]) - 2*pio2_1t;
+				return 2;
+			} else {
+				z = x + 2*pio2_1;
+				y[0] = z + 2*pio2_1t;
+				y[1] = (z-y[0]) + 2*pio2_1t;
+				return -2;
+			}
+		}
+	}
+	if (ix <= 0x401c463b) {  /* |x| ~<= 9pi/4 */
+		if (ix <= 0x4015fdbc) {  /* |x| ~<= 7pi/4 */
+			if (ix == 0x4012d97c)  /* |x| ~= 3pi/2 */
+				goto medium;
+			if (!sign) {
+				z = x - 3*pio2_1;
+				y[0] = z - 3*pio2_1t;
+				y[1] = (z-y[0]) - 3*pio2_1t;
+				return 3;
+			} else {
+				z = x + 3*pio2_1;
+				y[0] = z + 3*pio2_1t;
+				y[1] = (z-y[0]) + 3*pio2_1t;
+				return -3;
+			}
+		} else {
+			if (ix == 0x401921fb)  /* |x| ~= 4pi/2 */
+				goto medium;
+			if (!sign) {
+				z = x - 4*pio2_1;
+				y[0] = z - 4*pio2_1t;
+				y[1] = (z-y[0]) - 4*pio2_1t;
+				return 4;
+			} else {
+				z = x + 4*pio2_1;
+				y[0] = z + 4*pio2_1t;
+				y[1] = (z-y[0]) + 4*pio2_1t;
+				return -4;
+			}
+		}
+	}
+	if (ix < 0x413921fb) {  /* |x| ~< 2^20*(pi/2), medium size */
+medium:
+		/* rint(x/(pi/2)) */
+		fn = (double)x*invpio2 + toint - toint;
+		n = (int32_t)fn;
+		r = x - fn*pio2_1;
+		w = fn*pio2_1t;  /* 1st round, good to 85 bits */
+		/* Matters with directed rounding. */
+		if ((r - w < -pio4)) {
+			n--;
+			fn--;
+			r = x - fn*pio2_1;
+			w = fn*pio2_1t;
+		} else if ((r - w > pio4)) {
+			n++;
+			fn++;
+			r = x - fn*pio2_1;
+			w = fn*pio2_1t;
+		}
+		y[0] = r - w;
+		u.f = y[0];
+		ey = u.i>>52 & 0x7ff;
+		ex = ix>>20;
+		if (ex - ey > 16) { /* 2nd round, good to 118 bits */
+			t = r;
+			w = fn*pio2_2;
+			r = t - w;
+			w = fn*pio2_2t - ((t-r)-w);
+			y[0] = r - w;
+			u.f = y[0];
+			ey = u.i>>52 & 0x7ff;
+			if (ex - ey > 49) {  /* 3rd round, good to 151 bits, covers all cases */
+				t = r;
+				w = fn*pio2_3;
+				r = t - w;
+				w = fn*pio2_3t - ((t-r)-w);
+				y[0] = r - w;
+			}
+		}
+		y[1] = (r - y[0]) - w;
+		return n;
+	}
+
+    (void)tx;
+    (void)ty;
+    (void)i;
+    return 0;
+#if 0
+	/*
+	 * all other (large) arguments
+	 */
+	if (ix >= 0x7ff00000) {  /* x is inf or NaN */
+		y[0] = y[1] = x - x;
+		return 0;
+	}
+	/* set z = scalbn(|x|,-ilogb(x)+23) */
+	u.f = x;
+	u.i &= (uint64_t)-1>>12;
+	u.i |= (uint64_t)(0x3ff + 23)<<52;
+	z = u.f;
+	for (i=0; i < 2; i++) {
+		tx[i] = (double)(int32_t)z;
+		z     = (z-tx[i])*0x1p24;
+	}
+	tx[i] = z;
+	/* skip zero terms, first term is non-zero */
+	while (tx[i] == 0.0)
+		i--;
+	n = __rem_pio2_large(tx,ty,(int)(ix>>20)-(0x3ff+23),i+1,1);
+	if (sign) {
+		y[0] = -ty[0];
+		y[1] = -ty[1];
+		return -n;
+	}
+	y[0] = ty[0];
+	y[1] = ty[1];
+	return n;
+#endif
+}
+
+void dmath_sincos(double x, double *sin, double *cos) {
+	double y[2], s, c;
+	uint32_t ix;
+	unsigned n;
+
+	//GET_HIGH_WORD(ix, x);
+    union Float64Bits u = { .f = x };
+    ix = (uint32_t)(u.i >> 32);
+
+	ix &= 0x7fffffff;
+
+	/* |x| ~< pi/4 */
+	if (ix <= 0x3fe921fb) {
+		/* if |x| < 2**-27 * sqrt(2) */
+		if (ix < 0x3e46a09e) {
+			/* raise inexact if x!=0 and underflow if subnormal */
+
+			// FORCE_EVAL(ix < 0x00100000 ? x/0x1p120f : x+0x1p120f);
+            volatile double y_force_eval;
+            y_force_eval = ix < 0x00100000 ? x/0x1p120f : x+0x1p120f;
+            (void)y_force_eval;
+
+			*sin = x;
+			*cos = 1.0;
+			return;
+		}
+		*sin = __sin(x, 0.0, 0);
+		*cos = __cos(x, 0.0);
+		return;
+	}
+
+	/* sincos(Inf or NaN) is NaN */
+	if (ix >= 0x7ff00000) {
+		*sin = *cos = x - x;
+		return;
+	}
+
+	/* argument reduction needed */
+	n = __rem_pio2(x, y);
+	s = __sin(y[0], y[1], 1);
+	c = __cos(y[0], y[1]);
+	switch (n&3) {
+	case 0:
+		*sin = s;
+		*cos = c;
+		break;
+	case 1:
+		*sin = c;
+		*cos = -s;
+		break;
+	case 2:
+		*sin = -s;
+		*cos = -c;
+		break;
+	case 3:
+	default:
+		*sin = -c;
+		*cos = s;
+		break;
+	}
+}
+
+double dmath_sin(double x) {
+    double s, c;
+    dmath_sincos(x, &s, &c);
+    return s;
+}
+
+double dmath_cos(double x) {
+    double s, c;
+    dmath_sincos(x, &s, &c);
+    return c;
+}
+
+double dmath_tan(double x) {
+    double s, c;
+    dmath_sincos(x, &s, &c);
+    return s / c;
+}
+
+static double zig_r64(double z) {
+    const double pS0 = 1.66666666666666657415e-01;
+    const double pS1 = -3.25565818622400915405e-01;
+    const double pS2 = 2.01212532134862925881e-01;
+    const double pS3 = -4.00555345006794114027e-02;
+    const double pS4 = 7.91534994289814532176e-04;
+    const double pS5 = 3.47933107596021167570e-05;
+    const double qS1 = -2.40339491173441421878e+00;
+    const double qS2 = 2.02094576023350569471e+00;
+    const double qS3 = -6.88283971605453293030e-01;
+    const double qS4 = 7.70381505559019352791e-02;
+
+    double p = z * (pS0 + z * (pS1 + z * (pS2 + z * (pS3 + z * (pS4 + z * pS5)))));
+    double q = 1.0 + z * (qS1 + z * (qS2 + z * (qS3 + z * qS4)));
+    return p / q;
+}
+
+// https://github.com/ziglang/zig/blob/master/lib/std/math/asin.zig
+double dmath_asin(double x) {
+    const double pio2_hi = 1.57079632679489655800e+00;
+    const double pio2_lo = 6.12323399573676603587e-17;
+
+    union Float64Bits ux_union;
+    ux_union.f = x;
+    uint64_t ux = ux_union.i;
+    uint32_t hx = (uint32_t)(ux >> 32);
+    uint32_t ix = hx & 0x7FFFFFFF;
+
+    /* |x| >= 1 or nan */
+    if (ix >= 0x3FF00000) {
+        uint32_t lx = (uint32_t)(ux & 0xFFFFFFFF);
+
+        /* asin(1) = +-pi/2 with inexact */
+        if (((ix - 0x3FF00000) | lx) == 0) {
+            return x * pio2_hi + 0x1.0p-120;
+        } else {
+            return DMATH_NAN;
+        }
+    }
+
+    /* |x| < 0.5 */
+    if (ix < 0x3FE00000) {
+        /* if 0x1p-1022 <= |x| < 0x1p-26 avoid raising overflow */
+        if (ix < 0x3E500000 && ix >= 0x00100000) {
+            return x;
+        } else {
+            return x + x * zig_r64(x * x);
+        }
+    }
+
+    /* 1 > |x| >= 0.5 */
+    double z = (1 - dmath_fabs(x)) * 0.5;
+    double s = dmath_sqrt(z);
+    double r = zig_r64(z);
+    double fx;
+
+    /* |x| > 0.975 */
+    if (ix >= 0x3FEF3333) {
+        fx = pio2_hi - 2 * (s + s * r);
+    } else {
+        union Float64Bits jx_union = { .f = s };
+        uint64_t jx = jx_union.i;
+        union Float64Bits df_union = { .i = jx & 0xFFFFFFFF00000000ULL };
+        double df = df_union.f;
+        double c = (z - df * df) / (s + df);
+        fx = 0.5 * pio2_hi - (2 * s * r - (pio2_lo - 2 * c) - (0.5 * pio2_hi - 2 * df));
+    }
+
+    if (hx >> 31 != 0) {
+        return -fx;
+    } else {
+        return fx;
+    }
+}
+
+double dmath_acos(double x) {
+    return DMATH_PI / 2 - dmath_asin(x);
+}
+
+double dmath_atan(double x) {
+    return dmath_asin(x / dmath_sqrt(1 + x * x));
+}
+
+double dmath_atan2(double y, double x) {
+    if (x > 0) {
+        return dmath_atan(y / x);
+    } else if (x < 0 && y >= 0) {
+        return dmath_atan(y / x) + DMATH_PI;
+    } else if (x < 0 && y < 0) {
+        return dmath_atan(y / x) - DMATH_PI;
+    } else if (x == 0 && y > 0) {
+        return DMATH_PI / 2;
+    } else if (x == 0 && y < 0) {
+        return -DMATH_PI / 2;
+    } else {
+        return DMATH_NAN;
+    }
+}
+
+////////////////////////////////////////////////////////////////////
+
+int dmath_isinf(double x) {
+    union Float64Bits u = { .f = x };
+    return (u.i & -1ULL>>1) == 0x7ffULL<<52;
+}
+
+int dmath_isnan(double x) {
+    union Float64Bits u = { .f = x };
+    return (u.i & -1ULL>>1) > 0x7ffULL<<52;
+}
+
+int dmath_isnormal(double x) {
+    union Float64Bits u = { .f = x };
+    return ((u.i+(1ULL<<52)) & -1ULL>>1) >= 1ULL<<53;
+}
+
+int dmath_isfinite(double x) {
+    union Float64Bits u = { .f = x };
+    return (u.i & -1ULL>>1) < 0x7ffULL<<52;
+}
+
+// https://github.com/kraj/musl/blob/kraj/master/src/math/fmod.c
+double dmath_fmod(double x, double y) {
+	union Float64Bits ux = { .f = x }, uy = { .f = y };
+	int ex = ux.i>>52 & 0x7ff;
+	int ey = uy.i>>52 & 0x7ff;
+	int sx = ux.i>>63;
+	uint64_t i;
+
+	/* in the followings uxi should be ux.i, but then gcc wrongly adds */
+	/* float load/store to inner loops ruining performance and code size */
+	uint64_t uxi = ux.i;
+
+	if (uy.i<<1 == 0 || dmath_isnan(y) || ex == 0x7ff)
+		return (x*y)/(x*y);
+	if (uxi<<1 <= uy.i<<1) {
+		if (uxi<<1 == uy.i<<1)
+			return 0*x;
+		return x;
+	}
+
+	/* normalize x and y */
+	if (!ex) {
+		for (i = uxi<<12; i>>63 == 0; ex--, i <<= 1);
+		uxi <<= -ex + 1;
+	} else {
+		uxi &= -1ULL >> 12;
+		uxi |= 1ULL << 52;
+	}
+	if (!ey) {
+		for (i = uy.i<<12; i>>63 == 0; ey--, i <<= 1);
+		uy.i <<= -ey + 1;
+	} else {
+		uy.i &= -1ULL >> 12;
+		uy.i |= 1ULL << 52;
+	}
+
+	/* x mod y */
+	for (; ex > ey; ex--) {
+		i = uxi - uy.i;
+		if (i >> 63 == 0) {
+			if (i == 0)
+				return 0*x;
+			uxi = i;
+		}
+		uxi <<= 1;
+	}
+	i = uxi - uy.i;
+	if (i >> 63 == 0) {
+		if (i == 0)
+			return 0*x;
+		uxi = i;
+	}
+	for (; uxi>>52 == 0; uxi <<= 1, ex--);
+
+	/* scale result */
+	if (ex > 0) {
+		uxi -= 1ULL << 52;
+		uxi |= (uint64_t)ex << 52;
+	} else {
+		uxi >>= -ex + 1;
+	}
+	uxi |= (uint64_t)sx << 63;
+	ux.i = uxi;
+	return ux.f;
+}
+
+// https://github.com/kraj/musl/blob/kraj/master/src/math/copysign.c
+double dmath_copysign(double x, double y) {
+	union Float64Bits ux = { .f = x }, uy = { .f = y };
+	ux.i &= -1ULL/2;
+	ux.i |= uy.i & 1ULL<<63;
+	return ux.f;
+}
+
+double dmath_fabs(double x) {
+    return (x < 0) ? -x : x;
+}
+
+double dmath_ceil(double x) {
+	if(!dmath_isfinite(x)) return x;
+    int int_part = (int)x;
+    if (x > 0 && x != (double)int_part) {
+        return (double)(int_part + 1);
+    }
+    return (double)int_part;
+}
+
+double dmath_floor(double x) {
+	if(!dmath_isfinite(x)) return x;
+    int int_part = (int)x;
+    if (x < 0 && x != (double)int_part) {
+        return (double)(int_part - 1);
+    }
+    return (double)int_part;
+}
+
+double dmath_trunc(double x) {
+    return (double)((int)x);
+}
+
+// https://github.com/kraj/musl/blob/kraj/master/src/math/modf.c
+double dmath_modf(double x, double* iptr) {
+	union Float64Bits u = { .f = x };
+	uint64_t mask;
+	int e = (int)(u.i>>52 & 0x7ff) - 0x3ff;
+
+	/* no fractional part */
+	if (e >= 52) {
+		*iptr = x;
+		if (e == 0x400 && u.i<<12 != 0) /* nan */
+			return x;
+		u.i &= 1ULL<<63;
+		return u.f;
+	}
+
+	/* no integral part*/
+	if (e < 0) {
+		u.i &= 1ULL<<63;
+		*iptr = u.f;
+		return x;
+	}
+
+	mask = -1ULL>>12>>e;
+	if ((u.i & mask) == 0) {
+		*iptr = x;
+		u.i &= 1ULL<<63;
+		return u.f;
+	}
+	u.i &= ~mask;
+	*iptr = u.f;
+	return x - u.f;
+}
+
+double dmath_fmin(double x, double y) {
+    return (x < y) ? x : y;
+}
+
+double dmath_fmax(double x, double y) {
+    return (x > y) ? x : y;
+}

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/compiler/compiler.c

@@ -2801,6 +2801,8 @@ static Error* compile_stmt(Compiler* self) {
         case TK_WITH: {
             check(EXPR(self));  // [ <expr> ]
             Ctx__s_emit_top(ctx());
+            // Save context manager for later __exit__ call
+            Ctx__emit_(ctx(), OP_DUP_TOP, BC_NOARG, prev()->line);
             Ctx__enter_block(ctx(), CodeBlockType_WITH);
             NameExpr* as_name = NULL;
             if(match(TK_AS)) {
@@ -2809,17 +2811,33 @@ static Error* compile_stmt(Compiler* self) {
                 as_name = NameExpr__new(prev()->line, name, name_scope(self));
             }
             Ctx__emit_(ctx(), OP_WITH_ENTER, BC_NOARG, prev()->line);
-            // [ <expr> <expr>.__enter__() ]
             if(as_name) {
                 bool ok = vtemit_store((Expr*)as_name, ctx());
                 vtdelete((Expr*)as_name);
                 if(!ok) return SyntaxError(self, "invalid syntax");
             } else {
-                // discard `__enter__()`'s return value
                 Ctx__emit_(ctx(), OP_POP_TOP, BC_NOARG, BC_KEEPLINE);
             }
+            // Wrap body in try-except to ensure __exit__ is called even on exception
+            Ctx__enter_block(ctx(), CodeBlockType_TRY);
+            Ctx__emit_(ctx(), OP_BEGIN_TRY, BC_NOARG, prev()->line);
             check(compile_block_body(self));
+            Ctx__emit_(ctx(), OP_END_TRY, BC_NOARG, BC_KEEPLINE);
+            // Normal exit: call __exit__(None, None, None)
+            Ctx__emit_(ctx(), OP_LOAD_NONE, BC_NOARG, prev()->line);
+            Ctx__emit_(ctx(), OP_LOAD_NONE, BC_NOARG, prev()->line);
+            Ctx__emit_(ctx(), OP_LOAD_NONE, BC_NOARG, prev()->line);
             Ctx__emit_(ctx(), OP_WITH_EXIT, BC_NOARG, prev()->line);
+            int jump_patch = Ctx__emit_(ctx(), OP_JUMP_FORWARD, BC_NOARG, BC_KEEPLINE);
+            Ctx__exit_block(ctx());
+            // Exception handler: call __exit__ with exception info, then re-raise
+            Ctx__emit_(ctx(), OP_PUSH_EXCEPTION, BC_NOARG, BC_KEEPLINE);
+            Ctx__emit_(ctx(), OP_LOAD_NONE, BC_NOARG, BC_KEEPLINE);  // exc_type
+            Ctx__emit_(ctx(), OP_ROT_TWO, BC_NOARG, BC_KEEPLINE);     // reorder: [cm, None, exc]
+            Ctx__emit_(ctx(), OP_LOAD_NONE, BC_NOARG, BC_KEEPLINE);  // exc_tb
+            Ctx__emit_(ctx(), OP_WITH_EXIT, BC_NOARG, prev()->line);
+            Ctx__emit_(ctx(), OP_RE_RAISE, BC_NOARG, BC_KEEPLINE);
+            Ctx__patch_jump(ctx(), jump_patch);
             Ctx__exit_block(ctx());
         } break;
         /*************************************************/

src/interpreter/ceval.c

@@ -1122,14 +1122,35 @@ __NEXT_STEP:
             DISPATCH();
         }
         case OP_WITH_EXIT: {
-            // [expr]
+            // Stack: [cm, exc_type, exc_val, exc_tb]
-            py_push(TOP());
+            // Call cm.__exit__(exc_type, exc_val, exc_tb)
+            py_Ref exc_tb = TOP();
+            py_Ref exc_val = SECOND();
+            py_Ref exc_type = THIRD();
+            py_Ref cm = FOURTH();
+            
+            // Save all values from stack
+            py_TValue saved_cm = *cm;
+            py_TValue saved_exc_type = *exc_type;
+            py_TValue saved_exc_val = *exc_val;
+            py_TValue saved_exc_tb = *exc_tb;
+            self->stack.sp -= 4;
+            
+            // Push cm and get __exit__ method
+            py_push(&saved_cm);
             if(!py_pushmethod(__exit__)) {
-                TypeError("'%t' object does not support the context manager protocol", TOP()->type);
+                TypeError("'%t' object does not support the context manager protocol", saved_cm.type);
                 goto __ERROR;
             }
-            if(!py_vectorcall(0, 0)) goto __ERROR;
+            
-            POP();
+            // Push arguments: exc_type, exc_val, exc_tb
+            PUSH(&saved_exc_type);
+            PUSH(&saved_exc_val);
+            PUSH(&saved_exc_tb);
+            
+            // Call __exit__(exc_type, exc_val, exc_tb)
+            if(!py_vectorcall(3, 0)) goto __ERROR;
+            py_pop();  // discard return value
             DISPATCH();
         }
         ///////////

src/interpreter/vm.c

@@ -242,7 +242,7 @@ void VM__ctor(VM* self) {
     pk__add_module_stdc();
     pk__add_module_vmath();
     pk__add_module_array2d();
-    pk__add_module_colorcvt();
+    // pk__add_module_colorcvt();
 
     // add modules
     pk__add_module_os();

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/520_context.py

@@ -27,4 +27,29 @@ assert path == ['enter', 'in', 'exit']
 
 path.clear()
 
+# Test that __exit__ is called even when an exception occurs
+class B:
+    def __init__(self):
+        self.path = []
+    
+    def __enter__(self):
+        path.append('enter')
+        return self
+    
+    def __exit__(self, exc_type, exc_val, exc_tb):
+        path.append('exit')
+        if exc_type is not None:
+            path.append('exception')
+        return False  # propagate exception
+
+try:
+    with B():
+        path.append('before_raise')
+        raise ValueError('test')
+        path.append('after_raise')  # should not be reached
+except ValueError:
+    pass
+
+assert path == ['enter', 'before_raise', 'exit', 'exception'], f"Expected ['enter', 'before_raise', 'exit', 'exception'], got {path}"
+
 

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)

tests/903_colorcvt.py

@@ -1,3 +1,5 @@
+exit()
+
 import colorcvt
 from vmath import vec3
 

tests/930_deterministic_float.py

@@ -8,7 +8,7 @@ if config["PK_ENABLE_DETERMINISM"] == 0:
 def assertEqual(a, b):
     if a == b:
         return
-    print(f'{a} != {b}')
+    print(f'{a} != {b} ({a-b})')
     raise AssertionError
 
 # test constants
@@ -54,18 +54,18 @@ assertEqual(math.isnan(math.nan), True)
 # test exp
 assertEqual(math.exp(0), 1.0)
 assertEqual(math.exp(1), math.e)
-assertEqual(math.exp(1.5), 4.481689070338065 - 8.881784197001252e-16)
+assertEqual(math.exp(1.5), 4.48168907033806362960604019463) #4.481689070338065 - 8.881784197001252e-16)
-assertEqual(math.exp(3), 20.08553692318767 - 3.552713678800501e-15)
+assertEqual(math.exp(3), 20.0855369231876608182574273087) #20.08553692318767 - 3.552713678800501e-15)
-assertEqual(math.exp(-3), 0.04978706836786394 + 6.938893903907228e-18)
+assertEqual(math.exp(-3), 0.04978706836786396527916309651) #0.04978706836786394 + 6.938893903907228e-18)
 assertEqual(math.exp(-2.253647), 0.1050155336754953 - 1.387778780781446e-17)
-assertEqual(math.exp(4.729036), 113.1863980522005 - 4.263256414560601e-14)
+assertEqual(math.exp(4.729036), 113.186398052200445363268954679) #113.1863980522005 - 4.263256414560601e-14)
 
 # test log series
 assertEqual(math.log(0), -math.inf)
 assertEqual(math.log(1), 0.0)
 assertEqual(math.log(2), 0.69314718055994530942)
 assertEqual(math.log(math.e), 1.0)
-assertEqual(math.log(10), 2.30258509299404568402)
+assertEqual(math.log(10), 2.30258509299404545700440394284) #2.30258509299404568402)
 assertEqual(math.log(28.897124), 3.363742074595449)
 assertEqual(math.log2(math.e), 1.4426950408889634074)
 assertEqual(math.log2(78.781291), 6.299781153677818)
@@ -77,13 +77,13 @@ assertEqual(math.pow(2,2), 4.0)
 assertEqual(math.pow(1.41421356237309504880, 2), 2.0 + 4.440892098500626e-16)
 assertEqual(math.pow(0.70710678118654752440, 2), 0.5000000000000001)
 assertEqual(math.pow(-1.255782,-3), -0.5049603042167915)
-assertEqual(math.pow(6.127042, 4.071529), 1604.407546456745 + 2.273736754432321e-13)
+assertEqual(math.pow(6.127042, 4.071529), 1604.40754645674428502388764172) #1604.407546456745 + 2.273736754432321e-13)
 
 # test sqrt
-assertEqual(math.sqrt(2), 1.41421356237309504880)
+assertEqual(math.sqrt(2), 1.41421356237309492343001693370) #1.41421356237309504880)
 assertEqual(math.sqrt(math.pi), 1.772453850905516 - 2.220446049250313e-16)
 assertEqual(math.sqrt(125.872509), 11.21929182257062)
-assertEqual(math.sqrt(1225.296280), 35.00423231553579)
+assertEqual(math.sqrt(1225.296280), 35.0042323155358019448613049462) #35.00423231553579)
 
 # test cos, sin, tan
 assertEqual(math.cos(0), 1.0)
@@ -112,19 +112,19 @@ assertEqual(math.asin(1), 1.570796326794897 - 4.440892098500626e-16)
 assertEqual(math.asin(-0.225895), -0.2278616865773913 + 2.775557561562891e-17)
 assertEqual(math.asin(0.955658), 1.271886195819423 + 4.440892098500626e-16)
 assertEqual(math.atan(0), 0.0)
-assertEqual(math.atan(1), 0.7853981633974483)
+assertEqual(math.atan(1), 0.78539816339744839002179332965) #0.7853981633974483)
-assertEqual(math.atan(-3.758927), -1.310785284610617 - 4.440892098500626e-16)
+assertEqual(math.atan(-3.758927), -1.3107852846106160527028805518) #-1.310785284610617 - 4.440892098500626e-16)
-assertEqual(math.atan(35.789293), 1.542862277280122)
+assertEqual(math.atan(35.789293), 1.54286227728011748894232368911) #1.542862277280122)
 
 # test atan2
-assertEqual(math.atan2(math.pi/4, math.pi/4), 0.7853981633974483)
+assertEqual(math.atan2(math.pi/4, math.pi/4), 0.78539816339744839002179332965) #0.7853981633974483)
-assertEqual(math.atan2(-math.pi/4, math.pi/4), -0.7853981633974483)
+assertEqual(math.atan2(-math.pi/4, math.pi/4), -0.7853981633974483900217933296) #-0.7853981633974483)
 assertEqual(math.atan2(-math.pi/4, -math.pi/4), -2.356194490192345)
 assertEqual(math.atan2(math.pi/4, -math.pi/4), 2.356194490192345)
-assertEqual(math.atan2(1.573823, 0.685329), 1.160103682924653)
+assertEqual(math.atan2(1.573823, 0.685329), 1.16010368292465315676054160576) #1.160103682924653)
-assertEqual(math.atan2(-0.899663, 0.668972), -0.9314162757114095)
+assertEqual(math.atan2(-0.899663, 0.668972), -0.9314162757114096136135117376) #-0.9314162757114095)
-assertEqual(math.atan2(-0.762894, -0.126497), -1.735113347173296 - 4.440892098500626e-16)
+assertEqual(math.atan2(-0.762894, -0.126497), -1.7351133471732969049128314509) #-1.735113347173296 - 4.440892098500626e-16)
-assertEqual(math.atan2(0.468463, -0.992734), 2.700683410692374 - 4.440892098500626e-16)
+assertEqual(math.atan2(0.468463, -0.992734), 2.70068341069237316531825854326) #2.700683410692374 - 4.440892098500626e-16)
 
 # test fsum, sum
 fsum_sin = math.fsum([math.sin(i) for i in range(5000)])