Implemented multivariate polynomial interpolation algorithm.

CHANGELOG.md

@@ -7,6 +7,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
+### Added
+
+- Multivariate Newton interpolation algorithm, `multivariate_Newton_polynomial_interpolation`.
+
 ## [0.2.0] - 2024-01-02
 
 ### Added

pyadic/interpolation.py

@@ -2,6 +2,7 @@
 import numpy
 import random
 import sympy
+import inspect
 
 from collections.abc import Iterator
 
@@ -72,6 +73,58 @@ def fsubtracted(i, t):
     return tpoly.as_expr()
 
 
+def multivariate_Newton_polynomial_interpolation(f, prime, seed=0, depth=0, verbose=False):
+    """Recursive multivariate polynomial interpolation of f(ts), samples taken modulo prime."""
+
+    # End of recursion condition: univariate function
+    function_signature = inspect.signature(f)
+    num_args = len(function_signature.parameters)
+    if num_args == 1:
+        tpoly = Newton_polynomial_interpolation(f, prime, seed=seed, verbose=verbose)
+        tpoly = tpoly.subs({'t': 't1', })
+        return tpoly
+
+    t_sequence_generator = FFSequenceGenerator(prime, seed)
+    avals, tvals, subtracted = [], [], [f]
+
+    dynamic_fsubtracted_string = f"""def fsubtracted(i, {''.join([f't{i}, ' for i in range(1, num_args + 1)])}):
+    return (subtracted[i - 1]({''.join([f't{i}, ' for i in range(1, num_args + 1)])}) - \
+    ModP(int(avals[i - 1]({''.join([f'int(t{i}), ' for i in range(2, num_args + 1)])})), {prime})) / (t1 - tvals[i - 1])
+    """
+    local_dict = {'ModP': ModP, 'subtracted': subtracted, 'avals': avals, 'tvals': tvals}  # Include any necessary objects in the local dictionary
+    exec(dynamic_fsubtracted_string, local_dict)
+    fsubtracted = local_dict['fsubtracted']
+
+    MAX_SAMPLES = 1000  # sets iteration limit
+    for counter in range(MAX_SAMPLES):
+        if avals[-2:] == [0, 0]:
+            break
+        if verbose:
+            if depth != 0:
+                print()
+            print(f"@ depth: {depth} - samples: {len(avals)}, {(avals, tvals)}", end="\n")
+        tvals += [next(t_sequence_generator)]
+        local_dict = {'ModP': ModP, 'subtracted': subtracted}  # Include any necessary objects in the local dictionary
+        function_string = f"""def rest_function({''.join([f't{i}, ' for i in range(2, num_args + 1)])}):
+        return subtracted[-1](ModP('{tvals[-1]}'), {''.join([f't{i}, ' for i in range(2, num_args + 1)])})"""
+        exec(function_string, local_dict)
+        rest_function = local_dict['rest_function']
+        avals += [multivariate_Newton_polynomial_interpolation(rest_function, prime, seed=seed, depth=depth + 1, verbose=verbose)]
+        avals[-1] = avals[-1].subs({f't{i}': f't{i+1}' for i in range(1, num_args)}, simultaneous=True)
+        avals[-1] = sympy.poly(avals[-1], sympy.symbols([f't{i+1}' for i in range(1, num_args)]), modulus=prime)
+        fsubtracted_partial = functools.partial(fsubtracted, len(avals))
+        subtracted += [fsubtracted_partial]
+
+    if verbose:
+        print(f"\nFinished after {len(avals)} samples: {avals}.", end="\n")
+
+    FFGF = sympy.GF(prime).frac_field(*sympy.symbols([f't{i}' for i in range(1, num_args + 1)]))
+    tpoly = FFGF(0)
+    for aval, tval in zip(avals[:-2][::-1], tvals[:-2][::-1]):
+        tpoly = FFGF(aval.as_expr()) + (FFGF(sympy.symbols('t1')) - int(tval)) * tpoly
+    return tpoly.as_expr()
+
+
 class FFSequenceGenerator(Iterator):
     """Random number generator decoupled from global state."""
 

tests/test_interpolation.py

@@ -1,8 +1,10 @@
 import sympy
 
-from pyadic.interpolation import Newton_polynomial_interpolation, Thiele_rational_interpolation
+from pyadic.interpolation import Newton_polynomial_interpolation, Thiele_rational_interpolation, \
+    multivariate_Newton_polynomial_interpolation
 
 t = sympy.symbols('t')
+t1, t2, t3 = sympy.symbols('t1:4')
 
 
 def Ptest0(tval):
@@ -13,6 +15,10 @@ def Ptest1(tval):
     return (tval ** 20 + tval - 1)
 
 
+def Ptest2(t1, t2, t3):
+    return (t1 + 2 * t2) + 3 * (t1 * t2) + t3 ** 5
+
+
 def Rtest0(tval):
     return (1)
 
@@ -29,6 +35,10 @@ def test_Newton_polynomial_interpolation():
     assert Newton_polynomial_interpolation(Ptest1, 2 ** 31 - 1, verbose=True) == t ** 20 + t - 1
 
 
+def test_Newton_polynomial_interpolation_multivariate():
+    assert multivariate_Newton_polynomial_interpolation(Ptest2, 2 ** 31 - 1, verbose=True) == (t1 + 2 * t2) + 3 * (t1 * t2) + t3 ** 5
+
+
 def test_Thiele_rational_interpolation_trivial():
     assert Thiele_rational_interpolation(Rtest0, 2 ** 31 - 61, verbose=True) == 1