Skip to content
Permalink

Comparing changes

Choose two branches to see what’s changed or to start a new pull request. If you need to, you can also or learn more about diff comparisons.

Open a pull request

Create a new pull request by comparing changes across two branches. If you need to, you can also . Learn more about diff comparisons here.
base repository: crystal-lang/crystal
Failed to load repositories. Confirm that selected base ref is valid, then try again.
Loading
base: 0a375f504ca5
Choose a base ref
...
head repository: crystal-lang/crystal
Failed to load repositories. Confirm that selected head ref is valid, then try again.
Loading
compare: 32f583fbf091
Choose a head ref
  • 5 commits
  • 11 files changed
  • 1 contributor

Commits on Jun 4, 2017

  1. Add FloatPrinter based on Grisu3 

    This improves the speed of transforming floats to their string
representation. It is based on the 2004 paper "Printing Floating-Point
Numbers Quickly and Accurately with Integers" by Florian Loitsch[1].

Most of the code is a port from the BSD-licensed C++ project
"double-conversion"[2], which was extracted from the V8 engine.

The Grisu3 algorithm is fast because it deals only with fixed-sized
integer arithmetic. It takes advantage extra bits leftover from the
53-bit significand in a 64 bit number to help find the optimal string
representation. However this only works for 95.5% of floats and it
rejects the remaining 0.5%. Rejected numbers still need to be printed
with some other, slower method.

1: http://www.cs.tufts.edu/~nr/cs257/archive/florian-loitsch/printf.pdf
2: https://github.com/google/double-conversion
    @will @asterite
    will authored and asterite committed Jun 4, 2017

    Unverified

    This user has not yet uploaded their public signing key.
    GPG key ID: CEF2BBFE18BD0E67
    Learn about vigilant mode
    Copy the full SHA
    08dcf53 View commit details

  2. FloatPrinter: make docs and interface more consistent

    @will @asterite
    will authored and asterite committed Jun 4, 2017
    Copy the full SHA
    d8c8417 View commit details

  3. FloatPrinter: fix -NaN test, fix formatting

    @will @asterite
    will authored and asterite committed Jun 4, 2017
    Copy the full SHA
    5e1df56 View commit details

  4. FloatPrinter: add Float32 support

    @will @asterite
    will authored and asterite committed Jun 4, 2017
    Copy the full SHA
    576f815 View commit details

  5. FloatPrinter: add license to specs

    @will @asterite
    will authored and asterite committed Jun 4, 2017
    Copy the full SHA
    32f583f View commit details
178 changes: 178 additions & 0 deletions spec/std/float_printer/diy_fp_spec.cr
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,178 @@
# The example numbers in these specs are ported from the C++
# "double-conversions" library. The following is their license:
# Copyright 2012 the V8 project authors. All rights reserved.
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are
# met:
#
# * Redistributions of source code must retain the above copyright
# notice, this list of conditions and the following disclaimer.
# * 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.
# * Neither the name of Google Inc. 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
# OWNER 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.
require "spec"
require "float_printer/diy_fp"
include FloatPrinter

describe DiyFP do
it "multiply" do
fp1 = DiyFP.new(3_u64, 0)
fp2 = DiyFP.new(2_u64, 0)
prod = fp1 * fp2

prod.frac.should eq 0
prod.exp.should eq 64
end

it "multiply" do
fp1 = DiyFP.new(0x8000000000000000, 11)
fp2 = DiyFP.new(2_u64, 13)
prod = fp1 * fp2

prod.frac.should eq 1
prod.exp.should eq 11 + 13 + 64
end

it "multiply rounding" do
fp1 = DiyFP.new(0x8000000000000001_u64, 11)
fp2 = DiyFP.new(1_u64, 13)
prod = fp1 * fp2

prod.frac.should eq 1
prod.exp.should eq 11 + 13 + 64
end

it "multiply rounding" do
fp1 = DiyFP.new(0x7fffffffffffffff_u64, 11)
fp2 = DiyFP.new(1_u64, 13)
prod = fp1 * fp2

prod.frac.should eq 0
prod.exp.should eq 11 + 13 + 64
end

it "multiply big numbers" do
fp1 = DiyFP.new(0xffffffffffffffff_u64, 11)
fp2 = DiyFP.new(0xffffffffffffffff_u64, 13)
prod = fp1 * fp2

prod.frac.should eq 0xfffffffffffffffe_u64
prod.exp.should eq 11 + 13 + 64
end

it "converts ordered 64" do
ordered = 0x0123456789ABCDEF_u64
f = pointerof(ordered).as(Float64*).value
f.should eq 3512700564088504e-318 # ensure byte order

fp = DiyFP.from_f(f)

fp.exp.should eq 0x12 - 0x3FF - 52
# The 52 mantissa bits, plus the implicit 1 in bit 52 as a UINT64.
fp.frac.should eq 0x0013456789ABCDEF
end

it "converts ordered 32" do
ordered = 0x01234567_u32
f = pointerof(ordered).as(Float32*).value
f.should eq(2.9988165487136453e-38_f32)

fp = DiyFP.from_f(f)

fp.exp.should eq 0x2 - 0x7F - 23
# The 23 mantissa bits, plus the implicit 1 in bit 24 as a uint32.

fp.frac.should eq 0xA34567
end

it "converts min f64" do
min = 0x0000000000000001_u64
f = pointerof(min).as(Float64*).value
f.should eq 5e-324 # ensure byte order

fp = DiyFP.from_f(f)

fp.exp.should eq -0x3FF - 52 + 1
# This is denormal, so no hidden bit
fp.frac.should eq 1
end

it "converts min f32" do
min = 0x00000001_u32
f = pointerof(min).as(Float32*).value
fp = DiyFP.from_f(f)

fp.exp.should eq -0x7F - 23 + 1
# This is a denormal; so no hidden bit.
fp.frac.should eq 1
end

it "converts max f64" do
max = 0x7fefffffffffffff_u64
f = pointerof(max).as(Float64*).value
f.should eq 1.7976931348623157e308 # ensure byte order

fp = DiyFP.from_f(f)

fp.exp.should eq 0x7FE - 0x3FF - 52
fp.frac.should eq 0x001fffffffffffff_u64
end

it "converts max f32" do
max = 0x7f7fffff_u64
f = pointerof(max).as(Float32*).value
f.should eq 3.4028234e38_f32 # ensure byte order

fp = DiyFP.from_f(f)

fp.exp.should eq 0xFE - 0x7F - 23
fp.frac.should eq 0x00ffffff_u64
end

