# シンプソンの公式

# 被積分関数の定義
f <- function(x) {
  return(-x^4 - x^3 + 10 * x^2 + x - 2)
}

# 区間数
n <- 1000

# シンプソンの公式による [a,b] の数値積分
simpson <- function(a, b) {
  # 幅
  h <- (b - a) / n

  s1 <- f(a + h * 0.5)
  s2 <- 0
  for (i in 1 : (n - 1)) {
    x <- a + i * h
    s1 <- s1 + f(x + h * 0.5)
    s2 <- s2 + f(x)
  }
  s <- (f(a) + 4 * s1 + 2 * s2 + f(b)) * h / 6
  return(s)
}

x1 <- -3.63423216189423481381
x2 <- -0.49319311896925860372
x3 <- 0.41073330471801877684
x4 <- 2.71669197614524948747

s1 <- simpson(x1, x2)
s2 <- simpson(x2, x3)
s3 <- simpson(x3, x4)
print(sprintf("S1 = %f", s1), quote = FALSE)
print(sprintf("S2 = %f", s2), quote = FALSE)
print(sprintf("S3 = %f", s3), quote = FALSE)
print(sprintf("S = %f", abs(s1) + abs(s2) + abs(s3)), quote = FALSE)