diff --git a/.github/workflows/python-app.yml b/.github/workflows/python-app.yml
index 7d35ac5..4c8d28e 100644
--- a/.github/workflows/python-app.yml
+++ b/.github/workflows/python-app.yml
@@ -1,4 +1,4 @@
-# This workflow will install Python dependencies, run tests and lint with a single version of Python
+# This workflow will install Python dependencies, run tests and lint with the supported versions of Python
 # For more information see: https://docs.github.com/en/actions/automating-builds-and-tests/building-and-testing-python
 name: Python application
 on:
@@ -11,12 +11,15 @@ permissions:
 jobs:
   build:
     runs-on: ubuntu-latest
+    strategy:
+      matrix:
+        python-version: ["3.11", "3.14"]
     steps:
     - uses: actions/checkout@v4
-    - name: Set up Python 3.10
-      uses: actions/setup-python@v3
+    - name: Set up Python ${{ matrix.python-version }}
+      uses: actions/setup-python@v5
       with:
-        python-version: "3.10"
+        python-version: ${{ matrix.python-version }}
     - name: Install dependencies
       run: |
         python -m pip install --upgrade pip
diff --git a/.github/workflows/test-coverage.yml b/.github/workflows/test-coverage.yml
index b048fed..1954393 100644
--- a/.github/workflows/test-coverage.yml
+++ b/.github/workflows/test-coverage.yml
@@ -9,15 +9,18 @@ on:
 jobs:
   test:
     runs-on: ubuntu-latest
+    strategy:
+      matrix:
+        python-version: ["3.11", "3.14"]
 
     steps:
       - name: Checkout Repository
         uses: actions/checkout@v4
 
-      - name: Set up Python
-        uses: actions/setup-python@v4
+      - name: Set up Python ${{ matrix.python-version }}
+        uses: actions/setup-python@v5
         with:
-          python-version: "3.9"
+          python-version: ${{ matrix.python-version }}
 
       - name: Install Dependencies
         run: |
diff --git a/CHANGELOG.md b/CHANGELOG.md
new file mode 100644
index 0000000..97a449b
--- /dev/null
+++ b/CHANGELOG.md
@@ -0,0 +1,37 @@
+# Changelog
+
+## 1.10.0-4 (package version 1.14.0.post4)
+
+The Git tags of this fork (`1.10.0-x`) and the package version differ: tags `1.10.0-2` and `1.10.0-3` already declare version `1.14.0`, because upstream 1.14.0 was merged between `1.10.0-1` and `1.10.0-2`. Tag `1.10.0-4` declares `1.14.0.post4`, so the installed version keeps increasing (`pip show mathjson-solver`).
+
+| Git tag | Package version |
+|---|---|
+| 1.10.0-1 | 1.10.0 |
+| 1.10.0-2 | 1.14.0 |
+| 1.10.0-3 | 1.14.0 |
+| 1.10.0-4 | 1.14.0.post4 |
+
+### Breaking changes
+- **`["Log", x]` is now the base-10 logarithm** (it was the natural logarithm). This aligns the solver with the CortexJS Compute Engine used by the front end. Consumers must migrate stored `["Log", x]` expressions to `["Ln", x]` to keep the previous behaviour.
+- **Python >= 3.11 is required** (`math.cbrt`).
+- Out-of-domain and non-finite results now raise instead of returning a value (see below).
+
+### Added
+- `["Log", x, b]`: logarithm of `x` in base `b`.
+- `MathJSONDomainError` (subclass of `ValueError`, exported from `mathjson_solver`). It is the inner exception (`MathJSONException.e`) when arguments are outside the mathematical domain of an operator:
+  - `Log`, `Ln`, `Log2`, `Log10` with `x <= 0`; `Log` with base `<= 0` or base `== 1`.
+  - `Sqrt` of a negative number.
+  - `Power` with a negative base and a non-integer exponent. (`Power(0, negative)` still raises `ZeroDivisionError`.)
+  - `Root` of a negative number with an even or fractional index.
+- `Root` supports odd roots of negative numbers: `["Root", -8, 3]` is `-2` and `["Root", -32, 5]` is `-2`. `∛` uses `math.cbrt`, and an exact integer root is returned exactly on every platform: `["Root", 27, 3]` is `3.0`, although glibc's `cbrt(27)` is `3.0000000000000004`.
+- Overflow handling: `OverflowError` is converted to `MathJSONException`. A float result that is infinite raises `OverflowError`, and a NaN result raises `MathJSONDomainError`, e.g. `["Multiply", 1e300, 1e300]`, `["Power", 10, 400]`. The check applies to float results of any construct only; booleans, strings, lists and dates are unaffected.
+- `ExponentialE` evaluates to `math.e`. Parameters and `Constants` with the same name take precedence. `Pi` is intentionally not provided.
+- `create_solver(parameters, strict=False)`. With `strict=True`, an expression whose head is not a supported operator raises `MathJSONException` (caused by `NotImplementedError`) instead of being returned unchanged. The default behaviour is unchanged.
+- `extract_variables` knows `Ln` and `Rational`, and does not report `ExponentialE` as a variable.
+
+### Changed
+- `Power`, `Sqrt`, `Root` and `Square` use `math.pow` (never produce complex numbers) and return floats, e.g. `["Power", 2, 3]` is `8.0` (equal to `8`).
+
+### Packaging / CI
+- `setup.cfg` and `setup.py` version set to `1.14.0.post4` (tag `1.10.0-4`).
+- GitHub workflows test on Python 3.11 and 3.14 with `actions/setup-python@v5`.
diff --git a/README.md b/README.md
index 8173b50..4acda47 100644
--- a/README.md
+++ b/README.md
@@ -43,6 +43,10 @@ except MathJSONException:
 
 Left unhandled, the exception will look like `MathJSONException("Problem in Divide. ['Divide', 1, 0]. division by zero")`.
 