it "normalizes ordered" do
ordered = 0x0123456789ABCDEF_u64
f = pointerof(ordered).as(Float64*).value

fp = DiyFP.from_f_normalized(f)

fp.exp.should eq 0x12 - 0x3FF - 52 - 11
fp.frac.should eq 0x0013456789ABCDEF_u64 << 11
end

it "normalizes min f64" do
min = 0x0000000000000001_u64
f = pointerof(min).as(Float64*).value

fp = DiyFP.from_f_normalized(f)

fp.exp.should eq -0x3FF - 52 + 1 - 63
# This is a denormal; so no hidden bit
fp.frac.should eq 0x8000000000000000
end

it "normalizes max f64" do
max = 0x7fefffffffffffff_u64
f = pointerof(max).as(Float64*).value

fp = DiyFP.from_f_normalized(f)

fp.exp.should eq 0x7FE - 0x3FF - 52 - 11
fp.frac.should eq 0x001fffffffffffff << 11
end
end
193 changes: 193 additions & 0 deletions spec/std/float_printer/grisu3_spec.cr
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,193 @@
# The example numbers in these specs are ported from the C++
# "double-conversions" library. The following is their license:
# Copyright 2012 the V8 project authors. All rights reserved.
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are
# met:
#
# * Redistributions of source code must retain the above copyright
# notice, this list of conditions and the following disclaimer.
# * 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.
# * Neither the name of Google Inc. 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
# OWNER 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.
require "spec"
require "float_printer/grisu3"
include FloatPrinter

private def test_grisu(v : UInt64)
test_grisu pointerof(v).as(Float64*).value
end

private def test_grisu(v : UInt32)
test_grisu pointerof(v).as(Float32*).value
end

private def test_grisu(v : Float64 | Float32)
buffer = StaticArray(UInt8, 128).new(0_u8)
status, decimal_exponent, length = Grisu3.grisu3(v, buffer.to_unsafe)
point = decimal_exponent + length
return status, point, String.new(buffer.to_unsafe)
end

describe "grisu3" do
context "float64" do
it "min float" do
status, point, str = test_grisu 5e-324
status.should eq true
str.should eq "5"
point.should eq -323
end

it "max float" do
status, point, str = test_grisu 1.7976931348623157e308
status.should eq true
str.should eq "17976931348623157"
point.should eq 309
end

it "point at end" do
status, point, str = test_grisu 4294967272.0
status.should eq true
str.should eq "4294967272"
point.should eq 10
end

it "large number" do
status, point, str = test_grisu 4.1855804968213567e298
status.should eq true
str.should eq "4185580496821357"
point.should eq 299
end

it "small number" do
status, point, str = test_grisu 5.5626846462680035e-309
status.should eq true
str.should eq "5562684646268003"
point.should eq -308
end

it "another no point move" do
status, point, str = test_grisu 2147483648.0
status.should eq true
str.should eq "2147483648"
point.should eq 10
end

it "failure case" do
# grisu should not be able to do this number
# this number is reused to ensure the fallback works
status, point, str = test_grisu 3.5844466002796428e+298
status.should eq false
str.should_not eq "35844466002796428"
end

it "smallest normal" do
status, point, str = test_grisu 0x0010000000000000_u64
status.should eq true
str.should eq "22250738585072014"
point.should eq -307
end

it "largest denormal" do
status, point, str = test_grisu 0x000FFFFFFFFFFFFF_u64
status.should eq true
str.should eq "2225073858507201"
point.should eq -307
end
end

context "float32" do
it "min" do
status, point, str = test_grisu 1e-45_f32
status.should eq true
str.should eq "1"
point.should eq -44
end

it "max" do
status, point, str = test_grisu 3.4028234e38_f32
status.should eq true
str.should eq "34028235"
point.should eq 39
end

it "general whole number, rounding" do
status, point, str = test_grisu 4294967272.0_f32
status.should eq true
str.should eq "42949673"
point.should eq 10
end

it "general whole number, rounding" do
status, point, str = test_grisu 4294967272.0_f32
status.should eq true
str.should eq "42949673"
point.should eq 10
end

it "large number, rounding" do
status, point, str = test_grisu 3.32306998946228968226e35_f32
status.should eq true
str.should eq "332307"
point.should eq 36
end

it "small number" do
status, point, str = test_grisu 1.2341e-41_f32
status.should eq true
str.should eq "12341"
point.should eq -40
end

it "general no rounding" do
status, point, str = test_grisu 3.3554432e7_f32
status.should eq true
str.should eq "33554432"
point.should eq 8
end

it "general with rounding up" do
status, point, str = test_grisu 3.26494756798464e14_f32
status.should eq true
str.should eq "32649476"
point.should eq 15
end

it "general with rounding down" do
status, point, str = test_grisu 3.91132223637771935344e37_f32
status.should eq true
str.should eq "39113222"
point.should eq 38
end

it "smallest normal" do
status, point, str = test_grisu 0x00800000_u32
status.should eq true
str.should eq "11754944"
point.should eq -37
end

it "largest denormal" do
status, point, str = test_grisu 0x007FFFFF_u32
status.should eq true
str.should eq "11754942"
point.should eq -37
end
end
end
153 changes: 153 additions & 0 deletions spec/std/float_printer/ieee_spec.cr
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,153 @@
# The example numbers in these specs are ported from the C++
# "double-conversions" library. The following is their license:
# Copyright 2012 the V8 project authors. All rights reserved.
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are
# met:
#
# * Redistributions of source code must retain the above copyright
# notice, this list of conditions and the following disclaimer.
# * 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.
# * Neither the name of Google Inc. 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
# OWNER 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.
require "spec"
require "float_printer/diy_fp"
require "float_printer/ieee"
include FloatPrinter

private def gen_bound(v : UInt64)
f = pointerof(v).as(Float64*).value
gen_bound(f)
end

private def gen_bound(v : UInt32)
f = pointerof(v).as(Float32*).value
gen_bound(f)
end

private def gen_bound(v : Float64 | Float32)
fp = DiyFP.from_f_normalized(v)
b = IEEE.normalized_boundaries(v)
b[:minus].exp.should eq fp.exp
b[:plus].exp.should eq fp.exp

return fp.frac, b[:minus].frac, b[:plus].frac
end

describe "Float64 boundaires" do
it "boundaries 1.5" do
fp, mi, pl = gen_bound(1.5)
# 1.5 does not have a significand of the form 2^p (for some p).
# Therefore its boundaries are at the same distance.
(pl - fp).should eq(fp - mi)
(fp - mi).should eq(1 << 10)
end

