Storing high precision constants to use with fixed precision libraries (DoubleFloats.jl, Quadmath.jl, DecFP.jl, etc.)

General Usage
,
1

I am hoping to design algorithms that can be used efficiently with other fixed higher precision libraries such as DoubleFloats.jl, Quadmath.jl, DecFP.jl, MultiFloats.jl, etc… This usually involves polynomial approximations to a designated precision (128, 256). Is there a flexible way to store high precision coefficients to be used with all of these libraries?

For example, ideally something like this would work…

function test(x::T) where T
    # strings usually contain 30-50 digits
    P = (parse(T, "1.2"), parse(T, "1.3"), parse(T, "1.4"), parse(T, "1.5"))
    return evalpoly(x, P)
end

with benchmarks…

julia> @btime test(x) setup=(x=Float128(rand()))
  330.556 ns (0 allocations: 0 bytes)
1.36276565872626782310144123826112398e+00

julia> @btime test(x) setup=(x=Dec128(rand()))
  781.630 ns (0 allocations: 0 bytes)
3.119328959904789231759682793174765

# allocates for some reason
julia> @btime test(x) setup=(x=Double64(rand()))
  1.258 μs (20 allocations: 848 bytes)
2.2168204194286036

But we are paying the cost to parse the strings. A solution to avoid this is to just hard code the constants in a higher precision…

const P_f128 = (parse(Float128, "1.2"), parse(Float128, "1.3"), parse(Float128, "1.4"), parse(Float128, "1.5"))

function test2(x)
    return evalpoly(x, P_f128)
end

julia> @btime test2(x) setup=(x=Float128(rand()))
  71.869 ns (0 allocations: 0 bytes)
2.02825985806247608356097004488423765e+00

This approach works well and fast but a little less flexible than I was hoping for as it requires hard coding the constants for a specific package (and take it on as a dependency.) You could also cast any AbstractFloat to whatever internal precision you use but slightly less flexible.

Is there perhaps a good approach to be flexible enough to handle all the different higher precision packages in the current and future ecosystem but also avoid allocations and parsing of string constants during runtime?

2

If the numerator and denominator are small enough you could use rationals:

function test(x::T) where T
    P = T.((12//10, 13//10, 14//10, 15//10))
    return evalpoly(x, P)
end

The nice thing is that in many cases the compiler can do the conversion at compile time (as of JuliaLang/julia#41258).

1 Like
3

Thank you! There are times where we know the exact coefficients so this works well (as long as the Integer fits nicely with Int128). I like the idea using rationals, but sometimes terms look like this…


function test(x)
    return evalpoly(x, (1990919576929585163761247434580873723897677550573731281207275390625, -1094071724893410272201314846735800345180930752601602224440504638671875, 103359845691236916270005420886293061281368263471845448279855946394531250, -3993558675585831919973605001813276532326292885934464373884268722794718750, 83832389176247515253189196480455786065824463879054932596114704927571390625, -1101977288759423115916484463959767795703317649005495798822550381533967842875, 9865413224996821913093061266347613112205666995211786544733850490740903131000, -63509799534443193918054738326913581450607131662586232327271654588312820660616, 305080179205137176095342856930260393560550505786421683990077008259370104309410, -1122099959965657526344757894103766710826629020929459670948834078441128579791750, 3217185258543282976637373169047731813093578126603884759856059695249999306047500, -7276835259402589085458442196686990645669654847254510877738804991391058411412500, 13076651508858985557227319801010895818335803028865299751194038965212274849406250, -18719023753854332443081892087572754119280558462994056787647908433841920235468750, 21307327225910629020600045595173860611314592374884901403059587980149276271875000, -19156728115567236234371593788102013666866274144599851347429312225076676478125000, 13430107219321888006291381503763318985372222835071804983658005207838642173828125, -7186150593017474773099947110903785965879266917528290917207712073663895771484375, 2834031566467842832851805860262261718160136985649155528263587796394390332031250, -776308737033643856044315139790308583914612827783322139063211433790569824218750, 131897976398031751060811874321878385924211075364673046295410087596954345703125, -10468093364923154846096180501736379835254847251164527483762705364837646484375)) / 5456052820246562178600318839637653652351064473600000000
end

so we get the promotion to BigInt first. Of course, this is when our coefficients are known exactly but higher order terms usually require at least Int256 to store…

So I guess there are really two instances I’m looking at. 1. Where we know the higher order coefficients exactly as integers (usually fit within Int256) and 2. Where we just know the minimax floating point number in high precision (usually a Float128 number).

4

In this case you could write it as T(significand) * exp2(T(exponent)), using an Int128 significand, at the cost of an exp2 call.

3 Likes