Cubic Spline ASM

[MikeTrader]

Charles,
wanted to submit this for your consideration:
http://www.jose.it-berater.org/smfforum/index.php?topic=1383.new#new
to see if it is worthy of ASM treatment - if it provides fun for you of course.

[Charles Pegge]

This is an interesting algorithm Mike so I'll give it a go.

I thought it might be useful to show the stages we need to go through to reach an assemblerised version. As Lewis Carroll describes his poem 'The Hunting of the Snark'  as an Agony in Eight Fits. So here is Fit the First:

Code Select Expand




' Cubic Spline Interpolation (c)1986 Mike Trader - Ported to PB 2006

' Refactoring to prepare for Assembler translation Charles E V Pegge 22 Oct 2007


#COMPILE EXE "SplineDemo.exe"
#DIM ALL
     

'¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤'
FUNCTION CubicSpline( x() AS SINGLE, y() AS SINGLE, c() AS SINGLE ) AS LONG ' Cubic Spline

  #REGISTER NONE

  LOCAL z() AS SINGLE

  LOCAL s AS SINGLE
  LOCAL i, k, n AS LONG
                     
    n = UBOUND( x() ) ' Number of nodes specified ( >2 )
    DIM z(5*n) ' Derivative
   
    '------------------------------'

    LOCAL pc AS SINGLE PTR: pc=VARPTR(c(0))
    LOCAL px AS SINGLE PTR: px=VARPTR(z(0))
    LOCAL py AS SINGLE PTR: py=VARPTR(z(0))
    LOCAL pz AS SINGLE PTR: pz=VARPTR(z(0))
   
    ' trial coding to see how it translates
    ''------------------------------'
    '! mov esi,pz                   '
    '! mov edx,n                    '
    '! xor ecx,ecx                  '
    ''------------------------------'
    'lk1:                           '
    '! inc ecx                      ' k
    '! cmp ecx,edx                  ' n
    '! jge xit1                     '
    ''------------------------------'
    '! mov edi,px                   '
    '! mov edi,px                   '
    '! fld  dword ptr [edi+ecx*4+4] '
    '! fsub dword ptr [edi+ecx*4]   '
    '! mov  ebx,n                   '
    '! shl  ebx,1                   '
    '! add  ebx,n                   '
    '! add  ebx,ecx                 '
    '! fst dword ptr [esi+ebx*4]    '
    ''------------------------------'
    '! fld st(0)                    ' dup
    '! push dword 6                 '
    '! fdiv dword ptr [esp]         '
    '! add esp,4                    '
    '! push ecx                     '
    '! add ecx,n                    '
    '! fstp dword ptr [esi+ecx*4]   '
    '! pop   ecx                    '
    ''------------------------------'
    '! mov  edi,py                  '
    '! fld  dword ptr [edi+ecx*4+4] '
    '! fsub dword ptr [edi+ecx*4]   '
    '! fdiv st(0),st(1)             '
    '! mov ebx,n                    '
    '! shl ebx,1                    '
    '! add ebx,ecx                  '
    '! fstp dword ptr [esi+ebx*4]   '
    '! fcomp             '          ' pop only
    '! jmp lk1                      '
    'xit1:                          '
   
    LOCAL ss AS SINGLE
    LOCAL ii AS LONG
   
    k=0
    DO
     INCR k
     IF k>=n THEN EXIT DO
      ss = x(k+1) - x(k)
      ii=2*n+k
      z(ii+n)=ss
      z(n+k) = ss / 6
      z(ii) = (y(k+1) - y(k)) / ss
     LOOP
   
    'FOR k = 1 TO n-1
    '  z(3*n+k) = x(k+1) - x(k)
    '  z(  n+k) = z(3*n+k) / 6
    '  z(2*n+k) = (y(k+1) - y(k)) / z(3*n+k)
    'NEXT k

    k=1
    DO
     INCR k
     IF k>=n THEN EXIT DO
     ii=2*n+k
      z(ii+n+n) = z(ii) - z(ii-1)
    LOOP

    'FOR k = 2 TO n-1
    '  z(4*n+k) = z(2*n+k) - z(2*n-1+k)
    'NEXT k

    DIM m AS LONG
    m=n+1
    ii=4*n+2
    c(2+m)  = -1 - (z(3*n+1) / z(3*n+2))
    c(3+m)  = 1 - c(2+m)
    c(2+m*2)  = 2*(z(m) + z(m+1)) - (z(m) * c(2, 1))
    c(3+m*2)  = (z(m+1) - (z(m)*c(3+m))) / c(2+m*2)
    z(ii) = z(ii) / c(2+m*2)

    'c(2, 1)  = -1 - (z(3*n+1) / z(3*n+2))
    'c(3, 1)  = 1 - c(2, 1)
    'c(2, 2)  = 2*(z(n+1) + z(n+2)) - (z(n+1) * c(2, 1))
    'c(3, 2)  = (z(n+2) - (z(n+1)*c(3, 1))) / c(2, 2)
    'z(4*n+2) = z(4*n+2) / c(2, 2)

    DIM jj AS LONG
   
    k=2
    DO
     INCR k
     IF k>=n THEN EXIT DO
      ii=n-1+k
      jj=4*n+k
      c(2+k*m)  = 2*(z(ii) + z(ii+1)) - (z(ii) * c(3+ (k-1)*m))
      c(3+k*m)  = z(n+k) / c(2+k*m)
      z(jj) = (z(jj) - (z(ii)*z(jj-1))) / c(2+k*m)
    LOOP

    'FOR k = 3 TO n-1
    '  c(2, k)  = 2*(z(n-1+k) + z(n+k)) - (z(n-1+k) * c(3, k-1))
    '  c(3, k)  = z(n+k) / c(2, k)
    '  z(4*n+k) = (z(4*n+k) - (z(n-1+k)*z(4*n-1+k))) / c(2, k)
    'NEXT k
   
    DIM kk AS LONG

    ii=4*n
    jj=ii+n
    kk=n*m
    s         = z(ii-2) / z(ii-1)
    c(1+kk)  = 1 + s + c(3, n-2)
    c(2+kk)  = -s - (c(1+kk) * c(3+kk-m))
    z(jj) = z(jj-2) - (c(1+kk) * z(jj-1))
    z(n)     = z(jj) / c(2, n)

    's        = z(3*n-2+n) / z(3*n-1+n)
    'c(1, n)  = 1 + s + c(3, n-2)
    'c(2, n)  = -s - (c(1, n) * c(3, n-1))
    'z(4*n+n) = z(4*n-2+n) - (c(1, n) * z(4*n-1+n))
    'z(n)     = z(4*n+n) / c(2, n)

    k=0
    jj=n-1
    DO
      INCR k
      IF k>=jj THEN EXIT DO
      i  = n - k
      ii = 4*n+i
      z(i) = z(ii) - (c(3+i*m) * z(i+1))
    LOOP

    'FOR k = 1 TO n-2
    '  i    = n - k
    '  z(i) = z(4*n+i) - (c(3, i) * z(i+1))
    'NEXT k

    z(1) = (-c(2, 1)*z(2)) - (c(3, 1)*z(3))
    ' z(1) = (-c(2, 1)*z(2)) - (c(3, 1)*z(3))
 
    k=0
    DO
      INCR k
      IF k>=n THEN EXIT DO
      s       = 6*z(3*n+k)
      ii=k*m
      jj=3*n+k
      kk=n+k
      c(1+ii) = z(k)/s
      c(2+ii) = z(k+1)/s
      c(3+ii) = y(k)/z(jj)   - z(k)*z(kk)
      c(4+ii) = y(k+1)/z(jj) - z(k+1)*z(kk)
    LOOP

    'FOR k = 1 TO n-1
    '  s       = 6*z(3*n+k)
    '  c(1, k) = z(k)/s
    '  c(2, k) = z(k+1)/s
    '  c(3, k) = y(k)/z(3*n+k)   - z(k)*z(n+k)
    '  c(4, k) = y(k+1)/z(3*n+k) - z(k+1)*z(n+k)
    'NEXT k

  FUNCTION = 1

