Inline Float256 functions for speed

LICENSE.md

@@ -1,22 +1,40 @@
-The MIT License (MIT) <https://opensource.org/licenses/MIT>
-
-Copyright (c) 2018 Erik Schnetter
-
-Permission is hereby granted, free of charge, to any person obtaining
-a copy of this software and associated documentation files (the
-"Software"), to deal in the Software without restriction, including
-without limitation the rights to use, copy, modify, merge, publish,
-distribute, sublicense, and/or sell copies of the Software, and to
-permit persons to whom the Software is furnished to do so, subject to
-the following conditions:
-
-The above copyright notice and this permission notice shall be
-included in all copies or substantial portions of the Software.
-
-THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
-EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
-MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
-IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
-CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
-TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
-SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+The MIT License (MIT) <https://opensource.org/licenses/BSD-3-Clause>
+
+Copyright 2018 Erik Schnetter
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+1. Redistributions of source code must retain the above copyright
+   notice, this list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright
+   notice, this list of conditions and the following disclaimer in the
+   documentation and/or other materials provided with the
+   distribution.
+
+3. Neither the name of the copyright holder nor the names of its
+   contributors may be used to endorse or promote products derived
+   from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+
+
+This package contains source code based on the [QD
+package](http://crd-legacy.lbl.gov/~dhbailey/mpdist/) with its own set
+of authors, and which is under copyright by The Regents of the
+University of California, through Lawrence Berkeley National
+Laboratory. The original QD copyright notices are reproduced here in
+the files "QD-COPYING" and "QD-BSD-LBNL-License.doc".

Manifest.toml

@@ -1,6 +1,12 @@
 [[Base64]]
 uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 
+[[BenchmarkTools]]
+deps = ["JSON", "Printf", "Statistics", "Test"]
+git-tree-sha1 = "e686f1754227e4748259f400839b83a1e8773e02"
+uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+version = "0.4.1"
+
 [[BinaryProvider]]
 deps = ["Libdl", "Pkg", "SHA", "Test"]
 git-tree-sha1 = "055eb2690182ebc31087859c3dd8598371d3ef9e"
@@ -29,6 +35,12 @@ uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 deps = ["LinearAlgebra", "Markdown"]
 uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 
+[[JSON]]
+deps = ["Dates", "Distributed", "Mmap", "Sockets", "Test", "Unicode"]
+git-tree-sha1 = "1f7a25b53ec67f5e9422f1f551ee216503f4a0fa"
+uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
+version = "0.20.0"
+
 [[LibGit2]]
 uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
 

Project.toml

@@ -1,10 +1,12 @@
 name = "QD"
 uuid = "231d7462-ef32-11e8-24b1-7375831ffe71"
 authors = ["Erik Schnetter <schnetter@gmail.com>"]
-version = "0.1.0"
+version = "1.0.0"
 
 [deps]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 BinaryProvider = "b99e7846-7c00-51b0-8f62-c81ae34c0232"
 Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
 Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

QD-COPYING

@@ -0,0 +1,14 @@
+This work was supported by the Director, Office of Science, Division
+of Mathematical, Information, and Computational Sciences of the
+U.S. Department of Energy under contract numbers DE-AC03-76SF00098 and
+DE-AC02-05CH11231.
+
+Copyright (c) 2003-2009, The Regents of the University of California,
+through Lawrence Berkeley National Laboratory (subject to receipt of
+any required approvals from U.S. Dept. of Energy) All rights reserved. 
+
+By downloading or using this software you are agreeing to the modified
+BSD license that is in file "BSD-LBNL-License.doc" in the main ARPREC
+directory. If you wish to use the software for commercial purposes
+please contact the Technology Transfer Department at TTD@lbl.gov or
+call 510-286-6457."

README.md

@@ -55,28 +55,30 @@ julia> n=100;
 
 julia> A64 = rand(Float64, n, n);
 julia> B64 = @btime inv(A64);
-  317.306 μs (6 allocations: 129.27 KiB)
+  306.157 μs (6 allocations: 129.27 KiB)
 julia> norm(B64 * A64 - I)
-1.7329473539985953e-13
+1.4440108319436736e-12
 
 julia> A128 = rand(Float128, n, n);
 julia> B128 = @btime inv(A128);
-  5.170 ms (15 allocations: 474.42 KiB)
+  5.177 ms (15 allocations: 474.42 KiB)
 julia> norm(B128 * A128 - I)
-1.37702799848709929942060146535197e-28
+1.19496652320293845810089276691420e-28
 
 julia> A256 = rand(Float256, n, n);
 julia> B256 = @btime inv(A256);
-  153.904 ms (15 allocations: 946.14 KiB)
+  95.517 ms (15 allocations: 946.14 KiB)
 julia> norm(B256 * A256 - I)
-3.6861760395079082468609457967696988863273208819055135969958525024e-61
+3.4293792195278235762102730788487570191920941211322204339008314344e-61
 
 julia> Abig = rand(BigFloat, n, n);
 julia> Bbig = @btime inv(Abig);
-  268.576 ms (5334037 allocations: 285.10 MiB)
+  271.446 ms (5334037 allocations: 285.10 MiB)
 julia> norm(Bbig * Abig - I)
-2.866190771077340288694283208193327964200118337345917496588261302618102449593413e-74
+1.187140215116236991926392678171305368288858808728624740197177252518979524096639e-73
 ```
 
-(Times measured on a 2.8 GHz Intel Core i7.)
-
+(Times measured on a 2.8 GHz Intel Core i7.) This benchmarks inverts a
+`100 × 100` matrix. `Float128` is here about 50 times faster than
+`BigFloat` (but less than half as accurate), and `Float256` is about
+twice as fast and requires no memory allocations.

src/QD.jl

@@ -15,9 +15,226 @@ function __init__()
     check_deps()
 end
 
+
+
 # TODO: Check whether this is actually true
 const fma_is_fast = true
 
+
+
+# Basic (internal) Functions
+
+# Computes fl(a+b) and err(a+b).  Assumes |a| >= |b|.
+function quick_two_sum(a::T, b::T) where {T <: AbstractFloat}
+    s = a + b
+    err = b - (s - a)
+    s, err
+end
+
+# Computes fl(a-b) and err(a-b).  Assumes |a| >= |b|.
+function quick_two_diff(a::T, b::T) where {T <: AbstractFloat}
+    s = a - b
+    err = (a - s) - b
+    s, err
+end
+
+# Computes fl(a+b) and err(a+b).
+function two_sum(a::T, b::T) where {T <: AbstractFloat}
+    s = a + b
+    bb = s - a
+    err = (a - (s - bb)) + (b - bb)
+    s, err
+end
+
+# Computes fl(a-b) and err(a-b).
+function two_diff(a::T, b::T) where {T <: AbstractFloat}
+    s = a - b
+    bb = s - a
+    err = (a - (s - bb)) - (b + bb)
+    s, err
+end
+
+# Computes high word and lo word of a
+function split(a::T) where {T <: AbstractFloat}
+    _QD_SPLITTER = 134217729.0               # = 2^27 + 1
+    _QD_SPLIT_THRESH = 6.69692879491417e+299 # = 2^996
+    if a > _QD_SPLIT_THRESH || a < -_QD_SPLIT_THRESH
+        a *= 3.7252902984619140625e-09 # 2^-28
+        temp = _QD_SPLITTER * a
+        hi = temp - (temp - a)
+        lo = a - hi
+        hi *= 268435456.0       # 2^28
+        lo *= 268435456.0       # 2^28
+    else
+        temp = _QD_SPLITTER * a
+        hi = temp - (temp - a)
+        lo = a - hi
+    end
+    hi, lo
+end
+
+# Computes fl(a*b) and err(a*b).
+function two_prod(a::T, b::T) where {T <: AbstractFloat}
+    if fma_is_fast
+        p = a * b
+        err = fma(a, b, -p)
+    else
+        p = a * b
+        a_hi, a_lo = split(a)
+        b_hi, b_lo = split(b)
+        err = ((a_hi * b_hi - p) + a_hi * b_lo + a_lo * b_hi) + a_lo * b_lo
+    end
+    p, err
+end
+
+# Computes fl(a*a) and err(a*a).  Faster than the above method.
+function two_sqr(a::T) where {T <: AbstractFloat}
+    if fma_is_fast
+        p = a * a
+        err = fma(a, a, - p)
+        p, err
+    else
+        q = a * a
+        hi, lo = split(a)
+        err = ((hi * hi - q) + 2.0 * hi * lo) + lo * lo
+        q, err
+    end
+end
+
+
+
+# Additions
+
+function three_sum(a::T, b::T, c::T) where {T <: AbstractFloat}
+    t1, t2 = two_sum(a, b)
+    a , t3 = two_sum(c, t1)
+    b , c  = two_sum(t2, t3)
+    a, b, c
+end
+
+function three_sum2(a::T, b::T, c::T) where {T <: AbstractFloat}
+    t1, t2 = two_sum(a, b)
+    a , t3 = two_sum(c, t1)
+    b = t2 + t3
+    a, b, c
+end
+
+# Renormalization
+function quick_renorm(c0::T, c1::T, c2::T, c3::T, c4::T) where
+        {T <: AbstractFloat}
+    s , t3 = quick_two_sum(c3, c4)
+    s , t2 = quick_two_sum(c2, s )
+    s , t1 = quick_two_sum(c1, s )
+    c0, t0 = quick_two_sum(c0, s )
+
+    s , t2 = quick_two_sum(t2, t3)
+    s , t1 = quick_two_sum(t1, s )
+    c1, t0 = quick_two_sum(t0, s )
+
+    s , t1 = quick_two_sum(t1, t2)
+    c2, t0 = quick_two_sum(t0, s )
+
+    c3 = t0 + t1
+
+    c0, c1, c2, c3, c4
+end
+
+function renorm(c0::T, c1::T, c2::T, c3::T) where {T <: AbstractFloat}
+    s2 = 0.0
+    s3 = 0.0
+
+    isinf(c0) && return c0, c1, c2, c3
+
+    s0, c3 = quick_two_sum(c2, c3)
+    s0, c2 = quick_two_sum(c1, s0)
+    c0, c1 = quick_two_sum(c0, s0)
+
+    s0 = c0
+    s1 = c1
+    if s1 != 0.0
+        s1, s2 = quick_two_sum(s1, c2)
+        if s2 != 0.0
+            s2, s3 = quick_two_sum(s2, c3)
+        else
+            s1, s2 = quick_two_sum(s1, c3)
+        end
+    else
+        s0, s1 = quick_two_sum(s0, c2)
+        if s1 != 0.0
+            s1, s2 = quick_two_sum(s1, c3)
+        else
+            s0, s1 = quick_two_sum(s0, c3)
+        end
+    end
+
+    c0 = s0
+    c1 = s1
+    c2 = s2
+    c3 = s3
+
+    c0, c1, c2, c3
+end
+
+function renorm(c0::T, c1::T, c2::T, c3::T, c4::T) where {T <: AbstractFloat}
+    s2 = 0.0
+    s3 = 0.0
+
+    isinf(c0) && return c0, c1, c2, c3, c4
+
+    s0, c4 = quick_two_sum(c3, c4)
+    s0, c3 = quick_two_sum(c2, s0)
+    s0, c2 = quick_two_sum(c1, s0)
+    c0, c1 = quick_two_sum(c0, s0)
+
+    s0 = c0
+    s1 = c1
+
+    if s1 != 0.0
+        s1, s2 = quick_two_sum(s1, c2)
+        if s2 != 0.0
+            s2, s3 = quick_two_sum(s2, c3)
+            if s3 != 0.0
+                s3 += c4
+            else
+                s2, s3 = quick_two_sum(s2, c4)
+            end
+        else
+            s1, s2 = quick_two_sum(s1, c3)
+            if s2 != 0.0
+                s2, s3 = quick_two_sum(s2, c4)
+            else
+                s1, s2 = quick_two_sum(s1, c4)
+            end
+        end
+    else
+        s0, s1 = quick_two_sum(s0, c2)
+        if s1 != 0.0
+            s1, s2 = quick_two_sum(s1, c3)
+            if s2 != 0.0
+                s2, s3 = quick_two_sum(s2, c4)
+            else
+                s1, s2 = quick_two_sum(s1, c4)
+            end
+        else
+            s0, s1 = quick_two_sum(s0, c3)
+            if s1 != 0.0
+                s1, s2 = quick_two_sum(s1, c4)
+            else
+                s0, s1 = quick_two_sum(s0, c4)
+            end
+        end
+    end
+
+    c0 = s0
+    c1 = s1
+    c2 = s2
+    c3 = s3
+
+    c0, c1, c2, c3, c4
+end
+
+
+
 include("QD128.jl")
 include("QD256.jl")