class NewtonInterpolation {
  private x: number[];
  private y: number[];
  private coefficients: number[];

  constructor(x: number[], y: number[]) {
    if (x.length !== y.length) {
      throw new Error('x and y arrays must have the same length');
    }
    this.x = x;
    this.y = y;
    this.coefficients = this.computeCoefficients();
  }

  // 计算差商表
  private computeCoefficients(): number[] {
    const n = this.x.length;
    const f = Array.from({ length: n }, () => Array(n).fill(0));

    // 初始化第一列
    for (let i = 0; i < n; i++) {
      f[i][0] = this.y[i];
    }

    // 计算差商
    for (let j = 1; j < n; j++) {
      for (let i = 0; i < n - j; i++) {
        f[i][j] = (f[i + 1][j - 1] - f[i][j - 1]) / (this.x[i + j] - this.x[i]);
      }
    }

    // 提取对角线元素作为系数
    const coefficients: number[] = [];
    for (let i = 0; i < n; i++) {
      coefficients.push(f[0][i]);
    }

    return coefficients;
  }

  // 插值计算
  interpolate(value: number): number {
    let result = this.coefficients[0];
    let product = 1;

    for (let i = 1; i < this.coefficients.length; i++) {
      product *= value - this.x[i - 1];
      result += this.coefficients[i] * product;
    }

    return result;
  }

  // 获取系数
  getCoefficients(): number[] {
    return [...this.coefficients];
  }
}

// 使用示例
function newtonInterpolationExample() {
  const x = [1, 2, 3, 4, 5];
  const y = [1, 4, 9, 16, 25]; // y = x^2

  const interpolator = new NewtonInterpolation(x, y);
  console.log('Coefficients:', interpolator.getCoefficients());

  // 插值 x = 2.5
  const result = interpolator.interpolate(2.5);
  console.log('f(2.5) =', result); // 应该接近 6.25
}

// 龙格-库塔法求解常微分方程
function rungeKutta4(f: (t: number, y: number) => number, t0: number, y0: number, tEnd: number, h: number): { t: number[]; y: number[] } {
  const t: number[] = [];
  const y: number[] = [];
  let currentT = t0;
  let currentY = y0;

  t.push(currentT);
  y.push(currentY);

  while (currentT < tEnd) {
    // 确保不超过结束时间
    const step = Math.min(h, tEnd - currentT);

    // 计算四个斜率
    const k1 = f(currentT, currentY);
    const k2 = f(currentT + step / 2, currentY + (step * k1) / 2);
    const k3 = f(currentT + step / 2, currentY + (step * k2) / 2);
    const k4 = f(currentT + step, currentY + step * k3);

    // 更新
    currentY = currentY + (step / 6) * (k1 + 2 * k2 + 2 * k3 + k4);
    currentT = currentT + step;

    t.push(currentT);
    y.push(currentY);
  }

  return { t, y };
}

// 自适应步长的龙格-库塔法
function adaptiveRungeKutta(f: (t: number, y: number) => number, t0: number, y0: number, tEnd: number, initialH: number, tolerance: number = 1e-6): { t: number[]; y: number[] } {
  const t: number[] = [t0];
  const y: number[] = [y0];
  let currentT = t0;
  let currentY = y0;
  let h = initialH;

  while (currentT < tEnd) {
    // 确保不超过结束时间
    h = Math.min(h, tEnd - currentT);

    // 一步法
    const k1 = f(currentT, currentY);
    const k2 = f(currentT + h / 2, currentY + (h * k1) / 2);
    const k3 = f(currentT + h / 2, currentY + (h * k2) / 2);
    const k4 = f(currentT + h, currentY + h * k3);
    const y1 = currentY + (h / 6) * (k1 + 2 * k2 + 2 * k3 + k4);

    // 两步法
    const h2 = h / 2;
    const k1a = f(currentT, currentY);
    const k2a = f(currentT + h2 / 2, currentY + (h2 * k1a) / 2);
    const k3a = f(currentT + h2 / 2, currentY + (h2 * k2a) / 2);
    const k4a = f(currentT + h2, currentY + h2 * k3a);
    const y2a = currentY + (h2 / 6) * (k1a + 2 * k2a + 2 * k3a + k4a);

    const k1b = f(currentT + h2, y2a);
    const k2b = f(currentT + h2 + h2 / 2, y2a + (h2 * k1b) / 2);
    const k3b = f(currentT + h2 + h2 / 2, y2a + (h2 * k2b) / 2);
    const k4b = f(currentT + h, y2a + h2 * k3b);
    const y2 = y2a + (h2 / 6) * (k1b + 2 * k2b + 2 * k3b + k4b);

    // 计算误差
    const error = Math.abs(y2 - y1);

    if (error < tolerance) {
      // 接受这一步
      currentT += h;
      currentY = y2;
      t.push(currentT);
      y.push(currentY);

      // 增加步长
      h *= Math.min(2, 0.9 * Math.pow(tolerance / error, 1 / 4));
    } else {
      // 拒绝这一步，减少步长
      h *= Math.max(0.1, 0.9 * Math.pow(tolerance / error, 1 / 5));
    }
  }

  return { t, y };
}

