Wat

// реализация обобщенной составной ИКФ
trait IntegralMethod {
  // "шаг": ИКФ на заданном отрезке
  def methodStep(x1 : Double, x2: Double): Double

  // АСТ метода (3 для Ньютона-Котса, 5 для Гаусса)
  def m: Int

  // составная ИКФ с заданным шагом step на отрезке [x1; x2]
  def method(x1: Double, x2: Double, step: Double): Double = {
    var curX = x1
    var result = 0.0
    while(curX + step < x2) {
      result += methodStep(curX, curX + step)
      curX += step
    }
    result + methodStep(curX, x2)
  }

  // поиск оптимального шага на двух исходных шагах step1 и step2
  def optimalStep(epsilon: Double, step1: Double, step2: Double): Double = {
    val s1 = method(a, b, step1) // значение интеграла при шаге step1
    val s2 = method(a, b, step2) // при шаге step2
    val error = rungeError(s1, s2, step1 / step2)
    // вычисление оптимального шага по формуле из методички
    step2 * pow(epsilon / error, 1.0 / m)
  }

  // ошибка по Рунге (см. методичку)
  def rungeError(s1: Double, s2: Double, L: Double): Double = abs((s2 - s1) / pow(L, m))

  // вычисление интеграла с точностью epsilon, исходным шагом initStep и уменьшением шага в L раз
  // метод возвращает результат, оптимальный шаг и число итераций (последнее хз, зачем, это необязательно вроде)
  def apply(epsilon: Double, initStep: Double, L: Double): (Double, Double, Int) = {
    var curValue = method(a, b, initStep)
    var prevValue = 0.0
    var curStep = initStep / L
    var iterations = 0
    do {
      prevValue = curValue
      curValue = method(a, b, curStep)
      curStep /= L
      iterations += 1
    } while(rungeError(prevValue, curValue, L) > epsilon) // цикл продолжается, пока ошибка по Рунге больше epsilon
    (curValue, curStep * L, iterations)
  }
}


// Ньютон-Котс
object NewtonCatsMethod extends IntegralMethod {
  // АСТ == 3
  override val m: Int = 3

  def newtonMoments(x1: Double, x2: Double): (Double, Double, Double) = (newM0(x1, x2), newM1(x1, x2), newM2(x1, x2))

  override def methodStep(x1: Double, x2: Double): Double = {
    val (m0, m1, m2) = newtonMoments(x1, x2)
    val xMid = (x1 + x2) / 2
    val A1 = (m2 - m1 * (xMid + x2) + m0 * xMid * x2) / (xMid - x1) / (x2 - x1)
    val A2 = -(m2 - m1 * (x1 + x2) + m0 * x1 * x2) / (xMid - x1) / (x2 - xMid)
    val A3 = (m2 - m1 * (xMid + x1) + m0 * xMid * x1) / (x2 - xMid) / (x2 - x1)
    A1 * f(x1) + A2 * f(xMid) + A3 * f(x2)
  }
}

// Гаусс
object GaussMethod extends IntegralMethod {
  type GaussMoments = (Double, Double, Double, Double, Double, Double)
  val solver = new LaguerreSolver

  // АСТ == 6
  override val m: Int = 6

  def gaussMoments(x1: Double, x2: Double): GaussMoments =
    (newM0(x1, x2), newM1(x1, x2), newM2(x1, x2), newM3(x1, x2), newM4(x1, x2), newM5(x1, x2))

  // решение кубического уравнения (в numpy вроде есть функция типа np.roots или np.solve для этого)
  def solveCubicEquation(a0: Double, a1: Double, a2: Double, a3: Double): (Double, Double, Double) = {
    val coeffs = Array(a3,a2,a1,a0)
    val roots = solver.solveAllComplex(coeffs, 0)
    assert(roots.forall(_.getImaginary == 0))
    val Array(r1, r2, r3) = roots.map(_.getReal)
    (r1,r2,r3)
  }