END FUNCTION
                 

'¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤'
FUNCTION PBMAIN() AS LONG ' Spline
   
  #REGISTER NONE
                     
  LOCAL x(), y(), c() AS SINGLE

  LOCAL aa, bb, xval, yval, Lastx, Lasty AS SINGLE
  LOCAL k, nNodes, nPts AS LONG
  LOCAL hWin AS DWORD


                     
    GRAPHIC WINDOW "Spline Demo", 30, 30, 999, 600 TO hWin
    GRAPHIC ATTACH hWin, 0
    GRAPHIC COLOR -1, -1
    GRAPHIC SCALE (0, 600) - (999, 0)
               


    GRAPHIC COLOR %GRAY, -1
    GRAPHIC BOX (0, 0) - (999, 600), 0, -1, -1, 0
    GRAPHIC COLOR %RED, %BLUE


    nNodes = 4
    REDIM x(nNodes)
    REDIM y(nNodes)
    REDIM c(4, nNodes)
         
    'Define Nodes
    x(1) = 100 : y(1) = 580
    x(2) = 180 : y(2) = 200
    x(3) = 700 : y(3) = 500
    x(4) = 900 : y(4) = 200
     
    ' Draw Nodes
    FOR k = 1 TO nNodes
      xval = x(k)
      yval = y(k)
      GRAPHIC BOX (xval-4, yval-4) - (xval+4, yval+4), 20, %RED, RGB(191,191,191), 5
    NEXT k
             
    ' Calculate the Cubic Spline
    CALL CubicSpline( x(), y(), c() ) ' Returns 3rd order Polynominal for each Node in c()
                             

    ' Read off values along the spline
    nPts = 900
    Lastx = x(1)
    Lasty = y(1)
   
    xval=100

    LOCAL m  AS LONG: m=5 ' c stride of second dimension
    LOCAL ii AS LONG:     ' temp var

    DO
      IF xval>nPts THEN EXIT DO ' step 20
      k=0
      DO
        INCR k
        IF k>nNodes THEN EXIT DO
        IF xval < x(k) OR xval > x(k+1) THEN ITERATE ' Must be between nodes
        aa = x(k+1) - xval
        bb = xval - x(k)
        ii=k*m
        yval = c(1+ii)*aa*aa*aa + c(2+ii)*bb*bb*bb + c(3+ii)*aa + c(4+ii)*bb
        EXIT DO
      LOOP
      IF k = nNodes+1 THEN ITERATE
      GRAPHIC LINE (Lastx, Lasty) - (xval, yval)
      Lastx = xval
      Lasty = yval
      xval=xval+20
    LOOP
   
    'FOR xval = 100 TO nPts STEP 20
    '  FOR k = 1 TO nNodes-1
    '    IF xval < x(k) OR xval > x(k+1) THEN ITERATE ' Must be between nodes
    '    aa = x(k+1) - xval
    '    bb = xval - x(k)
    '    yval = c(1, k)*aa^3 + c(2, k)*bb^3 + c(3, k)*aa + c(4, k)*bb
    '    EXIT FOR
    '  NEXT k
    '  IF k = nNodes+1 THEN ITERATE
    '  GRAPHIC LINE (Lastx, Lasty) - (xval, yval)
    '  Lastx = xval
    '  Lasty = yval
    'NEXT

