```import java.util.*;
import java.awt.*;

// a cubic spline to be drawn through a vector of VPoint's
//
// Author:      Alexander Bogomolny, CTK Software, Inc.
// URL:         http://www.cut-the-knot.com
// Date:        November 30, 2000
//              Permission to use and modify the file is therefore granted
//              as long as this comment remains unchanged. Do this at your
//              own risk.
//
class CubicSpline
{
int N;
Vector P;           // N+1 points to interpolate through (screen coordinates)

public void SetMoveables(Vector mp)
{
N = mp.size()-1;
P = mp;
}

public Vector GetMoveables() { return P; }

double[] M;         // (N-1) second derivatives at "inner" points
double[] In;        // N integrals over N intervals

private int I = 0;  // points to the current interval where value is computed

public CubicSpline(Vector mp)
{
try
{
N = mp.size()-1;
P = mp;
}
catch (Exception e)
{
N = 0;
P = new Vector();
}
}

// solves for the vector of
// second derivatives
//
public void Resolve()
{
if (N < 2)
return;

double[] h = new double[N];
for (int i = 0; i < N; i++)
h[i] = ((VPoint)P.elementAt(i+1)).ScaledX() - ((VPoint)P.elementAt(i)).ScaledX();

double[] d = new double[N-1];
for (int i = 0; i < N-1; i++)
d[i] = (h[i] + h[i+1])/3;

double[] m = new double[N-2];
double[] p = new double[N-2];
for (int i = 0; i < N-2; i++)
m[i] = p[i] = h[i+1]/6;

Tridiagonal tr = new Tridiagonal(N-1, m, d, p);

double[] dy = new double[N];
for (int i = 0; i < N; i++)
dy[i] = ((VPoint)P.elementAt(i+1)).ScaledY() - ((VPoint)P.elementAt(i)).ScaledY();

double[] r = new double[N-1];
for (int i = 0; i < N-1; i++)
r[i] = dy[i+1]/h[i+1] - dy[i]/h[i];

M = tr.Solve(r);
}

public double Value(double x)
{
I = FindInterval(x);
return VPoint.ReverseY(Compute(I, VPoint.ScaledX(x)));
}

private int FindInterval(double x)
{
int Lo = 0;
int Hi = N;

do
{
int i = (Lo + Hi)/2;

double X = ((VPoint)P.elementAt(i)).x;

if (X > x)
Hi = i;
else
Lo = i;

}
while (Hi - Lo > 1);

return Lo;
}

// x comes in screen coordinates
//
private double Compute(int index, double x)
{
double h = ((VPoint)P.elementAt(index+1)).ScaledX() - ((VPoint)P.elementAt(index)).ScaledX();
double A = (((VPoint)P.elementAt(index+1)).ScaledX() - x)/h;
double B = (x - ((VPoint)P.elementAt(index)).ScaledX())/h;

double yA = ((VPoint)P.elementAt(index)).ScaledY();
double yB = ((VPoint)P.elementAt(index+1)).ScaledY();

double mA = index == 0 ? 0 : M[index-1];
double mB = index > N-2 ? 0 : M[index];

return A*yA + B*yB + ((A*A*A - A)*mA + (B*B*B - B)*mB)*h*h/6;
}

public double Derivative(double x)
{
I = FindInterval(x);
return VPoint.ReverseY(ComputeDerivative(I, VPoint.ScaledX(x)));
}

private double ComputeDerivative(int index, double x)
{
double h = ((VPoint)P.elementAt(index+1)).ScaledX() - ((VPoint)P.elementAt(index)).ScaledX();
double A = (((VPoint)P.elementAt(index+1)).ScaledX() - x)/h;
double B = (x - ((VPoint)P.elementAt(index)).ScaledX())/h;

double yA = ((VPoint)P.elementAt(index)).ScaledY();
double yB = ((VPoint)P.elementAt(index+1)).ScaledY();

double mA = index == 0 ? 0 : M[index-1];
double mB = index > N-2 ? 0 : M[index];

return (yB - yA)/h - (mB - mA)*h/6 + (mB*B*B - mA*A*A)*h/2;
}

public void Draw(Graphics g, int xFrom, int xTo, int dx)
{
int x = xFrom;
Point p = new Point(x, (int)Math.round(Value(x)));
int y = (int)Math.round(VPoint.ReverseY(0));

g.drawLine(x, y, x, p.y);

while (x < xTo)
{
x += dx;
Point q = new Point(x, (int)Math.round(Value(x)));

g.drawLine(p.x, p.y, q.x, q.y);

g.drawLine(x, y, x, q.y);

p = q;
}
}

public void DrawDerivative(Graphics g, int xFrom, int xTo, int dx)
{
int x = xFrom;
Point p = new Point(x, (int)Math.round(Derivative(x)));

while (x < xTo)
{
x += dx;
Point q = new Point(x, (int)Math.round(Derivative(x)));

g.drawLine(p.x, p.y, q.x, q.y);
p = q;
}
}

// saves time when computing integrals
// for a variable x
//
public void FillIntegralArray()
{
double[] h = new double[N];
for (int i = 0; i < N; i++)
h[i] = ((VPoint)P.elementAt(i+1)).ScaledX() - ((VPoint)P.elementAt(i)).ScaledX();

In = new double[N];

for (int i = 0; i < N; i++)
{
double mA = i == 0 ? 0 : M[i-1];
double mB = i > N-2 ? 0 : M[i];
double yA = ((VPoint)P.elementAt(i)).ScaledY();
double yB = ((VPoint)P.elementAt(i+1)).ScaledY();

double G = h[i]*(yB+yA)/2 - h[i]*h[i]*h[i]*(mB+mA)/24;

In[i] = i > 0 ? In[i-1] + G : G;
}
}

// definite integral from the first
// data point through x
//
public double Integral(double x)
{
I = FindInterval(x);
return VPoint.ReverseY(ComputeIntegral(I, VPoint.ScaledX(x)));
}

private double ComputeIntegral(int index, double x)
{
if (index >= N)
return In[N-1];

double h = ((VPoint)P.elementAt(index+1)).ScaledX() - ((VPoint)P.elementAt(index)).ScaledX();
double A = (((VPoint)P.elementAt(index+1)).ScaledX() - x)/h;
double B = (x - ((VPoint)P.elementAt(index)).ScaledX())/h;
double A2 = A*A;
double B2 = B*B;
double mA = index == 0 ? 0 : M[index-1];
double mB = index > N-2 ? 0 : M[index];
double yA = ((VPoint)P.elementAt(index)).ScaledY();
double yB = ((VPoint)P.elementAt(index+1)).ScaledY();

double f =  (index > 0 ? In[index-1] : 0) +
(mB*B2*B2 - mA*A2*A2)*h*h*h/24 +
(yB*B2 - yA*A2)*h/2 -
(mB*B2 - mA*A2)*h*h*h/12 +
h*yA/2 - h*h*h*mA/24;

return f;
}

public void DrawIntegral(Graphics g, int xFrom, int xTo, int dx)
{
int x = xFrom;
Point p = new Point(x, (int)Math.round(Integral(x)));

while (x < xTo)
{
x += dx;
Point q = new Point(x, (int)Math.round(Integral(x)));

g.drawLine(p.x, p.y, q.x, q.y);

p = q;
}
}

// Checks whether a mouse press has occured
// close to one of the interpolation points.
// If not creates such a point and adds it
// to the array
//
public Point CloseEnough(int x, int y)
{
int i = FindInterval(x);
double Y = VPoint.ReverseY(Compute(i, VPoint.ScaledX(x)));

if (Math.abs(Y - y) < 5)
{
VPoint p = new VPoint(x, y);

P.insertElementAt(p, i+1);
N++;

return p;
}

return null;
}

public double SecondDerivatve(double x)
{
I = FindInterval(x);
return VPoint.ReverseY(ComputeSecondDerivative(I, VPoint.ScaledX(x)));
}

private double ComputeSecondDerivative(int index, double x)
{
double h = ((VPoint)P.elementAt(index+1)).ScaledX() - ((VPoint)P.elementAt(index)).ScaledX();
double A = (((VPoint)P.elementAt(index+1)).ScaledX() - x)/h;
double B = (x - ((VPoint)P.elementAt(index)).ScaledX())/h;

double mA = index == 0 ? 0 : M[index-1];
double mB = index > N-2 ? 0 : M[index];

return (A*mA + B*mB)*h;
}
}
```