Merge pull request #73 from milankl/mk/restructure

Complete restructure

Project.toml

@@ -1,22 +1,24 @@
 name = "StochasticRounding"
 uuid = "3843c9a1-1f18-49ff-9d99-1b4c8a8e97ed"
-authors = ["Milan Kloewer"]
-version = "0.7"
+authors = ["Milan Kloewer and StochasticRounding.jl contributors"]
+version = "0.8.0"
 
 [deps]
 BFloat16s = "ab4f0b2a-ad5b-11e8-123f-65d77653426b"
 DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RandomNumbers = "e6cf234a-135c-5ec9-84dd-332b85af5143"
 
 [compat]
 BFloat16s = "0.1, 0.2, 0.3, 0.4"
-RandomNumbers = "^1.4"
 DoubleFloats = "1"
+Random = "1"
+RandomNumbers = "^1.4"
 julia = "1"
 
 [extras]
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test","DifferentialEquations"]
+test = ["Test", "DifferentialEquations"]

README.md

@@ -4,7 +4,7 @@
 
 Stochastic rounding for floating-point arithmetic.
 
-This package exports `Float64sr`, `Float32sr`,`Float16sr`, and `BFloat16sr`, three number formats that behave
+This package exports `Float64sr`, `Float32sr`,`Float16sr`, and `BFloat16sr`, four number formats that behave
 like their deterministic counterparts but with stochastic rounding that is proportional to the
 distance of the next representable numbers and therefore
 [exact in expectation](https://en.wikipedia.org/wiki/Rounding#Stochastic_rounding)
@@ -35,8 +35,25 @@ ask questions or suggest any changes or new features.
 
 ## Usage
 
-`Float64sr`, `Float32sr`, `Float16sr` and `BFloat16sr` are supposed to be drop-in replacements for their
-deterministically rounded counterparts. You can create data of those types as expected
+StochasticRounding.jl defines `stochastic_round` as a general interface to round stochastically to
+a floating-point number, e.g.
+```julia
+julia> stochastic_round(Float32,π)
+3.1415927f0
+
+julia> stochastic_round(Float32,π)
+3.1415925f0
+```
+Here we round π (which is first converted to Float64) to Float32 where it is not exactly representable,
+hence we have results that differ in the last bit. The chance here is 
+```julia
+julia> chance_roundup(Float32,π)
+0.6333222836256027
+```
+about 63% to obtain the first (which therefore would be round to nearest, the default rounding mode),
+and hence 37% to obtain the latter (round away from nearest). The `stochastic_round` function is wrapped
+into the floating-point formats `Float64sr`, `Float32sr`, `Float16sr` and `BFloat16sr` are supposed to be drop-in
+replacements for their deterministically rounded counterparts. You can create data of those types as expected
 (which is bitwise identical to the deterministic formats respectively) and the type
 will trigger stochastic rounding on every arithmetic operation.
 
@@ -50,7 +67,8 @@ BFloat16sr(0.33203125)
 ```
 As `1/3` is not exactly representable the rounding will be at 66.6% chance towards 0.33398438 
 and at 33.3% towards 0.33203125 such that in expectation the result is 0.33333... and therefore exact. 
-You can use `BFloat16_chance_roundup(x::Float32)` to get the chance that `x` will be rounded up.
+
+## Example
 
 Solving a linear equation system with LU decomposition and stochastic rounding:
 ```julia

src/StochasticRounding.jl

@@ -1,17 +1,7 @@
 module StochasticRounding
 
-    export BFloat16sr,BFloat16_stochastic_round,        # BFloat16 + SR
-        BFloat16_chance_roundup,NaNB16sr,InfB16sr,
-        Float16sr,Float16_stochastic_round,             # Float16 + SR
-        Float16_chance_roundup,NaN16sr,Inf16sr,
-        Float32sr,Float32_stochastic_round,             # Float32 + SR
-        Float32_chance_roundup,NaN32sr,Inf32sr,
-        Float64sr,Float64_stochastic_round,             # Float64 + SR
-        Float64_chance_roundup,
-        NaNsr,Infsr
-
     # use BFloat16 from BFloat16s.jl
-    import BFloat16s: BFloat16
+    import BFloat16s: BFloat16s, BFloat16
 
     # faster random number generator
     import RandomNumbers.Xorshifts.Xoroshiro128Plus
@@ -30,12 +20,12 @@ module StochasticRounding
         return nothing
     end
 
-    include("bfloat16sr.jl")
-    include("float16sr.jl")
-    include("float32sr.jl")
-    include("float64sr.jl")
+    # define abstract type
+    export AbstractStochasticFloat
+    abstract type AbstractStochasticFloat <: AbstractFloat end
 
-    include("general.jl")
-    include("promotion.jl")
-    include("conversions.jl")
+    include("types.jl")             # define concrete types
+    include("promotions.jl")        # their promotions
+    include("general.jl")           # and general functions
+    include("bfloat16.jl")          # some functions that aren't yet in BFloat16s.jl
 end

src/bfloat16.jl

@@ -0,0 +1,58 @@
+# from BFloat16s.jl#main (22/12/23) to be removed with next release of BFloat16s.jl
+function Base.nextfloat(f::BFloat16, d::Integer)
+
+    F = typeof(f)
+    fumax = reinterpret(Base.uinttype(F), F(Inf))
+    U = typeof(fumax)
+
+    isnan(f) && return f
+    fi = reinterpret(Int16, f)
+    fneg = fi < 0
+    fu = unsigned(fi & typemax(fi))
+
+    dneg = d < 0
+    da = Base.uabs(d)
+    if da > typemax(U)
+        fneg = dneg
+        fu = fumax
+    else
+        du = da % U
+        if fneg ⊻ dneg
+            if du > fu
+                fu = min(fumax, du - fu)
+                fneg = !fneg
+            else
+                fu = fu - du
+            end
+        else
+            if fumax - fu < du
+                fu = fumax
+            else
+                fu = fu + du
+            end
+        end
+    end
+    if fneg
+        fu |= Base.sign_mask(F)
+    end
+    reinterpret(F, fu)
+end
+
+function Base.decompose(x::BFloat16)::NTuple{3,Int}
+    isnan(x) && return 0, 0, 0
+    isinf(x) && return ifelse(x < 0, -1, 1), 0, 0
+    n = reinterpret(UInt16, x)
+    s = (n & 0x007f) % Int16
+    e = ((n & 0x7f80) >> 7) % Int
+    s |= Int16(e != 0) << 7
+    d = ifelse(signbit(x), -1, 1)
+    s, e - 134 + (e == 0), d
+end
+
+BFloat16s.BFloat16(x::BigFloat) = BFloat16(Float32(x))
+Base.BigFloat(x::BFloat16) = BigFloat(Float32(x))
+
+Base.maxintfloat(::Type{BFloat16}) = reinterpret(BFloat16,0x4380)
+
+Base.round(x::BFloat16, r::RoundingMode{:ToZero}) = BFloat16(trunc(Float32(x)))
+Base.round(x::BFloat16, r::RoundingMode{:Down}) = BFloat16(floor(Float32(x)))