SLEEP 2000 ' quick test for use during tranlation
'MSGBOX "All Done?",64,"Spline Demo"

GRAPHIC WINDOW END
EXIT FUNCTION


END FUNCTION

'¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤'


[Charles Pegge]

Fit the Second

In which the refactored basic is turned line by line into Assembler:
each block is tested as we go.

Sample Block

Code Select Expand


    ! mov ecx,0                     ' k=0
    ! mov edx,n                     '
    ! mov esi,px                    '
    ! mov edi,py                    '
    ! mov ebx,pz                    '
                                    '
    do01:                           ' do
     ! inc ecx                      '  incr k
     ! cmp ecx,edx                  '  if k>=n then exit do
     ! jge loop01                   '
     ! fld  dword ptr [esi+ecx*4+4] '  ss = x(k+1) - x(k)
     ! fsub dword ptr [esi+ecx*4  ] '
     ! mov eax,edx                  '  ii=2*n+k
     ! shl eax,1                    '
     ! add eax,ecx                  '
     ! push eax                     ' [hold ii]
     ! add eax,edx                  '
     ! fst dword ptr [ebx+eax*4]    '  z(ii+n)=ss
     ! fld st(0)                    '  [dup]
     ! push dword 6                 '  z(n+k) = ss / 6
     ! fidiv dword ptr [esp]        '
     ! add esp,4                    '
     ! mov eax,edx                  '
     ! add eax,ecx                  '
     ! fstp dword ptr [ebx+eax*4]   '
     ! fld dword ptr [edi+ecx*4+4]  '  z(ii) = (y(k+1) - y(k)) / ss
     ! fsub dword ptr [edi+ecx*4]   '
     ! fxch                         '
     ! fdivp st(1),st(0)            '
     ! pop eax                      ' [ii]
     ! fstp dword ptr [ebx+eax*4]   '
     ! jmp do01                     '
    loop01:                         'loop
   
    'FOR k = 1 TO n-1
    '  z(3*n+k) = x(k+1) - x(k)
    '  z(  n+k) = z(3*n+k) / 6
    '  z(2*n+k) = (y(k+1) - y(k)) / z(3*n+k)
    'NEXT k



[MikeTrader]

Makes sense. I wondered if you would do it this way or take an intermediate step of converting it to run with pointers.
I have not made any attempt to optimize it in that regard. perhaps I should?

[Charles Pegge]

The only pointers needed are the base addresses of each array.

Once these are loaded into the registers, everything else is done with indexes. Indexed addressing is strongly supported in the x86, especially with SIBs (Scaleable Index Base adressing), for example:

fst dword ptr [ebx+ecx*4+8]

where ebx contains the base address of the array
and ecx contains the index, and the 8 is used as an additional offset of 2 singles

The BASIC equivalent of this might be

x(i+2)=...
where x is an array of  singles

[Charles Pegge]

Fit the Third

These things always take longer than expected, but here it is. If I had a better understanding of the algorithm, I could probably make it much shorter, but this follows the original basic fairly closely, making economys on the way.

I have chopped up the refactored basic expressions and put them in the comment column. This is essential  to do for debugging, but also enables you to see how the translation was acheived.

Note that unlike the Bezier curve the control points of the cubic spline lie on the path of the curve, which makes it more useful in many situations

Thank you Mike!

Code Select Expand


