数值计算工具

#ifndef _NUMERICAL_H_
#define _NUMERICAL_H_

#include <algorithm>
#include <cmath>
#include <limits>
#include <vector>

namespace {

// 全局精度容差，用于判断收敛和误差控制
inline constexpr double accuracy = 1e-6;

// 辅助函数：计算混合绝对/相对容差，避免在接近零时过度敏感
inline double tolerance(double x) { return accuracy * (1.0 + std::abs(x)); }

}  // namespace

// ==================================================================
// 二分法求根（适用于已知区间内有符号变化的情况）
// ==================================================================

/// 在 [0, middle] 区间内搜索函数 fn 的根（假设 fn(middle) 有定义）
template <typename Fn>
double numeric_search_r1(double middle, Fn&& fn) {
  double min = 0.0;
  double max = middle;
  double ref = fn(middle);
  // 若参考点已是根或无效值，直接返回
  if (ref == 0.0 || !std::isfinite(ref)) return middle;

  double md = (min + max) * 0.5;  // 初始中点
  int cnt = 1000;                 // 最大迭代次数防死循环
  // 使用混合容差判断收敛，避免 md 接近 0 时除零
  while ((max - min) > tolerance(md) && --cnt > 0) {
    double v = fn(md);
    if (!std::isfinite(v)) break;  // 遇到 NaN/Inf 提前退出

    double ml = v * ref;  // 判断是否与参考点同号
    if (ml > 0) {
      max = md;  // 同号，根在左半区间
    } else if (ml < 0) {
      min = md;  // 异号，根在右半区间
    } else {
      return md;  // 找到精确根
    }
    md = (min + max) * 0.5;
  }
  return md;
}

/// 在 [middle, +∞) 区间内搜索根，通过变量替换 x = 1/t 转换为 (0, 1/middle] 区间
template <typename Fn>
double numeric_search_r2(double middle, Fn&& fn) {
  if (middle <= 0.0) return std::numeric_limits<double>::quiet_NaN();

  double min = 0.0;
  double max = 1.0 / middle;
  double ref = fn(middle);
  if (ref == 0.0 || !std::isfinite(ref)) return middle;

  double md = (min + max) * 0.5;
  int cnt = 1000;
  while (md > 0.0 && (max - min) > tolerance(md) && --cnt > 0) {
    double x = 1.0 / md;  // 变量替换回原始空间
    double v = fn(x);
    if (!std::isfinite(v)) break;

    double ml = v * ref;
    if (ml > 0) {
      max = md;
    } else if (ml < 0) {
      min = md;
    } else {
      return x;
    }
    md = (min + max) * 0.5;
  }
  // 防止除零
  return (md != 0.0) ? 1.0 / md : std::numeric_limits<double>::quiet_NaN();
}

// ==================================================================
// 数值积分：使用 Simpson 3/8 规则
// ==================================================================

/// 计算 ∫ₐᵇ fn(x) dx，使用 Simpson 3/8 复合公式
/**
 * @param fn 被积函数
 * @param a 积分下限
 * @param b 积分上限
 * @param steps 子区间数量（默认33，需为3的倍数效果更佳，但非强制）
 * @return 积分近似值
 */
template <typename Fn>
double numeric_integral(Fn&& fn, double a, double b, unsigned int steps = 33) {
  if (a == b || steps == 0) return 0.0;

  double h = (b - a) / static_cast<double>(steps);  // 步长
  double res = 0.0;
  double prev_fn = fn(a);  // 缓存前一个端点函数值，避免重复计算

  for (unsigned int i = 0; i < steps; ++i) {
    double ia = a + i * h;  // 当前子区间左端点
    double ib = ia + h;     // 当前子区间右端点
    double fnb = fn(ib);    // 右端点函数值

    // Simpson 3/8 需要两个内点
    double x1 = ia + h / 3.0;        // 1/3 点
    double x2 = ia + 2.0 * h / 3.0;  // 2/3 点
    double fnc = fn(x1);
    double fnd = fn(x2);

    // Simpson 3/8 公式：(b-a)/8 * [f(a) + 3f((2a+b)/3) + 3f((a+2b)/3) + f(b)]
    double segment = (ib - ia) * (prev_fn + 3.0 * fnc + 3.0 * fnd + fnb) / 8.0;
    res += segment;
    prev_fn = fnb;  // 为下一轮缓存
  }
  return res;
}

// ==================================================================
// 牛顿法求根（带保护机制）
// ==================================================================

