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

[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