  // поиск необходимых коэффициентов (см. методичку)
  def findPolyCoeffs(x1: Double, x2: Double): (Double, Double, Double) = {
    val (m0,m1,m2,m3,m4,m5) = gaussMoments(x1, x2)
    val col = new ColumnMatrix(Array(-m3,-m4,-m5))
    val matrix: SquareMatrix =
      m0 c m1 c m2 r
      m1 c m2 c m3 r
      m2 c m3 c m4
    val solution = SLAEUtils.solveNonDegSLAE(matrix, col)
    (solution(0), solution(1), solution(2))
  }

  // поиск коэффициентов, необходимых для вычисления интеграла (см. методичку)
  def findCoeffs(x1: Double, x2: Double): ((Double, Double, Double), (Double, Double, Double)) = {
    val (a0, a1, a2) = findPolyCoeffs(x1, x2)
    val (r1, r2, r3) = solveCubicEquation(1, a2, a1, a0)
    val (m0, m1, m2, _, _, _) = gaussMoments(x1, x2)
    val col = new ColumnMatrix(Array(m0,m1,m2))
    val matrix: SquareMatrix =
      1         c   1         c   1       r
      r1        c   r2        c   r3      r
      r1 * r1   c   r2 * r2   c   r3 * r3
    val sol = SLAEUtils.solveNonDegSLAE(matrix, col)
    ((sol(0), sol(1), sol(2)), (r1, r2, r3))
  }

  override def methodStep(x1: Double, x2: Double): Double = {
    val ((a1, a2, a3), (r1, r2, r3)) = findCoeffs(x1, x2)
    a1 * f(r1) + a2 * f(r2) + a3 * f(r3)
  }
}

// ...

// вспомогательные штучки
object TaskFunctions {
  def f(x: Double) = ...
  val a = ...
  val b = ...
  val alpha = ...
  val beta = ...

  val exactResult = ...

  // моменты порядков от 0 до 5 включительно, для случая, когда beta = 0, alpha != 0
  // если у тебя beta != 0, alpha = 0, можно либо вывести как-нибудь их аналогично,
  // либо посчитать интегралы руками,
  // либо своровать у кого-нибудь (:
  def newM0(x1: Double, x2: Double): Double =
    (pow(x2 - a, 1 - alpha) - pow(x1 - a, 1 - alpha)) / (1 - alpha)
  def newM1(x1: Double, x2: Double): Double =
    (pow(x2 - a, 2 - alpha) - pow(x1 - a, 2 - alpha)) / (2 - alpha) + a * newM0(x1, x2)
  def newM2(x1: Double, x2: Double): Double =
    (pow(x2 - a, 3 - alpha) - pow(x1 - a, 3 - alpha)) / (3 - alpha) + 2 * a * newM1(x1, x2) - a * a * newM0(x1, x2)
  def newM3(x1: Double, x2: Double): Double =
    (pow(x2 - a, 4 - alpha) - pow(x1 - a, 4 - alpha)) / (4 - alpha) + 3 * a * newM2(x1, x2) - 3 * a * a * newM1(x1, x2) + pow(a,3) * newM0(x1, x2)
  def newM4(x1: Double, x2: Double): Double =
    (pow(x2 - a, 5 - alpha) - pow(x1 - a, 5 - alpha)) / (5 - alpha) + 4 * a * newM3(x1, x2) - 6 * a * a * newM2(x1, x2) + 4 * pow(a,3) * newM1(x1, x2) - pow(a, 4) * newM0(x1, x2)
  def newM5(x1: Double, x2: Double): Double =
    (pow(x2 - a, 6 - alpha) - pow(x1 - a, 6 - alpha)) / (6 - alpha) +
    5 * a * newM4(x1, x2) - 10 * a * a * newM3(x1, x2) + 10 * a * a * a * newM2(x1, x2) -
    5 * a * a * a * a * newM1(x1, x2) + a * a * a * a * a * newM0(x1, x2)


  def newtonMoments(x1: Double, x2: Double): (Double, Double, Double) = (newM0(x1, x2), newM1(x1, x2), newM2(x1, x2))
}