/// 使用牛顿-拉弗森法求解 fn(x) = ofs 的根
/**
 * @param fn 原函数
 * @param dfn 导函数
 * @param ofs 目标值（通常为0）
 * @param initg 初始猜测值
 * @return 近似根
 */
template <typename Fn, typename dFn>
double newtonRoot(Fn&& fn, dFn&& dfn, double ofs, double initg) {
  double x = initg;
  double y = fn(x) - ofs;
  int max_iter = 100;

  // 相对误差判断收敛，避免 ofs=0 时除零
  while (std::abs(y) / (std::abs(ofs) + std::abs(y) + accuracy) > accuracy &&
         --max_iter > 0) {
    double dy = dfn(x);
    // 导数太小 → 牛顿法失效，退出
    if (std::abs(dy) < accuracy) break;
    x -= y / dy;
    y = fn(x) - ofs;
    // 若出现非有限值，提前终止
    if (!std::isfinite(x) || !std::isfinite(y)) break;
  }
  return x;
}

// ==================================================================
// 自适应积分表生成（递归实现）
// ==================================================================

/// 递归生成积分表（内部函数）
/**
 * 使用自适应策略细分区间，当局部误差小于阈值时停止递归。
 * @param fn 被积函数
 * @param a, b 当前区间
 * @param fa, fb 端点函数值（避免重复计算）
 * @param error 允许的最大局部误差
 * @param y 当前累积积分值
 * @param lev 递归深度（防止栈溢出）
 * @param out 输出回调函数 (x, cumulative_integral)
 */
template <typename Fn, typename Out>
static double generateIntTable2(Fn&& fn, double a, double b, double fa,
                                double fb, double error, double y, int lev,
                                Out&& out) {
  double w = b - a;
  double pa = w * fa;            // 左端点贡献（梯形近似）
  double pb = w * fb;            // 右端点贡献
  double e = std::abs(pa - pb);  // 粗略误差估计

  // 若误差过大且未达最大深度，则细分
  if (e > error && lev < 16) {
    double m = (a + b) * 0.5;
    double fm = fn(m);
    // 递归处理左半区间
    double sa = generateIntTable2(fn, a, m, fa, fm, error, y, lev + 1, out);
    y += sa;
    out(m, y);  // 输出中点积分值
    // 递归处理右半区间
    double sb = generateIntTable2(fn, m, b, fm, fb, error, y, lev + 1, out);
    return sa + sb;
  } else {
    // 否则使用梯形法则近似
    return (pa + pb) * 0.5;
  }
}

/// 对外接口：生成从 a 到 b 的积分表
template <typename Fn, typename Out>
static void generateIntTable(Fn&& fn, double a, double b, double error,
                             double y, Out&& out) {
  out(a, y);  // 输出起点
  double fa = fn(a);
  double fb = fn(b);
  double r = generateIntTable2(fn, a, b, fa, fb, error, y, 0, out);
  out(b, y + r);  // 输出终点
}

// ==================================================================
// 通用求根工具类
// ==================================================================

/// 数值计算工具类（支持多种求根策略）
template <unsigned int iterations = 30>
class Numerics {
 public:
  /// 在 [from, to] 区间内查找根（要求函数在端点处符号相反）
  template <typename Fn>
  static double find_root(double from, double to, Fn&& fn) {
    double f_from = fn(from);
    if (!std::isfinite(f_from)) return std::numeric_limits<double>::quiet_NaN();
    if (std::abs(f_from) < accuracy) return from;  // 起点已是根

    double f_to = fn(to);
    if (!std::isfinite(f_to)) return std::numeric_limits<double>::quiet_NaN();

    // 必须存在符号变化，否则无根
    if (f_from * f_to > 0) {
      return std::numeric_limits<double>::quiet_NaN();
    }

    double a = from, b = to;
    bool found = false;
    for (unsigned int i = 0; i < iterations; ++i) {
      double mid = (a + b) * 0.5;
      double v = fn(mid);
      if (!std::isfinite(v)) break;
      if (std::abs(v) < accuracy) return mid;  // 精确命中

      if (f_from * v < 0) {
        b = mid;  // 根在左半
        found = true;
      } else {
        a = mid;     // 根在右半
        f_from = v;  // 更新左端点函数值
      }
    }
    return found ? (a + b) * 0.5 : std::numeric_limits<double>::quiet_NaN();
  }

