Interactive PowerBasic Forum

IT-Consultant: Charles Pegge => OxygenBasic Examples => Topic started by: Charles Pegge on March 11, 2024, 05:40:28 PM

Title: Successive Approximation
Post by: Charles Pegge on March 11, 2024, 05:40:28 PM
Using the Newton Raphson method.

Make two estimates of the value, and apply as inputs to a function:

Example, calculating phi:

1/phi=phi-1  What is the value of phi?

First turn this into a null equation

(1/phi) - phi + 1 = 0

then make two different estimates, say 0.4 and 0.5

and feed the function and the estimates into the macro:

uses console

  -----------------------------
  macro approx(fun, est1, est2)
  =============================
  '
  double ans1,ans2,scale
  'est1               'first estimate of answer
  'est2               'second estimate of answer
  ans1=fun(est1) 'error from first estimate
  ans2=fun(est2) 'error from second estimate
  if ans1<>ans2
    scale=(est2-est1)/(ans2-ans1)
    est1=est2
    est2-=ans2*scale
  endif
  end macro

  double e1,e2
  '
  e1=0.4
  e2=0.5
  function fn(double a) as double
    'null expressions:
    '
    'return a*a-2       'square root of 2
    'return (1/(a*a))-2 'inverse square root of 2
    'return (1/a)-a-1   'little phi
     return (1/a)-a+1   'big phi
  end function
  '
  int i
  for i=1 to 20
    double e3=e2
    approx(fn,e1,e2)
    if e3=e2
      exit for
    endif
    print e2 cr
  next
  wait
Title: Re: Successive Approximation
Post by: Theo Gottwald on March 11, 2024, 10:21:48 PM
The given code is a macro in O2 Basic that uses the approximation method to find the root of a function. The function used in this case is `fn(double a) as double` which returns the value of `(1/a)-a+1`, representing the "big phi" function. The macro `approx` takes in three parameters: `fun` (the function to approximate), `est1` and `est2` (the first and second estimates of the answer).

Here's an optimized version of the code in x86 ASM using the FPU (Floating Point Unit) to perform the calculations:

```assembly
section .data
    e1 dq 0.4 ; initial estimate 1
    e2 dq 0.5 ; initial estimate 2
    fmt db "%f", 10, 0 ; print format

section .bss
    ans1 resq 1
    ans2 resq 1
    scale resq 1

section .text
    global main
    extern printf

approx: ; macro approx(fun, est1, est2)
    ; calculate error from first estimate
    movsd xmm0, [rdi] ; load est1
    call fn
    movsd [ans1], xmm0 ; store ans1

    ; calculate error from second estimate
    movsd xmm0, [rsi] ; load est2
    call fn
    movsd [ans2], xmm0 ; store ans2

    ; check if ans1 != ans2
    movsd xmm0, [ans1]
    ucomisd xmm0, [ans2]
    jp .end_approx

    ; calculate scale
    movsd xmm0, [rsi] ; load est2
    subsd xmm0, [rdi] ; subtract est1
    movsd xmm1, [ans2]
    subsd xmm1, [ans1]
    divsd xmm0, xmm1
    movsd [scale], xmm0

    ; update estimates
    movsd xmm0, [ans2]
    mulsd xmm0, [scale]
    movsd xmm1, [rsi] ; load est2
    subsd xmm1, xmm0
    movsd [rsi], xmm1 ; store new est2
    movsd [rdi], xmm1 ; store new est1

.end_approx:
    ret

fn: ; function fn(double a) as double
    ; return (1/a)-a+1
    movsd xmm0, [rsp + 8] ; load a
    movsd xmm1, xmm0
    movsd xmm2, xmm0
    mulsd xmm1, xmm1 ; a^2
    cvtsd2ss xmm0, xmm0 ; convert double to single
    cvtsi2ss xmm2, dword 1 ; load constant 1
    divss xmm0, xmm1 ; 1/a
    subss xmm0, xmm2 ; (1/a) - 1
    cvtsd2ss xmm1, xmm1 ; convert double to single
    subss xmm0, xmm1 ; (1/a) - a - 1
    cvtsi2sd xmm2, dword 2 ; load constant 2
    addsd xmm0, xmm2 ; (1/a) - a + 1
    ret

main:
    ; main loop
    xor rdi, rdi
.loop:
    inc rdi
    cmp rdi, 20
    jge .end_loop

    ; call approx(fn, e1, e2)
    mov rax, e1
    mov rdx, e2
    call approx

    ; check if e3 == e2
    movsd xmm0, [e2]
    ucomisd xmm0, [e2 + 8]
    jp .end_loop

    ; print e2
    movsd xmm0, [e2]
    mov rdi, fmt
    call printf

    jmp .loop

.end_loop:
    ; wait and exit
    xor eax, eax
    ret
```

