Solving difficult formulae efficiently Using Successive Approximation
v + v^2 + v^3 = 42
What is v?
This shows how to solve such headaches to 14 decimal places in 9 iterations, (7 if your estimate is very close).
This example shows how to do it in both BASIC and Inline Assembler. It uses scaled negative feedback to home into the answer.
PowerBasic
' Solving difficult formulae Efficiently Using Successive Approximation
' 1 using BASIC
' 2 using Assembler (saving about 40 clocks per loop)
' Charles E V Pegge
' 11 July 2007
' PowerBasic ver 8x
#COMPILE EXE
#DIM ALL
FUNCTION NastyFormula( BYVAL piv AS DOUBLE, BYVAL target AS DOUBLE ) AS DOUBLE
DIM pov AS DOUBLE, ov AS DOUBLE, ee AS DOUBLE, fbr AS DOUBLE
DIM i AS DWORD
' this is primed with 1 call only
ov=piv+piv*piv+piv*piv*piv-target
ee=1
DO
INCR i: IF (ABS(ee)<1e-14)OR(i>100) THEN FUNCTION=piv: EXIT LOOP
piv=piv+ee: pov=ov
ov=piv+piv*piv+piv*piv*piv-target
ee=ov*ee/(pov-ov)
LOOP
END FUNCTION
FUNCTION NastyFormulaASM( BYVAL piv AS DOUBLE, BYVAL target AS DOUBLE ) AS DOUBLE
#REGISTER NONE
DIM ee AS DOUBLE, pov AS DOUBLE, ov AS DOUBLE, precis AS DOUBLE
DIM nstate AS BYTE PTR
DIM nstates AS ASCIIZ*100
nstate=VARPTR(nstates)
ee=1
precis=1e-14
' uses all of the 8 fpu registers available.
' '
' initial '
'--------------------'
! mov ebx,nstate '
'! fsave [ebx] ' save fpu state
' '
! mov ecx,100 ' down counter to limit number of cycles
'
' initial stack level
! fld target ' 5
! fld precis ' 4
! fld piv ' 3
! fld ee ' 2
! fld ov ' 1
! fld pov ' 0
'
' make initial estimate to start off
'-------------------------
! call calcfun ' piv countains first guess, error will be
' returned in ov.
' ee contains an offset to the first guess. This
' provides a second estimate.
' both estimates together provide initial
' feedback scaling info.
' '
cycle: ' loop start
' '
'--------------------'
! dec ecx ' downcount
! jle xit ' exit when too many loops
'
' update piv '
'--------------------'
! fld st(2) ' load ee
' stack depth +1
! fadd st(0),st(4) ' add piv to ee
! fstp st(4) ' store result in piv
' stack depth +0
' update pov '
'! jmp xit
'--------------------'
! fld st(1) ' load ov
' stack depth +1
! fstp st(1) ' store in pov
' stack depth +0
! call calcfun ' do nul function: result in ov
' calc new ee '
'--------------------'
! fld st(1) ' ov
' stack depth +1
! fmul st(0),st(3) ' ee
! fld st(1) ' pov
' stack depth +2
! fsub st(0),st(3) ' ov
! fdivp st(1),st(0)' ee result in st(0)
' stack depth +1
' new ee is now in st(0)
' check limit of ee '
'--------------------'
! fst st(3) ' save to ee
! fabs ' make absolute
''! fcomip st(5) ' compare ee with precision threshhold
! db &hdf,&hf5 ' opcodes for fcomip st, st(5)
' stack depth +0
! jae cycle ' continue if above/eq to precision threshold
'
! jmp xit ' finish
'
'--------------------'
' calc function '
'--------------------'
calcfun: '
! fld st(3) ' piv
' stack level +1
! fmul st(0),st(4) ' square
! fld st(0) ' stack level +2
! fmul st(0),st(5) ' cube
! fadd st(0),st(5)
! faddp st(1),st(0)'
' stack level +1
! fsub st(0),st(6) ' target
! fstp st(2) ' save in ov
' stack level +0
! ret '
'--------------------'
xit: ' these come off the stack in reverse order
! fcomp st(0) ' pov
! fcomp st(0) ' ov
! fcomp st(0) ' ee
! fstp piv ' piv
! fcomp st(0) ' precis
! fcomp st(0) ' target
'--------------------'
'! frstor [ebx] ' restore state
'--------------------'
FUNCTION=piv
END FUNCTION
FUNCTION PBMAIN () AS LONG: DIM i AS LONG
#REGISTER NONE
DIM tick AS QUAD, tock AS QUAD
DIM c AS LONG,p AS LONG,q AS LONG
DIM v AS DOUBLE
p=VARPTR(tock):q=VARPTR(tick)
'---------------------------'
' CPU Clock count
'---------------------------'
' ' approx because it is not a serialised instruction
' ' it may execute before or after other instructions
' ' in the pipeline.
! mov ebx,p ' 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.
'---------------------------'
' Nasty Formula
' v+v^2+v^3=42
' what is v?
v=NastyFormulaASM(5,42) ' parameters: estimate then target
'---------------------------'
' CPU Clock count
'---------------------------'
' ' approx because it is not a serialised instruction
' ' it may execute before or after other instructions
' ' in the pipeline.
! mov ebx,q ' 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.
'---------------------------'
!
!
! mov c,ecx
MSGBOX "v ="+STR$(v)+" repeat loops ="+STR$(100-c)+" clocks ="+STR$(tick-tock)
END FUNCTION