BASIC Cubic Spline Interpolation


This is a fourth order polynominal for creating a cubic spline that can be used to interpolate between nodes. The math is here.

'Basic can be handy for creating graphic windows without large amounts of code. This example will compile and run as is:


(the red line is interpolated between the four nodes - the missing section is an artifact of screen capture)

The example also contains code for timing the implementation. Below is an ASM translation that is about 2x faster.


' Cubic Spline Interpolation (c)1986  - Ported to PB 2006 - 3800 Clks


#COMPILE EXE "SplineDemo.exe" 
#DIM ALL
   
GLOBAL hDbg AS LONG
                      

'¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤'
SUB time_stamp_count(tick AS QUAD) ' approx CPU Clock count

  '---------------------------'
  '                           ' approx because it is not a serialised instruction
  '                           ' it may execute before or after other instructions
  '                           ' in the pipeline.
  ! mov ebx,tick              ' var address where count is to be stored.
  ! db  &h0f,&h31             ' RDTSC read time-stamp counter into edx:eax hi lo.
  ! mov [ebx],eax             ' save low order 4 bytes.
  ! mov [ebx+4],edx           ' save high order 4 bytes.
  '---------------------------'

END SUB
        

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

  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) ' Derrivative

    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
 
    FOR k = 2 TO n-1
      z(4*n+k) = z(2*n+k) - z(2*n-1+k)
    NEXT k 
   
    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)
 
    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
  
    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)
 
    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))
 
    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 i, k, nNodes, nLoops AS LONG
  LOCAL hWin AS DWORD
  LOCAL cBeg, cEnd AS QUAD


    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
              

    nLoops = 1000 ' even out timing calc
 
    ' Calculate the Cubic Spline   
    time_stamp_count(cBeg) ' measure cpu clock cycles.
      FOR i = 1 TO nLoops
        CALL CubicSpline( x(), y(), c() ) ' Returns 3rd order Polynominal for each Node in c() 
      NEXT
    time_stamp_count(cEnd) ' measure cpu clock cycles.


    ' Read off values along the spline
    Lastx = x(1)
    Lasty = y(1)
    k = 1
    FOR xval = x(1) TO x(4) STEP 20                          
        IF xval > x(k+1) THEN INCR k ' Must be between nodes
        IF xval < x(k) OR k > nNodes THEN EXIT FOR ' Must be between nodes
                         
        aa = x(k+1) - xval
        bb = xval   - x(k) ' PRINT #hDbg, "k="+STR$(k) + ", xval="+STR$(xval)  
        yval = c(1, k)*aa*aa*aa + c(2, k)*bb*bb*bb + c(3, k)*aa + c(4, k)*bb 

        GRAPHIC LINE (Lastx, Lasty) - (xval, yval) ' Join points  
        Lastx = xval
        Lasty = yval
    NEXT

MSGBOX "Clks="+STR$((cEnd-cBeg)\nLoops),64,"Spline Demo"

GRAPHIC WINDOW END

END FUNCTION












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

'  Inline Assembler translation - 1700 Clks

#COMPILE EXE "SplineASMDemo.exe"
#DIM ALL
             

'¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤'
SUB time_stamp_count(tick AS QUAD) ' CPU Clock count    Charles E V Pegge

  '---------------------------'
  '                           ' approx because it is not a serialised instruction
  '                           ' it may execute before or after other instructions
  '                           ' in the pipeline.
  ! mov ebx,tick              ' var address where count is to be stored.
  ! db  &h0f,&h31             ' RDTSC read time-stamp counter into edx:eax hi lo.
  ! mov [ebx],eax             ' save low order 4 bytes.
  ! mov [ebx+4],edx           ' save high order 4 bytes.
  '---------------------------'

END SUB
        

'¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤'
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 i, k, nNodes, nPts, nLoops AS LONG
  LOCAL hWin AS DWORD
  LOCAL cBeg, cEnd AS QUAD


    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
            
 
    nLoops = 1000

    ' Calculate the Cubic Spline   
    time_stamp_count(cBeg) ' measurecpu clock cycles.
      FOR i = 1 TO nLoops
        CALL CubicSpline( x(), y(), c() ) ' Returns 3rd order Polynominal for each Node in c() 
      NEXT
    time_stamp_count(cEnd) ' measure cpu clock cycles.

    ' Read off values along the spline
    Lastx = x(1)
    Lasty = y(1)
    k = 1
    FOR xval = x(1) TO x(4) STEP 20                          
        IF xval > x(k+1) THEN INCR k ' Must be between nodes
        IF xval < x(k) OR k > nNodes THEN EXIT FOR ' Must be between nodes
                         
        aa = x(k+1) - xval
        bb = xval   - x(k) ' PRINT #hDbg, "k="+STR$(k) + ", xval="+STR$(xval)  
        yval = c(1, k)*aa*aa*aa + c(2, k)*bb*bb*bb + c(3, k)*aa + c(4, k)*bb 

        GRAPHIC LINE (Lastx, Lasty) - (xval, yval) ' Join points  
        Lastx = xval
        Lasty = yval
    NEXT


MSGBOX "Clks="+STR$((cEnd-cBeg)\nLoops),64,"Spline Demo"

GRAPHIC WINDOW END

END FUNCTION

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


powerbasic

Comments