it "boundaries 1.0" do
fp, mi, pl = gen_bound(1.0)
# 1.0 does have a significand of the form 2^p (for some p).
# Therefore its lower boundary is twice as close as the upper boundary.
(pl - fp).should be > fp - mi
(fp - mi).should eq 1 << 9
(pl - fp).should eq 1 << 10
end

it "boundaries min float64" do
fp, mi, pl = gen_bound(0x0000000000000001_u64)
# min-value does not have a significand of the form 2^p (for some p).
# Therefore its boundaries are at the same distance.
(pl - fp).should eq fp - mi
(fp - mi).should eq 1_u64 << 62
end

it "boundaries min normal f64" do
fp, mi, pl = gen_bound(0x0010000000000000_u64)
# Even though the significand is of the form 2^p (for some p), its boundaries
# are at the same distance. (This is the only exception).
(fp - mi).should eq(pl - fp)
(fp - mi).should eq(1 << 10)
end

it "boundaries max denormal f64" do
fp, mi, pl = gen_bound(0x000FFFFFFFFFFFFF_u64)

(fp - mi).should eq(pl - fp)
(fp - mi).should eq(1 << 11)
end

it "boundaries max f64" do
fp, mi, pl = gen_bound(0x7fEFFFFFFFFFFFFF_u64)
# max-value does not have a significand of the form 2^p (for some p).
# Therefore its boundaries are at the same distance.
(fp - mi).should eq(pl - fp)
(fp - mi).should eq(1 << 10)
end
end

describe "Float32 boundaires" do
it "boundaries 1.5" do
fp, mi, pl = gen_bound(1.5_f32)
# 1.5 does not have a significand of the form 2^p (for some p).
# Therefore its boundaries are at the same distance.
(pl - fp).should eq(fp - mi)
# Normalization shifts the significand by 8 bits. Add 32 bits for the bigger
# data-type, and remove 1 because boundaries are at half a ULP.
(fp - mi).should eq(1_u64 << 39)
end

it "boundaries 1.0" do
fp, mi, pl = gen_bound(1.0_f32)
# 1.0 does have a significand of the form 2^p (for some p).
# Therefore its lower boundary is twice as close as the upper boundary.
(pl - fp).should be > fp - mi
(fp - mi).should eq(1_u64 << 38)
(pl - fp).should eq(1_u64 << 39)
end

it "min Float32" do
fp, mi, pl = gen_bound(0x00000001_u32)
# min-value does not have a significand of the form 2^p (for some p).
# Therefore its boundaries are at the same distance.
(pl - fp).should eq(fp - mi)
# Denormals have their boundaries much closer.
(fp - mi).should eq(1_u64 << 62)
end

it "smallest normal 32" do
fp, mi, pl = gen_bound(0x00800000_u32)
# Even though the significand is of the form 2^p (for some p), its boundaries
# are at the same distance. (This is the only exception).
(pl - fp).should eq(fp - mi)
(fp - mi).should eq(1_u64 << 39)
end

it "largest denormal 32" do
fp, mi, pl = gen_bound(0x007FFFFF_u32)
(pl - fp).should eq(fp - mi)
(fp - mi).should eq(1_u64 << 40)
end

it "max Float32" do
fp, mi, pl = gen_bound(0x7F7FFFFF_u32)
# max-value does not have a significand of the form 2^p (for some p).
# Therefore its boundaries are at the same distance.
(pl - fp).should eq(fp - mi)
(fp - mi).should eq(1_u64 << 39)
end
end
160 changes: 160 additions & 0 deletions spec/std/float_printer_spec.cr
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,160 @@
# The example numbers in these specs are ported from the C++
# "double-conversions" library. The following is their license:
# Copyright 2012 the V8 project authors. All rights reserved.
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are
# met:
#
# * Redistributions of source code must retain the above copyright
# notice, this list of conditions and the following disclaimer.
# * 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.
# * Neither the name of Google Inc. 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
# OWNER 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.
require "spec"
require "float_printer"
include FloatPrinter

private def float_to_s(v)
String.build(22) do |buff|
FloatPrinter.print(v, buff)
end
end

private def test_str(s, file = __FILE__, line = __LINE__)
[s, "-#{s}"].each do |str|
test_pair str.to_f64, str, file, line
end
end

private def test_pair(v : UInt64, str, file = __FILE__, line = __LINE__)
d = pointerof(v).as(Float64*).value
test_pair(d, str, file, line)
end

private def test_pair(v : UInt32, str, file = __FILE__, line = __LINE__)
d = pointerof(v).as(Float32*).value
test_pair(d, str, file, line)
end

private def test_pair(v : Float64 | Float32, str, file = __FILE__, line = __LINE__)
float_to_s(v).should eq(str), file, line
end

describe "#print Float64" do
it { test_str "0.0" }

it { test_str "Infinity" }

it { test_pair 0x7ff8000000000000_u64, "NaN" }
it { test_pair 0xfff8000000000000_u64, "-NaN" }

it { test_str "0.01" }
it { test_str "0.1" }
it { test_str "1.0" }
it { test_str "1.2" }
it { test_str "123.456" }

it { test_str "1.0e+234" }
it { test_str "1.1e+234" }
it { test_str "1.0e-234" }

it { test_str "111111111111111.0" }
it { test_str "111111111111111.1" }

it { test_pair 0.001, "0.001" }
it { test_pair 0.0001, "0.0001" }
it { test_pair 0.00001, "1.0e-5" }
it { test_pair 0.000001, "1.0e-6" }
it { test_pair -0.0001, "-0.0001" }
it { test_pair -0.00001, "-1.0e-5" }
it { test_pair -12345e23, "-1.2345e+27" }

it { test_pair 10.0, "10.0" }
it { test_pair 1100.0, "1100.0" }

it { test_pair 100000000000000.0, "100000000000000.0" }
it { test_pair 1000000000000000.0, "1.0e+15" }
it { test_pair 1111111111111111.0, "1.111111111111111e+15" }

it "min float64" do
test_pair 5e-324, "5.0e-324"
end

it "max float64" do
test_pair 1.7976931348623157e308, "1.7976931348623157e+308"
end

it "large number, rounded" do
test_pair 4.1855804968213567e298, "4.185580496821357e+298"
end

it "small number, rounded" do
test_pair 5.5626846462680035e-309, "5.562684646268003e-309"
end

