### Author Topic: Cubic Spline ASM  (Read 7669 times)

0 Members and 1 Guest are viewing this topic.

• Guest
##### Cubic Spline ASM
« on: October 22, 2007, 06:51:10 AM »
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

• Global Moderator
• Hero Member
• Posts: 674
• User-Rate: +27/-1
##### Re: Cubic Spline ASM
« Reply #1 on: October 22, 2007, 02:24:09 PM »
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]
`' 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 = 1END 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    'NEXTSLEEP 2000 ' quick test for use during tranlation'MSGBOX "All Done?",64,"Spline Demo"GRAPHIC WINDOW ENDEXIT FUNCTIONEND FUNCTION'¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤'`

#### Charles Pegge

• Global Moderator
• Hero Member
• Posts: 674
• User-Rate: +27/-1
##### Re: Cubic Spline ASM: An Agony in n Fits
« Reply #2 on: October 22, 2007, 04:26:44 PM »
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]
`    ! 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`
« Last Edit: October 22, 2007, 04:28:53 PM by Charles Pegge »

• Guest
##### Re: Cubic Spline ASM
« Reply #3 on: October 22, 2007, 06:50:02 PM »
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

• Global Moderator
• Hero Member
• Posts: 674
• User-Rate: +27/-1
##### Re: Cubic Spline ASM
« Reply #4 on: October 22, 2007, 10:26:20 PM »
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

• Global Moderator
• Hero Member
• Posts: 674
• User-Rate: +27/-1
##### Re: Cubic Spline ASM: Fit the third (No Boojum encountered)
« Reply #5 on: October 23, 2007, 05:03:22 PM »

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]
`' 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 = 1END 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 tranlationMSGBOX "All Done?",64,"Spline Demo"GRAPHIC WINDOW ENDEXIT FUNCTIONEND FUNCTION'¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤'`