+The original error is available as `e.e`. Arguments outside an operator's domain (e.g. `["Log", 0]`, `["Sqrt", -4]`) give a `MathJSONDomainError` (a `ValueError` subclass) and overflowing results give an `OverflowError`.
+
+By default, expressions with an unknown operator are returned unchanged. Use `create_solver(parameters, strict=True)` to raise a `MathJSONException` (caused by `NotImplementedError`) instead.
+
 
 ## How to run tests
 Make sure you have `pytest` installed. Then `cd` into project directory and run:
diff --git a/docs/README.md b/docs/README.md
index 8049ffe..2857c30 100644
--- a/docs/README.md
+++ b/docs/README.md
@@ -16,6 +16,7 @@
 * [Root](#root-and-square-root)
 * [Sqrt](#root-and-square-root)
 * [Exp](#exponents-and-logarithms)
+* [Ln](#exponents-and-logarithms)
 * [Log](#exponents-and-logarithms)
 * [Log2](#exponents-and-logarithms)
 * [Log10](#exponents-and-logarithms)
@@ -139,27 +140,38 @@ Performs a division.
 ```
 
 #### Square and Power
-`Power` raises a number to given power. `Square` is a special case of `Power`.
+`Power` raises a number to given power. `Square` is a special case of `Power`. A negative base with a non-integer exponent is undefined and raises a `MathJSONDomainError`.
 ```python
 ["Power", 2, 3]                   # 2^3=8
+["Power", 10, -3]                 # 10^-3=0.001
 ["Square", 4]                     # 4^2=16
 ```
 
 #### Root and square root
+`Root` also supports odd roots of negative numbers. Even or fractional roots of negative numbers (and `Sqrt` of a negative number) raise a `MathJSONDomainError`.
 ```python
 ["Root", 9, 2]                    # √9=3.0
 ["Root", 8, 3]                    # ∛8=2.0
+["Root", -8, 3]                   # ∛-8=-2.0
+["Root", -16, 4]                  # MathJSONDomainError
 ["Sqrt", 9]                       # √9=3.0
 ```
 
 #### Exponents and logarithms
+`Log` is the base-10 logarithm (as in CortexJS Compute Engine); use `Ln` for the natural logarithm. `Log` accepts an optional base as the second argument. Arguments outside the domain (`x <= 0`, base `<= 0` or base `== 1`) raise a `MathJSONDomainError`. `ExponentialE` evaluates to _e_ unless a parameter or constant with that name is defined.
 ```python
 ["Exp", 2]                        # e^2≅7.389
-["Log", 2.7183]                   # ln(2.7183)≅1.0000
+["Ln", 2.7183]                    # ln(2.7183)≅1.0000
+["Ln", "ExponentialE"]            # ln(e)=1.0
+["Log", 1000]                     # log10(1000)=3.0
+["Log", 8, 2]                     # log2(8)=3.0
 ["Log2", 8]                       # log2(8)=3.0
 ["Log10", 1000]                   # log10(1000)=3.0
 ```
 
+#### Out-of-range results
+Floating-point results that overflow to infinity raise a `MathJSONException` caused by an `OverflowError`; undefined (NaN) results raise one caused by a `MathJSONDomainError`.
+
 ### Conditionality
 
 #### Equality, greater than, less than
diff --git a/setup.cfg b/setup.cfg
index d8ec57c..4a7cbec 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -1,6 +1,6 @@
 [metadata]
 name = mathjson-solver
-version = 1.14.0
+version = 1.14.0.post4
 author = Martins Mednis
 author_email = mrt@mednis.info
 description = Utilities for MathJSON evaluation
@@ -18,7 +18,7 @@ classifiers =
 package_dir = 
 	= src
 packages = find:
-python_requires = >=3.6
+python_requires = >=3.11
 
 [options.packages.find]
 where = src
diff --git a/setup.py b/setup.py
index 5e32614..84810dc 100644
--- a/setup.py
+++ b/setup.py
@@ -5,7 +5,7 @@
 
 setuptools.setup(
     name="mathjson-solver",
-    version="1.14.0",
+    version="1.14.0.post4",
     author="Martins Mednis",
     author_email="mrt@mednis.info",
     description="Utilities for MathJSON evaluation",
@@ -22,5 +22,5 @@
     ],
     package_dir={"": "src"},
     packages=setuptools.find_packages(where="src"),
-    python_requires=">=3.6",
+    python_requires=">=3.11",
 )
diff --git a/src/mathjson_solver/__init__.py b/src/mathjson_solver/__init__.py
index 97e433e..cb3285f 100644
--- a/src/mathjson_solver/__init__.py
+++ b/src/mathjson_solver/__init__.py
@@ -1,2 +1,2 @@
 from .__main__ import create_mathjson_solver as create_solver
-from .__main__ import MathJSONException, extract_variables
+from .__main__ import MathJSONException, MathJSONDomainError, extract_variables
diff --git a/src/mathjson_solver/__main__.py b/src/mathjson_solver/__main__.py
index e8ebdff..922d86f 100644
--- a/src/mathjson_solver/__main__.py
+++ b/src/mathjson_solver/__main__.py
@@ -29,6 +29,10 @@ def __str__(self):
         return f"Problem in {self.construct}. {self.expr}. {m}"
 
 
+class MathJSONDomainError(ValueError):
+    """The arguments are outside the mathematical domain of the operator"""
+
+
 def is_numeric(x):
     try:
         float(x)
@@ -149,7 +153,7 @@ def has_sublist2(
         return count >= required_match_count
 
 
-def create_mathjson_solver(solver_parameters):
+def create_mathjson_solver(solver_parameters, strict=False):
     def f(s, *args):
         if args:
             c = deepcopy(args[0])
@@ -529,6 +533,69 @@ def safe_float(val):
                 except TypeError:
                     return False
 
+            def domain_safe(fn, *args):
+                """Call `fn`; math domain errors become MathJSONDomainError"""
+                try:
+                    return fn(*args)
+                except ValueError as e:
+                    raise MathJSONDomainError(str(e))
+
+            def Power(s):
+                x = f(s[1], c)
+                y = f(s[2], c)
+                if x == 0 and y < 0:
+                    # Keep the division-by-zero error that pow() raised
+                    raise ZeroDivisionError("0 cannot be raised to a negative power")
+                return domain_safe(math.pow, x, y)
+
+            def positive_root(x, n, r):
+                y = math.cbrt(x) if n == 3 else domain_safe(math.pow, x, r)
+                # cbrt and pow(x, 1/n) can be one ulp off an exact root:
+                # give exact roots back exactly
+                if (
+                    math.isfinite(y)
+                    and isinstance(n, (int, float))
+                    and float(n).is_integer()
+                    and n > 0
+                ):
+                    k = round(y)
+                    if k ** int(n) == x:
+                        return float(k)
+                return y
+
+            def Root(s):
+                x = f(s[1], c)
+                n = f(s[2], c)
+                r = 1.0 / n
+                if x < 0:
+                    # Only odd roots of negative numbers are real
+                    is_odd = (
+                        isinstance(n, (int, float))
+                        and float(n).is_integer()
+                        and n % 2 == 1
+                    )
+                    if not is_odd:
+                        raise MathJSONDomainError(
+                            "Even or fractional root of a negative number"
+                        )
+                    return -positive_root(-x, n, r)
+                return positive_root(x, n, r)
+
+            def Log(s):
+                x = f(s[1], c)
+                if x <= 0:
+                    raise MathJSONDomainError(
+                        "Logarithm of a non-positive number"
+                    )
+                if len(s) < 3:
+                    return math.log10(x)
+                b = f(s[2], c)
+                if b <= 0 or b == 1:
+                    raise MathJSONDomainError(
+                        "Logarithm base must be positive and not equal to 1"
+                    )
+                return math.log(x, b)
+
             constructs = {
                 "Sum": Sum,
                 "Add": Add,
@@ -545,15 +612,15 @@ def safe_float(val):
                 "Divide": lambda s: f(s[1], c) / f(s[2], c),
                 "Rational": lambda s: f(s[1], c) / f(s[2], c),
                 "Negate": lambda s: -f(s[1], c),
-                "Power": lambda s: pow(f(s[1], c), f(s[2], c)),
-                "Root": lambda s: pow(f(s[1], c), 1.0 / f(s[2], c)),
-                "Sqrt": lambda s: pow(f(s[1], c), 1.0 / 2),
-                "Square": lambda s: pow(f(s[1], c), 2),
+                "Power": Power,
+                "Root": Root,
+                "Sqrt": lambda s: domain_safe(math.pow, f(s[1], c), 1.0 / 2),
+                "Square": lambda s: domain_safe(math.pow, f(s[1], c), 2),
                 "Exp": lambda s: math.exp(f(s[1], c)),
-                "Log": lambda s: math.log(f(s[1], c)),
-                "Ln": lambda s: math.log(f(s[1], c)),
-                "Log2": lambda s: math.log2(f(s[1], c)),
-                "Log10": lambda s: math.log10(f(s[1], c)),
+                "Log": Log,
+                "Ln": lambda s: domain_safe(math.log, f(s[1], c)),
+                "Log2": lambda s: domain_safe(math.log2, f(s[1], c)),
+                "Log10": lambda s: domain_safe(math.log10, f(s[1], c)),
                 "Less": lambda s: Compare(operator.lt, s, c),
                 "Greater": lambda s: Compare(operator.gt, s, c),
                 "Equal": lambda s: Compare(operator.eq, s, c),
@@ -603,7 +670,13 @@ def safe_float(val):
             }
             if s[0] in constructs:
                 try:
-                    return constructs[s[0]](s)
+                    res = constructs[s[0]](s)
+                    # Floats must be finite: inf is overflow, nan is undefined
+                    if isinstance(res, float) and not math.isfinite(res):
+                        if math.isnan(res):
+                            raise MathJSONDomainError("Result is not a number")
+                        raise OverflowError("Result is too large")
+                    return res
 
                 # except RecursionError:
                 #     return s[0]
@@ -616,10 +689,15 @@ def safe_float(val):
                     raise MathJSONException(e, s, mathjson_construct=s[0])
                 except ZeroDivisionError as e:
                     raise MathJSONException(e, s, mathjson_construct=s[0])
+                except OverflowError as e:
+                    raise MathJSONException(e, s, mathjson_construct=s[0])
             else:
-                # raise MathJSONException(
-                #     NotImplementedError(f"'{s[0]}' is not supported"), s
-                # )
+                if strict and isinstance(s[0], str):
+                    raise MathJSONException(
+                        NotImplementedError(f"'{s[0]}' is not supported"),
+                        s,
+                        mathjson_construct=s[0],
+                    )
                 return s
         elif s in solver_parameters:
             try:
@@ -628,6 +706,9 @@ def safe_float(val):
                 return solver_parameters[s]
         elif s in c:
             return f(c[s], c)
+        elif s == "ExponentialE":
+            # Parameters and context take precedence over the constant
+            return math.e
         else:
             # raise KeyError(f"Parameter '{s}' is not defined")
             return s
@@ -652,8 +733,10 @@ def extract_variables(
         "Root",
         "Sqrt",
         "Square",
+        "Rational",
         "Exp",
         "Log",
+        "Ln",
         "Log2",
         "Log10",
         "Equal",
@@ -696,6 +779,7 @@ def extract_variables(
         "TimeDeltaHours",
         "TimeDeltaMinutes",
         "TimeDeltaDays",
+        "ExponentialE",
     ]
     if isinstance(s, str):
         if s in ignore_list:
diff --git a/tests/test_with_pytest.py b/tests/test_with_pytest.py
index 4317921..e807cad 100644
--- a/tests/test_with_pytest.py
+++ b/tests/test_with_pytest.py
@@ -1,5 +1,6 @@
 import sys
 import os
+import math
 import pytest
 import re
 
@@ -7,7 +8,12 @@
 
 sys.path.append(os.path.join(os.path.dirname(__file__), "../src/"))
 
-from mathjson_solver import create_solver, MathJSONException, extract_variables
+from mathjson_solver import (
+    create_solver,
+    MathJSONException,
+    MathJSONDomainError,
+    extract_variables,
+)
 
 
 @pytest.mark.parametrize(
@@ -37,8 +43,8 @@
         ({}, ["Square", 4], 16),
         ({}, ["Round", ["Exp", 2], 3], 7.389),
         ({}, ["Divide", 10, ["Add", 2, 3]], 2),
-        ({}, ["Round", ["Log", 2.7183], 3], 1.0),
         ({}, ["Round", ["Ln", 2.7183], 3], 1.0),
+        ({}, ["Round", ["Log", 1000], 3], 3.0),
         ({}, ["Round", ["Ln", 2.7183, 12], 3], 1.0),
         ({}, ["Log2", 8], 3),
         ({}, ["Log10", 1000], 3),
@@ -606,6 +612,135 @@ def test_solver_simple(parameters, expression, expected_result):
     assert solver(expression) == expected_result
 
 
+@pytest.mark.parametrize(
+    "parameters, expression, expected_result",
+    [
+        ({}, ["Log", 100], pytest.approx(2)),
+        ({}, ["Log", 1], pytest.approx(0)),
+        ({}, ["Log", 0.001], pytest.approx(-3)),
+        ({}, ["Log", 8, 2], pytest.approx(3)),
+        ({}, ["Ln", "ExponentialE"], pytest.approx(1)),
+        ({"ExponentialE": 10}, ["Add", "ExponentialE", 1], 11),
+        ({}, ["Constants", ["ExponentialE", 3], "ExponentialE"], 3),
+        ({}, ["Power", 10, 2], 100),
+        ({}, ["Power", 10, 0], 1),
+        ({}, ["Power", 10, -3], pytest.approx(0.001)),
+        ({}, ["Power", -2, 3], -8),
+        ({}, ["Root", 27, 3], 3),
+        ({}, ["Root", -8, 3], -2),
+        ({}, ["Root", -27, 3], -3),
+        ({}, ["Root", -32, 5], -2),
+        ({}, ["Root", -8, 3.0], -2),
+    ],
+)
+def test_solver_log_power_root(parameters, expression, expected_result):
+    solver = create_solver(parameters)
+    assert solver(expression) == expected_result
+
+
+@pytest.mark.parametrize(
+    "expression, inner_exception",
+    [
+        (["Log", 0], MathJSONDomainError),
+        (["Log", -5], MathJSONDomainError),
+        (["Log", 8, 1], MathJSONDomainError),
+        (["Log", 8, 0], MathJSONDomainError),
+        (["Log", 8, -2], MathJSONDomainError),
+        (["Ln", 0], MathJSONDomainError),
+        (["Log2", -1], MathJSONDomainError),
+        (["Log10", 0], MathJSONDomainError),
+        (["Sqrt", -4], MathJSONDomainError),
+        (["Root", -16, 4], MathJSONDomainError),
+        (["Root", -8, 2.5], MathJSONDomainError),
+        (["Power", -8, 0.333], MathJSONDomainError),
+        (["Power", 10, 400], OverflowError),
+        (["Power", 0, -1], ZeroDivisionError),
+        (["Multiply", 1e300, 1e300], OverflowError),
+        (["Exp", 1000], OverflowError),
+    ],
+)
+def test_raises_domain_and_overflow(expression, inner_exception):
+    solver = create_solver({})
+    with pytest.raises(MathJSONException) as exc_info:
+        solver(expression)
+    assert isinstance(exc_info.value.e, inner_exception)
+
+
+@pytest.mark.parametrize(
+    "expression, expected_result",
+    [
+        (["Root", 27, 3], 3),
+        (["Root", -27, 3], -3),
+        (["Root", 26, 3], 3.0000000000000004),
+    ],
+)
+def test_root_is_exact_when_cbrt_is_one_ulp_off(
+    monkeypatch, expression, expected_result
+):
+    # glibc's cbrt, as on Linux, answers 3.0000000000000004 for 27
+    monkeypatch.setattr(math, "cbrt", lambda x: 3.0000000000000004)
+    assert create_solver({})(expression) == expected_result
+
+
+@pytest.mark.parametrize(
+    "expression, pow_result, expected_result",
+    [
+        (["Root", 1e10, 10], 10.000000000000002, 10),
+        (["Root", -32, 5], 2.0000000000000004, -2),
+        (["Root", 33, 5], 2.0000000000000004, 2.0000000000000004),
+    ],
+)
+def test_root_is_exact_when_pow_is_one_ulp_off(
+    monkeypatch, expression, pow_result, expected_result
+):
+    monkeypatch.setattr(math, "pow", lambda x, y: pow_result)
+    assert create_solver({})(expression) == expected_result
+
+
+@pytest.mark.parametrize(
+    "value, inner_exception",
+    [
+        (float("nan"), MathJSONDomainError),
+        (float("inf"), OverflowError),
+    ],
+)
+def test_root_of_a_non_finite_value(value, inner_exception):
+    with pytest.raises(MathJSONException) as exc_info:
+        create_solver({"x": value})(["Root", "x", 3])
+    assert isinstance(exc_info.value.e, inner_exception)
+
+
+def test_domain_error_is_value_error():
+    assert issubclass(MathJSONDomainError, ValueError)
+
+
+def test_unknown_head_strict():
+    assert create_solver({})(["Sin", 1]) == ["Sin", 1]
+    assert create_solver({}, strict=False)(["Sin", 1]) == ["Sin", 1]
+    solver = create_solver({}, strict=True)
+    with pytest.raises(MathJSONException) as exc_info:
+        solver(["Sin", 1])
+    assert isinstance(exc_info.value.e, NotImplementedError)
+    assert exc_info.value.construct == "Sin"
+    with pytest.raises(MathJSONException):
+        solver(["Add", 1, ["Sin", 1]])
+
+
+def test_known_heads_in_strict_mode():
+    solver = create_solver({"a": ["Array", 1, 2, 3]}, strict=True)
+    assert solver(["Map", ["Array", 1, 2, 3], ["Square"]]) == [
+        "Array",
+        1,
+        4,
+        9,
+    ]
+    assert solver(["StrictMap", "a", ["Power"], 2]) == ["Array", 1, 4, 9]
+    assert solver(["Filter", "a", ["Greater"], 1]) == ["Array", 2, 3]
+    assert solver(["Sum", "a"]) == 6
+    assert solver(["Max", "a"]) == 3
+    assert solver(["Round", ["Log", 1000], 3]) == 3.0
+
+
 def test_raises_divide():
     solver = create_solver({})
     with pytest.raises(