' Cubic Spline Interpolation (c)1986 Mike Trader - Ported to PB 2006

' PB Inline Assembler translation Charles E V Pegge 23 Oct 2007


#COMPILE EXE "SplineDemo.exe"
#DIM ALL


'¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤'
FUNCTION CubicSpline( x() AS SINGLE, y() AS SINGLE, c() AS SINGLE ) AS LONG ' Cubic Spline

  #REGISTER NONE

  LOCAL z() AS SINGLE

  LOCAL s AS SINGLE
  LOCAL i, k, n AS LONG

    n = UBOUND( x() ) ' Number of nodes specified ( >2 )
    DIM z(5*n) ' Derivative

    '------------------------------'

    LOCAL pc AS SINGLE PTR: pc=VARPTR(c(0))
    LOCAL px AS SINGLE PTR: px=VARPTR(x(0))
    LOCAL py AS SINGLE PTR: py=VARPTR(y(0))
    LOCAL pz AS SINGLE PTR: pz=VARPTR(z(0))


    LOCAL ss AS SINGLE
    LOCAL ii AS LONG
    LOCAL jj AS LONG
    LOCAL kk AS LONG
    LOCAL m  AS LONG


    ! mov edx,n                     '
    ! mov esi,px                    '
    ! mov edi,py                    '
    ! mov ebx,pz                    '
    ! mov ecx,0                     ' k=0
'-----------------------------------'                                    '
    do01:                           ' do
     ! inc ecx                      '  incr k
     ! cmp ecx,edx                  '  if k>=n then exit do
     ! jge loop01                   '
     '------------------------------'
     ! fld  dword ptr [esi+ecx*4+4] '  x(k+1)  ' ss = x(k+1) - x(k)
     ! fsub dword ptr [esi+ecx*4  ] '  x(k)
     '------------------------------'
     ! mov eax,edx                  '  n ii=2*n+k
     ! shl eax,1                    '  2*n
     ! add eax,ecx                  '  +k
     ! push eax                     '  [hold ii]
     ! add eax,edx                  '  [ii+n]
     '------------------------------'
     ! fst dword ptr [ebx+eax*4]    '  z(ii+n)=ss
     ! fld st(0)                    '  [dup]
     '------------------------------'
     ! push dword 6                 '  z(n+k) = ss / 6
     ! fidiv dword ptr [esp]        '  ss/6
     ! add esp,4                    '
     ! mov eax,edx                  '  n
     ! add eax,ecx                  '  n+k
     ! fstp dword ptr [ebx+eax*4]   '  z(n+k)=
     '..............................'
     ! fld dword ptr [edi+ecx*4+4]  '  y(k+1) ' z(ii) = (y(k+1) - y(k)) / ss
     ! fsub dword ptr [edi+ecx*4]   '  -k(k)
     ! fxch                         '
     ! fdivp st(1),st(0)            '  / ss
     ! pop eax                      '  [ii]
     ! fstp dword ptr [ebx+eax*4]   '  z(ii)
     ! jmp do01                     '
    loop01:                         'loop
'-----------------------------------'

    'k=0
    'DO
    ' INCR k
    ' IF k>=n THEN EXIT DO
    ' ss = x(k+1) - x(k)
    ' ii=2*n+k
    ' z(ii+n)=ss
    ' z(n+k) = ss / 6
    ' z(ii) = (y(k+1) - y(k)) / ss
    'LOOP
   
'-----------------------------------'
    ! mov ecx,1                     'k=1
    do02:                           'do
     ! inc ecx                      ' incr k
     ! cmp ecx,edx                  ' if k>=n then exit do
     ! jge loop02                   '
     ! mov eax,n                    ' ii=2*n+k
     ! shl eax,1                    '
     ! add eax,ecx                  '
     ! fld dword ptr [ebx+eax*4]    ' z(ii+n+n) = z(ii) - z(ii-1)
     ! fsub dword ptr [ebx+eax*4-4] '
     ! add eax,edx                  '
     ! add eax,edx                    '
     ! fstp dword ptr [ebx+eax*4]   '
     ! jmp do02                     '
    loop02:                         'loop
'-----------------------------------'


    'k=1
    'DO
    ' INCR k
    ' IF k>=n THEN EXIT DO
    ' ii=2*n+k
    '  z(ii+n+n) = z(ii) - z(ii-1)
    'LOOP

