Newton had a fun way to approximate general roots...

[Pete]

[Pete]

Cool, I think I figured it out. Just applied a couple of tweaks to the above, and...

Code: (Select All)

DIM AS DOUBLE a, n, root
DO
    a = 1
    INPUT "Number: "; n
    INPUT "Root:   "; root
    DO
        olda = a
        a = a - ((a) ^ root - n) / (root * a ^ (root - 1))
        PRINT a
    LOOP UNTIL olda = a
LOOP

Pete

[bplus]

LOOP UNTIL olda = a

Because it hard to get olda to be exactly = a

Loop Until abs(a-olda) < someTinyNumber ' Might work better

[Pete]

LOOP UNTIL olda = a

Because it hard to get olda to be exactly = a

Loop Until abs(a-olda) < someTinyNumber  ' Might work better

Good point. Many of these iteration functions eventually repeat. It might be an interesting challenge to find a number and root combination that does not. In any case, I would imagine when I convert this to string math, I can exit it by limiting the progressively longer output of digits. That would be preferable to trying to figure out what that "tiny number" would be with so many different decimal possibilities to consider. That's another algorithm to itself. The goal, of course, is to prevent an infinite loop situation. I could always use LOOP UNTIL STEVE = funny_looking, but that would cause a premature exit situation, in most cases.

Pete

[PhilOfPerth]

That all looks a bit hard for me...
I've been using "successive approximation", which I suspect is just another way of doing the same thing but is easier (for me) to understand:

DefDbl R, I, N
root = 1: IncRate = 1
Input "Number"; number
While root * root <> number
root = root + IncRate
If prevroot = root Then Exit While
If root * root > number Then
root = root - IncRate
IncRate = .9 * IncRate
End If
Print IncRate;
prevroot = root
Wend
Print "Root is"; root
Sleep

Of course, this only works within the limits of Double. The prevroot bit is to stop it if the same root is given twice.

[Jack]

Pete
for high precision sqr approximation using the Newton-Raphson method you can start with low precision and double the precision on each iteration
some pseudo code

Code: (Select All)

l=log(NUM_DIGITS*0.0625)*1.5 'appromimation to l=log(NUM_DIGITS/16)/log(2) ''+ a little extra
'get the first approximation using double arithmetic
'then
    limit&&=16
    for k=1 to l+1
        limit&&=2*limit&&
        tmp = sm_div(n, r, limit&&)
        r2 = sm_add(r, tmp, limit&&)
        r = sm_mul(r2, half, limit&&)
    next
    return r

you can apply the logic for n-root approximation
the first approximation will be more complicated if you want to go beyond the exponential range of double, you could try _Float for the first approximation but unless you use the CRT library it won't work because there's very limited function support for _Float

[Jack]

for completeness here's how you could use the CRT functions as first approximation

Code: (Select All)

Declare Library
    Function snprintf& (Dest As String, Byval l As Long, frmt As String, Byval x As _Float)
    Function powl## (ByVal x As _Float, Byval y As _Float)
End Declare

Dim As _Float x, y, z, p
Dim As String frmt, s
Dim As Long ok
frmt = "%Lf" + Chr$(0)
x = 3.1415926535897932384626433832795##
p = 0.33333333333333333333##
y = powl(x, p)
s = Space$(64)
frmt = "%.18Lg" + Chr$(0)
ok = snprintf(s, Len(s), frmt, y)
Print s
z = y * y * y
s = Space$(64)
ok = snprintf(s, Len(s), frmt, z)
Print s

[bplus]

LOOP UNTIL olda = a

Because it hard to get olda to be exactly = a

Loop Until abs(a-olda) < someTinyNumber  ' Might work better

Good point. Many of these iteration functions eventually repeat. It might be an interesting challenge to find a number and root combination that does not. In any case, I would imagine when I convert this to string math, I can exit it by limiting the progressively longer output of digits. That would be preferable to trying to figure out what that "tiny number" would be with so many different decimal possibilities to consider. That's another algorithm to itself. The goal, of course, is to prevent an infinite loop situation. I could always use LOOP UNTIL STEVE = funny_looking, but that would cause a premature exit situation, in most cases.

Pete

Yes, I was thinking ahead when you did start your String Math. Actually I was surprised that loops terminated on the few tests I made with your code using Double Type but the tests were neither extremely large nor small.

SomeTinyNumber was meant to be whatever precision you want so that the answer is "good enough" usu 10^-n n being number of decimals for your precision.

