#module NCS _n_seg, _t, _a, _b, _c, _d

/****************************************************************************
 *	グローバル関数
 ****************************************************************************/

/* NCS補間係数を作成する */
#modinit array t, array p, int n, \
	local u, local i, local m, local v, local h0, local h1, local u0, local u1
	assert n >= 1
	if n == 1 {
		_n_seg = 1
		ddim _t, 1
		ddim _a, 1
		ddim _b, 1
		ddim _c, 1
		ddim _d, 1
		_t(0) = 0.0
		_a(0) = 0.0
		_b(0) = 0.0
		_c(0) = 0.0
		_d(0) = p(0)
		return
	}
	_n_seg = n - 1
	ddim _t, _n_seg
	ddim _a, _n_seg
	ddim _b, _n_seg
	ddim _c, _n_seg
	ddim _d, _n_seg
	ddim u, _n_seg - 1
	if(_n_seg >= 2) {
		ddim m, (_n_seg - 1) * (_n_seg * 1)
		ddim v, _n_seg - 1
		for i, 0, _n_seg - 1
			h0 = t(i + 1) - t(i)
			h1 = t(i + 2) - t(i + 1)
			m(i * (_n_seg - 1) + i) = (h0 + h1) * 2.0
			if i < _n_seg - 2 {
				m(i * (_n_seg - 1) + (i + 1)) = h1
				m((i + 1) * (_n_seg - 1) + i) = h1
			}
			v(i) = (((p(i + 2) - p(i + 1)) / h1) - ((p(i + 1) - p(i)) / h0)) * 6.0
		next
		solve_simultaneous_equations@NCS m, v, u, _n_seg - 1
	}
	for i, 0, _n_seg
		if i == 0 {
			u0 = 0.0
		} else {
			u0 = u(i - 1)
		}
		if i == _n_seg - 1 {
			u1 = 0.0
		} else {
			u1 = u(i)
		}
		_t(i) = t(i)
		_a(i) = (u1 - u0) / (t(i + 1) - t(i)) / 6.0
		_b(i) = u0 / 2.0
		_c(i) = ((p(i + 1) - p(i)) / (t(i + 1) - t(i))) - ((t(i + 1) - t(i)) * (u1 + u0 * 2.0) / 6.0)
		_d(i) = p(i)
	next
	return

/* NCS補間を行った座標を取得する */
#modcfunc NCS_interpolate double t, \
	local i, local t0
	for i, 0, _n_seg - 1
		if _t(i + 1) > t {
			_break
		}
	next
	t0 = t - _t(i)
	return ((_a(i) * t0 + _b(i)) * t0 + _c(i)) * t0 + _d(i)

#modcfunc NCS_get_t int i
	return _t(i)

#modcfunc NCS_get_a int i
	return _a(i)

#modcfunc NCS_get_b int i
	return _b(i)

#modcfunc NCS_get_c int i
	return _c(i)

#modcfunc NCS_get_d int i
	return _d(i)

/****************************************************************************
 *	ローカル関数
 ****************************************************************************/

/* 連立方程式を解く */
#deffunc local solve_simultaneous_equations array A, array b, array x, int n, \
	local i, local j, local d
	d = determinant@NCS(A, n)
	for i, 0, n
		x(i) = 0.0
		for j, 0, n
			x(i) += (cofactor@NCS(A, n, j, i) / d) * b(j)
		next
	next
	return

/* 行列式を求める */
#defcfunc local determinant array m, int n, \
	local i, local d
	assert n >= 1
	d = 0.0
	for i, 0, n
		d += m((i * n) + 0) * cofactor@NCS(m, n, i, 0)
	next
	return d

/* 余因子を求める */
#defcfunc local cofactor array m, int n, int i, int j, \
	local k, local l, local d, local s, local t, local u
	assert n >= 1
	assert (i >= 0) and (i <= n - 1)
	assert (j >= 0) and (j <= n - 1)
	if n == 1 {
		return 1.0
	}
	ddim s, (n - 1) * (n - 1)
	t = 0
	u = 0
	for k, 0, n
		for l, 0, n
			if (k != i) and (l != j) {
				s(t) = m(u)
				t++
			}
			u++
		next
	next
	d = determinant@NCS(s, (n - 1))
	if (i ^ j) & 1 {
		d = -d
	}
	return d

#global