it "falure case" do
# grisu cannot do this number, so it should fall back to libc
test_pair 3.5844466002796428e+298, "3.5844466002796428e+298"
end

it "smallest normal" do
test_pair 0x0010000000000000_u64, "2.2250738585072014e-308"
end

it "largest denormal" do
test_pair 0x000FFFFFFFFFFFFF_u64, "2.225073858507201e-308"
end
end

describe "#print Float32" do
it { test_pair 0_f32, "0.0" }
it { test_pair -0_f32, "-0.0" }
it { test_pair Float32::INFINITY, "Infinity" }
it { test_pair -Float32::INFINITY, "-Infinity" }
it { test_pair 0x7fc00000_u32, "NaN" }
it { test_pair 0xffc80000_u32, "-NaN" }
it { test_pair 0.000001_f32, "1.0e-6" }
it { test_pair -0.0001_f32, "-0.0001" }
it { test_pair -0.00001_f32, "-1.0e-5" }
it { test_pair -12345e23_f32, "-1.2345e+27" }
it { test_pair 100000000000000.0_f32, "100000000000000.0" }
it { test_pair 1000000000000000.0_f32, "1.0e+15" }
it { test_pair 1111111111111111.0_f32, "1.1111111e+15" }
it { test_pair -3.9292015898194142585311918e-10_f32, "-3.9292017e-10" }

it "largest float" do
test_pair 3.4028234e38_f32, "3.4028235e+38"
end

it "largest normal" do
test_pair 0x7f7fffff_u32, "3.4028235e+38"
end

it "smallest positive normal" do
test_pair 0x00800000_u32, "1.1754944e-38"
end

it "largest denormal" do
test_pair 0x007fffff_u32, "1.1754942e-38"
end

it "smallest positive denormal" do
test_pair 0x00000001_u32, "1.0e-45"
end
end
6 changes: 3 additions & 3 deletions spec/std/float_spec.cr
View file
Original file line number Diff line number Diff line change
@@ -129,7 +129,7 @@ describe "Float" do
0.65000000000000002.to_s.should eq("0.65")
1.234001.to_s.should eq("1.234001")
1.23499.to_s.should eq("1.23499")
1.23499999999999.to_s.should eq("1.235")
1.23499999999999999.to_s.should eq("1.235")
1.2345.to_s.should eq("1.2345")
1.23456.to_s.should eq("1.23456")
1.234567.to_s.should eq("1.234567")
@@ -141,7 +141,7 @@ describe "Float" do
1.23456789123.to_s.should eq("1.23456789123")
9525365.25.to_s.should eq("9525365.25")
12.9999.to_s.should eq("12.9999")
12.999999999999.to_s.should eq("13.0")
12.9999999999999999.to_s.should eq("13.0")
1.0.to_s.should eq("1.0")
2e20.to_s.should eq("2.0e+20")
1e-10.to_s.should eq("1.0e-10")
@@ -164,7 +164,7 @@ describe "Float" do
65432.1234567891e20.to_s.should eq("6.54321234567891e+24")
(1.0/0.0).to_s.should eq("Infinity")
(-1.0/0.0).to_s.should eq("-Infinity")
(0.99999999999999989).to_s.should eq("1.0")
(0.999999999999999989).to_s.should eq("1.0")
end

it "does to_s for f32" do
206 changes: 7 additions & 199 deletions src/float.cr
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
require "c/stdio"
require "c/string"
require "float_printer"

# Float is the base type of all floating point numbers.
#
@@ -138,37 +139,13 @@ struct Float32
end

def to_s
String.new(22) do |buffer|
len = to_s_internal(buffer)
{len, len}
String.build(22) do |buffer|
FloatPrinter.print(self, buffer)
end
end

def to_s(io : IO)
chars = StaticArray(UInt8, 22).new(0_u8)
len = to_s_internal(chars.to_unsafe)
io.write_utf8 chars.to_slice[0, len]
end

private def to_s_internal(buffer)
LibC.snprintf(buffer, 22, "%g", to_f64)
len = LibC.strlen(buffer)

# If it's "inf", return "Infinity"
if buffer[0] === 'i'
buffer.copy_from("Infinity".to_unsafe, 8)
len = 8
return len
end

# If it's "-inf", return "-inf"
if len >= 2 && buffer[1] === 'i'
buffer.copy_from("-Infinity".to_unsafe, 9)
len = 9
return len
end

len
FloatPrinter.print(self, io)
end

def hash
@@ -221,182 +198,13 @@ struct Float64
end

def to_s
String.new(28) do |buffer|
len = to_s_internal(buffer)
{len, len}
String.build(22) do |buffer|
FloatPrinter.print(self, buffer)
end
end

def to_s(io : IO)
chars = StaticArray(UInt8, 28).new(0_u8)
len = to_s_internal(chars.to_unsafe)
io.write_utf8 chars.to_slice[0, len]
end

private def to_s_internal(buffer)
LibC.snprintf(buffer, 28, "%.17g", self)
len = LibC.strlen(buffer)
slice = Slice.new(buffer, len)

# If it's "inf", return "Infinity"
if buffer[0] === 'i'
buffer.copy_from("Infinity".to_unsafe, 8)
len = 8
return len
end

# If it's "-inf", return "-inf"
if len >= 2 && buffer[1] === 'i'
buffer.copy_from("-Infinity".to_unsafe, 9)
len = 9
return len
end

# Check if we have a run of zeros or nines after
# the decimal digit. If so, we remove them
# (rounding, if needed). This is a very simple
# (and probably inefficient) algorithm, but a good
# one is much longer and harder to do: we can probably
# do that later.
dot_index = slice.index('.'.ord.to_u8)

# If there's no dot add ".0" to it
unless dot_index
# Check if we have an 'e'
e_index = slice.index('e'.ord.to_u8)

# If there's an "e", we must move it to the right
if e_index
(buffer + e_index).move_to(buffer + e_index + 2, len - e_index)
buffer[e_index] = '.'.ord.to_u8
buffer[e_index + 1] = '0'.ord.to_u8
else
buffer[len] = '.'.ord.to_u8
buffer[len + 1] = '0'.ord.to_u8
end
len += 2
return len
end

original_len = len
index = dot_index

# Also return if the dot is the last char (shouldn't happen)
return len if index + 1 == len

# And also return if the length is less than 7
# (digit, dot plus at least 5 digits)
return len if len < 7

