DoubleFloats.jl -> Double128 | Quad64?

Specific Domains Numerics
41

yes

42

Ok thanks ! I’ll try that using some multiple dispatch on the type. Hopefuly it’ll reduce the runtie of bigfloats, i’ll tell you by how much as soon as possible.

43 
function xzyt(tmp::T, tmp2::T, x::T, y::T, z::T, t::T) where T<:BigFloat
  mul!(tmp, z, y)
  mul!(tmp2, tmp, t)
  add!(tmp, tmp2, x)
  return tmp
end

similarly for your other case with another multiply

44

Well @JeffreySarnoff thanks a lot for taking the time, but sadly i did not see any notable improvements from these. Maybe I misused them though.

What is your oppinon on exp/log for MultiFloats.jl ? I think this package is the way to go for me, regarding the speedup for 256ish bits of precision, although it seems like no updates have been made these past 5 months.

45

@lmrv I have no opinion about MultiFloats. exp and log have to be coded for the precision in use, unless variable precision algorithms are used (BigFloat, ArbNumerics). It is difficult to do correctly for new precisions.

Are you reusing the same temporary BigFloats? Doing that should speed things some.

46

Okay so I clearly wont be able to do it myself.

Yes I am :worried:

47

Arblib.jl supports exp and log along with arithmetic. And the default precision is 256 bits. Give it a try with your code, if you have not yet (ignore the warning).

1 Like
48

Arblib tests are a work in progress (see this issue)

49

I did not see that you needed expm1. You could try without it. Otherwise, I am out of directions for you.

50

You already helped a lot ! Thanks for all. I will try ArbLib as soon as i have some time, but for the moment it seems like MultiFloats is fastr than everything for these kind of precisions, even with the conversion to bigfloat for transcendental functions;

51

As of v0.2.1 Arblib already extends expm1 (it was under Arblib.expm1!, just not hooked to Base.expm1); Let me know when/if you need something more.

@JeffreySarnoff what was your rationale for defining eps as you did in ArbNumerics? In here: https://github.com/kalmarek/Arblib.jl/issues/101#issuecomment-738653420 I have three equally valid implementations (or at least that’s how I see it:)

52

@abulak because I defined both ulp and eps where
eps(x) = 2*ulp(x). Perhaps I should have swapped the names, perhaps not. There is a distinction in at least some of the literature.

note:

void arf_mag_set_ulp(mag_t res, const arf_t x, slong prec)
    Sets res to the magnitude of the unit in the last place
    (ulp) of x at precision prec.
``
53

for one thing: why 2*ulp?

eps(x::AbstractFloat)

Return the unit in last place (ulp) of x. This is the distance between consecutive representable floating point values at x. In most cases, if the distance on either side of x is different, then the larger of the two is taken

I thought Arblibs ulp takes care of it? Also you set eps(x::ArbFloat) to a small ball centered at eps(midpoint of x). Why this, and not a ball centered at 0 of radius eps(midpoint of x)?

54

As mentioned, swapping those names (eps, ulp) may be preferred.
I adopted an academic distinction rather than moving within the standard Julia terminology – and this is worth revisiting shortly.

you set eps(x::ArbFloat) to a small ball centered at eps(midpoint of x) . Why this, and not a ball centered at 0 of radius eps(midpoint of x) ?

using what I had called ulp(x) or, equivalently my eps(x)/2:

...
w = Mag()
# Sets w to the magnitude of the ulp of x at precision P
#  typeof(w) is Mag
 ccall(@libarb(arf_mag_set_ulp), Cvoid, (Ref{Mag}, Ref{ArbFloat}, Clong), w, x, P)
#   Converts w to an ArbFloat, `z`, so z has magnitude w
#   and, thus, z == the ulp of x 
z = ArbFloat{P}(w)
...

arf_mag_set_ulp is a direct call into the Arb C library. This is preferable to deriving the ulp more indirectly, or directly with more than one access to the library. Often, internals of the library compensate for corner cases and handle other subtle effects that are not easily fully available to the C programming interface. As a matter of course, I have learned to prefer this approach.

55

Ok, maybe now i understand your ulp/eps. To clarify the second question: there are (at least) three tentative definitions of ulp as an Arb (or ArbReal) i can come up with:

w = Mag()
res = Arb(prec=precision(x))
Arblib.set!(res, Arblib.set_ulp!(w, Arblib.midref(x), precision(x))
# res is a ball centered at ulp(x) with the appropriate radius

and

res = Arb(prec=precision(x))
Arblib.set_ulp!(Arblib.radref(res), Arblib.midref(x), precision(x))
# res is a ball centered at 0 with radius ulp(x)

and

res = Arb(prec=precision(x))
Arblib.set_ulp!(Arblib.radref(res), Arblib.midref(x), precision(x))
Arblib.set!(Arblib.midref(res), Arblib.radref(res))
Arblib.zero!(Arblib.radref(res))
# res is a ball centered at ulp(x) with radius 0

Was there any particular reason why did you go for the first one?