Modify the coding regulation of suffix array. Add adaptive simpson's rule algorithm and manacher algorithm. Add a new chapter named Numeric to the book.

Author: haomingbai

main.tex

@@ -299,4 +299,30 @@ \section{后缀数组}
 
 \lstinputlisting[language=C++, caption=suffix\_array.cpp, style=MyCStyle]{./str/suffix_array.cpp}
 
+\section{最长回文子串}
+
+最长回文子串问题是指: 在一个字符串中寻找其中最长的对称子串. 换句话说, 这个子串正着读和反着读完全相同.
+
+对于最长回文子串问题, 一种朴素的解法是枚举中心并向两边扩展, 或者使用其它动态规划的方法, 但是这些方式在字符串较长时往往会遇到效率瓶颈.
+
+Manacher算法是一种在线性时间内解决最长回文子串的算法. 它通过在原始字符串中插入分隔符并利用回文的对称性, 避免了重复计算, 将复杂度降低到$O(n)$. 这种高效性使得它不仅是理论研究的成果, 也在一些实际场景中得到了应用. 例如, 在DNA序列分析中寻找对称片段, 在自然语言处理里发现对称的语言结构, 以及在数据压缩中利用回文模式进行优化等.
+
+\lstinputlisting[language=C++, caption=manacher.cpp, style=MyCStyle]{./str/manacher.cpp}
+
+\chapter{数值方法}
+
+在生产实践中, 尤其是大型工业软件的开发过程中, 数值计算的速度和精度通常和软件的运行效率息息相关. 在数学建模, 仿真, 有限元分析等领域, 使用精度高, 效率高的算法通常可以使程序的性能, 比使用低效算法的相同程序提升十倍甚至百倍. 而在生产环境, 微不足道的浮点误差可以造成惨重的人员和财产损失. 由此观之, 在建设中国式现代化的过程中, 为了解决工业软件领域的 "卡脖子" 难题, 提升数值计算的效率是非常重要的一环. 而掌握正确的计算方法, 则是每一位立志成为 "总师型人才" 的计算机专业学生的必修课. 下面笔者将提供一些简单常用的数值计算方法的实现以供参考, 所有代码都配备详尽注释, 力求通俗易懂, 方便移植改造.
+
+\section{辛普森法}
+
+辛普森法是一种常用的数值积分方法, 它的原理是利用二次插值多项式在小区间上近似被积函数, 从而得到定积分的近似值.
+
+然而, 当积分区间较大或被积函数变化复杂时, 单一的辛普森公式往往难以保证精度. 为此, 可以采用复化辛普森法: 将积分区间划分为若干个小区间, 在每个小区间分别应用辛普森公式, 然后将结果加总.这样虽然提高了精度, 但由于分块数目固定, 可能会在函数比较平滑的区间上产生不必要的计算.
+
+自适应辛普森法则通过递归的方式在尽可能少分块的前提下保证精度. 它先在整个区间上应用一次辛普森公式, 再将区间一分为二分别应用辛普森公式, 并比较两者结果.如果二者差异超过误差容许范围, 就继续对子区间递归分裂; 否则接受结果.这样, 自适应辛普森法能够在函数平滑的区域少分块, 在函数变化剧烈的区域多分块, 从而兼顾效率和精度.
+
+在现实应用中, 辛普森法和自适应辛普森法被广泛用于工程与科学计算. 例如, 在航天工程中进行复杂气动力模型的积分估算, 在医学成像中计算复杂函数的积分, 或在金融工程中估算涉及不规则函数的积分问题.
+
+\lstinputlisting[language=C++, caption=asr.cpp, style=MyCStyle]{./numeric/asr.cpp}
+
 \end{document}

numeric/asr.cpp