// 使用示例：求解 dy/dt = -2t*y, y(0) = 1
function odeExample() {
  const f = (t: number, y: number) => -2 * t * y;
  const result = rungeKutta4(f, 0, 1, 2, 0.1);

  console.log('t values:', result.t);
  console.log('y values:', result.y);

  // 解析解：y = e^(-t^2)
  const analytical = result.t.map((t) => Math.exp(-t * t));
  console.log('Analytical:', analytical);
}

// 牛顿法求方程根
function newtonMethod(
  f: (x: number) => number,
  df: (x: number) => number,
  initialGuess: number,
  tolerance: number = 1e-6,
  maxIterations: number = 100,
): { root: number; iterations: number; converged: boolean } {
  let x = initialGuess;
  let iterations = 0;

  while (iterations < maxIterations) {
    const fx = f(x);
    const dfx = df(x);

    if (Math.abs(dfx) < 1e-10) {
      // 导数太小，可能陷入局部极值
      return { root: x, iterations, converged: false };
    }

    const xNew = x - fx / dfx;
    iterations++;

    if (Math.abs(xNew - x) < tolerance) {
      return { root: xNew, iterations, converged: true };
    }

    x = xNew;
  }

  return { root: x, iterations, converged: false };
}

// 简化牛顿法（数值导数）
function newtonMethodNumerical(
  f: (x: number) => number,
  initialGuess: number,
  h: number = 1e-5,
  tolerance: number = 1e-6,
  maxIterations: number = 100,
): { root: number; iterations: number; converged: boolean } {
  const df = (x: number) => (f(x + h) - f(x - h)) / (2 * h);
  return newtonMethod(f, df, initialGuess, tolerance, maxIterations);
}

// 二分法求根（作为对比）
function bisectionMethod(f: (x: number) => number, a: number, b: number, tolerance: number = 1e-6, maxIterations: number = 100): { root: number; iterations: number; converged: boolean } {
  if (f(a) * f(b) >= 0) {
    throw new Error('Function must have opposite signs at endpoints');
  }

  let iterations = 0;
  while (iterations < maxIterations) {
    const c = (a + b) / 2;
    const fc = f(c);

    if (Math.abs(fc) < tolerance || Math.abs(b - a) < tolerance) {
      return { root: c, iterations, converged: true };
    }

    if (f(a) * fc < 0) {
      b = c;
    } else {
      a = c;
    }

    iterations++;
  }

  return { root: (a + b) / 2, iterations, converged: false };
}

// 使用示例：求解 x^3 - 2x - 5 = 0
function rootFindingExample() {
  const f = (x: number) => x * x * x - 2 * x - 5;
  const df = (x: number) => 3 * x * x - 2;

  // 牛顿法
  const newtonResult = newtonMethod(f, df, 3);
  console.log('Newton method result:', newtonResult);

  // 二分法
  const bisectionResult = bisectionMethod(f, 2, 3);
  console.log('Bisection method result:', bisectionResult);
}

// 梯形法则
function trapezoidalRule(f: (x: number) => number, a: number, b: number, n: number): number {
  const h = (b - a) / n;
  let sum = (f(a) + f(b)) / 2;

  for (let i = 1; i < n; i++) {
    sum += f(a + i * h);
  }

  return sum * h;
}

// 辛普森法则
function simpsonRule(f: (x: number) => number, a: number, b: number, n: number): number {
  if (n % 2 !== 0) {
    throw new Error('n must be even for Simpson rule');
  }

  const h = (b - a) / n;
  let sum = f(a) + f(b);

  for (let i = 1; i < n; i++) {
    const x = a + i * h;
    if (i % 2 === 0) {
      sum += 2 * f(x);
    } else {
      sum += 4 * f(x);
    }
  }

  return (sum * h) / 3;
}

// 自适应积分
function adaptiveIntegration(f: (x: number) => number, a: number, b: number, tolerance: number = 1e-6): number {
  const h = b - a;
  const mid = (a + b) / 2;

  // 一步积分
  const singleStep = ((f(a) + f(b)) * h) / 2;

  // 两步积分
  const twoStep = ((f(a) + 4 * f(mid) + f(b)) * h) / 6;

  const error = Math.abs(twoStep - singleStep);

  if (error < tolerance) {
    return twoStep;
  } else {
    // 递归细分
    return adaptiveIntegration(f, a, mid, tolerance / 2) + adaptiveIntegration(f, mid, b, tolerance / 2);
  }
}