'-----------------------------------'
    !                               ' m=n+1
    ! mov eax,edx                   ' n            ' ii=4*n+2
    ! shl eax,2                     ' n*4
    ! add eax,2                     ' n*4+2
    ! mov ii,eax                    ' ii=n*4+2
    ! mov esi,pc                    ' [c array]
    ''-------------------------------'
    ! sub eax,edx                   ' 3*n+2         c(n+3)  = -1 - (z(3*n+1) / z(3*n+2))
    ! fld dword ptr [ebx+eax*4-4]   ' z(3*n+1)
    ! fdiv dword ptr [ebx+eax*4]    ' /z(3*n+2)
    ! fld1                          ' 1
    ! faddp st(1),st(0)             ' +()
    ! fchs                          ' -()
    ! fstp dword ptr [esi+edx*4+12] ' c(n+3)=
    '-------------------------------'
    ! fld1                          ' 1         ' c(n+4)  = 1 - c(2+m)
    ! fsub dword ptr [esi+edx*4+12] ' - c(2+m)
    ! fstp dword ptr [esi+edx*4+16] ' c(n+4)=
    '-------------------------------'
    ! fld1                          ' c(2+m*2)  = 2*(z(m) + z(m+1)) - (z(m) * c(2+m))
    ! fadd  st(0),st(0)             ' 2
    ! fld   dword ptr [ebx+edx*4+4] ' z(m)
    ! fadd  dword ptr [ebx+edx*4+8] ' z(m+1)
    ! fmulp st(1),st(0)             ' 2*()
    '...............................'
    ! fld   dword ptr [ebx+edx*4+4] ' z(m)
    ! fmul  dword ptr [esi+edx*4+12]' c(2+m)
    '...............................'
    ! fsubp st(1),st(0)             ' )-(
    '...............................'
    ! mov eax,edx                   ' n
    ! inc eax                       ' m
    ! shl eax,1                     ' [m*2]
    ! fstp dword ptr [esi+eax*4+8]  ' c(2+m*2)  =
    '-------------------------------'
    ! fld dword ptr [ebx+edx*4+8]   ' z(m+1)     c(3+m*2)  = (z(m+1) - (z(m)*c(3+m))) / c(2+m*2)
    ! fld dword ptr [ebx+edx*4+4]   ' z(m)
    ! fmul dword ptr [esi+edx*4+16] ' *c(3+m)
    ! fsubp st(1),st(0)             ' )-(
    '...............................'
    ! fdiv dword ptr [esi+eax*4+8]  ' )/(
    ! fstp dword ptr [esi+eax*4+12] ' c(3+m*2)  =
    '-------------------------------'
    ! mov eax,ii                    ' [ii]           z(ii) = z(ii) / c(2+m*2)
    ! fld dword ptr [ebx+eax*4]     ' z(ii)
    ! fdiv dword ptr [esi+edx*8+16] ' / c(2+m*2)
    ! fstp dword ptr [ebx+eax*4]    ' z(ii) =
    '-------------------------------'
nx1:

    m=n+1
    ii=4*n+2

    'c(n+3)  = -1 - (z(3*n+1) / z(3*n+2))
    'c(n+4)  = 1 - c(2+m)
    'c(2+m*2)  = 2*(z(m) + z(m+1)) - (z(m) * c(2, 1))
    'c(3+m*2)  = (z(m+1) - (z(m)*c(3+m))) / c(2+m*2)
    'z(ii) = z(ii) / c(2+m*2)

    ! mov edx,n
    ! mov esi,pc
    ! mov ebx,pz

    '-------------------------------'
    '! jmp loop03
    ! mov ecx,2                     ' k=2
   do03:                            ' do
    ! inc ecx                       '  incr k
    ! cmp ecx,edx                   '  IF k>=n THEN EXIT DO
    ! jge loop03                    '
    '-------------------------------'
    ! mov eax,edx                   '  ii=n-1+k
    ! dec eax                       '
    ! add eax,ecx                   '
    ! mov ii,eax                    '
    '-------------------------------'
    ! mov eax,edx                   '  jj=4*n+k
    ! shl eax,2                     '
    ! add eax,ecx                   '
    ! mov jj,eax                    '
    '-------------------------------'
    ! mov eax,edx                   '
    ! inc eax                       '
    ! push edx                      '
    ! mul ecx                       '
    ! pop edx                       '
    ! mov kk,eax                    ' kk=k*m
    '-------------------------------'
    ! fld1                          '  c(2+k*m)  = 2*(z(ii) + z(ii+1)) - (z(ii) * c(3+ (k-1)*m))
    ! fadd st(0),st(0)              '  [2]
    ! mov eax,ii                    '
    ! fld dword ptr [ebx+eax*4]     '  z(ii)
    ! fadd dword ptr [ebx+eax*4+4]  '  z(ii+1)
    ! fmulp st(1),st(0)             '  2*
    '...............................'
    ! fld dword ptr [ebx+eax*4]     '  z(ii)
    ! mov eax,kk                    '
    ! sub eax,edx                    '  kk - n -1 is (k-1)*m
    ! dec eax                       '
    ! fmul dword ptr [esi+eax*4+12] '  * c(3+ (k-1)*m)
    '...............................'
    ! fsubp st(1),st(0)             '  -
    ! mov eax,kk                    ' k*m
    ! fstp dword ptr [esi+eax*4+8]  '  c(2+k*m)
    '-------------------------------'
    ! mov eax,ecx                   '  c(3+k*m)  = z(n+k) / c(2+k*m)
    ! add eax,edx                   '  n+k
    ! fld dword ptr [ebx+eax*4]     '  z(n+k)
    '-------------------------------'
    ! push edx                      '  n
    '-------------------------------'
    ! mov eax,kk                    '  k*m
    ! fdiv dword ptr [esi+eax*4+08] '  c(2+k*m)
    ! fstp dword ptr [esi+eax*4+12] '  c(3+k*m)
    '-------------------------------'
    !                               '  z(jj) = (z(jj) - (z(ii)*z(jj-1))) / c(2+k*m)
    ! mov edx,jj                    '  jj
    ! fld dword ptr [ebx+edx*4]     '  z(jj)
    ! mov edx,ii                    '  ii
    ! fld dword ptr [ebx+edx*4]     '  z(ii)
    ! mov edx,jj                    '  jj
    ! fmul dword ptr [ebx+edx*4-4]  '  *z(jj-1)
    ! fsubp st(1),st(0)             '  -
    ! fdiv dword ptr [esi+eax*4+8]  '  /c(k*m+2)  [k*m in eax]
    ! fstp dword ptr [ebx+edx*4]    '  z(jj)
    '-------------------------------'
    ! pop edx                       '
    '-------------------------------'
    ! jmp do03                      ' loop
   loop03:                          '
