REM>Quad1
REM Alain Brobecker (aka baah/Arm's Tech)
REM 03-Oct-1998 first version
REM 21-Nov-1998 faster version, and inc1/inc2 mixed in a register
scrX%=640:scrY%=480     :REMscreen dimension
m%=5:n%=2^m%            :REMblocks' size, must be <=5
grdX%=scrX%/n%+2        :REMgrid size

DIMasm% 65536:FORpass%=0TO2STEP2:P%=asm%+n%*4*2:[opt pass%

.AGrid:     dcd Grid
.AGridMvts: dcd GridMvts
.AGridWave: dcd GridWave

REM---------------------------------
REMQuadratic interpolation in dim 1
REMr0=vram adress
.quadratic
  stmfd     r13!,{r0-r12,r14}
REMAt first we compute the position of the points in the original grid
  ldr       r1,AGrid
  mov       r14,r1
  ldr       r2,AGridMvts
  ldr       r3,AGridWave
  mov       r12,#grdX%
.MkGridOne
  ldr       r4,[r2]
  add       r4,r4,r4,lsl #16
  str       r4,[r2],#4
  ldrB      r4,[r3,r4,lsr #22]
  strB      r4,[r14],#1
 subS r12,r12,#1:bNE MkGridOne
REMNow go on for the quadratic interpolation over those values
  mov       r14,#&ff
  mov       r2,#grdX%-2
.QuadOneBlock
  mov       r3,#n%
  ldrB      r8,[r1,#2]              ;Load the original points
  ldrB      r7,[r1,#1]
  ldrB      r6,[r1],#1
  FNquadcoefs(6,7,8)
.QuadOneValue
  mov       r8,r6,lsr #1+2*m%       ;r8=value in the range [0;255]
  add       r8,r8,r8,lsl #2         ;r8=y*5
  strB      r14,[r0,r8,lsl #7]      ;save at vram+x+y*5*128
  add       r6,r6,r7,asr #16        ;val+=inc1
  add       r7,r7,r7,lsl #17        ;inc1+=inc2
  add       r0,r0,#1
 subS r3,r3,#1:bNE QuadOneValue
 subS r2,r2,#1:bNE QuadOneBlock
  ldmfd     r13!,{r0-r12,pc}^

]:PRINT"Pass ";pass%/2;" : ";P%-asm%;" bytes"

REMAdresses for the original grid and movements
Grid=P%:P%+=(grdX%ANDNOT255)+256
GridMvts=P%:P%+=((grdX%ANDNOT255)+256)*4
GridWave=P%:P%+=1024
NEXTpass%

REMFill the points movements and curve
p1=PI/3
FORc%=0TO1023
 p=c%/1024*2*PI
 GridWave?c%=95+48*(COSp+COS(2*p+p1))
NEXTc%
a%=GridMvts
FORc%=0TOgrdX%-1
  !a%=RND AND&ffff00ff:a%+=4
NEXTc%

REMScreen and memory initialisations
MODE21+128:OFF:ORIGIN0,550
SYS"OS_ReadDynamicArea",2TOvram%

b%=1
REPEAT
 WAIT:SYS"OS_Byte",&71,b%:b%=b%EOR3:SYS"OS_Byte",&70,b%:CLS
 A%=vram%+(b%>>1)*scrX%*scrY%
 CALLquadratic
UNTILFALSE

REM======== COMPUTE COEFS FOR QUADRATIC INTERPOLATION ========
REMIN  - m0=h0
REM      m1=h1
REM      m2=h2
REMOUT - m0=val*2=a0*2<<2m%
REM      m1=inc1*2=(a1<<m%+a2)*2 | inc2=2*a2
REMWhere a0..a2 are the coefs of the curve a2*t^2+a1*t+a0
REMa0=(h1+h0)/2
REMa1=h1-h0
REMa2=(h2-2h1+h0)/2
DEFFNquadcoefs(m0,m1,m2)
[opt pass%
  sub       m2,m2,m1,lsl #1
  addS      m2,m2,m0                ;m2=h2-2h1+h0=a2*2
  addMI     m2,m2,#&10000
  add       m0,m0,m1                ;m0=h1+h0=a0*2
  rsb       m1,m0,m1,lsl #1         ;m1=h1-h0=a1
  mov       m0,m0,lsl #2*m%         ;m0=a0*2<<2m%
  add       m2,m2,m2,lsl #16        ;m1=2*a1 | 2*a2
  add       m1,m2,m1,lsl #17+m%     ;m1=2*(a1*2^m%+a2) | 2*a2
]:=""