Jack is probably more experienced with this stuff good to have his input here!

[Jack]

Pete
here's my decfloat version for the n-th root, I hope that it will be helpful to you

Code: (Select All)

Sub fpnroot (result As decfloat, x As decfloat, p As Long, digits_in As Long)
    Dim As Long digits
    digits = digits_in
    If digits > NUM_DWORDS Then digits = NUM_DWORDS
    Dim As decfloat ry, tmp, tmp2
    Dim As Double t, t2
    Dim As Long i, ex, l, prec

    l = Log((NUM_DIGITS + 8) * 0.0625) * 1.5 'calculate the number of iterations needed
    ex = (x.exponent And &H7FFFFFFF) - BIAS - 1 'extract the exponent
    t = git(x.mantissa, 0) + git(x.mantissa, 1) / 100000000 'get 16 digits from x.mantissa
    'for example, if x = 3.1415926535897932384626433832795028842 then the above would give
    't = 31415926.53589793 because the mantissa doesn't have a decimal point, it's an integer
    'each element of the mantissa holds 8 digits
    'in this example ex = 0
    t = t / 10000000 'now t = 3.141592653589793
    t2 = Log(t) / p + Log(10) * ex / p 'log(t ^ (1/p)) = log(t)/p + Log(10) * ex / p
    'in this example since ex = 0, it becomes: log(t ^ (1/p)) = log(t)/p
    t2 = Exp(t2) 't2=t ^ (1/p)
    Call str2fp(ry, Str$(t2)) 'convert the double t2 to decfloat in ry
    prec = 3 '3 * 8 = 24 digits, prec here means number of 8 digit elements

    Call fpipow(tmp, ry, p - 1, prec) 'tmp = ry ^ (p-1)
    Call fpdiv(tmp, x, tmp, prec)     'tmp = x * tmp
    Call fpmul_si(tmp2, ry, p - 1, prec) 'tmp2 = ry * (p-1)
    Call fpadd(tmp2, tmp2, tmp, prec) 'tmp2 = tmp2 + tmp
    Call fpdiv_si(ry, tmp2, p, prec)  'ry = tmp2 / p
    For i = 1 To l
        prec = 2 * prec - 1
        Call fpipow(tmp, ry, p - 1, prec)    'tmp  = ry^(p-1)
        Call fpdiv(tmp, x, tmp, prec)        'tmp  = x/tmp
        Call fpmul_si(tmp2, ry, p - 1, prec) 'tmp2 = ry*(p-1)
        Call fpadd(tmp2, tmp2, tmp, prec)    'tmp2 = tmp2+tmp
        Call fpdiv_si(ry, tmp2, p, prec)     'ry   = tmp2/p
    Next
    result = ry
End Sub

[bplus]

@Pete,

Jack is doing Logs which I know you want to avoid, and your and Newton's scheme still needs to do roots (powers)! the very thing you need a function for so here is old idea I had some time ago.

It does the integer part of root by mult in loop (not very fancy for sure).
Then it converts the fraction to Binary 0's and 1's, it only uses Square Root but over and over again on 2.
I might have it explained better in comments

Code: (Select All)

_Title "Power Function 2 by bplus"
'QB64 X 64 version 1.2 20180228/86  from git b301f92

' started 2018-07-23  Naalaa has no power function (or operator), so I wrote a power function for it.
' ''Power pack the short version.txt
' ''written for Naalaa 6 by bplus posted 2018-07-23
' ''extracted from large test of fractions, random number functions, string... in Power test.txt
' 2018-07-24 Power Function 2 split is replaced with two much smaller and very handy string functions.

' OMG the crazy thing worked! It produced decent estimates of roots given the limitations of precision...
' Now I want to see how well it works with far greater precision available. So here we are, looking to see
' how this function compares to the regualar ^ operator in QB64.

'from Naalaa comments:
' The main purpose of this code: to demo a power function for real numbers,

' I had an idea for how real numbers to the power of real numbers might be done ie x ^ y = ?
' This means that not only can you take SQR of a number, you can get cube or cube root, quartic, 5th 6th... roots and any multiple

' It came from this simple idea
' 2 ^ 3.5 = 2 ^ 3 * 2 ^ .5 = 8 * Sqr(2)
' 3 ^ 3.5 = 3 ^ 3 * 3 ^ .5 = 27 * Sqr(3)

' so 2 ^ 3.25 = 2 ^ 3 * 2 ^ .25
' what is 2 ^ .25 ?  It is sqr(sqr(2)) !

