package integral

import std.math.*
import std.util.*
foreign func malloc(size: UIntNative) : CPointer<Unit>

public func fun(x: Float64): Float64 {
  var y: Float64
  y = sqrt(4.0 - x * x)
  return y
}

public func Trape(fun: (Float64) -> Float64, down: Float64, up: Float64, n: Int32): Float64 {
  var h: Float64
  var s: Float64
  h = (up - down) / Float64(n)
  s = 0.5 * ( fun(down) + fun(down + h) ) * h 
  for (i in 1..n) {
    s = s + 0.5 * ( fun(down + Float64(i) * h) + fun(down + Float64(i + 1) * h) ) *h
  }

  return s
}

public class Trapzd {
  var a: Float64
  var b: Float64
  var s: Float64
  var n: Int32
  var f: (Float64) -> Float64
  init (aa: Float64, bb: Float64, fun: (Float64) -> Float64) {
    a = aa
	b = bb
	s = 0.0
	n = 0
	f = fun
  }
  func next(): Float64 {
    var x: Float64
	var tnm: Float64
	var sum: Float64
	var del: Float64
	var it: Int32
	n = n + 1
	if (n == 1) {
	  s = 0.5 * (b - a) * ( f(a) + f(b) )
	  return s
	} else {
	  it = 1
	  for (j in 1..n-1) {
        it = it * 2
	  }
	  tnm = Float64(it)
	  del = (b - a) / tnm
	  x = a + 0.5 * del
	  sum = 0.0
	  for (j in 0..it) {
	    sum = sum + f(x)
	    x = x + del
	  }
	  s = 0.5 * (s + (b - a) * sum / tnm )
	  return s
    }
  }
}

public func qtrap(fun: (Float64) -> Float64, a: Float64, b: Float64, eps: Float64): Float64 {
  let JMAX: Int32 = 20
  var s: Float64 = 0.0
  var old: Float64 = 0.0
  var t = Trapzd(a, b, fun)
  for (j in 0..JMAX) {
    s = t.next()
	if (j > 5) {
	  if (abs(s - old) < eps * abs(old) || (s == 0.0 && old == 0.0)) {
	    return s
	  }
	}
    old = s
  }
  throw IllegalArgumentException("Too many steps in routine qtrap")
}

public func qsimp(fun: (Float64) -> Float64, a: Float64, b: Float64, eps: Float64): Float64 {
  let JMAX: Int32 = 20
  var s: Float64
  var st: Float64
  var ost: Float64 = 0.0
  var os: Float64 = 0.0
  var t = Trapzd(a, b, fun)
  for (j in 0..JMAX) {
    st = t.next()
	s = ( 4.0 * st - ost ) / 3.0
	if (j > 5) {
	  if (abs(s - os) < eps * abs(os) || (s == 0.0 && os == 0.0)) {
	    return s
	  }
	}
	os = s
	ost = st
  }
  throw IllegalArgumentException("Too many steps in routine qsimp")
}

func qgaus(fun: (Float64) -> Float64, a: Float64, b: Float64): Float64 {
  let x = Array<Float64>([0.1488743389816312, 0.4333953941292472,
						  0.6794095682990244, 0.8650633666889845, 0.9739065285171717])
  let w = Array<Float64>([0.2955242247147529, 0.2692667193099963,
					      0.2190863625159821, 0.1494513491505806, 0.0666713443086881])
  var xm: Float64 = 0.5 * (b + a)
  var xr: Float64 = 0.5 * (b - a)
  var s: Float64 = 0.0
  for (j in 0..5) {
    var dx = xr * x[j]
	var x_r = xm + dx
	var x_l = xm - dx
	s = s + w[j] * (fun(x_r) + fun(x_l))
  }
  return s * xr
}

public func testIntegral()  {
  print("被积函数为： y = sqrt(4 - x ^ 2)\n\n")


  print("区间为[-1, 1]\n")//3.826445909962073
  var down: Float64 = -1.0
  var up: Float64 = 1.0
  var n: Int32 = 10

  var sum1 = Trape(fun, down, up, 513)
  var s = Trapzd(down, up, fun)
  var sum2: Float64
  for (i in 1..n) {
    sum2 = s.next()
  }
  var sum3 = qtrap(fun, down, up, 1e-6)
  var sum4 = qsimp(fun, down, up, 1e-6)
  var sum5 = qgaus(fun, down, up)
  print("梯形公式积分值为： ${ sum1 }\n")
  print("拓展梯形公式积分值为： ${ sum2 }\n")
  print("拓展梯形公式积分值(eps=1e-6)为： ${ sum3 }\n")
  print("辛普森公式积分值(eps=1e-6)为： ${ sum4 }\n")
  print("高斯公式积分值为： ${ sum5 }\n\n")

  print("区间为[-2, 2]\n")//6.283185307179593
  down = -2.0
  up = 2.0
  n = 10

  sum1 = Trape(fun, down, up, 513)
  s = Trapzd(down, up, fun)
  sum2 = 0.0
  for (i in 1..n) {
    sum2 = s.next()
  }
  sum3 = qtrap(fun, down, up, 1e-6)
  sum4 = qsimp(fun, down, up, 1e-6)
  sum5 = qgaus(fun, down, up)
  print("梯形公式积分值为： ${ sum1 }\n")
  print("拓展梯形公式积分值为： ${ sum2 }\n")
  print("拓展梯形公式积分值(eps=1e-6)为： ${ sum3 }\n")
  print("辛普森公式积分值(eps=1e-6)为： ${ sum4 }\n")
  print("高斯公式积分值为： ${ sum5 }\n")
}