'-----------------------------------'
    'k=2                            '
    'DO
    ' INCR k
    ' IF k>=n THEN EXIT DO
    '  ii=n-1+k
    '  jj=4*n+k
    '  kk=m*k
    '  c(2+kk)  = 2*(z(ii) + z(ii+1)) - (z(ii) * c(3+ (k-1)*m))
    '  c(3+kk)  = z(n+k) / c(2+kk)
    '  z(jj) = (z(jj) - (z(ii)*z(jj-1))) / c(2+kk)
    'LOOP


'-----------------------------------'
    ! mov eax, edx                  '
    ! shl eax,2                     '
    ! mov ii,eax                    '
    ! add eax,edx                   '
    ! mov eax,jj                    '
    ! mov eax,ii                    '
    ! fld dword ptr [ebx+eax*4-8]   ' z(ii-2)
    ! fdiv dword ptr [ebx+eax*4-4]  ' /z(ii-1) as s
    '! fld st(0)                     ' dup s [problem]
    ! fst dword ptr s
    '-------------------------------'
    ! fld1                          ' c(1+kk)  = 1 + s + c(3+kk-m-m) ' c(3,n-2)
    ! faddp st(1),st(0)             ' 1+s
    ! mov eax,kk                    ' kk
    ! sub eax,edx                   ' -n
    ! sub eax,edx                   ' -n
    ! sub eax,2                     ' -2
    ! fadd dword ptr [esi+eax*4+12] ' c(kk-m-m+3)
    ! mov eax,kk                    ' kk as m*n
    ! fstp dword ptr [esi+eax*4+4]  ' c(kk+1) =
    '-------------------------------'
    ! fld dword ptr s               ' c(2+kk)  = -s - (c(1+kk) * c(3+kk-m))
    ! fchs                          ' -s
    ! mov eax,kk                    ' kk
    ! fld dword ptr [esi+eax*4+4]   ' c(1+kk)
    ! sub eax,edx                   ' kk-n
    ! dec eax                       ' kk-m
    ! fmul dword ptr [esi+eax*4+12] ' c(kk-m+3)
    ! fsubp st(1),st(0)             ' -s -
    ! mov eax,kk                    ' kk
    ! fstp dword ptr [esi+eax*4+8]  ' c(kk+2)
    '-------------------------------'
    ! jmp nx2
                                    ' 'z(jj) = z(jj-2) - (c(1+kk) * z(jj-1))  ' redundant?
    ! mov eax,jj                    ' jj
    ! fld dword ptr [ebx+eax*4-8]   ' z(jj-2)
    ! mov eax,kk                    ' kk
    ! fld dword ptr [esi+eax*4+4]   ' c(kk+1)
    ! mov eax,jj                    ' jj
    ! fmul dword ptr [ebx+eax*4-4]  ' z(jj-1)
    ! fsubp st(1),st(0)             ' )-(
    ! fstp dword ptr [ebx+eax*4]    ' z(jj)
    '-------------------------------'
                                    '  'z(n)     = z(jj) / c(kk+2) ' redundant ?
    ! mov eax,jj                    '  jj
    ! fld dword ptr [ebx+eax*4]     '  z(jj)
    ! mov eax,kk                    '  kk
    ! fdiv dword ptr [esi+eax*4+8]  '  / c(kk+2)
    ! fstp dword ptr [ebx+edx*4]    ' z(n)=
    '-------------------------------'
    nx2:
    'ii=4*n
    'jj=ii+n
    'kk=n*m
    's         = z(ii-2) / z(ii-1)
    'c(1+kk)  = 1 + s + c(3+kk-m-m) ' c(3,n-2)
    'c(2+kk)  = -s - (c(1+kk) * c(3+kk-m))
    'z(jj) = z(jj-2) - (c(kk+1) * z(jj-1))  ' redundant?
    'z(n)     = z(jj) / c(kk+2)             ' redundant?


    '-------------------------------'
    ! mov eax,edx                   ' n
    ! dec eax                       ' n-1
    ! mov jj,eax                    ' n-1 as jj
    ! xor ecx,ecx                   ' k=0
    '-------------------------------'
   do04:                            '
    ! inc ecx                       '
    ! cmp ecx,jj                    ' if k>=n-1
    ! jge loop04                    '
    '-------------------------------'
    ! mov eax,edx                   ' n
    ! sub eax,ecx                   ' n-k as i      '  i  = n - k
    '-------------------------------'
    ! push edx                      ' n
    '-------------------------------'
    ! shl edx,2                     ' 4*n
    ! add edx,eax                   ' 4*n+i as ii   '  ii = 4*n+i
    ! fld dword ptr [ebx+edx*4]     ' z(ii)         '  z(i) = z(ii) - (c(3+i*m) * z(i+1))
    '-------------------------------'
    ! pop edx                       ' n
    '-------------------------------'
    ! push eax                      ' i
    '-------------------------------'
    ! push edx                      ' n
    '-------------------------------'
    ! inc edx                       ' m
    ! mul edx                       ' i*m in eax
    '-------------------------------'
    ! pop edx                       ' n
    '-------------------------------'
    ! fld dword ptr [esi+eax*4+12]  ' c(i*m+3)
    '-------------------------------'
    ! pop eax                       ' i in eax
    '-------------------------------'
    ! fmul dword ptr [ebx+eax*4+4]  ' * z(i+1)
    ! fsubp st(1),st(0)             ' -
    ! fstp dword ptr [ebx+eax*4]    ' z(i)=
    '-------------------------------'
    ! jmp do04
   loop04:                          '
    '-------------------------------'

    'k=0
    'jj=n-1
    'DO
    '  INCR k
    '  IF k>=jj THEN EXIT DO
    '  i  = n - k
    '  ii = 4*n+i
    '  z(i) = z(ii) - (c(3+i*m) * z(i+1))
    'LOOP


    '-------------------------------'
    ! mov eax,edx                   ' n         z(1) = (-c(2, 1)*z(2)) - (c(3, 1)*z(3))
    ! inc eax                       ' m
    ! fld dword  ptr [esi+eax*4+8]  ' c(2, 1)
    ! fchs                          ' -c(2, 1)
    ! fmul dword ptr [ebx+8]        ' *z(2)
    ! fld dword ptr [esi+eax*4+12]  ' c(3, 1)
    ! fmul dword ptr [ebx+12]       ' *z(3)
    ! fsubp st(1),st(0)             ' )-(
    ! fstp dword ptr [ebx+4]        ' z(1) =    redundant?
                                    '
    '-------------------------------'

    'z(1) = (-c(2, 1)*z(2)) - (c(3, 1)*z(3)) ' redundant?


    '! push edi
    ! mov edi,py                    ' y()
    ! xor ecx,ecx                   ' k=0
    '-------------------------------'
   do05:                            '
    ! inc ecx                       '
    ! cmp ecx,edx                    ' if k>=n
    ! jge loop05                    '
    '-------------------------------'
    !                               ' s       = 6*z(3*n+k)
    ! push dword 6                  ' 6
    ! fild dword ptr [esp]          ' 6
    ! add esp,4                     ' 6
    ! mov eax,edx                   ' n
    ! shl eax,1                     ' n*2
    ! add eax,edx                   ' n*3
    ! add eax,ecx                   ' n*3+k
    ! fmul dword ptr [ebx+eax*4]    ' s
    ! fstp dword ptr s              ' save s
    ! mov jj,eax                    ' jj=n*3+k
    ! mov eax,edx                   ' n
    ! inc eax                       ' m
    ! push edx                      '
    ! mul ecx                       ' m*k
    ! pop edx                       '
    ! mov ii, eax                   ' ii=m*k
    ! mov eax,edx                   ' n
    ! add eax,ecx                   ' n+k
    ! mov kk,eax                    ' kk=n+k
    ! fld dword ptr [ebx+ecx*4]     ' z(k)    ' c(1+ii) = z(k)/s
    ! fdiv dword ptr s              ' /s
    ! mov eax,ii                    ' ii
    ! fstp dword ptr [esi+eax*4+4]  ' c(i+ii)=
    '-------------------------------'
    !                               ' c(2+ii) = z(k+1)/s
    ! fld dword ptr [ebx+ecx*4+4]   ' z(k+1)
    ! fdiv dword ptr s              ' /s
    ! mov eax,ii                    ' ii
    ! fstp dword ptr [esi+eax*4+8]  ' c(ii+2)
    '-------------------------------'
    !                               '       ' c(3+ii) = y(k)/z(jj)   - z(k)*z(kk)
    ! fld dword ptr [edi+ecx*4]     ' y(k)
    ! mov eax,jj                    ' jj
    ! fdiv dword ptr [ebx+eax*4]    ' /z(jj)
    ! fld dword ptr [ebx+ecx*4]     ' z(k)
    ! mov eax,kk                    ' kk
    ! fmul dword ptr [ebx+eax*4]    ' *z(kk)
    ! fsubp st(1),st(0)             ' )-(
    ! mov eax,ii                    ' ii
    ! fstp dword ptr [esi+eax*4+12] ' c(3+ii) =
    '-------------------------------'
                                    '       '  c(4+ii) = y(k+1)/z(jj) - z(k+1)*z(kk)
    ! fld dword ptr [edi+ecx*4+4]   ' y(k+1)
    ! mov eax,jj                    ' jj
    ! fdiv dword ptr [ebx+eax*4]    ' /z(jj)
    ! fld dword ptr [ebx+ecx*4+4]   ' z(k+1)
    ! mov eax,kk                    ' kk
    ! fmul dword ptr [ebx+eax*4]    ' *z(kk)
    ! fsubp st(1),st(0)             ' )-(
    ! mov eax,ii                    ' ii
    ! fstp dword ptr [esi+eax*4+16] ' c(4+ii) =
    '-------------------------------'
    ! jmp do05                      '
   loop05:                          '
    '-------------------------------'
    'k=0
    'DO
    '  INCR k
    '  IF k>=n THEN EXIT DO
    '  jj=3*n+k
    '  s       = 6*z(jj)
    '  ii=k*m
    '  kk=n+k
    '  c(1+ii) = z(k)/s
    '  c(2+ii) = z(k+1)/s
    '  c(3+ii) = y(k)/z(jj)   - z(k)*z(kk)
    '  c(4+ii) = y(k+1)/z(jj) - z(k+1)*z(kk)
    'LOOP


  FUNCTION = 1

