How To Implement Multi-Precision Functions?

General Usage
1

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?

2

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.

3

I would look at https://github.com/JuliaLang/julia/blob/master/base/mathconstants.jl

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: https://github.com/JuliaLang/julia/pull/41258)

4 Likes
4

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
5

how about


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

With one method for each type?

6

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?

7

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

8

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.

9

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