func (lsf *LeastSquaresFunction) Eval(x float64) float64 {
	slope := lsf.slope()
	result := lsf.intercept(slope) + slope*x
	return result
}
Eval evaluates the function for a given X value and
returns the Y value.
package numerics
import (
	"fmt"
	"math"
	"sort"
)
type Point struct {
	x float64
	y float64
}
func NewPoint(x, y float64) *Point {
	return &Point{x, y}
}
func (p Point) IsInf() bool {
	return math.IsInf(p.x, 0) || math.IsInf(p.y, 0)
}
func (p Point) IsNaN() bool {
	return math.IsNaN(p.x) || math.IsNaN(p.y)
}
func (p Point) X() float64 {
	return p.x
}
func (p Point) Y() float64 {
	return p.y
}
func (p Point) DistanceTo(op *Point) float64 {
	dx := p.x - op.x
	dy := p.y - op.y
	return math.Sqrt(dx*dx + dy*dy)
}
func (p Point) VectorTo(op *Point) *Vector {
	return NewVector(op.X()-p.x, op.Y()-p.y)
}
func (p Point) String() string {
	return fmt.Sprintf("(%f, %f)", p.x, p.y)
}
func MiddlePoint(a, b *Point) *Point {
	return NewPoint((a.x+b.x)/2, (a.y+b.y)/2)
}
func PointVector(a, b *Point) *Vector {
	return a.VectorTo(b)
}
type Points []*Point
func NewPoints(p ...*Point) Points {
	if len(p) > 0 {
		return p
	}
	return []*Point{}
}
func (ps Points) XAt(idx int) float64 {
	return ps[idx].X()
}
func (ps Points) YAt(idx int) float64 {
	return ps[idx].Y()
}
func (ps Points) XDifference(idxA, idxB int) float64 {
	return ps[idxA].X() - ps[idxB].X()
}
func (ps Points) YDifference(idxA, idxB int) float64 {
	return ps[idxA].Y() - ps[idxB].Y()
}
func (ps Points) XInRange(x float64) bool {
	minX := ps[0].X()
	maxX := ps[0].X()
	for _, p := range ps[1:] {
		if p.X() < minX {
			minX = p.X()
		}
		if p.X() > maxX {
			maxX = p.X()
		}
	}
	return minX <= x && x <= maxX
}
func (ps Points) SearchNextIndex(x float64) int {
	sf := func(i int) bool {
		return x < ps[i].X()
	}
	return sort.Search(len(ps), sf)
}
func (ps Points) Len() int {
	return len(ps)
}
func (ps Points) Less(i, j int) bool {
	switch {
	case ps[i].x < ps[j].x:
		return true
	case ps[i].x > ps[j].x:
		return false
	}
	switch {
	case ps[i].y < ps[j].y:
		return true
	case ps[i].y > ps[j].y:
		return false
	}
	return false
}
func (ps Points) Swap(i, j int) {
	ps[i], ps[j] = ps[j], ps[i]
}
func (ps Points) CubicSplineFunction() *CubicSplineFunction {
	return NewCubicSplineFunction(ps)
}
func (ps Points) LeastSquaresFunction() *LeastSquaresFunction {
	return NewLeastSquaresFunction(ps)
}
func (ps Points) String() string {
	pss := "{"
	for _, p := range ps {
		pss += p.String()
	}
	pss += "}"
	return pss
}
type Vector struct {
	x float64
	y float64
}
func NewVector(x, y float64) *Vector {
	return &Vector{x, y}
}
func (v Vector) X() float64 {
	return v.x
}
func (v Vector) Y() float64 {
	return v.y
}
func (v Vector) Len() float64 {
	return math.Sqrt(v.x*v.x + v.y*v.y)
}
func (v Vector) String() string {
	return fmt.Sprintf("<%f, %f>", v.x, v.y)
}
func AddVectors(a, b *Vector) *Vector {
	return NewVector(a.x+b.x, a.y+b.y)
}
func SubVectors(a, b *Vector) *Vector {
	return NewVector(a.x-b.x, a.y-b.y)
}
func ScaleVector(v *Vector, s float64) *Vector {
	return NewVector(v.x*s, v.y*s)
}
type Function interface {
	Eval(x float64) float64
	EvalPoint(x float64) *Point
	EvalPoints(fromX, toX float64, count int) Points
}
type PolynomialFunction struct {
	coefficients []float64
}
func NewPolynomialFunction(coefficients []float64) *PolynomialFunction {
	if len(coefficients) < 1 {
		return nil
	}
	pf := &PolynomialFunction{
		coefficients: coefficients,
	}
	return pf
}
func (pf PolynomialFunction) Eval(x float64) float64 {
	n := len(pf.coefficients)
	result := pf.coefficients[n-1]
	for i := n - 2; i >= 0; i-- {
		result = x*result + pf.coefficients[i]
	}
	return result
}
func (pf PolynomialFunction) EvalPoint(x float64) *Point {
	return NewPoint(x, pf.Eval(x))
}
func (pf PolynomialFunction) EvalPoints(fromX, toX float64, count int) Points {
	return evalPoints(pf, fromX, toX, count)
}
func (pf PolynomialFunction) Differentiate() *PolynomialFunction {
	n := len(pf.coefficients)
	if n == 1 {
		return NewPolynomialFunction([]float64{0.0})
	}
	newCoefficients := make([]float64, n-1)
	for i := n - 1; i > 0; i-- {
		newCoefficients[i-1] = float64(i) * pf.coefficients[i]
	}
	return NewPolynomialFunction(newCoefficients)
}
func (pf PolynomialFunction) String() string {
	pfs := "f(x) := "
	for i := len(pf.coefficients) - 1; i > 0; i-- {
		if pf.coefficients[i] != 0.0 {
			pfs += fmt.Sprintf("%vx", pf.coefficients[i])
			if i > 1 {
				pfs += fmt.Sprintf("^%v", i)
			}
			if pf.coefficients[i-1] > 0 {
				pfs += "+"
			}
		}
	}
	if pf.coefficients[0] != 0.0 {
		pfs += fmt.Sprintf("%v", pf.coefficients[0])
	}
	return pfs
}
type CubicSplineFunction struct {
	polynomials []*PolynomialFunction
	points      Points
}
func NewCubicSplineFunction(points Points) *CubicSplineFunction {
	if points.Len() < 3 {
		return nil
	}
	csf := &CubicSplineFunction{
		points: points,
	}
	intervals := points.Len() - 1
	differences := make([]float64, intervals)
	for i := 0; i < intervals; i++ {
		differences[i] = points[i+1].X() - points[i].X()
	}
	mu := make([]float64, intervals)
	z := make([]float64, points.Len())
	var g float64
	for i := 1; i < intervals; i++ {
		g = 2.0*points.XDifference(i+1, i-1) - differences[i-1]*mu[i-1]
		mu[i] = differences[i] / g
		z[i] = (3.0*(points.YAt(i+1)*differences[i-1]-points.YAt(i)*
			points.XDifference(i+1, i-1)+points.YAt(i-1)*differences[i])/
			(differences[i-1]*differences[i]) - differences[i-1]*z[i-1]) / g
	}
	b := make([]float64, intervals)
	c := make([]float64, points.Len())
	d := make([]float64, intervals)
	for i := intervals - 1; i >= 0; i-- {
		c[i] = z[i] - mu[i]*c[i+1]
		b[i] = points.YDifference(i+1, i)/differences[i] - differences[i]*(c[i+1]+2.0*c[i])/3.0
		d[i] = (c[i+1] - c[i]) / (3.0 * differences[i])
	}
	csf.polynomials = make([]*PolynomialFunction, intervals)
	coefficients := make([]float64, 4)
	for i := 0; i < intervals; i++ {
		coefficients[0] = points.YAt(i)
		coefficients[1] = b[i]
		coefficients[2] = c[i]
		coefficients[3] = d[i]
		csf.polynomials[i] = NewPolynomialFunction(coefficients)
	}
	return csf
}
func (csf *CubicSplineFunction) Eval(x float64) float64 {
	if !csf.points.XInRange(x) {
		panic("X out of range!")
	}
	idx := csf.points.SearchNextIndex(x)
	if idx >= len(csf.polynomials) {
		idx = len(csf.polynomials) - 1
	}
	return csf.polynomials[idx].Eval(x - csf.points.XAt(idx))
}
func (csf *CubicSplineFunction) EvalPoint(x float64) *Point {
	return NewPoint(x, csf.Eval(x))
}
func (csf *CubicSplineFunction) EvalPoints(fromX, toX float64, count int) Points {
	return evalPoints(csf, fromX, toX, count)
}
type LeastSquaresFunction struct {
	sumX, sumXX float64
	sumY, sumYY float64
	sumXY       float64
	xBar, yBar  float64
	count       int
}
func NewLeastSquaresFunction(points Points) *LeastSquaresFunction {
	lsf := new(LeastSquaresFunction)
	if points != nil {
		lsf.AppendPoints(points)
	}
	return lsf
}
func (lsf *LeastSquaresFunction) AppendPoint(x, y float64) {
	p := NewPoint(x, y)
	if lsf.count == 0 {
		lsf.xBar = p.X()
		lsf.yBar = p.Y()
	} else {
		dx := p.X() - lsf.xBar
		dy := p.Y() - lsf.yBar
		lsf.sumXX += dx * dx * float64(lsf.count) / float64(lsf.count+1.0)
		lsf.sumYY += dy * dy * float64(lsf.count) / float64(lsf.count+1.0)
		lsf.sumXY += dx * dy * float64(lsf.count) / float64(lsf.count+1.0)
		lsf.xBar += dx / float64(lsf.count+1.0)
		lsf.yBar += dy / float64(lsf.count+1.0)
	}
	lsf.sumX += p.X()
	lsf.sumY += p.Y()
	lsf.count++
}
func (lsf *LeastSquaresFunction) AppendPoints(points Points) {
	for _, p := range points {
		lsf.AppendPoint(p.X(), p.Y())
	}
}