  /// 在 [from, +∞) 区间查找根（通过变量替换）
  template <typename Fn>
  static double find_root_to_inf(double from, Fn&& fn) {
    if (from <= 0.0) return std::numeric_limits<double>::quiet_NaN();
    // 变量替换：x = 1/t，将 [from, ∞) 映射为 (0, 1/from]
    auto trnfn = [&](double t) -> double {
      if (t <= 0.0) return std::numeric_limits<double>::quiet_NaN();
      return fn(1.0 / t);
    };
    double res = find_root(1.0 / from, 0.0, trnfn);
    if (std::isnan(res) || res <= 0.0)
      return std::numeric_limits<double>::quiet_NaN();
    return 1.0 / res;
  }

  /// 在 (0, from] 区间查找根
  template <typename Fn>
  static double find_root_to_zero(double from, Fn&& fn) {
    if (from <= 0.0) return std::numeric_limits<double>::quiet_NaN();
    return find_root(from, 0.0, fn);
  }

  /// 根据方向自动选择搜索区间
  /**
   * @param from 起始点（必须 > 0）
   * @param dir 方向：>0 向 +∞ 搜索，<0 向 0 搜索
   */
  template <typename Fn>
  static double find_root_pos(double from, double dir, Fn&& fn) {
    if (dir > 0) {
      return find_root_to_inf(from, fn);
    } else if (dir < 0) {
      return find_root_to_zero(from, fn);
    } else {
      return std::numeric_limits<double>::quiet_NaN();
    }
  }
};

// ==================================================================
// 工具函数
// ==================================================================

/// 若 x 为 NaN 或无穷大，返回默认值 val；否则返回 x
inline double not_nan(double x, double val = 0.0) {
  return std::isfinite(x) ? x : val;
}

#endif

#ifndef _INTEGRAL_H_
#define _INTEGRAL_H_

#include <algorithm>
#include <functional>
#include <memory>
#include <type_traits>
#include <vector>

/**
 * @brief 数值积分实现类模板
 *
 * 该类通过递归梯形法则在区间 [a, b] 上对函数 fn 进行积分，
 * 并记录所有采样点 (x, ∫ₐˣ fn(t) dt)，最终构建一个可插值查询的积分函数。
 * 积分结果会整体减去 ∫ₐᵐ fn(t) dt，使得 F(m) = 0。
 *
 * @tparam X 自变量类型（如 double、float 等）
 * @tparam Y 函数值和积分值类型（通常与 X 相同或兼容）
 */
template <typename X, typename Y>
class NumericIntegralImplT {
 public:
  /**
   * @brief 生成积分表
   *
   * @param fn 被积函数（可调用对象）
   * @param a 积分下限
   * @param b 积分上限
   * @param m 偏移参考点（最终 F(m) = 0）
   * @param levels 递归层数，总采样点数 = 2^levels + 1
   */
  template <typename Fn>
  void generate(Fn&& fn, X a, X b, X m, unsigned int levels) {
    // 清空旧数据
    xvals.clear();
    yvals.clear();

    // 预分配内存：总点数 = 2^levels + 1
    std::size_t need_points = (static_cast<std::size_t>(1) << levels) + 1;
    xvals.reserve(need_points);
    yvals.reserve(need_points);

    // 起始点：x = a, 积分值 = 0
    xvals.push_back(a);
    yvals.push_back(Y{0});

    // 递归生成中间点并计算总积分 ∫ₐᵇ fn(t) dt
    Y total = do_generate(fn, a, b, fn(a), fn(b), levels);

    // 添加终点：x = b
    xvals.push_back(b);
    yvals.push_back(total);

    // 计算偏移量：F(m) = ∫ₐᵐ fn(t) dt（使用未偏移的原始积分表插值得到）
    Y offset = get_raw(m);

    // 对所有积分值减去 offset，使得 F(m) = 0
    for (auto& y : yvals) {
      y -= offset;
    }
  }

  /**
   * @brief 查询偏移后的积分值 F(x) = ∫ₐˣ fn(t) dt - ∫ₐᵐ fn(t) dt
   */
  Y operator()(X x) const { return get_shifted(x); }

 protected:
  std::vector<X> xvals;  ///< 采样点 x 坐标（严格递增）
  std::vector<Y> yvals;  ///< 对应的积分值（已偏移）