this_run = 0 # number of chars in this run
max_run = 0 # maximum consecutive chars of a run
run_byte = 0_u8 # the run character
last_run_start = -1 # where did the last run start
max_run_byte = 0_u8 # the byte of the last run
max_run_start = -1 # the index where the maximum run starts
max_run_end = -1 # the index where the maximum run ends
e_index = nil

while index < len
byte = slice.to_unsafe[index]

if byte == run_byte
this_run += 1
if this_run >= max_run
max_run = this_run
max_run_byte = byte
max_run_start = last_run_start
max_run_end = index
end
elsif byte === '0' || byte === '9'
run_byte = byte
last_run_byte = byte
last_run_start = index
this_run = 1
if this_run >= max_run
max_run = this_run
max_run_byte = byte
max_run_start = index
max_run_end = index
end
elsif byte === 'e'
e_index = index
break
else
run_byte = 0_u8
this_run = 0
end

index += 1
end

if e_index
# If we have an 'e', remove a sequence of 0s or 9s or
# any length, as long as they are near the end
tolerance = 1
else
# The more digits we have to the left of the dot,
# the less run digits we adjust
tolerance = case dot_index
when 0..8 then 5
when 9..11 then 4
when 12 then 3
when 13 then 2
else 1
end
end

if e_index
len = e_index
end

# If the maximum run ends one or two chars before
# the end of the string, we replace the run
if (len - 3 <= max_run_end < len) && max_run >= tolerance
case max_run_byte
when '0'
# Just trim
len = max_run_start
when '9'
# Need to add one and carry to the left
len = max_run_start
index = len - 1
while index >= 0
byte = slice.to_unsafe[index]
case byte
when '.'
# Nothing, continue
when '9'
# If this is the last char, remove it,
# otherwise turn into a zero
if index == len
len -= 1
else
slice.to_unsafe[index] = '0'.ord.to_u8
end
else
slice.to_unsafe[index] = byte + 1
break
end
index -= 1
end
end
end

# Add a zero if the last char is a dot
if slice.to_unsafe[len - 1] === '.'
slice.to_unsafe[len] = '0'.ord.to_u8
len += 1
end

# Add back the e and what's to the right of it
if e_index
e_len = original_len - e_index
(buffer + e_index).move_to(buffer + len, e_len)
len += e_len
end

len
FloatPrinter.print(self, io)
end

def hash
97 changes: 97 additions & 0 deletions src/float_printer.cr
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,97 @@
require "float_printer/*"

# FloatPrinter is based on Grisu3 algorithm described in the 2004 paper
# "Printing Floating-Point Numbers Quickly and Accurately with Integers" by
# Florian Loitsch.
module FloatPrinter
extend self
BUFFER_SIZE = 128

# Converts Float *v* to a string representation and prints it onto *io*
#
# It is used by `Float64#to_s` and it is probably not necessary to use
# this directly.
def print(v : Float64 | Float32, io : IO)
d = IEEE.to_uint(v)

if IEEE.sign(d) < 0
io << '-'
v = -v
end

if v == 0.0
io << "0.0"
elsif IEEE.special?(d)
if IEEE.inf?(d)
io << "Infinity"
else
io << "NaN"
end
else
internal(v, io)
end
end

private def internal(v : Float64 | Float32, io : IO)
buffer = StaticArray(UInt8, BUFFER_SIZE).new(0_u8)
success, decimal_exponent, length = Grisu3.grisu3(v, buffer.to_unsafe)

unless success
# grisu3 does not work for ~0.5% of floats
# when this happens, fallback to another, slower approach
if v.class == Float64
LibC.snprintf(buffer.to_unsafe, BUFFER_SIZE, "%.17g", v)
else
LibC.snprintf(buffer.to_unsafe, BUFFER_SIZE, "%g", v)
end
len = LibC.strlen(buffer)
io.write_utf8 buffer.to_slice[0, len]
return
end

point = decimal_exponent + length

exp = point
exp_mode = point > 15 || point < -3
point = 1 if exp_mode

# add leading zero
io << '0' if point < 1

i = 0

# add integer part digits
if decimal_exponent > 0 && !exp_mode
# whole number but not big enough to be exp form
io.write_utf8 buffer.to_slice[i, length - i]
i = length
(point - length).times { io << '0' }
elsif i < point
io.write_utf8 buffer.to_slice[i, point - i]
i = point
end

io << '.'

# add leading zeros after point
if point < 0
(-point).times { io << '0' }
end

# add fractional part digits
io.write_utf8 buffer.to_slice[i, length - i]
i = length

# print trailing 0 if whole number or exp notation of power of ten
if (decimal_exponent >= 0 && !exp_mode) || (exp != point && length == 1)
io << '0'
end

# exp notation
if exp != point
io << 'e'
io << '+' if exp > 0
(exp - 1).to_s(io)
end
end
end
174 changes: 174 additions & 0 deletions src/float_printer/cached_powers.cr
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,174 @@
# CachedPowers is ported from the C++ "double-conversions" library.
# The following is their license:
# Copyright 2006-2008 the V8 project authors. All rights reserved.
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are
# met:
#
# * Redistributions of source code must retain the above copyright
# notice, this list of conditions and the following disclaimer.
# * 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.
# * Neither the name of Google Inc. 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
# OWNER 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.

module FloatPrinter::CachedPowers
record Power, significand : UInt64, binary_exp : Int16, decimal_exp : Int16
# The minimal and maximal target exponent define the range of w's binary
# exponent, where 'w' is the result of multiplying the input by a cached power
# of ten.
#
# A different range might be chosen on a different platform, to optimize digit
# generation, but a smaller range requires more powers of ten to be cached.
MIN_TARGET_EXP = -60
MAX_TARGET_EXP = -32

CACHED_POWER_OFFSET = 348 # -1 * the first decimal_exp

# Not all powers of ten are cached. The decimal exponent of two neighboring
# cached numbers will differ by `CACHED_EXP_STEP`
CACHED_EXP_STEP = 8
MIN_CACHED_EXP = -348
MAX_CACHED_EXP = 340

D_1_LOG2_10 = 0.30102999566398114 # 1 / lg(10)

