juliamath / doubledouble.jl

Extended precision arithmetic for Julia (deprecated)

DoubleDouble.jl

Note: This package is no longer maintained. I suggest using DoubleFloats.jl instead.

DoubleDouble.jl is a Julia package for performing extended-precision arithmetic using pairs of floating-point numbers. This is commonly known as "double-double" arithmetic, as the most common format is a pair of C-doubles (Float64 in Julia), although DoubleDouble.jl will actually work for any floating-point type. Its aim is to provide accurate results without the overhead of BigFloat types.

The core routines are based on the ideas and algorithms of Dekker (1971).

Interface

The main type is Double, with two floating-point fields: hi, storing the leading bits, and lo storing the remainder. hi is stored to full precision and rounded to nearest; hence, for any Double x, we have abs(x.lo) <= 0.5 * eps(x.hi). Although these types can be created directly, the usual interface is the Double function:

julia> using DoubleDouble

julia> x = Double(pi)
Double{Float64}(3.141592653589793,1.2246467991473532e-16)

julia> eps(x.hi)
4.440892098500626e-16

The other type defined is Single, which is simply a wrapper for a floating-point type, but whose results will be promoted to Double.

Examples

Exact products and remainders

By exploiting this property, we can compute exact products of floating point numbers.

julia> u, v = 64 * rand(), 64 * rand()
(15.59263373822506,39.07676672446341)

julia> w = Single(u) * Single(v)
Double{Float64}(609.3097112086186, -5.3107663829696295e-14)

Note that the product of two Singles is a Double: the hi element of this double is equal to the usual rounded product, and the lo element contains the exact difference between the exact value and the rounded.

This can be used to get an accurate remainder

julia> r = rem(w, 1.0)
Double{Float64}(0.309711208618584, 1.6507898617445858e-17)

julia> Float64(r)
0.309711208618584

This is much more accurate than taking ordinary products, and gives the same answer as using BigFloats:

julia> rem(u*v, 1.0)
0.3097112086186371

julia> Float64(rem(big(u) * big(v), 1.0))
0.309711208618584

However, since the DoubleDouble version is carried out using ordinary floating-point operations, it is of the order of 1000x faster than the BigFloat version.

Correct rounding with non-exact floats

If a number cannot be exactly represented by a floating-point number, it may be rounded incorrectly when used later, e.g.

julia> pi * 0.1
0.3141592653589793

julia> Float64(big(pi) * 0.1)
0.31415926535897937

We can also do this computation using Doubles (note that the promotion rules mean that only one needs to be specified):

julia> Float64(Double(pi) * 0.1)
0.31415926535897937

julia> Float64(pi * Single(0.1))
0.31415926535897937

Emulated FMA

The fused multiply-add (FMA) operation is an intrinsic floating-point operation that allows the evaluation of a * b + c, with rounding occurring only at the last step. This operation is unavailable on 32-bit x86 architecture, and available only on the most recent x86_64 chips, but can be emulated via double-double arithmetic:

julia> 0.1 * 0.1 + 0.1
0.11000000000000001

julia> Float64(big(0.1) * 0.1 + 0.1)
0.11

julia> Base.fma(a::Float64,b::Float64,c::Float64) = Float64(Single(a) * Single(b) + Single(c))
fma (generic function with 1 method)

julia> fma(0.1, 0.1, 0.1)
0.11

doubledouble.jl's Issues

Tag initial version

@simonbyrne Could you please tag an initial version of DoubleDouble. I am using it in https://github.com/ivanslapnicar/Arrowhead.jl

DoubleDouble doesn't work with 0.4

Dear Simon:
Using Julia v0.4, I notice a few warning/error messages with your package; probably easy to fix, but I'm new to Julia...

