Const As Integer n1 = 32
Const As Integer n2 = (n1 * (n1 - 1) / 2)
Const As Double paso = 0.05
Const As Double INF  = 1e308
Const As Double NaN  = -(INF/INF)

Dim As Double xVal(n1), tSin(n1), tCos(n1), tTan(n1)
For i As Integer = 1 To n1
    xVal(i) = (i-1) * paso
    tSin(i) = Sin(xVal(i))
    tCos(i) = Cos(xVal(i))
    tTan(i) = Tan(xVal(i))
Next i

Enum
    rSin
    rCos
    rTan
    rTrig
End Enum

Dim Shared rhot(rTrig, n2) As Double
For i As Integer = 0 To rTrig
    For j As Integer = 0 To n2
        rhot(i, j) = NaN
    Next j
Next i

Function rho(x() As Double, y() As Double, Byval rdx As Integer, _
    Byval i As Integer, Byval n As Integer) As Double
    If n < 0 Then Return 0
    If n = 0 Then Return y(i+1)

    Dim As Integer idx = (n1 - 1 - n) * (n1 - n) / 2 + i + 1
    If rhot(rdx, idx) = NaN Then    'valor aún no calculado
        rhot(rdx, idx) = (x(i+1) - x(i+1 + n)) _
        / (rho(x(), y(), rdx, i, n-1) - rho(x(), y(), rdx, i+1, n-1)) _
        + rho(x(), y(), rdx, i+1, n-2)
    End If
    Return rhot(rdx, idx)
End Function

Function thieleInterpolator(x() As Double, y() As Double, _
    Byval rdx As Integer, Byval xin As Double, Byval n As Integer) As Double
    If n > n1-1 Then Return 1
    Return rho(x(), y(), rdx, 0, n) - rho(x(), y(), rdx, 0, n-2) _
    + (xin-x(n+1)) / thieleInterpolator(x(), y(), rdx, xin, n+1)
End Function

Print "                            PI : "; 4 * Atn(1)
Print "                 6*arcsin(0.5) : "; 6 * Asin(0.5)
Print "                 3*arccos(0.5) : "; 3 * Acos(0.5)
Print "                 4*arctan(1.0) : "; 4 * Atn(1.0)

Print "6*thiele(tSin,xVal,rSin,0.5,0) : "; 6 * thieleInterpolator(tSin(), xVal(), rSin, 0.5, 0)
Print "3*thiele(tCos,xVal,rCos,0.5,0) : "; 3 * thieleInterpolator(tCos(), xVal(), rCos, 0.5, 0)
Print "4*thiele(tTan,xVal,rTan,1.0,0) : "; 4 * thieleInterpolator(tTan(), xVal(), rTan, 1.0, 0)

Sleep