From 7cf701c7cfe9245461d5861842519e85108fcf65 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?I=C3=B1aki=20=C3=9Acar?= <iucar@fedoraproject.org>
Date: Fri, 17 Jul 2026 14:37:54 +0200
Subject: [PATCH 1/9] compute Taylor coefficients around an arbitrary constant
 in preallocated memory

---
 R/RcppExports.R     |  4 +--
 src/RcppExports.cpp |  9 ++++---
 src/nn2poly.cpp     |  6 +++--
 src/taylor.h        | 64 +++++++++++++++++++++------------------------
 4 files changed, 41 insertions(+), 42 deletions(-)

diff --git a/R/RcppExports.R b/R/RcppExports.R
index 65517b1..773bb23 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -5,8 +5,8 @@ obtain_taylor_vector <- function(taylor_orders, af_string_list) {
     .Call(`_nn2poly_obtain_taylor_vector`, taylor_orders, af_string_list)
 }
 
-obtain_derivatives_list <- function(taylor_orders, af_string_list) {
-    .Call(`_nn2poly_obtain_derivatives_list`, taylor_orders, af_string_list)
+obtain_derivatives_list <- function(taylor_orders, af_string_list, a = 0.0) {
+    .Call(`_nn2poly_obtain_derivatives_list`, taylor_orders, af_string_list, a)
 }
 
 alg_non_linear <- function(coeffs_input, labels_input, labels_output, previous_order, q_layer, g) {
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index 15f9a9c..48b5918 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -25,14 +25,15 @@ BEGIN_RCPP
 END_RCPP
 }
 // obtain_derivatives_list
-CoeffsList obtain_derivatives_list(const Term& taylor_orders, const Functions& af_string_list);
-RcppExport SEXP _nn2poly_obtain_derivatives_list(SEXP taylor_ordersSEXP, SEXP af_string_listSEXP) {
+CoeffsList obtain_derivatives_list(const Term& taylor_orders, const Functions& af_string_list, double a);
+RcppExport SEXP _nn2poly_obtain_derivatives_list(SEXP taylor_ordersSEXP, SEXP af_string_listSEXP, SEXP aSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
     Rcpp::traits::input_parameter< const Term& >::type taylor_orders(taylor_ordersSEXP);
     Rcpp::traits::input_parameter< const Functions& >::type af_string_list(af_string_listSEXP);
-    rcpp_result_gen = Rcpp::wrap(obtain_derivatives_list(taylor_orders, af_string_list));
+    Rcpp::traits::input_parameter< double >::type a(aSEXP);
+    rcpp_result_gen = Rcpp::wrap(obtain_derivatives_list(taylor_orders, af_string_list, a));
     return rcpp_result_gen;
 END_RCPP
 }
@@ -82,7 +83,7 @@ END_RCPP
 
 static const R_CallMethodDef CallEntries[] = {
     {"_nn2poly_obtain_taylor_vector", (DL_FUNC) &_nn2poly_obtain_taylor_vector, 2},
-    {"_nn2poly_obtain_derivatives_list", (DL_FUNC) &_nn2poly_obtain_derivatives_list, 2},
+    {"_nn2poly_obtain_derivatives_list", (DL_FUNC) &_nn2poly_obtain_derivatives_list, 3},
     {"_nn2poly_alg_non_linear", (DL_FUNC) &_nn2poly_alg_non_linear, 6},
     {"_nn2poly_nn2poly_algorithm", (DL_FUNC) &_nn2poly_nn2poly_algorithm, 5},
     {"_nn2poly_combinations_with_repetition", (DL_FUNC) &_nn2poly_combinations_with_repetition, 2},
diff --git a/src/nn2poly.cpp b/src/nn2poly.cpp
index 3634d55..a5b0abf 100644
--- a/src/nn2poly.cpp
+++ b/src/nn2poly.cpp
@@ -24,7 +24,8 @@ Term obtain_taylor_vector(const Term& taylor_orders,
 
 // [[Rcpp::export]]
 CoeffsList obtain_derivatives_list(const Term& taylor_orders,
-                                   const Functions& af_string_list) {
+                                   const Functions& af_string_list,
+                                   double a = 0.0) {
   if (taylor_orders.size() != af_string_list.size())
     throw std::invalid_argument(
       "`taylor_orders` length does not match provided number of layers");
@@ -35,7 +36,8 @@ CoeffsList obtain_derivatives_list(const Term& taylor_orders,
       throw std::invalid_argument("`taylor_orders` must be non-negative");
     // Obtain the vector with the derivatives of the activation function up to
     // the given degree centered at 0
-    out[i] = coeffs_taylor(af_string_list[i], taylor_orders[i]);
+    out[i].resize(static_cast<size_t>(taylor_orders[i] + 1));
+    coeffs_taylor(af_string_list[i], out[i], a);
   }
 
   return out;
diff --git a/src/taylor.h b/src/taylor.h
index 092da29..622d456 100644
--- a/src/taylor.h
+++ b/src/taylor.h
@@ -3,57 +3,53 @@
 
 #include "nn2poly_types.h"
 
-Coeffs coeffs_taylor_sigmoid(int order) {
-  Coeffs a(order + 1, 0.0);
-  a[0] = 0.5;
-  for (int n = 0; n < order; n++) {
+void coeffs_taylor_sigmoid(Coeffs& c, double a) {
+  c[0] = 1.0 / (1.0 + std::exp(-a));
+  for (size_t n = 0; n < c.size() - 1; n++) {
     double conv = 0.0;
-    for (int k = 0; k <= n; k++)
-      conv += a[k] * a[n - k];
-    a[n + 1] = (a[n] - conv) / static_cast<double>(n + 1);
+    for (size_t k = 0; k <= n; k++)
+      conv += c[k] * c[n - k];
+    c[n + 1] = (c[n] - conv) / static_cast<double>(n + 1);
   }
-  return a;
 }
 
-Coeffs coeffs_taylor_tanh(int order) {
-  Coeffs a(order + 1, 0.0);
-  for (int n = 0; n < order; n++) {
+void coeffs_taylor_tanh(Coeffs& c, double a) {
+  c[0] = std::tanh(a);
+  for (size_t n = 0; n < c.size() - 1; n++) {
     double conv = 0.0;
-    for (int k = 0; k <= n; k++)
-      conv += a[k] * a[n - k];
+    for (size_t k = 0; k <= n; k++)
+      conv += c[k] * c[n - k];
     const double rhs = (n == 0 ? 1.0 : 0.0) - conv;
-    a[n + 1] = rhs / static_cast<double>(n + 1);
+    c[n + 1] = rhs / static_cast<double>(n + 1);
   }
-  return a;
 }
 
-Coeffs coeffs_taylor_linear(int order) {
-  Coeffs a(order + 1, 0.0);
-  if (order >= 1)
-    a[1] = 1.0;
-  return a;
+void coeffs_taylor_linear(Coeffs& c, double a) {
+  c[0] = a;
+  if (c.size() > 1) {
+    c[1] = 1.0;
+    std::fill(c.begin() + 2, c.end(), 0.0);
+  }
 }
 
-Coeffs coeffs_taylor_softplus(int order) {
-  Coeffs a(order + 1, 0.0);
-  a[0] = std::log(2.0);
-  if (order == 0)
-    return a;
-  Coeffs sig = coeffs_taylor_sigmoid(order - 1);
-  for (int n = 1; n <= order; n++)
-    a[n] = sig[n - 1] / static_cast<double>(n);
-  return a;
+void coeffs_taylor_softplus(Coeffs& c, double a) {
+  coeffs_taylor_sigmoid(c, a);
+  for (size_t n = c.size() - 1; n > 0; n--)
+    c[n] = c[n - 1] / static_cast<double>(n);
+  if (a > 0.0)
+    c[0] = a + std::log(1.0 + std::exp(-a));
+  else c[0] = std::log(1.0 + std::exp(a));
 }
 
-Coeffs coeffs_taylor(const std::string& af, int order) {
+void coeffs_taylor(const std::string& af, Coeffs& c, double a = 0.0) {
   if (af == "sigmoid")
-    return coeffs_taylor_sigmoid(order);
+    return coeffs_taylor_sigmoid(c, a);
   if (af == "tanh")
-    return coeffs_taylor_tanh(order);
+    return coeffs_taylor_tanh(c, a);
   if (af == "softplus")
-    return coeffs_taylor_softplus(order);
+    return coeffs_taylor_softplus(c, a);
   if (af == "linear")
-    return coeffs_taylor_linear(order);
+    return coeffs_taylor_linear(c, a);
   throw std::invalid_argument("function '" + af + "' not supported");
 }
 

From 07f0eb01001b22a194a11337b2abbd6587cb2ece Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?I=C3=B1aki=20=C3=9Acar?= <iucar@fedoraproject.org>
Date: Fri, 17 Jul 2026 14:40:38 +0200
Subject: [PATCH 2/9] update smoke test for an arbitrary center

---
 tests/testthat/test-smoke.R | 439 ++++++++++++++++++++++++++++--------
 1 file changed, 339 insertions(+), 100 deletions(-)

diff --git a/tests/testthat/test-smoke.R b/tests/testthat/test-smoke.R
index a7f5f5e..d0e1930 100644
--- a/tests/testthat/test-smoke.R
+++ b/tests/testthat/test-smoke.R
@@ -12,7 +12,7 @@ smoke_test_data <- function() {
 # This script computes every coefficient explicitly. It intentionally does
 # not use nn2poly, symbolic algebra, loops, or a generic polynomial product.
 
-smoke_test_manual <- function(data) {
+smoke_test_manual <- function(data, taylor_center=0) {
 
   # Network weights ---------------------------------------------------------
   # The first value in each weight vector is the bias.
@@ -35,55 +35,169 @@ smoke_test_manual <- function(data) {
   w2 <- c(bias = b0, x1 = b1, x2 = b2, x3 = b3)
   w3 <- c(bias = w30, c = w31, d = w32)
 
-  # Taylor approximation and polynomial truncation -------------------------
+  # Fifth-order Taylor expansion of tanh at z = taylor_center ---------------
   #
-  #   T5(z) = z - z^3 / 3 + 2*z^5 / 15.
+  # tanh(z) is approximated by
   #
-  # Only monomials in x1, x2, x3 of total degree at most 4 are retained.
-  # The factors 3 and 6 below come from the expansion of p^3. The factors
-  # 5, 10, 20, 30 and 60 come from the expansion of p^5.
-
-  # Hidden neuron 1: c = T5(a), retaining degree <= 4 ----------------------
-
-  c0 <- a0 - a0^3 / 3 + (2 * a0^5) / 15
-
-  c1 <- a1 - (3 * a0^2 * a1) / 3 + (2 * 5 * a0^4 * a1) / 15
-  c2 <- a2 - (3 * a0^2 * a2) / 3 + (2 * 5 * a0^4 * a2) / 15
-  c3 <- a3 - (3 * a0^2 * a3) / 3 + (2 * 5 * a0^4 * a3) / 15
-
-  c11 <- -(3 * a0 * a1^2) / 3 + (2 * 10 * a0^3 * a1^2) / 15
-  c12 <- -(6 * a0 * a1 * a2) / 3 + (2 * 20 * a0^3 * a1 * a2) / 15
-  c13 <- -(6 * a0 * a1 * a3) / 3 + (2 * 20 * a0^3 * a1 * a3) / 15
-  c22 <- -(3 * a0 * a2^2) / 3 + (2 * 10 * a0^3 * a2^2) / 15
-  c23 <- -(6 * a0 * a2 * a3) / 3 + (2 * 20 * a0^3 * a2 * a3) / 15
-  c33 <- -(3 * a0 * a3^2) / 3 + (2 * 10 * a0^3 * a3^2) / 15
-
-  c111 <- -a1^3 / 3 + (2 * 10 * a0^2 * a1^3) / 15
-  c112 <- -(3 * a1^2 * a2) / 3 + (2 * 30 * a0^2 * a1^2 * a2) / 15
-  c113 <- -(3 * a1^2 * a3) / 3 + (2 * 30 * a0^2 * a1^2 * a3) / 15
-  c122 <- -(3 * a1 * a2^2) / 3 + (2 * 30 * a0^2 * a1 * a2^2) / 15
-  c123 <- -(6 * a1 * a2 * a3) / 3 + (2 * 60 * a0^2 * a1 * a2 * a3) / 15
-  c133 <- -(3 * a1 * a3^2) / 3 + (2 * 30 * a0^2 * a1 * a3^2) / 15
-  c222 <- -a2^3 / 3 + (2 * 10 * a0^2 * a2^3) / 15
-  c223 <- -(3 * a2^2 * a3) / 3 + (2 * 30 * a0^2 * a2^2 * a3) / 15
-  c233 <- -(3 * a2 * a3^2) / 3 + (2 * 30 * a0^2 * a2 * a3^2) / 15
-  c333 <- -a3^3 / 3 + (2 * 10 * a0^2 * a3^3) / 15
-
-  c1111 <- (2 * 5 * a0 * a1^4) / 15
-  c1112 <- (2 * 20 * a0 * a1^3 * a2) / 15
-  c1113 <- (2 * 20 * a0 * a1^3 * a3) / 15
-  c1122 <- (2 * 30 * a0 * a1^2 * a2^2) / 15
-  c1123 <- (2 * 60 * a0 * a1^2 * a2 * a3) / 15
-  c1133 <- (2 * 30 * a0 * a1^2 * a3^2) / 15
-  c1222 <- (2 * 20 * a0 * a1 * a2^3) / 15
-  c1223 <- (2 * 60 * a0 * a1 * a2^2 * a3) / 15
-  c1233 <- (2 * 60 * a0 * a1 * a2 * a3^2) / 15
-  c1333 <- (2 * 20 * a0 * a1 * a3^3) / 15
-  c2222 <- (2 * 5 * a0 * a2^4) / 15
-  c2223 <- (2 * 20 * a0 * a2^3 * a3) / 15
-  c2233 <- (2 * 30 * a0 * a2^2 * a3^2) / 15
-  c2333 <- (2 * 20 * a0 * a2 * a3^3) / 15
-  c3333 <- (2 * 5 * a0 * a3^4) / 15
+  #   tau0 + tau1*q + tau2*q^2 + tau3*q^3 + tau4*q^4 + tau5*q^5,
+  #
+  # where q = z - taylor_center and tau_k = tanh^(k)(taylor_center) / k!.
+
+  tanh_center <- tanh(taylor_center)
+
+  tau0 <- tanh_center
+  tau1 <- 1 - tanh_center^2
+  tau2 <- (-2 * tanh_center + 2 * tanh_center^3) / 2
+  tau3 <- (-2 + 8 * tanh_center^2 - 6 * tanh_center^4) / 6
+  tau4 <- (
+    16 * tanh_center - 40 * tanh_center^3 + 24 * tanh_center^5
+  ) / 24
+  tau5 <- (
+    16 - 136 * tanh_center^2 + 240 * tanh_center^4 -
+      120 * tanh_center^6
+  ) / 120
+
+  taylor_coefficients <- c(
+    order_0 = tau0,
+    order_1 = tau1,
+    order_2 = tau2,
+    order_3 = tau3,
+    order_4 = tau4,
+    order_5 = tau5
+  )
+
+  # Shifted intercepts used only in q = z - taylor_center.
+
+  a0_centered <- a0 - taylor_center
+  b0_centered <- b0 - taylor_center
+
+  # Hidden neuron 1: c -----------------------------------------------------
+
+  # Degree 0.
+
+  c0 <- tau0 +
+    tau1 * a0_centered +
+    tau2 * a0_centered^2 +
+    tau3 * a0_centered^3 +
+    tau4 * a0_centered^4 +
+    tau5 * a0_centered^5
+
+  # Degree 1. Multiplicities in q^1,...,q^5: 1, 2, 3, 4, 5.
+
+  c1 <- tau1 * a1 +
+    tau2 * (2 * a0_centered * a1) +
+    tau3 * (3 * a0_centered^2 * a1) +
+    tau4 * (4 * a0_centered^3 * a1) +
+    tau5 * (5 * a0_centered^4 * a1)
+  c2 <- tau1 * a2 +
+    tau2 * (2 * a0_centered * a2) +
+    tau3 * (3 * a0_centered^2 * a2) +
+    tau4 * (4 * a0_centered^3 * a2) +
+    tau5 * (5 * a0_centered^4 * a2)
+  c3 <- tau1 * a3 +
+    tau2 * (2 * a0_centered * a3) +
+    tau3 * (3 * a0_centered^2 * a3) +
+    tau4 * (4 * a0_centered^3 * a3) +
+    tau5 * (5 * a0_centered^4 * a3)
+
+  # Degree 2.
+  # Repeated variable (ii): 1, 3, 6, 10 in q^2,...,q^5.
+  # Distinct variables (ij): 2, 6, 12, 20 in q^2,...,q^5.
+
+  c11 <- tau2 * a1^2 +
+    tau3 * (3 * a0_centered * a1^2) +
+    tau4 * (6 * a0_centered^2 * a1^2) +
+    tau5 * (10 * a0_centered^3 * a1^2)
+  c12 <- tau2 * (2 * a1 * a2) +
+    tau3 * (6 * a0_centered * a1 * a2) +
+    tau4 * (12 * a0_centered^2 * a1 * a2) +
+    tau5 * (20 * a0_centered^3 * a1 * a2)
+  c13 <- tau2 * (2 * a1 * a3) +
+    tau3 * (6 * a0_centered * a1 * a3) +
+    tau4 * (12 * a0_centered^2 * a1 * a3) +
+    tau5 * (20 * a0_centered^3 * a1 * a3)
+  c22 <- tau2 * a2^2 +
+    tau3 * (3 * a0_centered * a2^2) +
+    tau4 * (6 * a0_centered^2 * a2^2) +
+    tau5 * (10 * a0_centered^3 * a2^2)
+  c23 <- tau2 * (2 * a2 * a3) +
+    tau3 * (6 * a0_centered * a2 * a3) +
+    tau4 * (12 * a0_centered^2 * a2 * a3) +
+    tau5 * (20 * a0_centered^3 * a2 * a3)
+  c33 <- tau2 * a3^2 +
+    tau3 * (3 * a0_centered * a3^2) +
+    tau4 * (6 * a0_centered^2 * a3^2) +
+    tau5 * (10 * a0_centered^3 * a3^2)
+
+  # Degree 3.
+  # Patterns iii, iij and ijk have multiplicities (1,4,10), (3,12,30)
+  # and (6,24,60), respectively, in q^3, q^4 and q^5.
+
+  c111 <- tau3 * a1^3 +
+    tau4 * (4 * a0_centered * a1^3) +
+    tau5 * (10 * a0_centered^2 * a1^3)
+  c112 <- tau3 * (3 * a1^2 * a2) +
+    tau4 * (12 * a0_centered * a1^2 * a2) +
+    tau5 * (30 * a0_centered^2 * a1^2 * a2)
+  c113 <- tau3 * (3 * a1^2 * a3) +
+    tau4 * (12 * a0_centered * a1^2 * a3) +
+    tau5 * (30 * a0_centered^2 * a1^2 * a3)
+  c122 <- tau3 * (3 * a1 * a2^2) +
+    tau4 * (12 * a0_centered * a1 * a2^2) +
+    tau5 * (30 * a0_centered^2 * a1 * a2^2)
+  c123 <- tau3 * (6 * a1 * a2 * a3) +
+    tau4 * (24 * a0_centered * a1 * a2 * a3) +
+    tau5 * (60 * a0_centered^2 * a1 * a2 * a3)
+  c133 <- tau3 * (3 * a1 * a3^2) +
+    tau4 * (12 * a0_centered * a1 * a3^2) +
+    tau5 * (30 * a0_centered^2 * a1 * a3^2)
+  c222 <- tau3 * a2^3 +
+    tau4 * (4 * a0_centered * a2^3) +
+    tau5 * (10 * a0_centered^2 * a2^3)
+  c223 <- tau3 * (3 * a2^2 * a3) +
+    tau4 * (12 * a0_centered * a2^2 * a3) +
+    tau5 * (30 * a0_centered^2 * a2^2 * a3)
+  c233 <- tau3 * (3 * a2 * a3^2) +
+    tau4 * (12 * a0_centered * a2 * a3^2) +
+    tau5 * (30 * a0_centered^2 * a2 * a3^2)
+  c333 <- tau3 * a3^3 +
+    tau4 * (4 * a0_centered * a3^3) +
+    tau5 * (10 * a0_centered^2 * a3^3)
+
+  # Degree 4.
+  # Patterns iiii, iiij, iijj and iijk have multiplicities (1,5), (4,20),
+  # (6,30) and (12,60), respectively, in q^4 and q^5.
+
+  c1111 <- tau4 * a1^4 +
+    tau5 * (5 * a0_centered * a1^4)
+  c1112 <- tau4 * (4 * a1^3 * a2) +
+    tau5 * (20 * a0_centered * a1^3 * a2)
+  c1113 <- tau4 * (4 * a1^3 * a3) +
+    tau5 * (20 * a0_centered * a1^3 * a3)
+  c1122 <- tau4 * (6 * a1^2 * a2^2) +
+    tau5 * (30 * a0_centered * a1^2 * a2^2)
+  c1123 <- tau4 * (12 * a1^2 * a2 * a3) +
+    tau5 * (60 * a0_centered * a1^2 * a2 * a3)
+  c1133 <- tau4 * (6 * a1^2 * a3^2) +
+    tau5 * (30 * a0_centered * a1^2 * a3^2)
+  c1222 <- tau4 * (4 * a1 * a2^3) +
+    tau5 * (20 * a0_centered * a1 * a2^3)
+  c1223 <- tau4 * (12 * a1 * a2^2 * a3) +
+    tau5 * (60 * a0_centered * a1 * a2^2 * a3)
+  c1233 <- tau4 * (12 * a1 * a2 * a3^2) +
+    tau5 * (60 * a0_centered * a1 * a2 * a3^2)
+  c1333 <- tau4 * (4 * a1 * a3^3) +
+    tau5 * (20 * a0_centered * a1 * a3^3)
+  c2222 <- tau4 * a2^4 +
+    tau5 * (5 * a0_centered * a2^4)
+  c2223 <- tau4 * (4 * a2^3 * a3) +
+    tau5 * (20 * a0_centered * a2^3 * a3)
+  c2233 <- tau4 * (6 * a2^2 * a3^2) +
+    tau5 * (30 * a0_centered * a2^2 * a3^2)
+  c2333 <- tau4 * (4 * a2 * a3^3) +
+    tau5 * (20 * a0_centered * a2 * a3^3)
+  c3333 <- tau4 * a3^4 +
+    tau5 * (5 * a0_centered * a3^4)
 
   coefficients_c <- c(
     "0" = c0,
@@ -101,47 +215,127 @@ smoke_test_manual <- function(data) {
     "2233" = c2233, "2333" = c2333, "3333" = c3333
   )
 
-  # Hidden neuron 2: d = T5(b), retaining degree <= 4 ----------------------
-
-  d0 <- b0 - b0^3 / 3 + (2 * b0^5) / 15
-
-  d1 <- b1 - (3 * b0^2 * b1) / 3 + (2 * 5 * b0^4 * b1) / 15
-  d2 <- b2 - (3 * b0^2 * b2) / 3 + (2 * 5 * b0^4 * b2) / 15
-  d3 <- b3 - (3 * b0^2 * b3) / 3 + (2 * 5 * b0^4 * b3) / 15
-
-  d11 <- -(3 * b0 * b1^2) / 3 + (2 * 10 * b0^3 * b1^2) / 15
-  d12 <- -(6 * b0 * b1 * b2) / 3 + (2 * 20 * b0^3 * b1 * b2) / 15
-  d13 <- -(6 * b0 * b1 * b3) / 3 + (2 * 20 * b0^3 * b1 * b3) / 15
-  d22 <- -(3 * b0 * b2^2) / 3 + (2 * 10 * b0^3 * b2^2) / 15
-  d23 <- -(6 * b0 * b2 * b3) / 3 + (2 * 20 * b0^3 * b2 * b3) / 15
-  d33 <- -(3 * b0 * b3^2) / 3 + (2 * 10 * b0^3 * b3^2) / 15
-
-  d111 <- -b1^3 / 3 + (2 * 10 * b0^2 * b1^3) / 15
-  d112 <- -(3 * b1^2 * b2) / 3 + (2 * 30 * b0^2 * b1^2 * b2) / 15
-  d113 <- -(3 * b1^2 * b3) / 3 + (2 * 30 * b0^2 * b1^2 * b3) / 15
-  d122 <- -(3 * b1 * b2^2) / 3 + (2 * 30 * b0^2 * b1 * b2^2) / 15
-  d123 <- -(6 * b1 * b2 * b3) / 3 + (2 * 60 * b0^2 * b1 * b2 * b3) / 15
-  d133 <- -(3 * b1 * b3^2) / 3 + (2 * 30 * b0^2 * b1 * b3^2) / 15
-  d222 <- -b2^3 / 3 + (2 * 10 * b0^2 * b2^3) / 15
-  d223 <- -(3 * b2^2 * b3) / 3 + (2 * 30 * b0^2 * b2^2 * b3) / 15
-  d233 <- -(3 * b2 * b3^2) / 3 + (2 * 30 * b0^2 * b2 * b3^2) / 15
-  d333 <- -b3^3 / 3 + (2 * 10 * b0^2 * b3^3) / 15
-
-  d1111 <- (2 * 5 * b0 * b1^4) / 15
-  d1112 <- (2 * 20 * b0 * b1^3 * b2) / 15
-  d1113 <- (2 * 20 * b0 * b1^3 * b3) / 15
-  d1122 <- (2 * 30 * b0 * b1^2 * b2^2) / 15
-  d1123 <- (2 * 60 * b0 * b1^2 * b2 * b3) / 15
-  d1133 <- (2 * 30 * b0 * b1^2 * b3^2) / 15
-  d1222 <- (2 * 20 * b0 * b1 * b2^3) / 15
-  d1223 <- (2 * 60 * b0 * b1 * b2^2 * b3) / 15
-  d1233 <- (2 * 60 * b0 * b1 * b2 * b3^2) / 15
-  d1333 <- (2 * 20 * b0 * b1 * b3^3) / 15
-  d2222 <- (2 * 5 * b0 * b2^4) / 15
-  d2223 <- (2 * 20 * b0 * b2^3 * b3) / 15
-  d2233 <- (2 * 30 * b0 * b2^2 * b3^2) / 15
-  d2333 <- (2 * 20 * b0 * b2 * b3^3) / 15
-  d3333 <- (2 * 5 * b0 * b3^4) / 15
+  # Hidden neuron 2: d -----------------------------------------------------
+
+  # Degree 0.
+
+  d0 <- tau0 +
+    tau1 * b0_centered +
+    tau2 * b0_centered^2 +
+    tau3 * b0_centered^3 +
+    tau4 * b0_centered^4 +
+    tau5 * b0_centered^5
+
+  # Degree 1. Multiplicities in q^1,...,q^5: 1, 2, 3, 4, 5.
+
+  d1 <- tau1 * b1 +
+    tau2 * (2 * b0_centered * b1) +
+    tau3 * (3 * b0_centered^2 * b1) +
+    tau4 * (4 * b0_centered^3 * b1) +
+    tau5 * (5 * b0_centered^4 * b1)
+  d2 <- tau1 * b2 +
+    tau2 * (2 * b0_centered * b2) +
+    tau3 * (3 * b0_centered^2 * b2) +
+    tau4 * (4 * b0_centered^3 * b2) +
+    tau5 * (5 * b0_centered^4 * b2)
+  d3 <- tau1 * b3 +
+    tau2 * (2 * b0_centered * b3) +
+    tau3 * (3 * b0_centered^2 * b3) +
+    tau4 * (4 * b0_centered^3 * b3) +
+    tau5 * (5 * b0_centered^4 * b3)
+
+  # Degree 2.
+
+  d11 <- tau2 * b1^2 +
+    tau3 * (3 * b0_centered * b1^2) +
+    tau4 * (6 * b0_centered^2 * b1^2) +
+    tau5 * (10 * b0_centered^3 * b1^2)
+  d12 <- tau2 * (2 * b1 * b2) +
+    tau3 * (6 * b0_centered * b1 * b2) +
+    tau4 * (12 * b0_centered^2 * b1 * b2) +
+    tau5 * (20 * b0_centered^3 * b1 * b2)
+  d13 <- tau2 * (2 * b1 * b3) +
+    tau3 * (6 * b0_centered * b1 * b3) +
+    tau4 * (12 * b0_centered^2 * b1 * b3) +
+    tau5 * (20 * b0_centered^3 * b1 * b3)
+  d22 <- tau2 * b2^2 +
+    tau3 * (3 * b0_centered * b2^2) +
+    tau4 * (6 * b0_centered^2 * b2^2) +
+    tau5 * (10 * b0_centered^3 * b2^2)
+  d23 <- tau2 * (2 * b2 * b3) +
+    tau3 * (6 * b0_centered * b2 * b3) +
+    tau4 * (12 * b0_centered^2 * b2 * b3) +
+    tau5 * (20 * b0_centered^3 * b2 * b3)
+  d33 <- tau2 * b3^2 +
+    tau3 * (3 * b0_centered * b3^2) +
+    tau4 * (6 * b0_centered^2 * b3^2) +
+    tau5 * (10 * b0_centered^3 * b3^2)
+
+  # Degree 3.
+
+  d111 <- tau3 * b1^3 +
+    tau4 * (4 * b0_centered * b1^3) +
+    tau5 * (10 * b0_centered^2 * b1^3)
+  d112 <- tau3 * (3 * b1^2 * b2) +
+    tau4 * (12 * b0_centered * b1^2 * b2) +
+    tau5 * (30 * b0_centered^2 * b1^2 * b2)
+  d113 <- tau3 * (3 * b1^2 * b3) +
+    tau4 * (12 * b0_centered * b1^2 * b3) +
+    tau5 * (30 * b0_centered^2 * b1^2 * b3)
+  d122 <- tau3 * (3 * b1 * b2^2) +
+    tau4 * (12 * b0_centered * b1 * b2^2) +
+    tau5 * (30 * b0_centered^2 * b1 * b2^2)
+  d123 <- tau3 * (6 * b1 * b2 * b3) +
+    tau4 * (24 * b0_centered * b1 * b2 * b3) +
+    tau5 * (60 * b0_centered^2 * b1 * b2 * b3)
+  d133 <- tau3 * (3 * b1 * b3^2) +
+    tau4 * (12 * b0_centered * b1 * b3^2) +
+    tau5 * (30 * b0_centered^2 * b1 * b3^2)
+  d222 <- tau3 * b2^3 +
+    tau4 * (4 * b0_centered * b2^3) +
+    tau5 * (10 * b0_centered^2 * b2^3)
+  d223 <- tau3 * (3 * b2^2 * b3) +
+    tau4 * (12 * b0_centered * b2^2 * b3) +
+    tau5 * (30 * b0_centered^2 * b2^2 * b3)
+  d233 <- tau3 * (3 * b2 * b3^2) +
+    tau4 * (12 * b0_centered * b2 * b3^2) +
+    tau5 * (30 * b0_centered^2 * b2 * b3^2)
+  d333 <- tau3 * b3^3 +
+    tau4 * (4 * b0_centered * b3^3) +
+    tau5 * (10 * b0_centered^2 * b3^3)
+
+  # Degree 4.
+
+  d1111 <- tau4 * b1^4 +
+    tau5 * (5 * b0_centered * b1^4)
+  d1112 <- tau4 * (4 * b1^3 * b2) +
+    tau5 * (20 * b0_centered * b1^3 * b2)
+  d1113 <- tau4 * (4 * b1^3 * b3) +
+    tau5 * (20 * b0_centered * b1^3 * b3)
+  d1122 <- tau4 * (6 * b1^2 * b2^2) +
+    tau5 * (30 * b0_centered * b1^2 * b2^2)
+  d1123 <- tau4 * (12 * b1^2 * b2 * b3) +
+    tau5 * (60 * b0_centered * b1^2 * b2 * b3)
+  d1133 <- tau4 * (6 * b1^2 * b3^2) +
+    tau5 * (30 * b0_centered * b1^2 * b3^2)
+  d1222 <- tau4 * (4 * b1 * b2^3) +
+    tau5 * (20 * b0_centered * b1 * b2^3)
+  d1223 <- tau4 * (12 * b1 * b2^2 * b3) +
+    tau5 * (60 * b0_centered * b1 * b2^2 * b3)
+  d1233 <- tau4 * (12 * b1 * b2 * b3^2) +
+    tau5 * (60 * b0_centered * b1 * b2 * b3^2)
+  d1333 <- tau4 * (4 * b1 * b3^3) +
+    tau5 * (20 * b0_centered * b1 * b3^3)
+  d2222 <- tau4 * b2^4 +
+    tau5 * (5 * b0_centered * b2^4)
+  d2223 <- tau4 * (4 * b2^3 * b3) +
+    tau5 * (20 * b0_centered * b2^3 * b3)
+  d2233 <- tau4 * (6 * b2^2 * b3^2) +
+    tau5 * (30 * b0_centered * b2^2 * b3^2)
+  d2333 <- tau4 * (4 * b2 * b3^3) +
+    tau5 * (20 * b0_centered * b2 * b3^3)
+  d3333 <- tau4 * b3^4 +
+    tau5 * (5 * b0_centered * b3^4)
 
   coefficients_d <- c(
     "0" = d0,
@@ -159,7 +353,7 @@ smoke_test_manual <- function(data) {
     "2233" = d2233, "2333" = d2333, "3333" = d3333
   )
 
-  # Linear output: e = w30 + w31*c + w32*d -------------------------------
+  # Linear output: e -------------------------------------------------------
 
   e0 <- w30 + w31 * c0 + w32 * d0
 
@@ -217,27 +411,72 @@ smoke_test_manual <- function(data) {
     "2233" = e2233, "2333" = e2333, "3333" = e3333
   )
 
-  # Final results -----------------------------------------------------------
+  # NN2Poly keep_layers structure -----------------------------------------
+  # Labels are integer vectors, matching the representation used by nn2poly.
+
+  input_labels <- list(0L, 1L, 2L, 3L)
+
+  polynomial_labels <- list(
+    0L,
+    1L, 2L, 3L,
+    c(1L, 1L), c(1L, 2L), c(1L, 3L),
+    c(2L, 2L), c(2L, 3L), c(3L, 3L),
+    c(1L, 1L, 1L), c(1L, 1L, 2L), c(1L, 1L, 3L),
+    c(1L, 2L, 2L), c(1L, 2L, 3L), c(1L, 3L, 3L),
+    c(2L, 2L, 2L), c(2L, 2L, 3L), c(2L, 3L, 3L),
+    c(3L, 3L, 3L),
+    c(1L, 1L, 1L, 1L), c(1L, 1L, 1L, 2L),
+    c(1L, 1L, 1L, 3L), c(1L, 1L, 2L, 2L),
+    c(1L, 1L, 2L, 3L), c(1L, 1L, 3L, 3L),
+    c(1L, 2L, 2L, 2L), c(1L, 2L, 2L, 3L),
+    c(1L, 2L, 3L, 3L), c(1L, 3L, 3L, 3L),
+    c(2L, 2L, 2L, 2L), c(2L, 2L, 2L, 3L),
+    c(2L, 2L, 3L, 3L), c(2L, 3L, 3L, 3L),
+    c(3L, 3L, 3L, 3L)
+  )
+
+  layer_1_input_values <- matrix(
+    c(a0, a1, a2, a3, b0, b1, b2, b3),
+    nrow = 4L,
+    ncol = 2L,
+    dimnames = NULL
+  )
+
+  layer_1_output_values <- matrix(
+    c(unname(coefficients_c), unname(coefficients_d)),
+    nrow = length(polynomial_labels),
+    ncol = 2L,
+    dimnames = NULL
+  )
+
+  layer_2_values <- matrix(
+    unname(coefficients_e),
+    nrow = length(polynomial_labels),
+    ncol = 1L,
+    dimnames = NULL
+  )
+
+  layer_2_polynomial <- list(
+    labels = polynomial_labels,
+    values = layer_2_values
+  )
 
   obj <- list(
     layer_1 = list(
       input = list(
-        labels = list(0, 1, 2, 3),
-        values = data$tanh
+        labels = input_labels,
+        values = layer_1_input_values
       ),
       output = list(
-        labels = sapply(strsplit(names(coefficients_c), ""), as.numeric),
-        values = cbind(unname(coefficients_c), unname(coefficients_d))
+        labels = polynomial_labels,
+        values = layer_1_output_values
       )
     ),
     layer_2 = list(
-      input = list(
-        labels = sapply(strsplit(names(coefficients_e), ""), as.numeric),
-        values = cbind(unname(coefficients_e))
-      )
+      input = layer_2_polynomial,
+      output = layer_2_polynomial
     )
   )
-  obj$layer_2$output <- obj$layer_2$input
   class(obj) <- "nn2poly"
   obj
 }

From 187df2b425c10be13b8b2fdaf0b0d74671400bbe Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?I=C3=B1aki=20=C3=9Acar?= <iucar@fedoraproject.org>
Date: Fri, 17 Jul 2026 14:41:49 +0200
Subject: [PATCH 3/9] implement and expose Taylor center

---
 R/RcppExports.R             |  4 ++--
 R/nn2poly.R                 |  7 ++++++-
 man/nn2poly.Rd              | 12 +++++++++++-
 src/RcppExports.cpp         |  9 +++++----
 src/linalg_arma.h           |  4 ++++
 src/linalg_eigen.h          |  4 ++++
 src/linalg_torch.h          |  4 ++++
 src/nn2poly.cpp             |  7 +++++--
 tests/testthat/test-smoke.R |  5 +++++
 9 files changed, 46 insertions(+), 10 deletions(-)

diff --git a/R/RcppExports.R b/R/RcppExports.R
index 773bb23..1c5caf6 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -13,8 +13,8 @@ alg_non_linear <- function(coeffs_input, labels_input, labels_output, previous_o
     .Call(`_nn2poly_alg_non_linear`, coeffs_input, labels_input, labels_output, previous_order, q_layer, g)
 }
 
-nn2poly_algorithm <- function(layers, af_list, max_order, keep_layers, taylor_orders) {
-    .Call(`_nn2poly_nn2poly_algorithm`, layers, af_list, max_order, keep_layers, taylor_orders)
+nn2poly_algorithm <- function(layers, af_list, max_order, keep_layers, taylor_orders, a = 0.0) {
+    .Call(`_nn2poly_nn2poly_algorithm`, layers, af_list, max_order, keep_layers, taylor_orders, a)
 }
 
 combinations_with_repetition <- function(n, k) {
diff --git a/R/nn2poly.R b/R/nn2poly.R
index b92a54e..e260d3f 100644
--- a/R/nn2poly.R
+++ b/R/nn2poly.R
@@ -45,6 +45,9 @@ NULL
 #  layers. If a vector is used, each value corresponds to the Taylor order used
 #' at each layer activation function. Default set to \code{8}.
 #'
+#' @param a \code{numeric} value that sets the point at which the Taylor expansion
+#' is computed. Default set to \code{0}.
+#'
 #' @param ... Ignored.
 #'
 #' @return Returns an object of class `nn2poly`.
@@ -114,6 +117,7 @@ nn2poly <- function(object,
                     max_order = 2,
                     keep_layers = FALSE,
                     taylor_orders = 8,
+                    a = 0.0,
                     ...) {
   UseMethod("nn2poly")
 }
@@ -123,9 +127,10 @@ nn2poly.list <- function(object,
                          max_order = 2,
                          keep_layers = FALSE,
                          taylor_orders = 8,
+                         a = 0.0,
                          ...) {
   result <- nn2poly_algorithm(object, names(object),
-                              max_order, keep_layers, taylor_orders)
+                              max_order, keep_layers, taylor_orders, a)
   class(result) <- "nn2poly"
   result
 }
diff --git a/man/nn2poly.Rd b/man/nn2poly.Rd
index 776c105..cf4c93f 100644
--- a/man/nn2poly.Rd
+++ b/man/nn2poly.Rd
@@ -4,7 +4,14 @@
 \alias{nn2poly}
 \title{Obtain polynomial representation}
 \usage{
-nn2poly(object, max_order = 2, keep_layers = FALSE, taylor_orders = 8, ...)
+nn2poly(
+  object,
+  max_order = 2,
+  keep_layers = FALSE,
+  taylor_orders = 8,
+  a = 0,
+  ...
+)
 }
 \arguments{
 \item{object}{An object for which the computation of the NN2Poly algorithm is
@@ -43,6 +50,9 @@ degree at which Taylor expansion is truncated at each layer. If a single
 value is used, that value is set for each non linear layer and 1 for linear
 at each layer activation function. Default set to \code{8}.}
 
+\item{a}{\code{numeric} value that sets the point at which the Taylor expansion
+is computed. Default set to \code{0}.}
+
 \item{...}{Ignored.}
 }
 \value{
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index 48b5918..1227a3a 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -54,8 +54,8 @@ BEGIN_RCPP
 END_RCPP
 }
 // nn2poly_algorithm
-List nn2poly_algorithm(const Layers& layers, const Functions& af_list, int max_order, bool keep_layers, const Term& taylor_orders);
-RcppExport SEXP _nn2poly_nn2poly_algorithm(SEXP layersSEXP, SEXP af_listSEXP, SEXP max_orderSEXP, SEXP keep_layersSEXP, SEXP taylor_ordersSEXP) {
+List nn2poly_algorithm(const Layers& layers, const Functions& af_list, int max_order, bool keep_layers, const Term& taylor_orders, double a);
+RcppExport SEXP _nn2poly_nn2poly_algorithm(SEXP layersSEXP, SEXP af_listSEXP, SEXP max_orderSEXP, SEXP keep_layersSEXP, SEXP taylor_ordersSEXP, SEXP aSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
@@ -64,7 +64,8 @@ BEGIN_RCPP
     Rcpp::traits::input_parameter< int >::type max_order(max_orderSEXP);
     Rcpp::traits::input_parameter< bool >::type keep_layers(keep_layersSEXP);
     Rcpp::traits::input_parameter< const Term& >::type taylor_orders(taylor_ordersSEXP);
-    rcpp_result_gen = Rcpp::wrap(nn2poly_algorithm(layers, af_list, max_order, keep_layers, taylor_orders));
+    Rcpp::traits::input_parameter< double >::type a(aSEXP);
+    rcpp_result_gen = Rcpp::wrap(nn2poly_algorithm(layers, af_list, max_order, keep_layers, taylor_orders, a));
     return rcpp_result_gen;
 END_RCPP
 }
@@ -85,7 +86,7 @@ static const R_CallMethodDef CallEntries[] = {
     {"_nn2poly_obtain_taylor_vector", (DL_FUNC) &_nn2poly_obtain_taylor_vector, 2},
     {"_nn2poly_obtain_derivatives_list", (DL_FUNC) &_nn2poly_obtain_derivatives_list, 3},
     {"_nn2poly_alg_non_linear", (DL_FUNC) &_nn2poly_alg_non_linear, 6},
-    {"_nn2poly_nn2poly_algorithm", (DL_FUNC) &_nn2poly_nn2poly_algorithm, 5},
+    {"_nn2poly_nn2poly_algorithm", (DL_FUNC) &_nn2poly_nn2poly_algorithm, 6},
     {"_nn2poly_combinations_with_repetition", (DL_FUNC) &_nn2poly_combinations_with_repetition, 2},
     {NULL, NULL, 0}
 };
diff --git a/src/linalg_arma.h b/src/linalg_arma.h
index 9247d13..7eb6714 100644
--- a/src/linalg_arma.h
+++ b/src/linalg_arma.h
@@ -28,6 +28,10 @@ inline Weights alg_linear(Weights& coeffs_list, const Weights& layer) {
   return arma::trans(layer) * arma::join_cols(intercept, coeffs_list);
 }
 
+inline void sub_scalar(Weights& mat, int i, double scalar) {
+  mat.col(i) -= scalar;
+}
+
 inline void add_partition(Weights& mat, int i, double scalar, const Vector& vec) {
   mat.col(i) += scalar * vec;
 }
diff --git a/src/linalg_eigen.h b/src/linalg_eigen.h
index 6d3ca39..f9d5b06 100644
--- a/src/linalg_eigen.h
+++ b/src/linalg_eigen.h
@@ -30,6 +30,10 @@ inline Weights alg_linear(const Weights& coeffs_list, const Weights& layer) {
   return layer.transpose() * joined;
 }
 
+inline void sub_scalar(Weights& mat, int i, double scalar) {
+  mat.col(i) -= scalar;
+}
+
 inline void add_partition(Weights& mat, int i, double scalar, const Vector& vec) {
   mat.col(i) += scalar * vec;
 }
diff --git a/src/linalg_torch.h b/src/linalg_torch.h
index 4175a4e..37c539a 100644
--- a/src/linalg_torch.h
+++ b/src/linalg_torch.h
@@ -31,6 +31,10 @@ inline Weights alg_linear(Weights& coeffs_list, const Weights& layer) {
   return torch::matmul(layer.t(), joined);
 }
 
+inline void sub_scalar(Weights& mat, int i, double scalar) {
+  mat.select(1, i).sub_(scalar);
+}
+
 inline void add_partition(Weights& mat, int i, double scalar, const Vector& vec) {
   mat.select(1, i).add_(vec * scalar);
 }
diff --git a/src/nn2poly.cpp b/src/nn2poly.cpp
index a5b0abf..8f2bff9 100644
--- a/src/nn2poly.cpp
+++ b/src/nn2poly.cpp
@@ -137,7 +137,7 @@ inline void check_weights_dimensions(const Layers& layers) {
 // [[Rcpp::export]]
 List nn2poly_algorithm(const Layers& layers, const Functions& af_list,
                        int max_order, bool keep_layers,
-                       const Term& taylor_orders) {
+                       const Term& taylor_orders, double a = 0.0) {
   if (layers.empty())
     throw std::invalid_argument("`layers` is empty");
   if (af_list.empty())
@@ -156,7 +156,7 @@ List nn2poly_algorithm(const Layers& layers, const Functions& af_list,
   // depending on the last layer being linear or not
   WeightsLists results(static_cast<size_t>(last_linear ? 2 * L - 1 : 2 * L));
   const Term taylor = obtain_taylor_vector(taylor_orders, af_list);
-  const CoeffsList af_dlist = obtain_derivatives_list(taylor, af_list);
+  const CoeffsList af_dlist = obtain_derivatives_list(taylor, af_list, a);
 
   // Starting point for the algorithm: Set weights as coefficients
   // of an order 1 polynomial.
@@ -224,6 +224,9 @@ List nn2poly_algorithm(const Layers& layers, const Functions& af_list,
 
     NN2POLY_DEBUG_LOG(2, "[layer", l, "]",
       DTAG(previous_order), DTAG(q_layer), DTAG(labels_list.size()));
+
+    // Subtract the center from the intercept
+    if (a) sub_scalar(coeffs_list, 0, a);
     // Apply non-linear algorithm
     coeffs_list = alg_non_linear_impl(
       coeffs_list,
diff --git a/tests/testthat/test-smoke.R b/tests/testthat/test-smoke.R
index d0e1930..7711365 100644
--- a/tests/testthat/test-smoke.R
+++ b/tests/testthat/test-smoke.R
@@ -483,7 +483,12 @@ smoke_test_manual <- function(data, taylor_center=0) {
 
 test_that("manual example passes", {
   data <- smoke_test_data()
+
   x <- smoke_test_manual(data)
   y <- nn2poly(data, max_order = 4, taylor = 5, keep_layers = TRUE)
   expect_equal(x, y)
+
+  x <- smoke_test_manual(data, 0.5)
+  y <- nn2poly(data, max_order = 4, taylor = 5, keep_layers = TRUE, a = 0.5)
+  expect_equal(x, y)
 })

From d39cac05c5ebd6e3d23ff099ea3b7731e31fc63e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?I=C3=B1aki=20=C3=9Acar?= <iucar@fedoraproject.org>
Date: Fri, 17 Jul 2026 14:57:02 +0200
Subject: [PATCH 4/9] set devel version

---
 DESCRIPTION | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index b7e939b..106ee7c 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: nn2poly
 Title: Neural Network Weights Transformation into Polynomial Coefficients
-Version: 1.0.0
+Version: 1.0.0.9000
 Authors@R: c(
     person(given = "Pablo", family = "Morala", role = c("aut", "cre"),
            email = "moralapablo@gmail.com",

From 0ac51f752a5765131b6a795bd8ffa71c971065c4 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?I=C3=B1aki=20=C3=9Acar?= <iucar@fedoraproject.org>
Date: Fri, 17 Jul 2026 15:20:02 +0200
Subject: [PATCH 5/9] fix comment

---
 src/nn2poly.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/nn2poly.cpp b/src/nn2poly.cpp
index 8f2bff9..aa62448 100644
--- a/src/nn2poly.cpp
+++ b/src/nn2poly.cpp
@@ -35,7 +35,7 @@ CoeffsList obtain_derivatives_list(const Term& taylor_orders,
     if (taylor_orders[i] < 0)
       throw std::invalid_argument("`taylor_orders` must be non-negative");
     // Obtain the vector with the derivatives of the activation function up to
-    // the given degree centered at 0
+    // the given degree centered at a
     out[i].resize(static_cast<size_t>(taylor_orders[i] + 1));
     coeffs_taylor(af_string_list[i], out[i], a);
   }

From 7ebe28741195fce5a83ee271b587953f04b63dac Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?I=C3=B1aki=20=C3=9Acar?= <iucar@fedoraproject.org>
Date: Sat, 18 Jul 2026 13:57:27 +0200
Subject: [PATCH 6/9] refactor approximation functions

---
 src/approx.h        | 14 ++++++++++
 src/approx_taylor.h | 64 +++++++++++++++++++++++++++++++++++++++++++++
 src/nn2poly.cpp     |  4 +--
 src/taylor.h        | 56 ---------------------------------------
 4 files changed, 80 insertions(+), 58 deletions(-)
 create mode 100644 src/approx.h
 create mode 100644 src/approx_taylor.h
 delete mode 100644 src/taylor.h

diff --git a/src/approx.h b/src/approx.h
new file mode 100644
index 0000000..1cec78e
--- /dev/null
+++ b/src/approx.h
@@ -0,0 +1,14 @@
+#ifndef nn2poly__approx_h
+#define nn2poly__approx_h
+
+#include "nn2poly_types.h"
+#include "approx_taylor.h"
+
+template <typename... Args>
+void coeffs(const std::string& type, const std::string& af, Coeffs& c, Args&&... args) {
+  if (type == "taylor") if constexpr (sizeof...(Args) < 2)
+    return nn2poly::detail::coeffs_taylor(af, c, std::forward<Args>(args)...);
+  throw std::invalid_argument("approximation type '" + type + "' not supported");
+}
+
+#endif
diff --git a/src/approx_taylor.h b/src/approx_taylor.h
new file mode 100644
index 0000000..5ae418d
--- /dev/null
+++ b/src/approx_taylor.h
@@ -0,0 +1,64 @@
+#ifndef nn2poly__approx_taylor_h
+#define nn2poly__approx_taylor_h
+
+namespace nn2poly {
+namespace detail {
+
+// Helper to handle the shared Taylor convolution loop structure
+template<typename Func>
+void compute_taylor_conv(Coeffs& c, Func step, double c0) {
+  c[0] = c0;
+  for (size_t n = 0; n < c.size() - 1; n++) {
+    double conv = 0.0;
+    for (size_t k = 0; k <= n; k++)
+      conv += c[k] * c[n - k];
+    c[n + 1] = step(n, conv, c) / static_cast<double>(n + 1);
+  }
+}
+
+inline void coeffs_taylor_sigmoid(Coeffs& c, double a) {
+  compute_taylor_conv(c, [](size_t n, double conv, const Coeffs& c) {
+    return (c[n] - conv);
+  }, 1.0 / (1.0 + std::exp(-a)));
+}
+
+inline void coeffs_taylor_tanh(Coeffs& c, double a) {
+  compute_taylor_conv(c, [](size_t n, double conv, const Coeffs&) {
+    return (n == 0 ? 1.0 : 0.0) - conv;
+  }, std::tanh(a));
+}
+
+inline void coeffs_taylor_linear(Coeffs& c, double a) {
+  c[0] = a;
+  if (c.size() > 1) {
+    c[1] = 1.0;
+    std::fill(c.begin() + 2, c.end(), 0.0);
+  }
+}
+
+inline void coeffs_taylor_softplus(Coeffs& c, double a) {
+  coeffs_taylor_sigmoid(c, a);
+  for (size_t n = c.size() - 1; n > 0; n--)
+    c[n] = c[n - 1] / static_cast<double>(n);
+  if (a > 0.0)
+    c[0] = a + std::log(1.0 + std::exp(-a));
+  else c[0] = std::log(1.0 + std::exp(a));
+}
+
+inline void coeffs_taylor(const std::string& af, Coeffs& c, double a = 0.0) {
+  if (c.empty()) return;
+  if (af == "sigmoid")
+    return coeffs_taylor_sigmoid(c, a);
+  if (af == "tanh")
+    return coeffs_taylor_tanh(c, a);
+  if (af == "softplus")
+    return coeffs_taylor_softplus(c, a);
+  if (af == "linear")
+    return coeffs_taylor_linear(c, a);
+  throw std::invalid_argument("function '" + af + "' not supported");
+}
+
+} // namespace detail
+} // namespace nn2poly
+
+#endif
diff --git a/src/nn2poly.cpp b/src/nn2poly.cpp
index aa62448..82bf17f 100644
--- a/src/nn2poly.cpp
+++ b/src/nn2poly.cpp
@@ -1,6 +1,6 @@
 #include "partitions.h"
 #include "multiset.h"
-#include "taylor.h"
+#include "approx.h"
 using namespace nn2poly::linalg;
 
 // [[Rcpp::export]]
@@ -37,7 +37,7 @@ CoeffsList obtain_derivatives_list(const Term& taylor_orders,
     // Obtain the vector with the derivatives of the activation function up to
     // the given degree centered at a
     out[i].resize(static_cast<size_t>(taylor_orders[i] + 1));
-    coeffs_taylor(af_string_list[i], out[i], a);
+    coeffs("taylor", af_string_list[i], out[i], a);
   }
 
   return out;
diff --git a/src/taylor.h b/src/taylor.h
deleted file mode 100644
index 622d456..0000000
--- a/src/taylor.h
+++ /dev/null
@@ -1,56 +0,0 @@
-#ifndef nn2poly__taylor_h
-#define nn2poly__taylor_h
-
-#include "nn2poly_types.h"
-
-void coeffs_taylor_sigmoid(Coeffs& c, double a) {
-  c[0] = 1.0 / (1.0 + std::exp(-a));
-  for (size_t n = 0; n < c.size() - 1; n++) {
-    double conv = 0.0;
-    for (size_t k = 0; k <= n; k++)
-      conv += c[k] * c[n - k];
-    c[n + 1] = (c[n] - conv) / static_cast<double>(n + 1);
-  }
-}
-
-void coeffs_taylor_tanh(Coeffs& c, double a) {
-  c[0] = std::tanh(a);
-  for (size_t n = 0; n < c.size() - 1; n++) {
-    double conv = 0.0;
-    for (size_t k = 0; k <= n; k++)
-      conv += c[k] * c[n - k];
-    const double rhs = (n == 0 ? 1.0 : 0.0) - conv;
-    c[n + 1] = rhs / static_cast<double>(n + 1);
-  }
-}
-
-void coeffs_taylor_linear(Coeffs& c, double a) {
-  c[0] = a;
-  if (c.size() > 1) {
-    c[1] = 1.0;
-    std::fill(c.begin() + 2, c.end(), 0.0);
-  }
-}
-
-void coeffs_taylor_softplus(Coeffs& c, double a) {
-  coeffs_taylor_sigmoid(c, a);
-  for (size_t n = c.size() - 1; n > 0; n--)
-    c[n] = c[n - 1] / static_cast<double>(n);
-  if (a > 0.0)
-    c[0] = a + std::log(1.0 + std::exp(-a));
-  else c[0] = std::log(1.0 + std::exp(a));
-}
-
-void coeffs_taylor(const std::string& af, Coeffs& c, double a = 0.0) {
-  if (af == "sigmoid")
-    return coeffs_taylor_sigmoid(c, a);
-  if (af == "tanh")
-    return coeffs_taylor_tanh(c, a);
-  if (af == "softplus")
-    return coeffs_taylor_softplus(c, a);
-  if (af == "linear")
-    return coeffs_taylor_linear(c, a);
-  throw std::invalid_argument("function '" + af + "' not supported");
-}
-
-#endif

From a254bfb0271ec1fa64414963784f10a886b0e426 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?I=C3=B1aki=20=C3=9Acar?= <iucar@fedoraproject.org>
Date: Wed, 22 Jul 2026 15:27:00 +0200
Subject: [PATCH 7/9] rename new variable, update NEWS

---
 NEWS.md                     |  4 ++++
 R/RcppExports.R             |  8 ++++----
 R/nn2poly.R                 | 10 +++++-----
 man/nn2poly.Rd              |  6 +++---
 src/RcppExports.cpp         | 16 ++++++++--------
 src/nn2poly.cpp             | 12 ++++++------
 tests/testthat/test-smoke.R |  5 +++--
 7 files changed, 33 insertions(+), 28 deletions(-)

diff --git a/NEWS.md b/NEWS.md
index f09c8f4..c045bcf 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,7 @@
+# nn2poly 1.0.1.9000
+
+- Implement Taylor expansion at any given point (#85 addressing #82).
+
 # nn2poly 1.0.1
 
 - Fix UB in hash implementation (#84).
diff --git a/R/RcppExports.R b/R/RcppExports.R
index 1c5caf6..21649b0 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -5,16 +5,16 @@ obtain_taylor_vector <- function(taylor_orders, af_string_list) {
     .Call(`_nn2poly_obtain_taylor_vector`, taylor_orders, af_string_list)
 }
 
-obtain_derivatives_list <- function(taylor_orders, af_string_list, a = 0.0) {
-    .Call(`_nn2poly_obtain_derivatives_list`, taylor_orders, af_string_list, a)
+obtain_derivatives_list <- function(taylor_orders, af_string_list, taylor_center = 0.0) {
+    .Call(`_nn2poly_obtain_derivatives_list`, taylor_orders, af_string_list, taylor_center)
 }
 
 alg_non_linear <- function(coeffs_input, labels_input, labels_output, previous_order, q_layer, g) {
     .Call(`_nn2poly_alg_non_linear`, coeffs_input, labels_input, labels_output, previous_order, q_layer, g)
 }
 
-nn2poly_algorithm <- function(layers, af_list, max_order, keep_layers, taylor_orders, a = 0.0) {
-    .Call(`_nn2poly_nn2poly_algorithm`, layers, af_list, max_order, keep_layers, taylor_orders, a)
+nn2poly_algorithm <- function(layers, af_list, max_order, keep_layers, taylor_orders, taylor_center = 0.0) {
+    .Call(`_nn2poly_nn2poly_algorithm`, layers, af_list, max_order, keep_layers, taylor_orders, taylor_center)
 }
 
 combinations_with_repetition <- function(n, k) {
diff --git a/R/nn2poly.R b/R/nn2poly.R
index e260d3f..04aac03 100644
--- a/R/nn2poly.R
+++ b/R/nn2poly.R
@@ -45,8 +45,8 @@ NULL
 #  layers. If a vector is used, each value corresponds to the Taylor order used
 #' at each layer activation function. Default set to \code{8}.
 #'
-#' @param a \code{numeric} value that sets the point at which the Taylor expansion
-#' is computed. Default set to \code{0}.
+#' @param taylor_center \code{numeric} value that sets the point at which the
+#' Taylor expansion is computed. Default set to \code{0}.
 #'
 #' @param ... Ignored.
 #'
@@ -127,10 +127,10 @@ nn2poly.list <- function(object,
                          max_order = 2,
                          keep_layers = FALSE,
                          taylor_orders = 8,
-                         a = 0.0,
+                         taylor_center = 0.0,
                          ...) {
-  result <- nn2poly_algorithm(object, names(object),
-                              max_order, keep_layers, taylor_orders, a)
+  result <- nn2poly_algorithm(object, names(object), max_order, keep_layers,
+                              taylor_orders, taylor_center)
   class(result) <- "nn2poly"
   result
 }
diff --git a/man/nn2poly.Rd b/man/nn2poly.Rd
index cf4c93f..f60bb7f 100644
--- a/man/nn2poly.Rd
+++ b/man/nn2poly.Rd
@@ -50,10 +50,10 @@ degree at which Taylor expansion is truncated at each layer. If a single
 value is used, that value is set for each non linear layer and 1 for linear
 at each layer activation function. Default set to \code{8}.}
 
-\item{a}{\code{numeric} value that sets the point at which the Taylor expansion
-is computed. Default set to \code{0}.}
-
 \item{...}{Ignored.}
+
+\item{taylor_center}{\code{numeric} value that sets the point at which the
+Taylor expansion is computed. Default set to \code{0}.}
 }
 \value{
 Returns an object of class \code{nn2poly}.
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index 1227a3a..330e3b9 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -25,15 +25,15 @@ BEGIN_RCPP
 END_RCPP
 }
 // obtain_derivatives_list
-CoeffsList obtain_derivatives_list(const Term& taylor_orders, const Functions& af_string_list, double a);
-RcppExport SEXP _nn2poly_obtain_derivatives_list(SEXP taylor_ordersSEXP, SEXP af_string_listSEXP, SEXP aSEXP) {
+CoeffsList obtain_derivatives_list(const Term& taylor_orders, const Functions& af_string_list, double taylor_center);
+RcppExport SEXP _nn2poly_obtain_derivatives_list(SEXP taylor_ordersSEXP, SEXP af_string_listSEXP, SEXP taylor_centerSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
     Rcpp::traits::input_parameter< const Term& >::type taylor_orders(taylor_ordersSEXP);
     Rcpp::traits::input_parameter< const Functions& >::type af_string_list(af_string_listSEXP);
-    Rcpp::traits::input_parameter< double >::type a(aSEXP);
-    rcpp_result_gen = Rcpp::wrap(obtain_derivatives_list(taylor_orders, af_string_list, a));
+    Rcpp::traits::input_parameter< double >::type taylor_center(taylor_centerSEXP);
+    rcpp_result_gen = Rcpp::wrap(obtain_derivatives_list(taylor_orders, af_string_list, taylor_center));
     return rcpp_result_gen;
 END_RCPP
 }
@@ -54,8 +54,8 @@ BEGIN_RCPP
 END_RCPP
 }
 // nn2poly_algorithm
-List nn2poly_algorithm(const Layers& layers, const Functions& af_list, int max_order, bool keep_layers, const Term& taylor_orders, double a);
-RcppExport SEXP _nn2poly_nn2poly_algorithm(SEXP layersSEXP, SEXP af_listSEXP, SEXP max_orderSEXP, SEXP keep_layersSEXP, SEXP taylor_ordersSEXP, SEXP aSEXP) {
+List nn2poly_algorithm(const Layers& layers, const Functions& af_list, int max_order, bool keep_layers, const Term& taylor_orders, double taylor_center);
+RcppExport SEXP _nn2poly_nn2poly_algorithm(SEXP layersSEXP, SEXP af_listSEXP, SEXP max_orderSEXP, SEXP keep_layersSEXP, SEXP taylor_ordersSEXP, SEXP taylor_centerSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
@@ -64,8 +64,8 @@ BEGIN_RCPP
     Rcpp::traits::input_parameter< int >::type max_order(max_orderSEXP);
     Rcpp::traits::input_parameter< bool >::type keep_layers(keep_layersSEXP);
     Rcpp::traits::input_parameter< const Term& >::type taylor_orders(taylor_ordersSEXP);
-    Rcpp::traits::input_parameter< double >::type a(aSEXP);
-    rcpp_result_gen = Rcpp::wrap(nn2poly_algorithm(layers, af_list, max_order, keep_layers, taylor_orders, a));
+    Rcpp::traits::input_parameter< double >::type taylor_center(taylor_centerSEXP);
+    rcpp_result_gen = Rcpp::wrap(nn2poly_algorithm(layers, af_list, max_order, keep_layers, taylor_orders, taylor_center));
     return rcpp_result_gen;
 END_RCPP
 }
diff --git a/src/nn2poly.cpp b/src/nn2poly.cpp
index 82bf17f..1747068 100644
--- a/src/nn2poly.cpp
+++ b/src/nn2poly.cpp
@@ -25,7 +25,7 @@ Term obtain_taylor_vector(const Term& taylor_orders,
 // [[Rcpp::export]]
 CoeffsList obtain_derivatives_list(const Term& taylor_orders,
                                    const Functions& af_string_list,
-                                   double a = 0.0) {
+                                   double taylor_center = 0.0) {
   if (taylor_orders.size() != af_string_list.size())
     throw std::invalid_argument(
       "`taylor_orders` length does not match provided number of layers");
@@ -35,9 +35,9 @@ CoeffsList obtain_derivatives_list(const Term& taylor_orders,
     if (taylor_orders[i] < 0)
       throw std::invalid_argument("`taylor_orders` must be non-negative");
     // Obtain the vector with the derivatives of the activation function up to
-    // the given degree centered at a
+    // the given degree centered at taylor_center
     out[i].resize(static_cast<size_t>(taylor_orders[i] + 1));
-    coeffs("taylor", af_string_list[i], out[i], a);
+    coeffs("taylor", af_string_list[i], out[i], taylor_center);
   }
 
   return out;
@@ -137,7 +137,7 @@ inline void check_weights_dimensions(const Layers& layers) {
 // [[Rcpp::export]]
 List nn2poly_algorithm(const Layers& layers, const Functions& af_list,
                        int max_order, bool keep_layers,
-                       const Term& taylor_orders, double a = 0.0) {
+                       const Term& taylor_orders, double taylor_center = 0.0) {
   if (layers.empty())
     throw std::invalid_argument("`layers` is empty");
   if (af_list.empty())
@@ -156,7 +156,7 @@ List nn2poly_algorithm(const Layers& layers, const Functions& af_list,
   // depending on the last layer being linear or not
   WeightsLists results(static_cast<size_t>(last_linear ? 2 * L - 1 : 2 * L));
   const Term taylor = obtain_taylor_vector(taylor_orders, af_list);
-  const CoeffsList af_dlist = obtain_derivatives_list(taylor, af_list, a);
+  const CoeffsList af_dlist = obtain_derivatives_list(taylor, af_list, taylor_center);
 
   // Starting point for the algorithm: Set weights as coefficients
   // of an order 1 polynomial.
@@ -226,7 +226,7 @@ List nn2poly_algorithm(const Layers& layers, const Functions& af_list,
       DTAG(previous_order), DTAG(q_layer), DTAG(labels_list.size()));
 
     // Subtract the center from the intercept
-    if (a) sub_scalar(coeffs_list, 0, a);
+    if (taylor_center) sub_scalar(coeffs_list, 0, taylor_center);
     // Apply non-linear algorithm
     coeffs_list = alg_non_linear_impl(
       coeffs_list,
diff --git a/tests/testthat/test-smoke.R b/tests/testthat/test-smoke.R
index 7711365..98095e1 100644
--- a/tests/testthat/test-smoke.R
+++ b/tests/testthat/test-smoke.R
@@ -485,10 +485,11 @@ test_that("manual example passes", {
   data <- smoke_test_data()
 
   x <- smoke_test_manual(data)
-  y <- nn2poly(data, max_order = 4, taylor = 5, keep_layers = TRUE)
+  y <- nn2poly(data, max_order = 4, keep_layers = TRUE, taylor_orders = 5)
   expect_equal(x, y)
 
   x <- smoke_test_manual(data, 0.5)
-  y <- nn2poly(data, max_order = 4, taylor = 5, keep_layers = TRUE, a = 0.5)
+  y <- nn2poly(data, max_order = 4, keep_layers = TRUE,
+               taylor_orders = 5, taylor_center = 0.5)
   expect_equal(x, y)
 })

From eee716384a464916cea8788c307ec92ff1262970 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?I=C3=B1aki=20=C3=9Acar?= <iucar@fedoraproject.org>
Date: Wed, 22 Jul 2026 15:28:26 +0200
Subject: [PATCH 8/9] update NEWS

---
 NEWS.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/NEWS.md b/NEWS.md
index c045bcf..192920b 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,6 +1,6 @@
 # nn2poly 1.0.1.9000
 
-- Implement Taylor expansion at any given point (#85 addressing #82).
+- Implement Taylor expansion at any given `taylor_center` (#85 addressing #82).
 
 # nn2poly 1.0.1
 

From 694be25288e391bba2bb5e2e613fffa79951f244 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?I=C3=B1aki=20=C3=9Acar?= <iucar@fedoraproject.org>
Date: Wed, 22 Jul 2026 15:33:18 +0200
Subject: [PATCH 9/9] forgot one instance

---
 R/nn2poly.R    | 2 +-
 man/nn2poly.Rd | 6 +++---
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/R/nn2poly.R b/R/nn2poly.R
index 04aac03..22405c8 100644
--- a/R/nn2poly.R
+++ b/R/nn2poly.R
@@ -117,7 +117,7 @@ nn2poly <- function(object,
                     max_order = 2,
                     keep_layers = FALSE,
                     taylor_orders = 8,
-                    a = 0.0,
+                    taylor_center = 0.0,
                     ...) {
   UseMethod("nn2poly")
 }
diff --git a/man/nn2poly.Rd b/man/nn2poly.Rd
index f60bb7f..6aa2cd0 100644
--- a/man/nn2poly.Rd
+++ b/man/nn2poly.Rd
@@ -9,7 +9,7 @@ nn2poly(
   max_order = 2,
   keep_layers = FALSE,
   taylor_orders = 8,
-  a = 0,
+  taylor_center = 0,
   ...
 )
 }
@@ -50,10 +50,10 @@ degree at which Taylor expansion is truncated at each layer. If a single
 value is used, that value is set for each non linear layer and 1 for linear
 at each layer activation function. Default set to \code{8}.}
 
-\item{...}{Ignored.}
-
 \item{taylor_center}{\code{numeric} value that sets the point at which the
 Taylor expansion is computed. Default set to \code{0}.}
+
+\item{...}{Ignored.}
 }
 \value{
 Returns an object of class \code{nn2poly}.