Version 0.4.0-dev+5497 (2015-06-21 17:56 UTC)
julia> using DoubleDouble
Warning: New definition
call(Type{DoubleDouble.Double{#T<:Union{Float64, Float32}}}, Any, Any)
is ambiguous with:
call(Type{#T<:FloatingPoint}, Real, Base.Rounding.RoundingMode) at rounding.jl:68.
To fix, define
call(Type{DoubleDouble.Double{#T<:Union{Float64, Float32}}}, Real, Base.Rounding.RoundingMode)
before the new definition.

julia> x = double(pi)
Error showing value of type DoubleDouble.Double{Float64}:
ERROR: no promotion exists for DoubleDouble.Double{Float64} and Int64
in _show at grisu.jl:62
...

Basic functionality seems OK:
julia> x = double(pi);
julia> float64(x)
WARNING: float64(x) is deprecated, use Float64(x) instead.
...
julia> Float64(x)
3.141592653589793

Dekker multiplication isn't exact

I was surprised to discover this:

x = Float16(0.992)
y = Float16(6.0e-8)   # subnormal

julia> d = Single(x)*Single(y)
Double(6.0e-8, 0.0)
 - value: 5.9604644775390625e-08

julia> bits(widen(d.hi) + widen(d.lo))
"00110011100000000000000000000000"

julia> bits(widen(x)*widen(y))
"00110011011111100000000000000000"

It can be fixed this way:

julia> ys, ye = frexp(y)
(Float16(0.5), -23)

julia> d = Single(x)*Single(ys)
Double(0.496, 0.0)
 - value: 0.49609375

julia> bits(ldexp(widen(d.hi) + widen(d.lo), ye))
"00110011011111100000000000000000"

But given that many treatises say "it's exact unless the splitting overflows," this was a surprising discovery to me.

Move to JuliaMath org?

Maybe it would make sense to move this to the JuliaMath org.

DoubleDouble Cannot Multiply Ints

Example:

julia> Double(2)*2
ERROR: MethodError: Cannot `convert` an object of type Int64 to an object of type DoubleDouble.Double{Float64}
This may have arisen from a call to the constructor DoubleDouble.Double{Float64}(...),
since type constructors fall back to convert methods.
 in *(::DoubleDouble.Double{Float64}, ::Int64) at ./promotion.jl:191
 in eval_user_input(::Any, ::Base.REPL.REPLBackend) at ./REPL.jl:62
 in macro expansion at ./REPL.jl:92 [inlined]
 in (::Base.REPL.##3#4{Base.REPL.REPLBackend})() at ./event.jl:46

Correctly-rounded operations

It would be great to have correctly-rounded basic operations, i.e. be able to fix the rounding mode to RoundDown or RoundUp and do a + b for a::Double{T} and b::Double{T}.

Error tagging new release

The URL of this package does not match that stored in METADATA.jl.
cc: @dpsanders

Enable Travis

@simonbyrne I have made a runtests.jl file on the tests branch. Could you please enable Travis on this repo?

ERROR: LoadError: - not defined for DoubleDouble.Single{Float64}

I'm trying to use this library but I get the following:

julia> Single(1.0) - Single(0.0)
ERROR: - not defined for DoubleDouble.Single{Float64}
Stacktrace:
 [1] -(::DoubleDouble.Single{Float64}, ::DoubleDouble.Single{Float64}) at ./promotion.jl:337

Basic math functions

Are there any plans to make this a drop-in replacement from Float64 or BigFloat? At this stage this is impossible because most of the missing basic math functions like sin, exp, etc.

[PkgEval] Your package doesn't have a test/runtests.jl file

PackageEvaluator.jl is a script that runs nightly.
It attempts to load all Julia packages and run their tests (if available) on both the stable
version of Julia (0.3) and the nightly build of the unstable version (0.4).
The results of this script are used to generate a package listing
enhanced with testing results. This service also benefits package developers by notifying them if
their package breaks for some reason (caused by e.g. changes in Julia, changes in dependencies,
or broken binary dependencies.)

Currently PackageEvaluator attempts to find your test scripts using a heuristic, preferring the
standarized test/runtests.jl whenever present. Using test/runtests.jl allows people to test
your package using simply Pkg.test("DoubleDouble"), with any testing-only dependencies being
installed by looking at test/REQUIRE.

Your package doesn't appear to have a test/runtests.jl file. PackageEvaluator is going to move
away from auto-detecting tests and will instead only test packages with a test/runtests.jl
file. This change will take place in about a month.

You can:

• Add the file and tag a new version. You may in fact have already added this file but not
  tagged a new version. PackageEvaluator only tests your latest tagged verison, so you must tag
  for the file to be detected.
• Chose to do nothing. PackageEvaluator will stop attempting to test your package, and the testing
  status will be reported as "not possible".

If you'd like help or more information, please just reply to this issue.

Moving things forward

Following the lively discussion on discourse, what do you think is now the best way going forward?
I think we all have the same goal: Getting a solid and fast double-double and quad-double precision library.
But the question is: How do we want to achieve this?

Here are some more or less unstructured thoughts from my side:

• I don't care whether my package survives or not. But I need something usable.
• Should the development happen in DoubleDouble.jl or maybe in a new repo under the JuliaMath organization? I think the release of Julia 1.0 would be a good point to move to get rid of "legacy" packages, i.e. I would prefer a clean slate.
• I think only using FMA instructions instead of things like Veltkamp splitting is the way to go
• I am not convinced that a parameterization with an T<:Union{Float64, Float32, Float16} is worth the added complexity.
• Who is actually motivated enough / has the time to contribute in a meaningful way to such a project?

cc @simonbyrne @JeffreySarnoff @dpsanders

Set up Attobot

Could you please enable Attobot on this repo, @simonbyrne?

Info about upcoming removal of packages in the General registry

As described in https://discourse.julialang.org/t/ann-plans-for-removing-packages-that-do-not-yet-support-1-0-from-the-general-registry/ we are planning on removing packages that do not support 1.0 from the General registry. This package has been detected to not support 1.0 and is thus slated to be removed. The removal of packages from the registry will happen approximately a month after this issue is open.

To transition to the new Pkg system using Project.toml, see https://github.com/JuliaRegistries/Registrator.jl#transitioning-from-require-to-projecttoml.
To then tag a new version of the package, see https://github.com/JuliaRegistries/Registrator.jl#via-the-github-app.

If you believe this package has erroneously been detected as not supporting 1.0 or have any other questions, don't hesitate to discuss it here or in the thread linked at the top of this post.

Package compatibility caps

Ref: https://discourse.julialang.org/t/package-compatibility-caps/15301

`split` clashes with Base function

This is likely to cause problems since it clashes with the Base function of the same name and rather different functionality.

No sqrt(Double(0.0)) and sign(Double(0.0))

Simon,

presently sqrt(Double(0.0)) and sign(Double(0.0)) throw an error.
Possible remedy for sqrt is to change

# Dekker sqrt2
function sqrt{T}(x::Double{T})
    if x.hi <= 0
        throw(DomainError("sqrt will only return a complex result if called with a complex argument."))
    end
    c = sqrt(x.hi)
    u = Single(c)*Single(c)
    cc = (((x.hi - u.hi) - u.lo) + x.lo)*map(typeof(x.hi),0.5)/c
    double(c,cc)
end

into

# Dekker sqrt2
function sqrt{T}(x::Double{T})
    if x.hi = 0
         return Double(zero(T))
    end
    if x.hi < 0
        throw(DomainError("sqrt will only return a complex result if called with a complex argument."))
    end
    c = sqrt(x.hi)
    u = Single(c)*Single(c)
    cc = (((x.hi - u.hi) - u.lo) + x.lo)*map(typeof(x.hi),0.5)/c
    double(c,cc)
end

Cheers, Ivan

ERROR: <= not defined for DoubleDouble.Double{Float64}

I thought <= might have been inferred from < and == but apparently, that's not the case:

julia> a = Double(1.0);

julia> b = Double(2);

julia> a < b
true

julia> a == b
false

julia> a ≤ b
ERROR: <= not defined for DoubleDouble.Double{Float64}
Stacktrace:
 [1] <=(::DoubleDouble.Double{Float64}, ::DoubleDouble.Double{Float64}) at ./promotion.jl:351

It seems to be a matter of

julia> import Base.≤

julia> ≤{T}(x::Double{T}, y::Double{T}) = x < y || x == y
<= (generic function with 58 methods)

Tag new version

Is it possible to tag a new version?

I wasn't aware that PR #18 had already fixed issue #17
Also I think PR #19 was unnecessary as zero and one will work after PR #18
This may have come about because the updates were not in released version.
If you like I can create a PR to remove explicit methods for one and zero (but of course leave the tests)
In any case, releasing a new version would be good.

compare with Bailey package

The most popular implementation of this sort of thing seems to be Bailey's QD package (which also includes double-double arithmetic).

Would be good to do some comparison with that package in terms of performance, accuracy, etcetera, to see how this Julia package compares, in order to have some confidence in this implementation.

See also this paper for a good review of the state of the art.

Replace double by Double

I tried to do this but haven't got it to work yet.

Constructor of Double for non-Floats

One can easily create a Double from a Real using Double(x), but you can't use Double{Float64}(x).

I understand Double is not yet a drop-in replacement datatype, but the latter could be useful when changing datatypes throughout some section of code.

I don't fully understand what would be required, but could one of the following work?

Double{T<:AbstractFloat}(x::Real) = convert(Double{T}, T(x))

Or maybe its better to widen x argument to include all Reals rather than just AbstractFloats:

function convert{T<:AbstractFloat}(::Type{Double{T}}, x::Real) 
    z = convert(T, x) 
    Double(z, convert(T, x-z)) 
end

Or maybe this doesn't make sense at all.