END FUNCTION


'¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤'
FUNCTION PBMAIN() AS LONG ' Spline

  #REGISTER NONE

  LOCAL x(), y(), c() AS SINGLE

  LOCAL aa, bb, xval, yval, Lastx, Lasty AS SINGLE
  LOCAL k, nNodes, nPts AS LONG
  LOCAL hWin AS DWORD



    GRAPHIC WINDOW "Spline Demo", 30, 30, 999, 600 TO hWin
    GRAPHIC ATTACH hWin, 0
    GRAPHIC COLOR -1, -1
    GRAPHIC SCALE (0, 600) - (999, 0)



    GRAPHIC COLOR %WHITE, -1
    GRAPHIC BOX (0, 0) - (999, 600), 0, -1, -1, 0
    GRAPHIC COLOR %RED, %BLUE


    nNodes = 4
    REDIM x(nNodes)
    REDIM y(nNodes)
    REDIM c(4, nNodes)

    'Define Nodes
    x(1) = 100 : y(1) = 580
    x(2) = 180 : y(2) = 200
    x(3) = 700 : y(3) = 500
    x(4) = 900 : y(4) = 200

    ' Draw Nodes
    FOR k = 1 TO nNodes
      xval = x(k)
      yval = y(k)
      GRAPHIC BOX (xval-4, yval-4) - (xval+4, yval+4), 20, %RED, RGB(191,191,191), 5
    NEXT k

    ' Calculate the Cubic Spline
    CALL CubicSpline( x(), y(), c() ) ' Returns 3rd order Polynominal for each Node in c()


    ' Read off values along the spline
    nPts = 900
    Lastx = x(1)
    Lasty = y(1)

    xval=100

    LOCAL m  AS LONG: m=5 ' c stride of second dimension
    LOCAL ii AS LONG:     ' temp var

    DO
      IF xval>nPts THEN EXIT DO ' step 20
      k=0
      DO
        INCR k
        IF k>nNodes THEN EXIT DO
        IF xval < x(k) OR xval > x(k+1) THEN ITERATE ' Must be between nodes
        aa = x(k+1) - xval
        bb = xval - x(k)
        ii=k*m
        yval = c(1+ii)*aa*aa*aa + c(2+ii)*bb*bb*bb + c(3+ii)*aa + c(4+ii)*bb
        EXIT DO
      LOOP
      IF k = nNodes+1 THEN ITERATE
      GRAPHIC LINE (Lastx, Lasty) - (xval, yval)
      Lastx = xval
      Lasty = yval
      xval=xval+20
    LOOP


'SLEEP 2000 ' quick test for use during tranlation
MSGBOX "All Done?",64,"Spline Demo"

GRAPHIC WINDOW END
EXIT FUNCTION


END FUNCTION

'¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤'


[MikeTrader]

Wow! that is substantially faster.
I get approx avg of 1700 Clks for your ASM version, compared to 3800 Clks before.

I forgot the underlying math. It has been a few years....

Thx Charles.