This x86 ASM code attempts to optimize the given O2 Basic code by using SIMD instructions and minimizing memory accesses. It utilizes the XMM registers for floating-point calculations and follows the System V AMD64 ABI for function calls.

MISTRAL-AI (One-Shot) (https://chat.mistral.ai/chat)
Title: Re: Successive Approximation
Post by: Charles Pegge on March 11, 2024, 11:17:30 PM
Looks interesting, Theo :).  I also have some assembler that uses the FPU, and thus calculates with extended precision. But the presentation here in Basic is for ease of understanding. I also want to extend it for solving multiple solutions/roots reliably. As you will see, the phi equation also produces negative little-phi when given negative estimates.
Title: Re: Successive Approximation
Post by: Theo Gottwald on March 11, 2024, 11:25:13 PM
It was meant as a motivation for you, to give MISTRAL AI a chance.
First ... its free. Second, it will not shorten the code like GPT-4.
Title: Re: Successive Approximation
Post by: Charles Pegge on March 11, 2024, 11:54:59 PM
I have 2 projects in mind:

The Marsenne Twister: Not a French gangster, but another high quality Randomizer.

Generating Voronoi diagrams: Computationally tricky to do efficiently.
Title: Re: Successive Approximation
Post by: Charles Pegge on March 13, 2024, 10:34:20 AM
This is an older example using the FPU to calculate the estimates:


'SUCCESSIVE APPROXIMATION USING NEGATIVE FEEDBACK
'================================================

  use console
  '
  macro ff(a)
  '==========
   '
   'NULL EXPRESSIONS
   '
   '0
   '1
   'a^5-2               'roots of 2
   '(a*a*a) + (a*a) -3 'compound roots
   'a*a*a-2            'cube root of 2
   abs (tan(a/4)-1)   'pi
   'abs((1/a)-a+1)     'phi
   'tan(rad(a))-2      'spherical angle of icosahedron facet edge
  end macro

  function approx(float e1,e2) as extended
  ========================================
  dim as double a,a1,a2,b1,b2,fbk
  dim as long i
  a1=e1      ' first estimate of answer
  a2=e2      ' second estimate of answer
  b1=ff(a1)  ' first result
  b2=ff(a2)  ' second result
  print i tab a2 cr
  do
    if i>=1000 then exit do 'restrict iterations
    if a1=a2 then exit do   'input precsion limit
    if b1=b2 then exit do   'feedback precision limit
    '
    'fbk=(a2-a1)/(b2-b1)    'feedback
    'a1=a2                  'new a1 input
    'a2=a2-b2*fbk           'new a2 input
    'b1=b2                  'new b1 output
      fld  qword a2
      fld  st0
      fld  st0
      fsub qword a1
      fld  qword b2
      fld st0
      fld st0
      fsub qword b1
      fxch st1
      fstp qword b1
      fdivp          'fbk
      fmulp          'qword b2
      fsubp          'qword a2
      fstp qword a2
      fstp qword a1
    b2=ff(a2)               'new b2 output
    i++
    '
    print i tab a2 cr
  end do
  return a2
  end function

'extended r=approx(-2,-1) ' lower est, upper est
extended r=approx(3.1414,3.1415) ' lower est, upper est
'print str(r)
waitkey