@@ -0,0 +1,135 @@
+/**
+ * @file asr.cpp
+ * @brief
+ * @author Haoming Bai <haomingbai@hotmail.com>
+ * @date   2025-08-20
+ *
+ * Copyright © 2025 Haoming Bai
+ * SPDX-License-Identifier: MIT
+ *
+ * @details
+ */
+
+#include <cmath>
+#include <concepts>
+#include <cstdlib>
+
+/**
+ * @brief Simpson's Rule for numerical integration.
+ *
+ * Simpson's Rule is a method for approximating the definite integral of a
+ * function by fitting a quadratic polynomial to the function at three points
+ * (the endpoints and the midpoint) and integrating that polynomial over the
+ * interval [a, b]. The formula for the approximation is: (b - a)/6 * [f(a) +
+ * 4f((a+b)/2) + f(b)].
+ *
+ * This method provides exact results for polynomials up to degree 3 and
+ * generally offers better accuracy than the trapezoidal rule for smooth
+ * functions.
+ *
+ * @tparam Func Callable type. Must be invocable with a floating-point argument
+ * and return a value convertible to a floating-point type.
+ * @tparam F1   Floating-point type for the lower bound.
+ * @tparam F2   Floating-point type for the upper bound.
+ *
+ * @param f The function to be integrated. Should be continuous on [a, b].
+ * @param a Lower limit of integration.
+ * @param b Upper limit of integration.
+ *
+ * @return Approximation of the integral ∫_{a}^{b} f(x) dx.
+ *         The return type is deduced as the common type of F1 and F2.
+ *
+ * @note This is a non-adaptive implementation that provides a single
+ * approximation over the entire interval. For higher accuracy over complex
+ * functions, consider using an adaptive method like Adaptive Simpson's Rule
+ * (ASR).
+ * @note Requires that `f` be defined and continuous on [a, b].
+ *
+ * Example usage:
+ * @code{.cpp}
+ * auto result = Simpson([](double x) { return std::sin(x); }, 0.0, 3.14159);
+ * @endcode
+ */
+template <typename Func, std::floating_point F1, std::floating_point F2>
+auto Simpson(Func f, F1 a, F2 b)
+    -> decltype(std::declval<F1>() + std::declval<F2>()) {
+  using ReturnType = decltype(std::declval<F1>() + std::declval<F2>());
+
+  ReturnType mid = a + (b - a) / 2.0;
+  ReturnType result = f(a) + f(b) + f(mid) * 4.0;
+  result *= (b - a);
+  result /= 6.0;
+
+  return result;
+}
+
+/**
+ * @brief Adaptive Simpson's Rule (ASR) for numerical integration.
+ *
+ * Adaptive Simpson's Rule is a recursive method for approximating the definite
+ * integral of a function. It dynamically refines the integration interval by
+ * subdividing it into smaller segments where the function exhibits more
+ * variation, and uses coarse segments where the function is smooth. This allows
+ * for efficient computation by minimizing the number of function evaluations
+ * while achieving the desired accuracy.
+ *
+ * The method compares the Simpson's rule estimates on the current segment and
+ * its two halves. If the difference between the composite estimate and the
+ * whole segment estimate is within the error tolerance, the result is accepted.
+ * Otherwise, the algorithm recurses on both halves.
+ *
+ * @tparam Func Callable type. Must be invocable with a floating-point argument
+ * and return a value convertible to a floating-point type.
+ * @tparam F1   Floating-point type for the lower bound.
+ * @tparam F2   Floating-point type for the upper bound.
+ * @tparam Ferr Floating-point type for the error tolerance.
+ *
+ * @param f     The function to be integrated. Must be continuous on [a, b].
+ * @param a     Lower limit of integration.
+ * @param b     Upper limit of integration.
+ * @param error Maximum allowable error tolerance for the integral
+ * approximation.
+ *
+ * @return Approximation of the integral ∫_{a}^{b} f(x) dx.
+ *         The return type is deduced as the common type of F1 and F2.
+ *
+ * @note The function uses recursion and may cause stack overflow for extremely
+ * tight tolerances or highly oscillatory functions over large intervals.
+ * @note Requires that `f` be defined and continuous on [a, b].
+ *
+ * Example usage:
+ * @code{.cpp}
+ * auto result = ASR([](double x) { return std::sin(x); }, 0.0, 3.14159, 1e-6);
+ * @endcode
+ */
+template <typename Func, std::floating_point F1, std::floating_point F2,
+          std::floating_point Ferr>
+auto ASR(Func f, F1 a, F2 b, Ferr error)
+    -> decltype(std::declval<F1>() + std::declval<F2>()) {
+  using ReturnType = decltype(std::declval<F1>() + std::declval<F2>());
+  ReturnType mid = a + (b - a) / 2.0;
+
+  // 先对整体使用辛普森法
+  auto S0 = Simpson(f, a, b);
+  // 再将区间分成两块求解
+  auto S1 = Simpson(f, a, mid);
+  auto S2 = Simpson(f, mid, b);
+
+  // 求解此时的误差
+  Ferr curr_error = std::abs(S0 - (S1 + S2));
+
+  // 误差如果超过最大允许范围,
+  // (因为在返回过程中, 差值会被除以15,
+  // 所以curr_error是15倍的当前误差)
+  // 那么就递归求解两侧的面积,
+  // 同时要求两侧的误差不超过最大误差的一半.
+  // 否则就直接返回结果.
+  if (curr_error >= 15.0 * error) {
+    ReturnType result =
+        ASR(f, a, mid, error / 2.0) + ASR(f, mid, b, error / 2.0);
+    return result;
+  } else {
+    ReturnType result = (S1 + S2) + (S1 + S2 - S0) / 15.0;
+    return result;
+  }
+}

str/manacher.cpp

