fix tests and implement cbrt and fourthroot

Project.toml

@@ -13,7 +13,10 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LogarithmicNumbersForwardDiffExt = "ForwardDiff"
 
 [compat]
+Aqua = "0.8"
 ForwardDiff = "0.10"
+Random = "1"
+Test = "1"
 julia = "1.6"
 
 [extras]

README.md

@@ -44,7 +44,7 @@ julia> log(x)
 999.8545865421312
 ```
 
-## Documentation
+## Features
 
 ### Exported types
 * `ULogarithmic{T}` represents a non-negative real number by its logarithm of type `T`.
@@ -59,7 +59,7 @@ julia> log(x)
   you already know the logarithm `logx` of your number `x`.
 
 ### Functions in Base
-* **Arithmetic:** `+`, `-`, `*`, `/`, `^`, `inv`, `prod`, `sum`.
+* **Arithmetic:** `+`, `-`, `*`, `/`, `^`, `inv`, `prod`, `sum`, `sqrt`, `cbrt`, `fourthroot`.
 * **Ordering:** `==`, `<`, `≤`, `cmp`, `isless`, `isequal`, `sign`, `signbit`, `abs`.
 * **Logarithm:** `log`, `log2`, `log10`, `log1p`. These are returned as the base (non-logarithmic) type.
 * **Conversion:** `float`, `unsigned`, `signed`, `widen`, `big`. These also operate on types.
@@ -84,7 +84,7 @@ possibilities for `func`:
 - [SpecialFunctions.jl](https://github.com/JuliaMath/SpecialFunctions.jl):
   `gamma`, `factorial`, `beta`, `erfc`, `erfcx`.
 
-## Autodiff
+### Autodiff
 
 On Julia >= 1.9, if you load [ForwardDiff.jl](https://github.com/JuliaDiff/ForwardDiff.jl), you should be allowed to compute
 

docs/src/index.md

@@ -1,35 +1,89 @@
 # LogarithmicNumbers.jl
 
-A [logarithmic number system](https://en.wikipedia.org/wiki/Logarithmic_number_system) for Julia.
+A [**logarithmic number system**](https://en.wikipedia.org/wiki/Logarithmic_number_system)
+for Julia.
 
-Provides the [`ULogarithmic`](@ref) and [`Logarithmic`](@ref) subtypes of `Real` for representing real numbers in log-space.
+Provides the signed `Logarithmic` and unsigned `ULogarithmic` types for representing real
+numbers on a logarithmic scale.
 
-## Features
+This is useful when numbers are too big or small to fit accurately into a `Float64` and you
+only really care about magnitude.
 
-* `ULogarithmic(x)` and `Logarithmic(x)` represent the number `x`.
-* `exp(ULogarithmic, x)` and `exp(Logarithmic, x)` represent `exp(x)` (and `x` can be huge).
-* Arithmetic: plus, minus, times, divide, power, `sqrt`, `inv`, `log`, `prod`, `sum`.
-* Comparisons: equality and ordering.
-* Random: `rand(ULogarithmic)` and `rand(Logarithmic)` produces a random number in the unit interval.
-* Other functions: `float`, `big`, `unsigned` (converts `ULogarithmic` to `Logarithmic`), `signed` (vice versa), `widen`, `typemin`, `typemax`, `zero`, `one`, `iszero`, `isone`, `isinf`, `isfinite`, `isnan`, `sign`, `signbit`, `abs`, `nextfloat`, `prevfloat`, `write`, `read`.
+For example, it can be useful to represent probabilities in this form, and you don't need to
+worry about getting zero when multiplying many of them together.
+
+## Installation
+
+```
+pkg> add LogarithmicNumbers
+```
 
 ## Example
+
 ```julia
 julia> using LogarithmicNumbers
 
 julia> ULogarithmic(2.7)
 exp(0.9932517730102834)
 
-julia> float(ans) # convert to Float64
+julia> float(ans)
 2.7
 
 julia> x = exp(ULogarithmic, 1000) - exp(ULogarithmic, 998)
 exp(999.8545865421312)
 
