<< Statistics`DescriptiveStatistics`
<< LinearAlgebra`MatrixManipulation`


 (*    Hermite spline and complete spline interpolation of f (x)  =  
        x^(2  k)        *) 
 (*        on [0, 1] with 2  k - 
        4 interior knots                                    *) 
 (*        generated as order statistics of 2  k - 
        4 i . i . d .  uniform (0, 
            1)            *) 
 (*        random variables                                                      *) 
 (*    with  the knot configuration                                              *) 
 (*    in the form                                                               *) 
 (*    {hermite error,  
        complete error,  {knot_  1 ,  knot_  2 ,   ...  ,  
          knot_ (2  k - 
                4)}}   *) 
 (*        Restriction :   k  ≥  
          3                                          nn       *) 

Print["k :  (2k-1) degree spline interpolation"]
k = 7
Print["d = number of digits of accuracy"]
d = 
      100


B1[x_, j_, l_] := ((j!)/(j - l)!)*
        x^(j - l)
B2[x_, y_, k_, 
        l_] := ((2*k - 1)!/(2*k - 1 - l)!)*((x - y)*
              UnitStep[x - y])^(2*k - 1 - l)
hd[x_, k_, l_] := ((2*k)!/(2*k - l)!)*
        x^(2*k - l)

Basis[x_, vec_List, k_] := 
      Join[{1}, Table[B1[x, j, 0], {j, 1, 2*k - 1}], 
        Table[B2[x, vec[[j]], k, 0], {j, 1, 
            2*k - 4}]]

Vect[k_, l_, vec_List] :=  
      Join[Table[0, {m, 1, l}], Table[B1[1, i, l], {i, l, 2*k - 1}], 
        Table[B2[1, vec[[i]], k, l], {i, 1, 
            2*k - 4}]]

HermiteInterp[x_, coeff_List, vec_List, k_] := 
      coeff . Basis[x, vec, k]
CompInterp[x_, coeff_List, vec_List, k_] := 
      coeff . Basis[x, vec, k]

mxlistH[coeff_List, vec_List,  k_] := 
      Table[Abs[
          HermiteInterp[i/1000, coeff, vec, k] - hd[i/1000,  k, 0]], {i, 1, 
          1000}]
mxlistC[coeff_List, vec_List, k_] := 
      Table[Abs[CompInterp[i/1000, coeff, vec, k] - hd[i/1000,  k, 0]], {i, 
          1, 1000}]


Print["interior knots, tau_1, ..., tau_{2k-4}"]
VectauUnif = 
      Sort[Table[
          Random[Real, {0, 1}, d], {j, 1, 
            2*k - 4}]]
Print["all the knots"]
Vectau = Append[Prepend[VectauUnif, 0], 1]
MMHermite = 
        Join[
Table[
            Join[{1}, 
Table[
                B1[Vectau[[j]], i, 0], {i, 1, 2*k - 1}], 
              Table[B2[Vectau[[j]], Vectau[[i]], k, 0], {i, 2, 
                  2*k - 3}]], {j, 1, 2*k - 2}], 
Table[
            Join[{0, 1}, 
Table[
                B1[Vectau[[j]], i, 1], {i, 2, 2*k - 1}], 
              Table[
                B2[Vectau[[j]], Vectau[[i]], k, 1], {i, 2, 
                  2*k - 3}]], {j, 1, 2*k - 2}]];
VecHermite = 
        Join[Table[hd[Vectau[[j]], k, 0], {j, 1, 2*k - 2}], 
          Table[hd[Vectau[[j]], k, 1], {j, 1, 
              2*k - 2}]];
HermiteCoeff = 
        LinearSolve[MMHermite, VecHermite];
Hermitemxl = 
        mxlistH[HermiteCoeff, VectauUnif,  k];
Print["Maximum Hermite Error (on the grid j/1000, \
j=1,...,1000)"]
Hmx = Max[Hermitemxl]
MM1 =  
        AppendRows[DiagonalMatrix[Table[j!, {j, 0, (k - 1)}]], 
          Table[0, {i, 1, k}, {j, 1, 3*k - 4}]];
MM2 =  
        Table[Join[{1}, 
            Table[B1[VectauUnif[[j]], i, 0], {i, 1, 2*k - 1}], 
            Table[B2[VectauUnif[[j]], VectauUnif[[i]], k, 0], {i, 1, 
                2*k - 4}]], {j, 1, 2*k - 4}];
Vec3 =  
        Join[Table[1, {m, 1, 2*k}],  
          Table[B2[1, VectauUnif[[i]], k, 0], {i, 1, 
              2*k - 4}]];
MM3 =  
        Join[{Vec3}, 
          Table[Vect[k, l, VectauUnif], {l, 1, 
              k - 1}]];
MMComp =  Join[MM1, MM2, MM3];
VecC = 
        Join[Table[hd[0, k, l], {l, 0, k - 1}],  
          Table[hd[VectauUnif[[j]], k, 0], {j, 1, 2*k - 4}],  
          Table[hd[1, k, l], {l, 0, k - 1}]];
CompCoeff = LinearSolve[MMComp, VecC];
Compmxl = 
        mxlistC[CompCoeff, VectauUnif, k];
Print["Maximum Complete Interpolation error (on the grid j/1000, \
j=1,...,1000)"]
Compmx =  Max[Compmxl]
Print["Max Hermite Error, Max Complete Error, knots"]
mx = {Hmx, Compmx, Vectau}
Print["Monomial f(x) = x^(2k) being interpolated"]
P1 = 
      Plot[hd[x, k, 0], {x, 0, 1}, 
        PlotRange -> All]
Print["Hermite Interpolant"]
P2 = 
      Plot[HermiteInterp[x,  HermiteCoeff, VectauUnif, k], {x, 0, 1}, 
        PlotRange -> All]
Print["Complete Interpolant"]
P3 = 
      Plot[CompInterp[x, CompCoeff, VectauUnif, k], {x, 0, 1}, 
        PlotRange -> All]
Print["Monomial with Hermite and Complete Interpolants"]
Show[P1, P2, 
      P3]