PowCache = [
{0xfa8fd5a0081c0288_u64, -1220_i16, -348_i16},
{0xbaaee17fa23ebf76_u64, -1193_i16, -340_i16},
{0x8b16fb203055ac76_u64, -1166_i16, -332_i16},
{0xcf42894a5dce35ea_u64, -1140_i16, -324_i16},
{0x9a6bb0aa55653b2d_u64, -1113_i16, -316_i16},
{0xe61acf033d1a45df_u64, -1087_i16, -308_i16},
{0xab70fe17c79ac6ca_u64, -1060_i16, -300_i16},
{0xff77b1fcbebcdc4f_u64, -1034_i16, -292_i16},
{0xbe5691ef416bd60c_u64, -1007_i16, -284_i16},
{0x8dd01fad907ffc3c_u64, -980_i16, -276_i16},
{0xd3515c2831559a83_u64, -954_i16, -268_i16},
{0x9d71ac8fada6c9b5_u64, -927_i16, -260_i16},
{0xea9c227723ee8bcb_u64, -901_i16, -252_i16},
{0xaecc49914078536d_u64, -874_i16, -244_i16},
{0x823c12795db6ce57_u64, -847_i16, -236_i16},
{0xc21094364dfb5637_u64, -821_i16, -228_i16},
{0x9096ea6f3848984f_u64, -794_i16, -220_i16},
{0xd77485cb25823ac7_u64, -768_i16, -212_i16},
{0xa086cfcd97bf97f4_u64, -741_i16, -204_i16},
{0xef340a98172aace5_u64, -715_i16, -196_i16},
{0xb23867fb2a35b28e_u64, -688_i16, -188_i16},
{0x84c8d4dfd2c63f3b_u64, -661_i16, -180_i16},
{0xc5dd44271ad3cdba_u64, -635_i16, -172_i16},
{0x936b9fcebb25c996_u64, -608_i16, -164_i16},
{0xdbac6c247d62a584_u64, -582_i16, -156_i16},
{0xa3ab66580d5fdaf6_u64, -555_i16, -148_i16},
{0xf3e2f893dec3f126_u64, -529_i16, -140_i16},
{0xb5b5ada8aaff80b8_u64, -502_i16, -132_i16},
{0x87625f056c7c4a8b_u64, -475_i16, -124_i16},
{0xc9bcff6034c13053_u64, -449_i16, -116_i16},
{0x964e858c91ba2655_u64, -422_i16, -108_i16},
{0xdff9772470297ebd_u64, -396_i16, -100_i16},
{0xa6dfbd9fb8e5b88f_u64, -369_i16, -92_i16},
{0xf8a95fcf88747d94_u64, -343_i16, -84_i16},
{0xb94470938fa89bcf_u64, -316_i16, -76_i16},
{0x8a08f0f8bf0f156b_u64, -289_i16, -68_i16},
{0xcdb02555653131b6_u64, -263_i16, -60_i16},
{0x993fe2c6d07b7fac_u64, -236_i16, -52_i16},
{0xe45c10c42a2b3b06_u64, -210_i16, -44_i16},
{0xaa242499697392d3_u64, -183_i16, -36_i16},
{0xfd87b5f28300ca0e_u64, -157_i16, -28_i16},
{0xbce5086492111aeb_u64, -130_i16, -20_i16},
{0x8cbccc096f5088cc_u64, -103_i16, -12_i16},
{0xd1b71758e219652c_u64, -77_i16, -4_i16},
{0x9c40000000000000_u64, -50_i16, 4_i16},
{0xe8d4a51000000000_u64, -24_i16, 12_i16},
{0xad78ebc5ac620000_u64, 3_i16, 20_i16},
{0x813f3978f8940984_u64, 30_i16, 28_i16},
{0xc097ce7bc90715b3_u64, 56_i16, 36_i16},
{0x8f7e32ce7bea5c70_u64, 83_i16, 44_i16},
{0xd5d238a4abe98068_u64, 109_i16, 52_i16},
{0x9f4f2726179a2245_u64, 136_i16, 60_i16},
{0xed63a231d4c4fb27_u64, 162_i16, 68_i16},
{0xb0de65388cc8ada8_u64, 189_i16, 76_i16},
{0x83c7088e1aab65db_u64, 216_i16, 84_i16},
{0xc45d1df942711d9a_u64, 242_i16, 92_i16},
{0x924d692ca61be758_u64, 269_i16, 100_i16},
{0xda01ee641a708dea_u64, 295_i16, 108_i16},
{0xa26da3999aef774a_u64, 322_i16, 116_i16},
{0xf209787bb47d6b85_u64, 348_i16, 124_i16},
{0xb454e4a179dd1877_u64, 375_i16, 132_i16},
{0x865b86925b9bc5c2_u64, 402_i16, 140_i16},
{0xc83553c5c8965d3d_u64, 428_i16, 148_i16},
{0x952ab45cfa97a0b3_u64, 455_i16, 156_i16},
{0xde469fbd99a05fe3_u64, 481_i16, 164_i16},
{0xa59bc234db398c25_u64, 508_i16, 172_i16},
{0xf6c69a72a3989f5c_u64, 534_i16, 180_i16},
{0xb7dcbf5354e9bece_u64, 561_i16, 188_i16},
{0x88fcf317f22241e2_u64, 588_i16, 196_i16},
{0xcc20ce9bd35c78a5_u64, 614_i16, 204_i16},
{0x98165af37b2153df_u64, 641_i16, 212_i16},
{0xe2a0b5dc971f303a_u64, 667_i16, 220_i16},
{0xa8d9d1535ce3b396_u64, 694_i16, 228_i16},
{0xfb9b7cd9a4a7443c_u64, 720_i16, 236_i16},
{0xbb764c4ca7a44410_u64, 747_i16, 244_i16},
{0x8bab8eefb6409c1a_u64, 774_i16, 252_i16},
{0xd01fef10a657842c_u64, 800_i16, 260_i16},
{0x9b10a4e5e9913129_u64, 827_i16, 268_i16},
{0xe7109bfba19c0c9d_u64, 853_i16, 276_i16},
{0xac2820d9623bf429_u64, 880_i16, 284_i16},
{0x80444b5e7aa7cf85_u64, 907_i16, 292_i16},
{0xbf21e44003acdd2d_u64, 933_i16, 300_i16},
{0x8e679c2f5e44ff8f_u64, 960_i16, 308_i16},
{0xd433179d9c8cb841_u64, 986_i16, 316_i16},
{0x9e19db92b4e31ba9_u64, 1013_i16, 324_i16},
{0xeb96bf6ebadf77d9_u64, 1039_i16, 332_i16},
{0xaf87023b9bf0ee6b_u64, 1066_i16, 340_i16},
].map { |t| Power.new t[0], t[1], t[2] }

Pow10Cache = {0, 1, 10, 100, 1000, 10000, 100000, 1000000, 10000000, 100000000, 1000000000}