  /**
   * @brief 使用原始（未偏移）积分表进行线性插值
   *
   * 用于计算偏移量 offset = F(m)，避免使用已偏移的数据造成循环依赖。
   */
  Y get_raw(X x) const {
    auto it = std::lower_bound(xvals.begin(), xvals.end(), x);
    // 若 x 小于所有点，返回第一个值
    if (it == xvals.begin()) {
      return yvals.front();
    }
    // 若 x 大于所有点，返回最后一个值
    if (it == xvals.end()) {
      return yvals.back();
    }

    // 找到包围 x 的两个点 [xb, xe]
    std::size_t idx_e =
        static_cast<std::size_t>(std::distance(xvals.begin(), it));
    std::size_t idx_b = idx_e - 1;

    X xb = xvals[idx_b], xe = xvals[idx_e];
    Y yb = yvals[idx_b], ye = yvals[idx_e];

    // 防止除零（理论上不会发生，因 xvals 严格递增）
    if (xe == xb) return yb;

    // 线性插值：F(x) = yb + (x - xb)/(xe - xb) * (ye - yb)
    X t = (x - xb) / (xe - xb);
    return yb + t * (ye - yb);
  }

  /**
   * @brief 使用偏移后的积分表进行线性插值（对外接口）
   */
  Y get_shifted(X x) const {
    auto it = std::lower_bound(xvals.begin(), xvals.end(), x);
    if (it == xvals.begin()) {
      return yvals.front();
    }
    if (it == xvals.end()) {
      return yvals.back();
    }

    std::size_t idx_e =
        static_cast<std::size_t>(std::distance(xvals.begin(), it));
    std::size_t idx_b = idx_e - 1;

    X xb = xvals[idx_b], xe = xvals[idx_e];
    Y yb = yvals[idx_b], ye = yvals[idx_e];

    if (xe == xb) return yb;

    X t = (x - xb) / (xe - xb);
    return yb + t * (ye - yb);
  }

  /**
   * @brief 递归生成积分点（使用梯形法则）
   *
   * 在区间 [a, b] 上递归细分，每层中点插入一次。
   * 返回该区间的积分值 ∫ₐᵇ fn(t) dt。
   *
   * @param fn 被积函数
   * @param a 区间左端点
   * @param b 区间右端点
   * @param fa fn(a)
   * @param fb fn(b)
   * @param levels 剩余递归层数
   * @return Y 区间 [a, b] 上的积分近似值
   */
  template <typename Fn>
  Y do_generate(Fn&& fn, X a, X b, Y fa, Y fb, unsigned int levels) {
    if (levels == 0) {
      // 最底层：使用梯形法则 ∫ ≈ (b - a) * (f(a) + f(b)) / 2
      X width = b - a;
      return width * (fa + fb) / static_cast<Y>(2);
    }

    // 计算中点
    X mid = (a + b) / static_cast<X>(2);
    Y fmid = fn(mid);

    // 递归处理左半区间 [a, mid]
    Y left_integral = do_generate(fn, a, mid, fa, fmid, levels - 1);
    // 插入中点：x = mid, y = ∫ₐ^mid fn(t) dt
    xvals.push_back(mid);
    yvals.push_back(left_integral);

    // 递归处理右半区间 [mid, b]
    Y right_integral = do_generate(fn, mid, b, fmid, fb, levels - 1);

    // 返回整个区间的积分值
    return left_integral + right_integral;
  }
};

/**
 * @brief 数值积分接口类（轻量包装）
 *
 * 持有一个实现类的共享指针，提供函数调用操作符。
 */
template <typename X, typename Y>
class NumericIntegralT {
 public:
  NumericIntegralT() = default;

  explicit NumericIntegralT(std::shared_ptr<NumericIntegralImplT<X, Y>> ptr)
      : ptr_(std::move(ptr)) {}

  /**
   * @brief 查询偏移后的积分值 F(x)
   */
  Y operator()(X x) const { return (*ptr_)(x); }

 private:
  std::shared_ptr<NumericIntegralImplT<X, Y>> ptr_;
};

/**
 * @brief 辅助命名空间：用于从可调用对象中推导参数和返回类型
 */
