The answer is 42 but what is the question?

Started by Charles Pegge, July 11, 2007, 02:28:52 PM

Previous topic - Next topic

0 Members and 1 Guest are viewing this topic.

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



' 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