def self.largest_pow10(n, n_bits)
# 1233/4096 is approximately 1/lg(10).
# We increment to skip over the first entry in the powers cache.
guess = ((n_bits + 1) * 1233 >> 12) + 1

# We don't have any guarantees that 2^number_bits <= number.<Paste>
guess -= 1 if n < Pow10Cache[guess]

return Pow10Cache[guess], guess
end

# Returns a cached power-of-ten with a binary exponent in the range
# around *exp* (boundaries included).
def self.get_cached_power_for_binary_exponent(exp) : {DiyFP, Int32}
min_exp = MIN_TARGET_EXP - (exp + DiyFP::SIGNIFICAND_SIZE)
max_exp = MAX_TARGET_EXP - (exp + DiyFP::SIGNIFICAND_SIZE)
k = ((min_exp + DiyFP::SIGNIFICAND_SIZE - 1) * D_1_LOG2_10).ceil
index = ((CACHED_POWER_OFFSET + k.to_i - 1) / CACHED_EXP_STEP) + 1
pow = PowCache[index]
_invariant min_exp <= pow.binary_exp
_invariant pow.binary_exp <= max_exp
return DiyFP.new(pow.significand, pow.binary_exp), pow.decimal_exp.to_i
end

private macro _invariant(exp, file = __FILE__, line = __LINE__)
{% if !flag?(:release) %}
unless {{exp}}
raise "Assertion Failed #{{{file}}}:#{{{line}}}"
end
{% end %}
end
end
150 changes: 150 additions & 0 deletions src/float_printer/diy_fp.cr
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,150 @@
# DiyFP is ported from the C++ "double-conversions" library.
# The following is their license:
# Copyright 2010 the V8 project authors. All rights reserved.
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are
# met:
#
# * Redistributions of source code must retain the above copyright
# notice, this list of conditions and the following disclaimer.
# * 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.
# * Neither the name of Google Inc. 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
# OWNER 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.

require "./ieee"

# This "Do It Yourself Floating Point" struct implements a floating-point number
# with a `UIht64` significand and an `Int32` exponent. Normalized DiyFP numbers will
# have the most significant bit of the significand set.
# Multiplication and Subtraction do not normalize their results.
# DiyFP is not designed to contain special Floats (NaN and Infinity).
struct FloatPrinter::DiyFP
SIGNIFICAND_SIZE = 64
# Also known as the significand
property frac : UInt64
# exponent
property exp : Int32

def initialize(@frac, @exp)
end

def initialize(@frac, exp : Int16)
@exp = exp.to_i32
end

def new(frac : Int32, exp)
new frac.to_u64, exp
end

# Returns a new `DiyFP` caculated as self - *other*.
#
# The exponents of both numbers must be the same and the frac of self must be
# greater than the other.
#
# This result is not normalized.
def -(other : DiyFP)
_invariant self.exp == other.exp && frac >= other.frac
self.class.new(frac - other.frac, exp)
end

MASK32 = 0xFFFFFFFF_u32

# Returns a new `DiyFP` caculated as self * *other*.
#
# Simply "emulates" a 128 bit multiplication.
# However: the resulting number only contains 64 bits. The least
# significant 64 bits are only used for rounding the most significant 64
# bits.
#
# This result is not normalized.
def *(other : DiyFP)
a = frac >> 32
b = frac & MASK32
c = other.frac >> 32
d = other.frac & MASK32
ac = a*c
bc = b*c
ad = a*d
bd = b*d
tmp = (bd >> 32) + (ad & MASK32) + (bc & MASK32)
# By adding 1U << 31 to tmp we round the final result.
# Halfway cases will be round up.
tmp += 1_u32 << 31
f = ac + (ad >> 32) + (bc >> 32) + (tmp >> 32)
e = exp + other.exp + 64

self.class.new(f, e)
end

def normalize
_invariant frac != 0
f = frac
e = exp

# This method is mainly called for normalizing boundaries. In general
# boundaries need to be shifted by 10 bits. We thus optimize for this case.
k10MSBits = 0xFFC0000000000000_u64
kUint64MSB = 0x8000000000000000_u64
while (f & k10MSBits) == 0
# puts " sig: #{f}"
# puts " exp: #{e}"
f <<= 10_u64
e -= 10
end
while (f & kUint64MSB) == 0
# puts " sig: #{f}"
# puts " exp: #{e}"
f <<= 1_u64
e -= 1
end
DiyFP.new(f, e)
end

def self.from_f(d : Float64 | Float32)
_invariant d > 0
frac, exp = IEEE.frac_and_exp(d)
new(frac, exp)
end

# Normalize such that the most signficiant bit of frac is set
def self.from_f_normalized(v : Float64 | Float32)
pre_normalized = from_f(v)
f = pre_normalized.frac
e = pre_normalized.exp

# could be a denormal
while (f & IEEE::HIDDEN_BIT_64) == 0
f <<= 1
e -= 1
end

# do the final shifts in one go
f <<= DiyFP::SIGNIFICAND_SIZE - IEEE::SIGNIFICAND_SIZE_64
e -= DiyFP::SIGNIFICAND_SIZE - IEEE::SIGNIFICAND_SIZE_64
DiyFP.new(f, e)
end

private macro _invariant(exp, file = __FILE__, line = __LINE__)
{% if !flag?(:release) %}
unless {{exp}}
raise "Assertion Failed #{{{file}}}:#{{{line}}}"
end
{% end %}
end
end
363 changes: 363 additions & 0 deletions src/float_printer/grisu3.cr
View file

Large diffs are not rendered by default.

202 changes: 202 additions & 0 deletions src/float_printer/ieee.cr
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,202 @@
# IEEE is ported from the C++ "double-conversions" library.
# The following is their license:
# Copyright 2012 the V8 project authors. All rights reserved.
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are
# met:
#
# * Redistributions of source code must retain the above copyright
# notice, this list of conditions and the following disclaimer.
# * 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.
# * Neither the name of Google Inc. 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
# OWNER 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.

module FloatPrinter::IEEE
extend self

EXPONENT_MASK_64 = 0x7FF0000000000000_u64
SIGNIFICAND_MASK_64 = 0x000FFFFFFFFFFFFF_u64
HIDDEN_BIT_64 = 0x0010000000000000_u64
PHYSICAL_SIGNIFICAND_SIZE_64 = 52 # Excludes the hidden bit
SIGNIFICAND_SIZE_64 = 53
EXPONENT_BIAS_64 = 0x3FF + PHYSICAL_SIGNIFICAND_SIZE_64
DENORMAL_EXPONENT_64 = -EXPONENT_BIAS_64 + 1
SIGN_MASK_64 = 0x8000000000000000_u64