// 使用示例：计算 ∫(0 to 1) x^2 dx = 1/3
function integrationExample() {
  const f = (x: number) => x * x;

  const trapezoidal = trapezoidalRule(f, 0, 1, 100);
  console.log('Trapezoidal rule:', trapezoidal);

  const simpson = simpsonRule(f, 0, 1, 100);
  console.log('Simpson rule:', simpson);

  const adaptive = adaptiveIntegration(f, 0, 1);
  console.log('Adaptive integration:', adaptive);

  console.log('Exact value: 1/3 =', 1 / 3);
}

class Matrix {
  private data: number[][];
  public rows: number;
  public cols: number;

  constructor(data: number[][]) {
    this.data = data;
    this.rows = data.length;
    this.cols = data[0].length;
  }

  // 矩阵加法
  add(other: Matrix): Matrix {
    if (this.rows !== other.rows || this.cols !== other.cols) {
      throw new Error('Matrix dimensions must match');
    }

    const result = [];
    for (let i = 0; i < this.rows; i++) {
      result[i] = [];
      for (let j = 0; j < this.cols; j++) {
        result[i][j] = this.data[i][j] + other.data[i][j];
      }
    }

    return new Matrix(result);
  }

  // 矩阵乘法
  multiply(other: Matrix | number): Matrix {
    if (typeof other === 'number') {
      // 标量乘法
      const result = [];
      for (let i = 0; i < this.rows; i++) {
        result[i] = [];
        for (let j = 0; j < this.cols; j++) {
          result[i][j] = this.data[i][j] * other;
        }
      }
      return new Matrix(result);
    } else {
      // 矩阵乘法
      if (this.cols !== other.rows) {
        throw new Error('Matrix dimensions incompatible for multiplication');
      }

      const result = [];
      for (let i = 0; i < this.rows; i++) {
        result[i] = [];
        for (let j = 0; j < other.cols; j++) {
          let sum = 0;
          for (let k = 0; k < this.cols; k++) {
            sum += this.data[i][k] * other.data[k][j];
          }
          result[i][j] = sum;
        }
      }

      return new Matrix(result);
    }
  }

  // 高斯消元法求逆矩阵
  inverse(): Matrix {
    if (this.rows !== this.cols) {
      throw new Error('Matrix must be square');
    }

    const n = this.rows;
    const augmented = [];

    // 创建增广矩阵 [A | I]
    for (let i = 0; i < n; i++) {
      augmented[i] = [...this.data[i]];
      for (let j = 0; j < n; j++) {
        augmented[i].push(i === j ? 1 : 0);
      }
    }

    // 高斯消元
    for (let i = 0; i < n; i++) {
      // 找到主元
      let maxRow = i;
      for (let k = i + 1; k < n; k++) {
        if (Math.abs(augmented[k][i]) > Math.abs(augmented[maxRow][i])) {
          maxRow = k;
        }
      }

      // 交换行
      [augmented[i], augmented[maxRow]] = [augmented[maxRow], augmented[i]];

      // 消元
      for (let k = i + 1; k < n; k++) {
        const factor = augmented[k][i] / augmented[i][i];
        for (let j = i; j < 2 * n; j++) {
          augmented[k][j] -= factor * augmented[i][j];
        }
      }
    }

    // 回代
    for (let i = n - 1; i >= 0; i--) {
      for (let k = i - 1; k >= 0; k--) {
        const factor = augmented[k][i] / augmented[i][i];
        for (let j = i; j < 2 * n; j++) {
          augmented[k][j] -= factor * augmented[i][j];
        }
      }
    }

    // 归一化主对角线
    for (let i = 0; i < n; i++) {
      const divisor = augmented[i][i];
      for (let j = i; j < 2 * n; j++) {
        augmented[i][j] /= divisor;
      }
    }

    // 提取逆矩阵
    const inverseData = [];
    for (let i = 0; i < n; i++) {
      inverseData[i] = augmented[i].slice(n);
    }

    return new Matrix(inverseData);
  }

  // 矩阵转置
  transpose(): Matrix {
    const result = [];
    for (let j = 0; j < this.cols; j++) {
      result[j] = [];
      for (let i = 0; i < this.rows; i++) {
        result[j][i] = this.data[i][j];
      }
    }
    return new Matrix(result);
  }

  toString(): string {
    return this.data.map((row) => row.join(' ')).join('\n');
  }
}