# How To Implement Multi-Precision Functions?

Say you’re writing a mathematical function for which you want to support multiple precisions – Float16, Float32, Float64, and BigFloat (up to 256 bits).

What is the best way to program the constants that are needed in the calculation?

Is below an appropriate strategy? It defines the constant at the highest precision, then let’s Julia reduce it as needed:

``````function multip(x)

println(eltype(x))

c = big"0.29200483948595689514283538207783029688471938696299037821234419190370162457800791"

T=eltype(x)

q = T(x*c)

return q

end
``````
``````julia> y=Float16(1.0)
Float16(1.0)

julia> multip(y)
Float16
Float16(0.292)

julia> y=Float32(1.0)
1.0f0

julia> multip(y)
Float32
0.29200485f0

julia> y=Float64(1.0)
1.0

julia> multip(y)
Float64
0.2920048394859569

julia> y=BigFloat(1.0)
1.0

julia> multip(y)
BigFloat
0.292004839485956895142835382077830296884719386962990378212344191903701624578007
``````

Or is there a better way?

I figure different precisions require different constants and using a BigFloat for all of them is going to be wasteful and super slow when dealing with anything other than BigFloat as the input. I’d define multiple specialized methods for each type and let dispatch handle it. That way, you’ll get fast code in either case.

Another option would be mimicking what Base is doing for Irrationals like Pi and e. You’ll have to check their implementation of e.g. multiplication.

I would look at julia/mathconstants.jl at master · JuliaLang/julia · GitHub

e.g.

``````Base.@irrational c 0.2920048394859569 big"0.29200483948595689514283538207783029688471938696299037821234419190370162457800791"
``````

More generally, it is better to replace the hard-coded `big"..."` with an expression that will compute the desired constant in the current `BigFloat` precision, whatever that is. This way you truly have an arbitrary precision constant.

Using `Base.@irrational` has the advantage of being much faster when working in double or single precision, because in that case `c` will use the precomputed `Float64` value directly without a costly runtime conversion from `BigFloat`. It also handles promotion for you automatically: `x*c` will automatically use `c` in the precision of `x`.

(For rational constants, you can just use exact rationals, e.g. `c = 2//3`. If `x` is a floating-point value, `c*x` will automatically evaluate `c` in the desired precision. This very recently got a performance boost so that the conversion of `2//3` to floating point can happen at compile time, but that optimization won’t be available until Julia 1.8: improve constant prop in gcd by simeonschaub · Pull Request #41258 · JuliaLang/julia · GitHub)

4 Likes

Well, that looked promising, but `Base.@irrational` defines global constants. It doesn’t work for local constants inside a function. It also doesn’t like to be fed `big"0.292"` constants.

``````function _ctest(x)
println(eltype(x))
Base.@irrational λ 0.29200483948595689514283538207783029688471938696299037821234419190370162457800791 c1
λ, const c1 = λ
y = λ * x;
println(eltype(y))

return y
end
``````

result:

``````ERROR: syntax: unsupported `const` declaration on local variable around irrationals.jl:187
``````

``````
inline function computemycconstant(Type)
#... do stuff
end
``````

With one method for each type?

Blockquote figure different precisions require different constants and using a BigFloat for all of them is going to be wasteful and super slow when dealing with anything other than BigFloat as the input. I’d define multiple specialized methods for each type and let dispatch handle it.

If I want to support Float16, Float32, Float64, BigFloat 128 bits, BigFloat 256 bits, that would mean copying the function 5 times, replacing the constants as appropriate. That would work, and would be fast, but I wonder if there’s a more elegant way to do it.

What I created a another function that just served up the constants?
`give_me_constants(Float64)`
and it returns a vector of the Float64 constants, and so forth.
Do you think that would be much slower than the “5 copies of the function” solution?

yes this is basically my solution. make it an inline function, and voila. I don’t think it’d be slower at all.

Yes, that’d be a possibility as well if the functions are otherwise identical. Using the type of what you need as a trait for the constants. If the number of constants per type is small, using a tuple will probably be faster & inline better.

You can just put it in a local `module` block to make its scope visible only to the function:

``````module _Foo
Base.@irrational c 0.2920048394859569 big"0.29200483948595689514283538207783029688471938696299037821234419190370162457800791"
function foo(x)
return x*c
end
end
using ._Foo: foo # import foo(x) to surrounding scope
``````

However, it’s not completely hidden, since it defines methods for the globally visible type `Irrational{:c}`. In order to hide it entirely (e.g. so it doesn’t conflict with `c` irrationals defined elsewhere, you’d want to define a local subtype of `AbstractIrrational`, or maybe generate a unique local symbol

``````using UUIDs
let c = Symbol(uuid4()) # unique local symbol
@eval module _Foo
Base.@irrational \$c 0.2920048394859569 big"0.29200483948595689514283538207783029688471938696299037821234419190370162457800791"
function foo(x)
return x*\$c
end
end
end
using ._Foo: foo
``````

This is a bit awkward, though — it would be nice to have a package that smooths out the process of defining a local `AbstractIrrational` constant.

1 Like