EXPONENT_MASK_32 = 0x7F800000_u32
SIGNIFICAND_MASK_32 = 0x007FFFFF_u32
HIDDEN_BIT_32 = 0x00800000_u32
PHYSICAL_SIGNIFICAND_SIZE_32 = 23 # Excludes the hidden bit
SIGNIFICAND_SIZE_32 = 24
EXPONENT_BIAS_32 = 0x7F + PHYSICAL_SIGNIFICAND_SIZE_32
DENORMAL_EXPONENT_32 = -EXPONENT_BIAS_32 + 1
SIGN_MASK_32 = 0x80000000_u32

def to_uint(v : Float64)
pointerof(v).as(UInt64*).value
end

def to_uint(v : Float32)
pointerof(v).as(UInt32*).value
end

def sign(d64 : UInt64)
(d64 & SIGN_MASK_64) == 0 ? 1 : -1
end

def sign(d32 : UInt32)
(d32 & SIGN_MASK_32) == 0 ? 1 : -1
end

def special?(d64 : UInt64)
(d64 & EXPONENT_MASK_64) == EXPONENT_MASK_64
end

def special?(d32 : UInt32)
(d32 & EXPONENT_MASK_32) == EXPONENT_MASK_32
end

def inf?(d64 : UInt64)
special?(d64) && (d64 & SIGNIFICAND_MASK_64 == 0)
end

def inf?(d32 : UInt32)
special?(d32) && (d32 & SIGNIFICAND_MASK_32 == 0)
end

def nan?(d64 : UInt64)
special?(d64) && (d64 & SIGNIFICAND_MASK_64 != 0)
end

def nan?(d32 : UInt32)
special?(d32) && (d32 & SIGNIFICAND_MASK_32 != 0)
end

# Computes the two boundaries of *v*.
# The bigger boundary (m_plus) is normalized. The lower boundary has the same
# exponent as m_plus.
# Precondition: the value encoded by this Flaot must be greater than 0.
def normalized_boundaries(v : Float64)
_invariant v > 0
w = DiyFP.from_f(v)
m_plus = DiyFP.new((w.frac << 1) + 1, w.exp - 1).normalize

d64 = to_uint(v)

# The boundary is closer if the significand is of the form f == 2^p-1 then
# the lower boundary is closer.
# Think of v = 1000e10 and v- = 9999e9.
# Then the boundary (== (v - v-)/2) is not just at a distance of 1e9 but
# at a distance of 1e8.
# The only exception is for the smallest normal: the largest denormal is
# at the same distance as its successor.
# Note: denormals have the same exponent as the smallest normals.
physical_significand_is_zero = (d64 & SIGNIFICAND_MASK_64) == 0

lower_bound_closer = physical_significand_is_zero && (exponent(d64) != DENORMAL_EXPONENT_64)
calcualted_exp = exponent(d64)
calc_denormal = denormal?(d64)
f, e = if lower_bound_closer
{(w.frac << 2) - 1, w.exp - 2}
else
{(w.frac << 1) - 1, w.exp - 1}
end
m_minus = DiyFP.new(f << (e - m_plus.exp), m_plus.exp)
return {minus: m_minus, plus: m_plus}
end

def normalized_boundaries(v : Float32)
_invariant v > 0
w = DiyFP.from_f(v)
m_plus = DiyFP.new((w.frac << 1) + 1, w.exp - 1).normalize

d32 = to_uint(v)

physical_significand_is_zero = (d32 & SIGNIFICAND_MASK_32) == 0

lower_bound_closer = physical_significand_is_zero && (exponent(d32) != DENORMAL_EXPONENT_32)
calcualted_exp = exponent(d32)
calc_denormal = denormal?(d32)
f, e = if lower_bound_closer
{(w.frac << 2) - 1, w.exp - 2}
else
{(w.frac << 1) - 1, w.exp - 1}
end
m_minus = DiyFP.new(f << (e - m_plus.exp), m_plus.exp)
return {minus: m_minus, plus: m_plus}
end

def frac_and_exp(v : Float64)
d64 = to_uint(v)
_invariant (d64 & EXPONENT_MASK_64) != EXPONENT_MASK_64

if (d64 & EXPONENT_MASK_64) == 0 # denormal float
frac = d64 & SIGNIFICAND_MASK_64
exp = 1 - EXPONENT_BIAS_64
else
frac = (d64 & SIGNIFICAND_MASK_64) + HIDDEN_BIT_64
exp = (((d64 & EXPONENT_MASK_64) >> PHYSICAL_SIGNIFICAND_SIZE_64) - EXPONENT_BIAS_64).to_i
end

{frac, exp}
end

def frac_and_exp(v : Float32)
d32 = to_uint(v)
_invariant (d32 & EXPONENT_MASK_32) != EXPONENT_MASK_32

if (d32 & EXPONENT_MASK_32) == 0 # denormal float
frac = d32 & SIGNIFICAND_MASK_32
exp = 1 - EXPONENT_BIAS_32
else
frac = (d32 & SIGNIFICAND_MASK_32) + HIDDEN_BIT_32
exp = (((d32 & EXPONENT_MASK_32) >> PHYSICAL_SIGNIFICAND_SIZE_32) - EXPONENT_BIAS_32).to_i
end

{frac.to_u64, exp}
end

private def denormal?(d64 : UInt64) : Bool
(d64 & EXPONENT_MASK_64) == 0
end

private def denormal?(d32 : UInt32) : Bool
(d32 & EXPONENT_MASK_32) == 0
end

private def exponent(d64 : UInt64)
return DENORMAL_EXPONENT_64 if denormal?(d64)
baised_e = ((d64 & EXPONENT_MASK_64) >> PHYSICAL_SIGNIFICAND_SIZE_64).to_i
baised_e - EXPONENT_BIAS_64
end

private def exponent(d32 : UInt32)
return DENORMAL_EXPONENT_32 if denormal?(d32)
baised_e = ((d32 & EXPONENT_MASK_32) >> PHYSICAL_SIGNIFICAND_SIZE_32).to_i
baised_e - EXPONENT_BIAS_32
end

private macro _invariant(exp, file = __FILE__, line = __LINE__)
{% if !flag?(:release) %}
unless {{exp}}
raise "Assertion Failed #{{{file}}}:#{{{line}}}"
end
{% end %}
end
end