namespace _details {

/**
 * @brief 通用推导：适用于 lambda、函数对象等
 */
template <typename T>
struct DeduceArg {
  using Type = typename DeduceArg<decltype(&T::operator())>::Type;
  using RetV = typename DeduceArg<decltype(&T::operator())>::RetV;
};

/**
 * @brief 普通函数指针特化
 */
template <typename Y, typename X>
struct DeduceArg<Y (*)(X)> {
  using Type = X;
  using RetV = Y;
};

/**
 * @brief 成员函数指针特化（非 const）
 */
template <typename Class, typename Ret, typename Arg>
struct DeduceArg<Ret (Class::*)(Arg)> {
  using Type = Arg;
  using RetV = Ret;
};

/**
 * @brief 成员函数指针特化（const）
 */
template <typename Class, typename Ret, typename Arg>
struct DeduceArg<Ret (Class::*)(Arg) const> {
  using Type = Arg;
  using RetV = Ret;
};

/**
 * @brief noexcept 成员函数指针特化（非 const）
 */
template <typename Class, typename Ret, typename Arg>
struct DeduceArg<Ret (Class::*)(Arg) noexcept> {
  using Type = Arg;
  using RetV = Ret;
};

/**
 * @brief noexcept 成员函数指针特化（const）
 */
template <typename Class, typename Ret, typename Arg>
struct DeduceArg<Ret (Class::*)(Arg) const noexcept> {
  using Type = Arg;
  using RetV = Ret;
};

}  // namespace _details

/**
 * @brief 创建数值积分器的工厂函数
 *
 * 自动推导函数参数类型 X 和返回类型 Y，构建 NumericIntegralT 对象。
 *
 * @param fn 被积函数
 * @param a 积分下限
 * @param b 积分上限
 * @param m 偏移参考点（使得 F(m) = 0）
 * @param levels 递归层数（控制精度和采样点数）
 * @return NumericIntegralT<X, Y> 可调用的积分函数对象
 */
template <typename Fn,
          typename X = typename _details::DeduceArg<std::decay_t<Fn>>::Type>
auto numeric_integral(Fn&& fn, X a, X b, X m, unsigned int levels) {
  using RetVal = std::invoke_result_t<Fn, X>;
  auto impl = std::make_shared<NumericIntegralImplT<X, RetVal>>();
  impl->generate(std::forward<Fn>(fn), a, b, m, levels);
  return NumericIntegralT<X, RetVal>(std::move(impl));
}

#endif

#ifndef _SGN_H_
#define _SGN_H_

#include <cmath>

/**
 * @brief 计算数值的符号（sign function）
 *
 * 对任意类型 T（需支持与 0 比较），返回：
 * - 1 如果 val > 0
 * - 0 如果 val == 0
 * -1 如果 val < 0
 *
 * 利用布尔表达式转整数的技巧：(true) - (false) = 1 - 0 = 1
 *
 * @tparam T 数值类型（如 int, float, double 等）
 * @param val 输入值
 * @return int 符号值（-1, 0, 或 1）
 */
template <typename T>
inline int sgn(T val) {
  return (T(0) < val) - (val < T(0));
}

/**
 * @brief 计算 double 类型数值的平方
 *
 * 避免调用 std::pow(x, 2) 带来的函数调用开销，
 * 直接使用乘法，性能更高。
 *
 * @param x 输入值
 * @return double x 的平方
 */
inline double pow2(double x) { return x * x; }

/**
 * @brief 判断两个 double 浮点数是否“近似相等”
 *
 * 由于浮点数存在精度误差，不能直接用 == 比较。
 * 本函数采用“混合容差”策略：
 * - 当两个数的绝对值都较小时（< 1.0），使用**绝对误差**判断；
 * - 否则，使用**相对误差**判断，以避免大数时的精度问题。
 *
 * 公式：
 * 若 max(|a|, |b|) < 1.0：
 * |a - b| <= e
 * 否则：
 * |a - b| <= e * max(|a|, |b|)
 *
 * 此方法兼顾小数值的绝对精度和大数值的相对精度。
 *
 * @param a 第一个浮点数
 * @param b 第二个浮点数
 * @param e 容差阈值（默认 1e-8）
 * @return true 如果 a 和 b 在容差范围内相等；否则 false
 */
inline bool similar(double a, double b, double e = 1e-8) {
  // 快速路径：完全相等（包括 +inf, -inf 等特殊情况）
  if (a == b) {
    return true;
  }

  const double absA = std::abs(a);
  const double absB = std::abs(b);
  const double diff = std::abs(a - b);

  // 当两个数的绝对值都小于 1.0 时，使用绝对误差判断
  // 这能有效处理接近零的数值（如 1e-10 和 2e-10）
  if (absA < 1.0 && absB < 1.0) {
    return diff <= e;
  }

  // 否则使用相对误差：允许的误差与数值本身的大小成比例
  // 使用 max(|a|, |b|) 作为尺度，避免除法，提升性能
  return diff <= e * std::max(absA, absB);
}

#endif