' likewise 2 ^ 3.125 = 2 ^ 3 * 2 ^ 1/8
' what is 2 ^ 1/8 ? It is sqr(sqr(sqr(2))) !

' any decimal can be written as a sum of fraction powers of 2 ie 1/2^n, as any integer can be written in powers of 2.
' in binary expansions
' 1/2 = .1       or .5 base 10
' 1/4 = .01      or .25 base 10
' 1/8 = .001     or .125 base 10
' 1/16 = .0001   or .0625 base 10

' So with binary expansion of decimal, we can SQR and multiply our way to an estimate
' of any real number to the power of another real number using binary expansion of
' the decimal parts as long as we stay in Naalaa integer limits and are mindful of precision.

Const wW = 800
Const wH = 600
Screen _NewImage(wW, wH, 32)
_ScreenMove 360, 60
_Define A-Z As _FLOAT

Do
    Print "Testing the power(x, pow) function:"
    Input "(nothing) quits, Please enter a real number to raise to some power. x = "; x
    If x = 0 Then Exit Do
    Input "(nothing) quits, Please enter a real number for the power. pow = ", pw
    If pw = 0 Then Exit Do
    result = power(x, pw)
    Print result; " is what we estimate for"; x; " raised to power of"; pw
    Print x ^ pw; " is what the ^ operator gives us."
    Print
Loop
Print
Print "This is matching the ^ operator very well! This code is clear proof of concept!"
Print " OMG, it worked!!!"
Sleep


' A power function for real numbers (small ones but still!)
' x to the power of pow
Function power## (x As _Float, pow As _Float)
    'this sub needs: bExpand60$, leftOf$, rightOf$

    Dim build As _Float
    r$ = "0" + Str$(pow) 'in case pow starts with decimal
    integer$ = leftOf$(r$, ".")
    build = 1.0
    If integer$ <> "0" Then
        p = Val(integer$)
        For i = 1 To p
            build = build * x
        Next
    End If
    'that takes care of integer part

    n$ = rightOf$(r$, ".")
    If n$ = "" Then power = build: Exit Function

    'remove 0's to right of main digits
    ld = Len(n$)
    While Right$(n$, 1) = "0"
        n$ = Left$(n$, ld - 1)
        ld = Len(n$)
    Wend

    'note: we are pretending that the ^ operator is not available, so this is hand made integer power
    denom& = 10
    For i = 2 To ld
        denom& = denom& * 10
    Next

    'OK for bExpand60$ don't have to simplify fraction and that saves us having to extract n and d again from n/d
    bs$ = bExpand60$(Val(n$), denom&)

    'at moment we haven't taken any sqr of x
    runningXSQR = x

    'run through all the 0's and 1's in the bianry expansion, bs$, the fraction part of the power float
    For i = 1 To Len(bs$)
        'this is the matching sqr of the sqr of the sqr... of x
        runningXSQR = Sqr(runningXSQR)
        'for every 1 in the expansion, multiple our build with the running sqr of ... sqr of x
        If Mid$(bs$, i, 1) = "1" Then build = build * runningXSQR
    Next

    'our build should now be an estimate or x to power of pow
    power = build
End Function

'write a series of 1s and 0s that represent the decimal fraction n/d in binary 60 places long
Function bExpand60$ (nOver&, d&)
    Dim b As _Float, r As _Float
    ' b for base
    b = 0.5
    ' r for remainder
    r = nOver& / d&
    ' s for string$ 0's and 1's that we will build and return for function value
    s$ = ""
    ' f for flag to stop
    f% = 0
    ' c for count to track how far we are, don't want to go past 20
    c% = 0
    While f% = 0
        If r < b Then
            s$ = s$ + "0"
        Else
            s$ = s$ + "1"
            If r > b Then
                r = r - b
            Else
                f% = 1
            End If
        End If
        c% = c% + 1
        If c% >= 60 Then f% = 1
        b = b * 0.5
    Wend
    bExpand60$ = s$
End Function

Function leftOf$ (source$, of$)
    posOf = InStr(source$, of$)
    If posOf > 0 Then leftOf$ = Mid$(source$, 1, posOf - 1)
End Function

Function rightOf$ (source$, of$)
    posOf = InStr(source$, of$)
    If posOf > 0 Then rightOf$ = Mid$(source$, posOf + Len(of$))
End Function

Actually this idea might go better here:
https://qb64phoenix.com/forum/showthread.php?tid=881

The approach is similar until we convert the fraction part to Binary.