-julia> float(x) # overflows Float64
+julia> float(x) # overflows
 Inf
 
 julia> log(x)
 999.8545865421312
 ```
 
+## Documentation
+
+### Exported types
+* `ULogarithmic{T}` represents a non-negative real number by its logarithm of type `T`.
+* `Logarithmic{T}` represents a real number by its absolute value as a `ULogarithmic{T}` and
+  a sign bit.
+* `LogFloat64` is an alias for `Logarithmic{Float64}`. There are also `ULogFloat16`,
+  `ULogFloat32`, `ULogFloat64`, `LogFloat16`, and `LogFloat32`.
+
+### Constructors
+* `ULogarithmic(x)` and `Logarithmic(x)` represent the number `x`.
+* `exp(ULogarithmic, logx)` represents `exp(logx)`, and `logx` can be huge. Use this when
+  you already know the logarithm `logx` of your number `x`.
+
+### Functions in Base
+* **Arithmetic:** `+`, `-`, `*`, `/`, `^`, `inv`, `prod`, `sum`, `sqrt`, `cbrt`, `fourthroot`.
+* **Ordering:** `==`, `<`, `≤`, `cmp`, `isless`, `isequal`, `sign`, `signbit`, `abs`.
+* **Logarithm:** `log`, `log2`, `log10`, `log1p`. These are returned as the base (non-logarithmic) type.
+* **Conversion:** `float`, `unsigned`, `signed`, `widen`, `big`. These also operate on types.
+* **Special values:** `zero`, `one`, `typemin`, `typemax`.
+* **Predicates:** `iszero`, `isone`, `isinf`, `isfinite`, `isnan`.
+* **IO:** `show`, `write`, `read`.
+* **Random:** `rand(ULogarithmic)` is a random number in the unit interval.
+* **Misc:** `nextfloat`, `prevfloat`, `hash`.
+* **Note:** Any functions not mentioned here might be inaccurate.
+
+### Interoperability with other packages
+
+It is natural to use this package in conjunction with other packages which return
+logarithms. The general pattern is that you can use `exp(ULogarithmic, logfunc(args...))`
+instead of `func(args...)` to get the answer as a logarithmic number. Here are some
+possibilities for `func`:
+
+- [StatsFuns.jl](https://github.com/JuliaStats/StatsFuns.jl):
+  `normpdf`, `normcdf`, `normccdf`, plus equivalents for other distributions.
+- [Distributions.jl](https://github.com/JuliaStats/Distributions.jl):
+  `pdf`, `cdf`, `ccdf`.
+- [SpecialFunctions.jl](https://github.com/JuliaMath/SpecialFunctions.jl):
+  `gamma`, `factorial`, `beta`, `erfc`, `erfcx`.
+
+### Autodiff
+
+On Julia >= 1.9, if you load [ForwardDiff.jl](https://github.com/JuliaDiff/ForwardDiff.jl), you should be allowed to compute
+
+- derivatives of functions involving `exp(Logarithmic, x)`
+- derivatives of functions evaluated at `Logarithmic(x)`
+
+This functionality is experimental, please report any bug or unexpected behavior.

src/LogarithmicNumbers.jl

@@ -605,8 +605,39 @@ function Base.exp(x::Logarithmic)
     x.signbit ? inv(exp(x.abs)) : exp(x.abs)
 end
 
-function Base.sqrt(x::AnyLogarithmic)
-    uexp(x.log/2)
+function Base.sqrt(x::ULogarithmic)
+    uexp(x.log / oftype(x.log, 2))
+end
+
+function Base.sqrt(x::Logarithmic)
+    if !x.signbit || iszero(x)
+        Logarithmic(sqrt(x.abs))
+    else
+        throw(DomainError(x))
+    end
+end
+
+function Base.cbrt(x::ULogarithmic)
+    uexp(x.log / oftype(x.log, 3))
+end
+
+function Base.cbrt(x::Logarithmic)
+    Logarithmic(cbrt(x.abs), x.signbit)
+end
+
+if hasproperty(Base, :fourthroot)
+    # fourthroot was introduced in julia 1.10
+    function Base.fourthroot(x::ULogarithmic)
+        uexp(x.log / oftype(x.log, 4))
+    end
+
+    function Base.fourthroot(x::Logarithmic)
+        if !x.signbit || iszero(x)
+            Logarithmic(fourthroot(x.abs))
+        else
+            throw(DomainError(x))
+        end
+    end
 end
 
 ### Hash

test/runtests.jl

@@ -15,6 +15,8 @@ _mul(args...) = @inferred *(args...)
 _div(args...) = @inferred /(args...)
 _pow(args...) = @inferred ^(args...)
 _sqrt(args...) = @inferred sqrt(args...)
+_cbrt(args...) = @inferred cbrt(args...)
+_fourthroot(args...) = @inferred fourthroot(args...)
 _float(args...) = @inferred float(args...)
 _inv(args...) = @inferred inv(args...)
 _prod(args...) = @inferred prod(args...)
@@ -39,7 +41,7 @@ atypes2 = (ULogarithmic, ULogFloat32, Logarithmic, LogFloat32)
 @testset verbose=true "LogarithmicNumbers" begin
 
     @testset verbose=true "Aqua" begin
-        Aqua.test_all(LogarithmicNumbers, project_toml_formatting=(VERSION >= v"1.7"))
+        Aqua.test_all(LogarithmicNumbers)
     end
 
     @testset "types" begin
@@ -461,8 +463,32 @@ atypes2 = (ULogarithmic, ULogFloat32, Logarithmic, LogFloat32)
 
         @testset "sqrt" begin
             for A in atypes, x in vals
-                x < 0 && continue
-                @test _approx(_float(_sqrt(A(x))), sqrt(float(x)))
+                if x < 0
+                    A <: ULogarithmic && continue
+                    @test_throws DomainError _sqrt(A(x))
+                else
+                    @test _approx(_float(_sqrt(A(x))), sqrt(float(x)))
+                end
+            end
+        end
+
+        @testset "cbrt" begin
+            for A in atypes, x in vals
+                x < 0 && A <: ULogarithmic && continue
+                @test _approx(_float(_cbrt(A(x))), cbrt(float(x)))
+            end
+        end
+
+        @testset "fourthroot" begin
+            if hasproperty(Base, :fourthroot)
+                for A in atypes, x in vals
+                    if x < 0
+                        A <: ULogarithmic && continue
+                        @test_throws DomainError _fourthroot(A(x))
+                    else
+                        @test _approx(_float(_fourthroot(A(x))), fourthroot(float(x)))
+                    end
+                end
             end
         end