What causes this loss of precision for increasing number of roots found?

Question

I have a while loop that calculates the roots of my equation. It works perfectly when I am working with low numbers, for example 100, the precision is a little more than 250 decimal places (I need 250 decimal places of precision).

The problem is when I increase the amount of my loop. The higher the amount, after 100 roots are found, the precision gets lower and lower. When I make it to a thousand roots, the precision of the last root is only 6 decimal places.

This is my code:

setprecision(BigFloat, 250; base=10)

Mp = BigFloat("250.0")
L = BigFloat("1.0")
rp = BigFloat("1.0")


exp = Mp - 10
pk = 50
h = 1/10
ts = -2*pk
te = -h
v0 = -pk
vmax = 0
u0 = 0
umax = pk
k1 = (vmax -v0)/h
k2 = (umax -u0)/h
 
max = 100.0
function arctanh(x)
    if abs2(x) < 1
        res = BigFloat(log( ( 1 + x ) / (1 - x)  ))
    else
        res = BigFloat(sign(x) * log1p( ( 2 * abs(x) ) / (1 - abs2(x)) ))
    end
    return 0.5*res
end

function tor_2(r)
    a = BigFloat( rp / r )
    b = BigFloat(arctanh(a))
    res = BigFloat(L^2 * (-b / rp))
    return res
end

function find_root(f, a, b, precision)
    oldm = big(a)
    n = big(a)
    p = big(b)
    m = (n + p) / big(2)
    epsilon = big(10)^(-precision)

    while abs(f(m)) > epsilon && oldm != m
        oldm = m
        val = f(m)
        if val > 0
            p = m
        elseif val < 0
            n = m
        else
            break
        end
        m = (n + p) / BigFloat(2.0)
    end
    return m
 
end
t = @elapsed begin
    tort = [find_root(r -> tor_2(r) - test, rp + 10.0^-exp, max, Mp) for test in ts:h:te]
end

where tort[1]= 1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000276779305347347506129736291395816937080609516467895441878785070622487206090198597561760179491324750990568982322199122021923608490714617064493293276254845275543133336, that is correct.

But, tort[1000]:10.033311132253989056880435010549127161920612410233012969744805223401820078850954716601912612057050747701277915837004398085186404284469580227204776685726634389654992756053090443362119682617587346589647165675974434032301224533470334376859854777986199219 is wrong.

The value was supposed to be: 10.03331113225398961014527049285149912673425410710831002117675758378789785273643535498548789545079535319530347609313690968614200882388184338974490564270302639426889010992841078245555898544027132450710912524700458841030045198039559947676313086937548943

Answer 1

I recommmend reading What Every Computer Scientist Should Know About Floating-Point Arithmetic, by David Goldberg and having a look at this answer.

I do not have time to go through the code in detail, but the following already gets you to a precision of 10.0333111322539896101452704928514991267342541071083100211767575837878978527364353549854878954507953531953034760931369096861420088238818433897449056427030263942688901099284107824555589854402713245071091252470045884103004519803955994767631308693754895186 which is correct up to the last two digits 10.0333111322539896101452704928514991267342541071083100211767575837878978527364353549854878954507953531953034760931369096861420088238818433897449056427030263942688901099284107824555589854402713245071091252470045884103004519803955994767631308693754895186 (it displays two more).

setprecision(BigFloat, 250; base=10)
    
    Mp = 250.0
    L = 1.0
    rp = 1
    
    
    exp = Mp - 10
    pk = 50
    h = 1//10
    ts = -2*pk
    te = -h
    v0 = -pk
    vmax = 0
    u0 = 0
    umax = pk
    k1 = (vmax -v0)//h
    k2 = (umax -u0)//h
     
    max = 100.0
    function arctanh(x)
        if abs2(x) < 1
            res = BigFloat(log( ( 1 + x ) / (1 - x)  ))
        else
            res = BigFloat(sign(x) * log1p( ( 2 * abs(x) ) / (1 - abs2(x)) ))
        end
        return 0.5*res
    end
    
    function tor_2(r)
        a = BigFloat( rp / r )
        b = BigFloat(arctanh(a))
        res = BigFloat(L^2 * (-b / rp))
        return res
    end
    
    function find_root(f, a, b, precision)
        oldm = big(a)
        n = big(a)
        p = big(b)
        m = (n + p) / big(2)
        epsilon = big(10)^(-precision)
    
        while abs(f(m)) > epsilon && oldm != m
            oldm = m
            val = f(m)
            if val > 0
                p = m
            elseif val < 0
                n = m
            else
                break
            end
            m = (n + p) / BigFloat(2.0)
        end
        return m
     
    end
    t = @elapsed begin
        tort = [find_root(r -> tor_2(r) - test, rp + 10.0^-exp, max, Mp) for test in ts:h:te]
    end

Quick explanation:

if you look at [test for test in ts:h:te] you see that you frequently divide by 5 or 10. However, these floating point numbers result in an inexact representation of rational numbers. In other words, 1/5 == 1//5 is false.