The answer is 42 but what is the question?

[Charles Pegge]

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

Code Select Expand



' 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