@@ -0,0 +1,146 @@
+/**
+ * @file manacher.cpp
+ * @brief
+ * @author Haoming Bai <haomingbai@hotmail.com>
+ * @date   2025-08-21
+ *
+ * Copyright © 2025 Haoming Bai
+ * SPDX-License-Identifier: MIT
+ *
+ * @details
+ */
+
+#include <algorithm>
+#include <array>
+#include <concepts>
+#include <vector>
+
+#include "../concepts.cpp"
+
+struct ManacherResult {
+  long start, len;
+};
+
+template <std::equality_comparable E, RandomStdContainer<E> Container>
+ManacherResult Manacher(Container &&str) {
+  if (str.empty()) {
+    return {.start = -1, .len = 0};
+  }
+
+  if (str.size() == 1) {
+    return {.start = 0, .len = 1};
+  }
+
+  enum Type : bool { ODD, EVEN };
+
+  std::array<std::vector<long>, 2> palindrome_lens{
+      std::vector<long>(str.size(), 1), std::vector<long>(str.size() - 1, 0)};
+
+  // 先讨论奇数长度的最长回文串.
+  {
+    // 这两个变量表示当前最靠右的回文串的左右边界.
+    // 初始状态下, 因为单字符是回文串,
+    // 所以二者下标为0, 因为第0个字符为中心的最长回文串就是str[0].
+    long L = 0, R = 0;
+
+    for (long i = 1; i < str.size(); i++) {
+      // 这个curr_len表示从这个长度开始扩张.
+      long curr_len;
+
+      // 分两种情况讨论, 如果被遍历到的i
+      // 在当前最右的回文串内部, 那么可以从对称位置开始扩张.
+      // 如果不在, 那么直接用朴素算法开始扩张.
+      if (i <= R) {
+        // 这个时候满足:
+        // R - i = mirror_i - L
+        auto mirror_i = R - i + L;
+
+        curr_len = std::min(palindrome_lens[ODD][mirror_i], (R - i) * 2 + 1);
+      } else {
+        curr_len = 1;
+      }
+
+      auto curr_stretch = curr_len / 2;
+
+      while (i > curr_stretch && i + curr_stretch + 1 < str.size()) {
+        // 尝试把边界延伸1
+        // 如果新加入的两个字符相等, 则还能延伸1
+        if (str[i + curr_stretch + 1] == str[i - curr_stretch - 1]) {
+          curr_stretch++;
+        } else {
+          break;
+        }
+      }
+
+      auto curr_l = i - curr_stretch, curr_r = i + curr_stretch;
+      // 如果延伸出去了, 就更新L和R
+      if (curr_r > R) {
+        L = curr_l, R = curr_r;
+      }
+
+      // 更新长度.
+      curr_len = curr_stretch * 2 + 1;
+      palindrome_lens[ODD][i] = curr_len;
+    }
+  }
+
+  // 偶数长度的串, 情况要复杂很多,
+  // 这里朴素算法只能默认从0长度开始.
+  {
+    long L = 0, R = -1;
+
+    for (long i = 0; i < str.size() - 1; i++) {
+      // 这个curr_len表示从这个长度开始扩张.
+      long curr_len;
+
+      if (i + 1 <= R) {
+        // R - i = mirror_i + 1 - L
+        auto mirror_i = R - i + L - 1;
+        curr_len = std::min(palindrome_lens[EVEN][mirror_i], (R - i) * 2);
+      } else {
+        curr_len = 0;
+      }
+
+      auto curr_stretch = curr_len / 2 - 1;
+      while (i - curr_stretch > 0 && i + curr_stretch + 1 < str.size() - 1) {
+        if (str[i + 1 + curr_stretch + 1] == str[i - curr_stretch - 1]) {
+          curr_stretch++;
+        } else {
+          break;
+        }
+      }
+
+      auto curr_l = i + curr_stretch + 1, curr_r = i - curr_stretch;
+      if (curr_r > R) {
+        L = curr_l, R = curr_r;
+      }
+
+      curr_len = (curr_stretch + 1) * 2;
+      palindrome_lens[EVEN][i] = curr_len;
+    }
+  }
+
+  ManacherResult result;
+
+  {
+    long start = 0, len = 1;
+    for (long i = 0; i < str.size(); i++) {
+      if (palindrome_lens[ODD][i] > len) {
+        auto curr_stretch = palindrome_lens[ODD][i] / 2;
+        start = i - curr_stretch;
+        len = palindrome_lens[ODD][i];
+      }
+    }
+    for (long i = 0; i < str.size() - 1; i++) {
+      if (palindrome_lens[EVEN][i] > len) {
+        auto curr_stretch = palindrome_lens[EVEN][i] / 2 - 1;
+        start = i - curr_stretch;
+        len = palindrome_lens[EVEN][i];
+      }
+    }
+
+    result = {.start = start, .len = len};
+  }
+
+  return result;
+}