From 8fef0a7d52605d9aa0a4de5833e6648e9646e2bf Mon Sep 17 00:00:00 2001 From: KlausC Date: Tue, 31 Oct 2023 15:46:49 +0100 Subject: [PATCH 01/15] run for ltr, release , nightly on all os --- .github/workflows/CI.yml | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index acf42c98..e1af050f 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -13,12 +13,13 @@ jobs: fail-fast: false matrix: version: - - '1.8.4' -# - 'nightly' + - '1.6' + - '1' + - 'nightly' os: - ubuntu-latest - # - macOS-latest - # - windows-latest + - macOS-latest + - windows-latest arch: - x64 steps: From 1ce8140633025702fdebd5171e21007bc03dd981 Mon Sep 17 00:00:00 2001 From: KlausC Date: Wed, 1 Nov 2023 16:12:57 +0100 Subject: [PATCH 02/15] Pack of tests and trivial fixes --- .gitignore | 1 + README.md | 86 ++++++++++++++++++-------------------- src/float/prearith.jl | 32 ++------------ src/libarb/ArbComplex.jl | 21 ---------- src/libarb/ArbFloat.jl | 23 +++++----- src/values/constructors.jl | 18 ++++---- test/complex.jl | 73 ++++++++++++++++++++++++++++++-- 7 files changed, 133 insertions(+), 121 deletions(-) diff --git a/.gitignore b/.gitignore index f7635ba5..305731c4 100644 --- a/.gitignore +++ b/.gitignore @@ -2,6 +2,7 @@ *.jl.*.cov *.jl.mem *.jl.*.mem +lcov.info Manifest.toml settings.json deps/deps.jl diff --git a/README.md b/README.md index 2adeb516..24648226 100644 --- a/README.md +++ b/README.md @@ -1,21 +1,19 @@ -# + # ArbNumerics.jl +**Copyright © 2015-2023 by Jeffrey Sarnoff.** -#### Copyright © 2015-2023 by Jeffrey Sarnoff. -#### This work is released under The MIT License. +**This work is released under The MIT License.** For multiprecision numerical computing using values with 25..2,500 digits. With arithmetic and higher level mathematics, this package offers you the best balance of performance and accuracy. This package uses the [Arb C Library](http://arblib.org/index.html), and adapts some C library interface work from [Nemo](https://github.com/wbhart/Nemo.jl) (see [_below_](https://github.com/JeffreySarnoff/ArbNumerics.jl/blob/master/README.md#acknowledgements)). Here is a presentation by the designer and architect of the [Arb C library](https://fredrikj.net/math/oxford2019.pdf). - ----- [![Travis Build Status](https://travis-ci.org/JeffreySarnoff/ArbNumerics.jl.svg?branch=master)](https://travis-ci.org/JeffreySarnoff/ArbNumerics.jl)   [![Docs](https://img.shields.io/badge/docs-stable-blue.svg)](http://jeffreysarnoff.github.io/ArbNumerics.jl/stable/)   [![Docs](https://img.shields.io/badge/docs-dev-blue.svg)](http://jeffreysarnoff.github.io/ArbNumerics.jl/dev/) - ----- +----- ## Introduction @@ -23,7 +21,6 @@ ArbNumerics exports three types: `ArbFloat`, `ArbReal`, `ArbComplex`. `ArbFloat While the bounds of an `ArbReal` or `ArbComplex` are available, the default is to show these values as digit sequences which almost assuredly are accurate, in a round to nearest sense, to the precision displayed. Math with `ArbFloat` does not provide the assurance one gets using `ArbReal`, as an `ArbFloat` is a point value. While some effort has been taken to provide you with more reliable results from math with `ArbFloat` values than would be the case using the underlying library itself, `ArbReal` or `ArbComplex` are suggested for work that is important to you. `ArbFloat` is appropriate when exactness is not required during development, or with applications that are approximating something at increasing precisions. - ## Installation ```julia @@ -38,7 +35,7 @@ When updating ArbNumerics, do `pkg> gc` to prevent accruing a great deal of unus ## StartUp `using ArbNumerics` -or, if you installed Readables, +or, if you installed Readables, `using ArbNumerics, Readables` ## Precision @@ -49,23 +46,23 @@ Otherwise, some extra bits are used to assist with printing values rounded to th You can set the internal working precision (which is the same as the displayed precision with `setextrabits(0)`) to a given number of bits or a given number of decimal digits: -`setworkingprecision(ArbFloat, bits=250)`, `setworkingprecision(ArbReal, digits=100)` +`setworkingprecision(ArbFloat, bits=250)`, `setworkingprecision(ArbReal, digits=100)` The type can be any of `ArbFloat`, `ArbReal`, `ArbComplex`. All types share the same precision so interconversion makes sense. You can set the external displayed precision (which is the the same as the internal precision with `setextrabits(0)`) to a given number of bits or a given number of decimal digits: -`setprecision(ArbFloat, bits=250)`, `setworkingprecision(ArbReal, digits=100)` +`setprecision(ArbFloat, bits=250)`, `setworkingprecision(ArbReal, digits=100)` The type can be any of `ArbFloat`, `ArbReal`, `ArbComplex`. All types share the same precision so interconversion makes sense. - -## Using ArbFloat +## Using ArbFloat Reading the sections that follow gives you a good platform from which to develop. - consider `using ArbNumerics, Readables` + ```julia julia> ArbFloat(pi, digits=30, base=10) 3.14159265358979323846264338328 @@ -82,7 +79,8 @@ Initially, the default precision is set to 106 bits. All ArbNumeric types use t The precision in use may be set globally, as with BigFloats, or it may be given with the constructor. For most purposes, you should work with a type at one, two, or three precisions. It is helps clarity to convert precisions explicitly, however, it is not necessary. -#### Constructors using the default precision +### Constructors using the default precision + ```julia julia> a = ArbFloat(3) 3.0000000000000000000000000000000 @@ -95,42 +93,39 @@ julia> c = one(ArbComplex) ``` #### Constructors using a specified precision + ```julia julia> BITS = 53; -julia> a = sqrt(ArbFloat(2, BITS)) +julia> a = sqrt(ArbFloat(2, bits=BITS)) 1.414213562373095 -julia> b = ArbReal(pi, BITS) +julia> b = ArbReal(pi, bits=BITS) 3.141592653589793 -julia> c = ArbComplex(a, b, BITS) +julia> c = ArbComplex(a, b, bits=BITS) 1.414213562373095 + 3.141592653589793im ``` + ```julia julia> DIGITS = 78; -julia> ArbFloat(pi, bits4digits(DIGITS)) +julia> ArbFloat(pi, digits=DIGITS) 3.14159265358979323846264338327950288419716939937510582097494459230781640628621 julia> DIGITS == length(string(ans)) - 1 # (-1 for the decimal point) true ``` + ### changing precision ```julia -julia> a = ArbFloat(2, 25) -2.000000 -julia> a = ArbFloat(a, 50) -2.00000000000000 - -julia> precision = 25 -julia> a = ArbFloat(2, precision) +julia> a = ArbFloat(2, bits=25) 2.000000 -julia> precision = 50 -julia> a = ArbFloat(a, precision) +julia> a = ArbFloat(a, bits=50) 2.00000000000000 ``` + ### interconversion ```julia @@ -153,11 +148,11 @@ julia> Float16(c) Float16(1.414) ``` ----- +----- Consider using ArbReals instead of ArbFloats if you want your results to be rock solid. That way you can examine the enclosures for your results with `radius(value)` or `bounds(value)`. This is strongly suggested when working with precisions that you are increasing dynamically. ----- +----- ### Math @@ -165,7 +160,7 @@ Consider using ArbReals instead of ArbFloats if you want your results to be rock - `+`,`-`, `*`, `/` - `square`, `cube`, `sqrt`, `cbrt`, `hypot` -- `pow(x,i)`, `root(x,i)` _where i is an integer > 0_ +- `pow(x, i)`, `root(x, i)` _where i is an integer > 0_ - `factorial`, `doublefactorial`, `risingfactorial` - `binomial` @@ -208,6 +203,7 @@ Consider using ArbReals instead of ArbFloats if you want your results to be rock - `elliptic_zeta`, `elliptic_sigma` ##### elliptic integrals of squared modulus + - `elliptic_e2`, `elliptic_k2` - `elliptic_p2`, `elliptic_pi2` - `elliptic_zeta2`, `elliptic_sigma2` @@ -218,10 +214,10 @@ Consider using ArbReals instead of ArbFloats if you want your results to be rock - `weierstrass_zeta`, `weierstrass_sigma` #### hypergeometric functions - + - `hypgeom0f1`, `hypgeom1f1`, `hypgeom1f2` - `hypgeom0f1reg`, `hypgeom1f1reg`, `hypgeom1f2reg` (regularized) - + #### other special functions - `ei`, `si`, `ci` @@ -241,13 +237,13 @@ Consider using ArbReals instead of ArbFloats if you want your results to be rock - `dft`, `inverse_dft` - see docs for use -## Intervals +### Intervals #### parts -- midpoint, radius -- upperbound, lowerbound, bounds -- upperbound_abs, lowerbound_abs, bounds_abs +- `midpoint`, `radius` +- `upperbound`, `lowerbound`, `bounds` +- `upperbound_abs`, `lowerbound_abs`, `bounds_abs` #### construction @@ -267,26 +263,26 @@ The radii are kept using an Arb C library internal structure that has a 30 bit u When constructing intervals , you should scale the radius to be as small as possible while preserving enclosure. ----- +----- ### a caution for BigFloat ```julia -julia> p=64;setprecision(BigFloat,p); +julia> p = 64; setprecision(BigFloat, p); -julia> ArbFloat(pi,p+8) +julia> ArbFloat(pi, bits=p+8) 3.14159265358979323846 -julia> ArbFloat(pi,p),BigFloat(pi) +julia> ArbFloat(pi, bits=p), BigFloat(pi) (3.141592653589793238, 3.14159265358979323851) -julia> [ArbFloat(pi,p), BigFloat(pi)] -2-element Array{ArbFloat{88},1}: +julia> [ArbFloat(pi, bits=p), BigFloat(pi)] +2-element Vector{ArbFloat{88}}: 3.141592653589793238 3.141592653589793239 ``` ----- +----- ## The Arb C library @@ -300,18 +296,18 @@ julia> [ArbFloat(pi,p), BigFloat(pi)] - The code is thread-safe, portable, and extensively tested. The library outperforms others. - ## Acknowledgements This work develops parts the Arb C library within Julia. It is entirely dependent on Arb by Fredrik Johansson and would not exist without the good work of William Hart, Tommy Hofmann and the Nemo.jl team. The libraries for `Arb` and `Flint`, and build file are theirs, used with permission. ----- +----- ## Alternatives For a numeric types like `Float64` and `ComplexF64` with about twice their precision, [Quadmath.jl](https://github.com/JuliaMath/Quadmath.jl) exports `Float128` and `ComplexF128`. For almost as much precision with better performance, [DoubleFloats.jl](https://github.com/JuliaMath/DoubleFloats.jl) exports `Double64` and `ComplexDF64`. ValidatedNumerics.jl and other packages available at [JuliaIntervals](https://github.com/JuliaIntervals) provide an alternative approach to developing correctly contained results. Those packages are very good and worthwhile when you do not require multiprecision numerics. ----- +----- + ## notes - To propose internal changes, please use pull requests. diff --git a/src/float/prearith.jl b/src/float/prearith.jl index 8531ba99..3ef0405a 100644 --- a/src/float/prearith.jl +++ b/src/float/prearith.jl @@ -91,36 +91,10 @@ abs2(x::ArbFloat{P}) where {P} = square( abs(x) ) abs2(x::ArbReal{P}) where {P} = square( abs(x) ) abs2(x::ArbComplex{P}) where {P} = square( abs(x) ) +# needed for GenericLinearAlgebra - -flipsign(x::ArbFloat{P}, y::U) where {P, U<:Unsigned} = +x -copysign(x::ArbFloat{P}, y::U) where {P, U<:Unsigned} = +x - -flipsign(x::ArbFloat{P}, y::S) where {P, S<:Signed} = signbit(y) ? -x : x -copysign(x::ArbFloat{P}, y::S) where {P, S<:Signed} = signbit(y) ? -abs(x) : abs(x) - -flipsign(x::ArbFloat{P}, y::F) where {P, F<:AbstractFloat} = signbit(y) ? -x : x -copysign(x::ArbFloat{P}, y::F) where {P, F<:AbstractFloat} = signbit(y) ? -abs(x) : abs(x) - -flipsign(x::ArbReal{P}, y::U) where {P, U<:Unsigned} = +x -copysign(x::ArbReal{P}, y::U) where {P, U<:Unsigned} = +x - -flipsign(x::ArbReal{P}, y::S) where {P, S<:Signed} = signbit(y) ? -x : x -copysign(x::ArbReal{P}, y::S) where {P, S<:Signed} = signbit(y) ? -abs(x) : abs(x) - -flipsign(x::ArbReal{P}, y::F) where {P, F<:AbstractFloat} = signbit(y) ? -x : x -copysign(x::ArbReal{P}, y::F) where {P, F<:AbstractFloat} = signbit(y) ? -abs(x) : abs(x) - -flipsign(x::ArbComplex{P}, y::U) where {P, U<:Unsigned} = +x -copysign(x::ArbComplex{P}, y::U) where {P, U<:Unsigned} = +x - -flipsign(x::ArbComplex{P}, y::S) where {P, S<:Signed} = signbit(y) ? -x : x -copysign(x::ArbComplex{P}, y::S) where {P, S<:Signed} = signbit(y) ? -abs(x) : abs(x) - -flipsign(x::ArbComplex{P}, y::F) where {P, F<:AbstractFloat} = signbit(y) ? -x : x -copysign(x::ArbComplex{P}, y::F) where {P, F<:AbstractFloat} = signbit(y) ? (signbit(x.re) ? x : -x) : x - - +flipsign(x::ArbNumber, y::Real) = signbit(y) ? -x : x +copysign(x::ArbNumber, y::Real) = (signbit(y) == signbit(x)) ? x : -x inv(x::ArbFloat{P}) where {P} = ArbFloat{P}( inv(ArbReal{P}(x)) ) diff --git a/src/libarb/ArbComplex.jl b/src/libarb/ArbComplex.jl index b83b4feb..b5850ca8 100644 --- a/src/libarb/ArbComplex.jl +++ b/src/libarb/ArbComplex.jl @@ -368,27 +368,6 @@ function Base.angle(x::ArbComplex{P}) where {P} !(signbit(a) || signbit(T(pi) - a)) ? a : (signbit(a) ? zero(T) : T(pi)) end -# needed for GenericLinearAlgebra - -flipsign(x::ArbComplex{P}, y::T) where {P, T<:Base.IEEEFloat} = - signbit(y) ? -x : x -flipsign(x::ArbComplex{P}, y::T) where {P, T<:Real} = - signbit(y) ? -x : x -flipsign(x::ArbComplex{P}, y::T) where {P, T<:ArbFloat} = - signbit(y) ? -x : x -flipsign(x::ArbComplex{P}, y::T) where {P, T<:ArbReal} = - signbit(y) ? -x : x - -copysign(x::ArbComplex{P}, y::T) where {P, T<:Base.IEEEFloat} = - signbit(y) ? (signbit(x) ? x : -x) : (signbit(x) ? -x : x) -copysign(x::ArbComplex{P}, y::T) where {P, T<:Real} = - signbit(y) ? (signbit(x) ? x : -x) : (signbit(x) ? -x : x) -copysign(x::ArbComplex{P}, y::T) where {P, T<:ArbFloat} = - signbit(y) ? (signbit(x) ? x : -x) : (signbit(x) ? -x : x) -copysign(x::ArbComplex{P}, y::T) where {P, T<:ArbReal} = - signbit(y) ? (signbit(x) ? x : -x) : (signbit(x) ? -x : x) - - # a type specific hash function helps the type to 'just work' const hash_arbcomplex_lo = (UInt === UInt64) ? 0x76143ad985246e79 : 0x5b6a64dc const hash_0_arbcomplex_lo = hash(zero(UInt), hash_arbcomplex_lo) diff --git a/src/libarb/ArbFloat.jl b/src/libarb/ArbFloat.jl index 325bf3cf..67580ac5 100644 --- a/src/libarb/ArbFloat.jl +++ b/src/libarb/ArbFloat.jl @@ -126,17 +126,14 @@ end Int32(x::ArbFloat{P}, roundingmode::RoundingMode) where {P} = Int32(Int64(x), roundingmode) Int16(x::ArbFloat{P}, roundingmode::RoundingMode) where {P} = Int16(Int64(x), roundingmode) -BigFloat(x::ArbFloat{P}) where {P} = BigFloat(x, RoundNearest) -function BigFloat(x::ArbFloat{P}, roundingmode::RoundingMode) where {P} - rounding = match_rounding_mode(roundingmode) - z = BigFloat(0, workingprecision(x)) - roundingdir = ccall(@libarb(arf_get_mpfr), Cint, (Ref{BigFloat}, Ref{ArbFloat}, Cint), z, x, rounding) - return z -end -BigFloat(x::ArbFloat{P}, bitprecision::Int) where {P} = BigFloat(x, bitprecision, RoundNearest) -function BigFloat(x::ArbFloat{P}, bitprecision::Int, roundingmode::RoundingMode) where {P} +""" + BigFloat(::ArbFloat; [precision=workingprecision(x), roundingmode=RoundNearest]) + +Construct a `BigFloat`from an `ArbFloat`. +""" +function BigFloat(x::ArbFloat{P}; precision::Int=workingprecision(x), roundingmode::RoundingMode=RoundNearest) where {P} rounding = match_rounding_mode(roundingmode) - z = BigFloat(0, bitprecision) + z = BigFloat(0; precision) roundingdir = ccall(@libarb(arf_get_mpfr), Cint, (Ref{BigFloat}, Ref{ArbFloat}, Cint), z, x, rounding) return z end @@ -192,7 +189,7 @@ function divrem(x::ArbFloat{P}, y::ArbFloat{P}) where {P} end fld(x::ArbFloat{P}, y::ArbFloat{P}) where {P} = - floor(x / y) + floor(x / y) mod(x::ArbFloat{P}, y::ArbFloat{P}) where {P} = x - (fld(x,y) * y) @@ -204,9 +201,9 @@ function fldmod(x::ArbFloat{P}, y::ArbFloat{P}) where {P} end cld(x::ArbFloat{P}, y::ArbFloat{P}) where {P} = - ceil(x / y) + ceil(x / y) + - # a type specific hash function helps the type to 'just work' const hash_arbfloat_lo = (UInt === UInt64) ? 0x37e642589da3416a : 0x5d46a6b4 const hash_0_arbfloat_lo = hash(zero(UInt), hash_arbfloat_lo) diff --git a/src/values/constructors.jl b/src/values/constructors.jl index 0a30d1fa..9fb62620 100644 --- a/src/values/constructors.jl +++ b/src/values/constructors.jl @@ -103,7 +103,7 @@ Complex(x::ArbComplex{P}) where {P} = Complex{Float64}(x, RoundNearest) for I in (:Int8, :Int16, :Int32, :Int64, :Int128) @eval begin - + function $I(x::ArbFloat{P}) where {P} !isinteger(x) && throw(InexactError("$(x) is not an integer")) bi = BigInt(BigFloat(x)) @@ -111,20 +111,20 @@ for I in (:Int8, :Int16, :Int32, :Int64, :Int128) return $I(bi) end - function $I(x::ArbReal{P}) where {P} + function $I(x::ArbReal{P}) where {P} (!isexact(x) | !isinteger(x)) && throw(InexactError("$(x) is not an integer")) bi = BigInt(BigFloat(x)) !(typemin($I) <= bi <= typemax($I)) && throw(InexactError("$(x)")) return $I(bi) end - function $I(x::ArbComplex{P}) where {P} + function $I(x::ArbComplex{P}) where {P} (!isexact(x) | !isinteger(x) | !iszero(imag(x))) && throw(InexactError("$(x) is not an integer")) bi = BigInt(BigFloat(x)) !(typemin($I) <= bi <= typemax($I)) && throw(InexactError("$(x)")) return $I(bi) - end - + end + end end @@ -140,20 +140,20 @@ ArbComplex(x::T, y::T) where {S, T<:Rational{S}} = ArbComplex(ArbReal(x), ArbRea function ArbReal{P}(x::Irrational{S}) where {P,S} mid = ArbFloat{P}(x) rad = ulp(mid) - return setball(mid, rad) + return setball(mid, rad) end function ArbReal(x::Irrational{S}) where {S} P = workingprecision(ArbReal) mid = ArbFloat{P}(x) rad = ulp(mid) - return setball(mid, rad) + return setball(mid, rad) end ArbComplex(x::Irrational{S}) where {S} = ArbComplex(ArbReal(x), ArbReal(0)) ArbComplex{P}(x::Irrational{S}) where {P,S} = ArbComplex{P}(ArbReal{P}(x), ArbReal{P}(0)) -ArbComplex(x::Irrational{S}, y::Real) where {S} = ArbComplex(ArbReal(x), ArbReal(0)) -ArbComplex{P}(x::Irrational{S}, y::Real) where {P,S} = ArbComplex{P}(ArbReal{P}(x), ArbReal{P}(0)) +ArbComplex(x::Irrational{S}, y::Real) where {S} = ArbComplex(ArbReal(x), ArbReal(y)) +ArbComplex{P}(x::Irrational{S}, y::Real) where {P,S} = ArbComplex{P}(ArbReal{P}(x), ArbReal{P}(y)) # fallback diff --git a/test/complex.jl b/test/complex.jl index 2edd99c1..7d6aaa8e 100644 --- a/test/complex.jl +++ b/test/complex.jl @@ -1,20 +1,85 @@ + +import ArbNumerics: sign_bit, sign_bits + @testset "ArbComplex" begin @testset "complex" begin @test (1 + im) / ArbNumerics.ArbComplex(1, 1) == 1 + @test one(ArbComplex) == 1 + BITS = 53 + a = sqrt(ArbFloat(2, bits=BITS)) + b = ArbReal(pi, bits=BITS) + c = ArbComplex(a, b, bits=BITS) + @test precision(c) == BITS end # issue 34 and similar Cint-related stuff @testset "ArbComplex constructors" begin - @test_throws ErrorException ArbComplex(5.0, base = 3) - @test_throws ErrorException ArbComplex(5.0,1.0, base = 3) - @test ArbComplex(5, digits = 53) == ArbComplex{53}(5) - @test ArbComplex(5, 0, digits = 53) == ArbComplex{53}(5, 0) + @test ismissing(ArbComplex{100}(missing)) + @test_throws ErrorException ArbComplex(5.0, base=3) + @test_throws ErrorException ArbComplex(5.0, 1.0, base=3) + @test ArbComplex{128}(1) == 1 + ac = ArbComplex{64}(1.0 + 2im) + rac = real(ac) + fac = float(rac) + @test abs(ac) ≈ sqrt(abs2(fac) + abs2(imag(ac))) + @test ArbComplex{64}(ac) === ac + @test ArbComplex{64}(fac) == ArbComplex(fac) + @test ArbComplex(5, digits=53) == 5 + @test ArbComplex(5, 0, digits=53) == 5 @test ArbComplex{53}(1 + 2im) == 1 + 2im + @test sign_bit(ac) == false + @test sign_bits(ac) == (false, false) + @test ArbComplex(fac) == ArbComplex(rac) + @test copy(ac) == deepcopy(ac) + @test ArbComplex(rac, fac) == Complex(rac, fac) == complex(rac, fac) + @test ArbComplex(rac) == Complex(rac) == complex(rac) + + @test ArbComplex(rac, Int128(1)) isa ArbComplex{64} + @test ArbComplex(fac, UInt8(1)) isa ArbComplex{128} # inconsistent for fac::ArbFloat + @test ArbComplex(Float32(1.0), fac) isa ArbComplex{64} + @test ArbComplex{25}(1, fac) isa ArbComplex{25} + @test ArbComplex(π) isa ArbComplex{128} + @test ArbComplex(Int32(-1)) == ArbComplex(-1.0) @test ArbComplex(Int32(-1), Int32(-2)) == ArbComplex(-1.0, -2.0) @test ArbComplex(typemax(Int32) + 1) == ArbComplex(float(typemax(Int32) + 1)) @test hash(ArbComplex(-2, -2)) isa UInt + @test_broken ArbComplex{256}(1.0 + im) == ArbComplex(1 + im, bits=256) # TODO missing + @test ArbComplex{256}(1.0 + im) == ArbComplex(1, 1, bits=256) + + @test midpoint(ac) == ac + @test 1e-38 <= abs(radius(ArbComplex(π))) <= 2e-38 + api = ArbComplex(pi, -pi) + @test trunc(api) == ArbComplex(3 - 3im) + @test floor(api) == ArbComplex(3 - 4im) + @test ceil(api) == ArbComplex(4 - 3im) + @test trunc(Int, api) == (3, -3) + @test floor(Int8, api) == (3, -4) + @test ceil(Int16, api) == (4, -3) + @test modf(ArbComplex(3.5, -2.5)) == ((0.5, 3.0), (-0.5, -2.0)) + @test fmod(modf(api)...) == api + @test_broken angle(api) ≈ angle(Complex(api)) + + @test flipsign(api, -1) == -api + @test flipsign(api, 0.0) == api + @test flipsign(api, ArbFloat(-1)) == -api + @test flipsign(api, ArbReal(0.0)) == api + + @test copysign(-api, -1) == -api + @test copysign(-api, 0.0) == api + @test copysign(-api, ArbFloat(-1)) == -api + @test copysign(-api, ArbReal(0.0)) == api + + @test sign(api) == 1 + @test signs(api) == (1, -1) + @test signs(ArbComplex(0)) == (0, 0) + + @test signbit(api) == false + @test signbits(api) == (false, true) + + @test abs2(api) ≈ abs2(real(api)) + abs2(imag(api)) + end end From 20847080575494c5b5b2a41f71220996e133286e Mon Sep 17 00:00:00 2001 From: KlausC Date: Wed, 1 Nov 2023 17:57:38 +0100 Subject: [PATCH 03/15] Introduce `Slong` for C-interface --- src/libarb/ArbComplex.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/libarb/ArbComplex.jl b/src/libarb/ArbComplex.jl index b5850ca8..08821719 100644 --- a/src/libarb/ArbComplex.jl +++ b/src/libarb/ArbComplex.jl @@ -60,6 +60,7 @@ ArbComplex{P}(x::T) where {P, T<:Number} = ArbComplex{P}(real(x), imag(x)) const PtrToArbComplex = Ptr{ArbComplex} # acb_ptr const PtrToPtrToArbComplex = Ptr{Ptr{ArbComplex}} # acb_ptr* +const Slong = Int # to accomodate windows clear_acb(x::ArbComplex{P}) where {P} = ccall(@libarb(acb_clear), Cvoid, (Ref{ArbComplex},), x) @@ -114,11 +115,11 @@ function ArbComplex{P}(rea::Float64) where {P} return z end -const ArbInts = Union{Int64,Int32,Int16,Int8} +const ArbInts = Union{Int,Int32,Int16,Int8} # Int is Int32 on some windows enviroment function ArbComplex{P}(rea::ArbInts) where {P} z = ArbComplex{P}() - ccall(@libarb(acb_set_si), Cvoid, (Ref{ArbComplex}, Clong), z, rea) + ccall(@libarb(acb_set_si), Cvoid, (Ref{ArbComplex}, Slong), z, rea) return z end @@ -148,7 +149,7 @@ end function ArbComplex{P}(x::ArbInts, y::ArbInts) where {P} z = ArbComplex{P}() - ccall(@libarb(acb_set_si_si), Cvoid, (Ref{ArbComplex}, Clong, Clong, Clong), z, x, y, P) + ccall(@libarb(acb_set_si_si), Cvoid, (Ref{ArbComplex}, Slong, Slong, Slong), z, x, y, P) return z end From 286036c7c3618719182efce582b31c8d2d155379 Mon Sep 17 00:00:00 2001 From: KlausC Date: Wed, 1 Nov 2023 23:14:49 +0100 Subject: [PATCH 04/15] CI: windows only on nightly (load packages fails) --- .github/workflows/CI.yml | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index e1af050f..1b268e21 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -22,6 +22,11 @@ jobs: - windows-latest arch: - x64 + exclude: + - version: '1.6' + os: windows-latest + - version: '1' + os: windows-latest steps: - uses: actions/checkout@v2 - uses: julia-actions/setup-julia@v1 From 0f0e9b511b6e7e47a00c8b875da3d3604f4b450d Mon Sep 17 00:00:00 2001 From: KlausC Date: Wed, 1 Nov 2023 23:34:46 +0100 Subject: [PATCH 05/15] CI: fixed bug in CompatHelper.yml --- .github/workflows/CompatHelper.yml | 2 +- .gitignore | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/.github/workflows/CompatHelper.yml b/.github/workflows/CompatHelper.yml index 6dcb434c..1aa15a15 100644 --- a/.github/workflows/CompatHelper.yml +++ b/.github/workflows/CompatHelper.yml @@ -2,7 +2,7 @@ name: CompatHelper on: schedule: - - chron 8 3 * * 1,3,5 + - cron: 8 3 * * 1,3,5 workflow_dispatch: permissions: contents: write diff --git a/.gitignore b/.gitignore index 305731c4..0d4b7fea 100644 --- a/.gitignore +++ b/.gitignore @@ -3,6 +3,7 @@ *.jl.mem *.jl.*.mem lcov.info +.vscode Manifest.toml settings.json deps/deps.jl From dd46766b4c790211b097ea610443583e75a7dbd0 Mon Sep 17 00:00:00 2001 From: KlausC Date: Thu, 2 Nov 2023 18:26:23 +0100 Subject: [PATCH 06/15] removed ambiguities and added test cases --- src/ArbNumerics.jl | 2 +- src/float/dot.jl | 18 -------- src/libarb/ArbComplex.jl | 88 ++++++++++++++++---------------------- src/libarb/ArbFloat.jl | 2 + src/values/compare.jl | 11 +---- src/values/constructors.jl | 42 +----------------- src/values/random.jl | 8 ++-- src/values/rounding.jl | 27 +++++++----- test/complex.jl | 79 ++++++++++++++++++++++++++++------ test/runtests.jl | 17 +++++--- 10 files changed, 140 insertions(+), 154 deletions(-) diff --git a/src/ArbNumerics.jl b/src/ArbNumerics.jl index e40f2a05..e1e0dde7 100644 --- a/src/ArbNumerics.jl +++ b/src/ArbNumerics.jl @@ -169,7 +169,7 @@ include("support/ArblibVector.jl") include("support/random.jl") -const ArbNumber = Union{ArbFloat, ArbReal, ArbComplex} +const ArbNumber{P} = Union{ArbFloat{P}, ArbReal{P}, ArbComplex{P}} include("libarb/ArbMatrix.jl") # must preceed ArbRealMatrix include("libarb/ArbRealMatrix.jl") # must preceed ArbFloatMatrix diff --git a/src/float/dot.jl b/src/float/dot.jl index 990bf7a3..0692ee85 100644 --- a/src/float/dot.jl +++ b/src/float/dot.jl @@ -18,21 +18,3 @@ for (TT,dot_f) in [(:ArbReal, @libarb(arb_dot)),(:ArbComplex, @libarb(acb_dot))] end end =# - -function LinearAlgebra.dot(x::AbstractVector{T}, y::AbstractVector{T}) where {T<:ArbNumber} - length(x) == length(y) || throw(DimensionMismatch("x and y must have the same lengths")) - xy = x .* y - return sum(xy) - - #= does not work properly - xv = ArblibVector(x) - yv = ArblibVector(y) - - d = dot(xv, yv) - - free!(xv) - free!(yv) - - T(d) - =# -end diff --git a/src/libarb/ArbComplex.jl b/src/libarb/ArbComplex.jl index 08821719..7a81c808 100644 --- a/src/libarb/ArbComplex.jl +++ b/src/libarb/ArbComplex.jl @@ -23,85 +23,79 @@ Good practice is to set the precision substantively greater than the precision of the resultant width required, and to check the radius of results. See also: [`ArbFloat`](@ref), [`ArbReal`](@ref), [`ball`](@ref), [`setball`](@ref), [`midpoint`](@ref), [`radius`](@ref) -""" ArbComplex +""" +ArbComplex # ArbComplex structs hold complex balls given as ArbReal pairs -mutable struct ArbComplex{P} <: Number # P is the precision in bits - # real midpoint +mutable struct ArbComplex{P} <: Number # P is the precision in bits + # real midpoint real_mid_exp::Int # fmpz exponent of 2 (2^exp) real_mid_size::UInt # mp_size_t nwords and sign (lsb holds sign of significand) real_mid_d1::UInt # significand unsigned, immediate value or the initial span real_mid_d2::UInt # (d1, d2) the final part indicating the significand, or 0 - # real radius + # real radius real_rad_exp::Int # fmpz exponent of 2 (2^exp) real_rad_man::UInt # mp_limb_t radius, unsigned by definition - # imaginary midpoint + # imaginary midpoint imag_mid_exp::Int # fmpz exponent of 2 (2^exp) imag_mid_size::UInt # mp_size_t nwords and sign (lsb holds sign of significand) imag_mid_d1::UInt # significand unsigned, immediate value or the initial span imag_mid_d2::UInt # (d1, d2) the final part indicating the significand, or 0 - # imaginary radius + # imaginary radius imag_rad_exp::Int # fmpz exponent of 2 (2^exp) imag_rad_man::UInt # mp_limb_t radius, unsigned by definition function ArbComplex{P}() where {P} - z = new{P}(0,0,0,0,0,0,0,0,0,0,0,0) + z = new{P}(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) ccall(@libarb(acb_init), Cvoid, (Ref{ArbComplex},), z) finalizer(clear_acb, z) return z end end -ArbComplex{P}(x::T) where {P, T<:Number} = ArbComplex{P}(real(x), imag(x)) +ArbComplex{P}(x::T) where {P,T<:Number} = ArbComplex{P}(real(x), imag(x)) # for use within a struct, eg. ArbComplexMatrix const PtrToArbComplex = Ptr{ArbComplex} # acb_ptr const PtrToPtrToArbComplex = Ptr{Ptr{ArbComplex}} # acb_ptr* const Slong = Int # to accomodate windows +const ArbFloatReal{P} = Union{ArbFloat{P},ArbReal{P}} +const ArbNumber1{P} = Union{ArbFloatReal{P},ArbComplex{P}} +# finalizer: clear_acb(x::ArbComplex{P}) where {P} = ccall(@libarb(acb_clear), Cvoid, (Ref{ArbComplex},), x) ArbComplex{P}(x::ArbComplex{P}) where {P} = x ArbComplex(x::ArbComplex{P}) where {P} = x -ArbComplex{P}(x::Missing) where {P} = missing -ArbComplex(x::Missing) = missing - +ArbComplex(::Missing, ::Union{Missing,Real}...; kw...) = missing +ArbComplex(::Real, ::Missing; kw...) = missing +ArbComplex{P}(::Missing, ::Union{Missing,Real}...) where {P} = missing +ArbComplex{P}(::Real, ::Missing) where {P} = missing -function ArbComplex(x::T; bits::Int=0, digits::Int=0, base::Int=iszero(bits) ? 10 : 2) where {T<:Real} - bits > 0 && bits < MINIMUM_PRECISION_BASE2 && throw(DomainError("bit precision $bits < $MINIMUM_PRECISION_BASE2")) - digits > 0 && digits < MINIMUM_PRECISION_BASE10 && throw(DomainError("digit precision $digits < $MINIMUM_PRECISION_BASE10")) - if base === 10 - bits = digits > 0 ? bits4digits(digits)+extrabits() : (bits > 0 ? bits+extrabits() : DEFAULT_PRECISION.x) - elseif base === 2 - bits = bits > 0 ? bits+extrabits() : (digits > 0 ? digits+extrabits() : DEFAULT_PRECISION.x) - else - throw(ErrorException("base expects 2 or 10")) - end - ArbComplex{bits}(x, zero(T)) -end +ArbComplex(x::T; kw...) where {T<:Number} = ArbComplex(real(x), imag(x); kw...) -function ArbComplex(x::T, y::T; bits::Int=0, digits::Int=0, base::Int=iszero(bits) ? 10 : 2) where {T<:Real} +function ArbComplex(x::T, y::Real=zero(x); bits::Int=0, digits::Int=0, base::Int=iszero(bits) ? 10 : 2) where {T<:Real} bits > 0 && bits < MINIMUM_PRECISION_BASE2 && throw(DomainError("bit precision $bits < $MINIMUM_PRECISION_BASE2")) digits > 0 && digits < MINIMUM_PRECISION_BASE10 && throw(DomainError("digit precision $digits < $MINIMUM_PRECISION_BASE10")) if base === 10 - bits = digits > 0 ? bits4digits(digits)+extrabits() : (bits > 0 ? bits+extrabits() : DEFAULT_PRECISION.x) + bits = digits > 0 ? bits4digits(digits) + extrabits() : (bits > 0 ? bits + extrabits() : DEFAULT_PRECISION.x) elseif base === 2 - bits = bits > 0 ? bits+extrabits() : (digits > 0 ? digits+extrabits() : DEFAULT_PRECISION.x) + bits = bits > 0 ? bits + extrabits() : (digits > 0 ? digits + extrabits() : DEFAULT_PRECISION.x) else throw(ErrorException("base expects 2 or 10")) end ArbComplex{bits}(x, y) end -function ArbComplex(x::T1, y::T2; bits::Int=0, digits::Int=0, base::Int=iszero(bits) ? 10 : 2) where {T1<:Number, T2<:Number} - ArbComplex(promote(x, y)...,) -end +ArbComplex(x::ArbNumber1{P}) where {P} = ArbComplex{P}(x) -ArbComplex(re::T, im::T) where T<:AbstractFloat = ArbComplex(ArbFloat(re), ArbFloat(im)) +ArbComplex(x::ArbFloatReal{P}, y::ArbFloatReal{S}) where {P,S} = ArbComplex{min(P, S)}(x, y) +ArbComplex(x::ArbFloatReal{P}, y::Real) where {P} = ArbComplex{P}(x, y) +ArbComplex(x::Real, y::ArbFloatReal{P}) where {P} = ArbComplex{P}(x, y) @inline sign_bit(x::ArbComplex{P}) where {P} = isodd(x.real_mid_size) @inline sign_bits(x::ArbComplex{P}) where {P} = isodd(x.real_mid_size), isodd(x.imag_mid_size) @@ -115,7 +109,7 @@ function ArbComplex{P}(rea::Float64) where {P} return z end -const ArbInts = Union{Int,Int32,Int16,Int8} # Int is Int32 on some windows enviroment +const ArbInts = Union{Int,Int32,Int16,Int8} # Int may be Int32 function ArbComplex{P}(rea::ArbInts) where {P} z = ArbComplex{P}() @@ -162,13 +156,6 @@ end deepcopy(x::ArbComplex{P}) where {P} = copy(x) -ArbComplex(x::ArbFloat{P}) where {P} = ArbComplex{P}(x) -ArbComplex(x::ArbReal{P}) where {P} = ArbComplex{P}(x) -ArbComplex(x::ArbFloat{P}, y::ArbFloat{P}) where {P} = ArbComplex{P}(x,y) -ArbComplex(x::ArbReal{P}, y::ArbReal{P}) where {P} = ArbComplex{P}(x,y) -ArbComplex(x::ArbFloat{P}, y::ArbReal{P}) where {P} = ArbComplex{P}(ArbReal{P}(x),y) -ArbComplex(x::ArbReal{P}, y::ArbFloat{P}) where {P} = ArbComplex{P}(x,ArbReal{P}(y)) - Base.Complex(x::ArbFloat{P}) where {P} = ArbComplex(x) Base.Complex(x::ArbReal{P}) where {P} = ArbComplex(x) Base.Complex(re::ArbFloat{P}, im::ArbFloat{P}) where {P} = ArbComplex(re, im) @@ -196,14 +183,14 @@ end function ArbComplex{P}(x::T1, y::T2) where {P,T1<:Union{Integer,AbstractFloat},T2<:Union{Integer,AbstractFloat}} x1 = ArbReal{P}(x) y1 = ArbReal{P}(y) - z = ArbComplex{P}(x1, y1) + z = ArbComplex{P}(x1, y1) return z end function ArbComplex{P}(x::Irrational{S1}, y::Irrational{S2}) where {P,S1,S2} x1 = ArbReal{P}(x) y1 = ArbReal{P}(y) - z = ArbComplex{P}(x1, y1) + z = ArbComplex{P}(x1, y1) return z end @@ -228,10 +215,8 @@ function ArbComplex{P}(x::T, y::ArbFloat{P}) where {P,T<:Union{Integer,AbstractF return z end -ArbFloat{P}(x::ArbComplex{P}, bits::Int) where P = ArbFloat(real(x), bits=bits) -ArbReal{P}(x::ArbComplex{P}, bits::Int) where P = ArbReal(real(x), bits=bits) -ArbComplex{P}(x::ArbReal{P}, bits::Int) where P = ArbComplex(x, bits=bits) -ArbComplex{P}(x::ArbFloat{P}, bits::Int) where P = ArbComplex(x, bits=bits) +ArbFloat{P}(x::ArbComplex{P}, bits::Int) where {P} = ArbFloat(real(x), bits=bits) +ArbReal{P}(x::ArbComplex{P}, bits::Int) where {P} = ArbReal(real(x), bits=bits) function ArbComplex{P}(x::ArbFloat{P}, y::T) where {P,T<:Union{Integer,AbstractFloat}} x1 = ArbReal{P}(x) @@ -340,13 +325,13 @@ end trunc(x::ArbComplex{P}) where {P} = ArbComplex{P}(trunc(real(x)), trunc(imag(x))) -trunc(::Type{T}, x::ArbComplex{P}) where {P, T<:Integer} = T(trunc(real(x))), T(trunc(imag(x))) +trunc(::Type{T}, x::ArbComplex{P}) where {P,T<:Integer} = T(trunc(real(x))), T(trunc(imag(x))) floor(x::ArbComplex{P}) where {P} = ArbComplex{P}(floor(real(x)), floor(imag(x))) -floor(::Type{T}, x::ArbComplex{P}) where {P, T<:Integer} = T(floor(real(x))), T(floor(imag(x))) +floor(::Type{T}, x::ArbComplex{P}) where {P,T<:Integer} = T(floor(real(x))), T(floor(imag(x))) ceil(x::ArbComplex{P}) where {P} = ArbComplex{P}(ceil(real(x)), ceil(imag(x))) -ceil(::Type{T}, x::ArbComplex{P}) where {P, T<:Integer} = T(ceil(real(x))), T(ceil(imag(x))) +ceil(::Type{T}, x::ArbComplex{P}) where {P,T<:Integer} = T(ceil(real(x))), T(ceil(imag(x))) function modf(x::ArbComplex{P}) where {P} repart = modf(real(x)) @@ -354,7 +339,7 @@ function modf(x::ArbComplex{P}) where {P} return (repart, impart) end -fmod(repart::Tuple{ArbReal{P1}, ArbReal{P1}}, impart::Tuple{ArbReal{P2},ArbReal{P2}}) where {P1, P2} = +fmod(repart::Tuple{ArbReal{P1},ArbReal{P1}}, impart::Tuple{ArbReal{P2},ArbReal{P2}}) where {P1,P2} = ArbComplex(fmod(repart), fmod(impart)) # phase angle @@ -364,7 +349,7 @@ function Base.angle(x::ArbComplex{P}) where {P} y = hypot(rea - 1.0, ima) x = hypot(rea + 1.0, ima) - a = 2 * atan(y , x) + a = 2 * atan(y, x) T = ArbFloat{P} !(signbit(a) || signbit(T(pi) - a)) ? a : (signbit(a) ? zero(T) : T(pi)) end @@ -374,6 +359,7 @@ const hash_arbcomplex_lo = (UInt === UInt64) ? 0x76143ad985246e79 : 0x5b6a64dc const hash_0_arbcomplex_lo = hash(zero(UInt), hash_arbcomplex_lo) function Base.hash(z::ArbComplex{P}, h::UInt) where {P} hash(z.real_mid_d1 ⊻ z.real_rad_exp ⊻ z.imag_mid_d1 ⊻ z.imag_rad_exp, - h ⊻ hash(z.real_mid_d2 ⊻ (~reinterpret(UInt,P)), hash_arbcomplex_lo) - ⊻ hash_0_arbcomplex_lo) + h ⊻ hash(z.real_mid_d2 ⊻ (~reinterpret(UInt, P)), hash_arbcomplex_lo) + ⊻ + hash_0_arbcomplex_lo) end diff --git a/src/libarb/ArbFloat.jl b/src/libarb/ArbFloat.jl index 67580ac5..708a0d3f 100644 --- a/src/libarb/ArbFloat.jl +++ b/src/libarb/ArbFloat.jl @@ -155,11 +155,13 @@ for (F,A) in ((:floor, :arf_floor), (:ceil, :arf_ceil)) ccall(@libarb($A), Cvoid, (Ref{ArbFloat}, Ref{ArbFloat}), z, x) return z end + $F(::Type{Bool}, x::ArbFloat{P}) where {P} = Bool($F(x)) $F(::Type{T}, x::ArbFloat{P}) where {P, T<:Integer} = T($F(x)) end end trunc(x::ArbFloat{P}) where {P} = signbit(x) ? ceil(x) : floor(x) +trunc(::Type{Bool}, x::ArbFloat{P}) where {P} = Bool(trunc(x)) trunc(::Type{T}, x::ArbFloat{P}) where {P, T<:Integer} = T(trunc(x)) midpoint(x::ArbFloat{P}) where {P} = x diff --git a/src/values/compare.jl b/src/values/compare.jl index 008e2bf4..f778ce95 100644 --- a/src/values/compare.jl +++ b/src/values/compare.jl @@ -81,12 +81,12 @@ function (>)(x::ArbReal{P}, y::ArbReal{P}) where {P} end function (<=)(x::ArbReal{P}, y::ArbReal{P}) where {P} - x < y || + x < y || 0 != ccall(@libarb(arb_contains), Cint, (Ref{ArbReal}, Ref{ArbReal}), x, y) end function (>=)(x::ArbReal{P}, y::ArbReal{P}) where {P} - x > y || + x > y || 0 != ccall(@libarb(arb_contains), Cint, (Ref{ArbReal}, Ref{ArbReal}), x, y) end @@ -214,10 +214,3 @@ isapprox(x::ArbReal{P}, y::F) where {P, F<:IEEEFloat} = isapprox(F(x), y) isapprox(x::F, y::ArbReal{P}) where {P, F<:IEEEFloat} = isapprox(x, F(y)) isapprox(x::ArbComplex{P}, y::F) where {P, F<:IEEEFloat} = isapprox(F(real(x)), y) isapprox(x::F, y::ArbComplex{P}) where {P, F<:IEEEFloat} = isapprox(x, F(real(y))) - -isapprox(x::ArbFloat{P}, y::Real) where {P} = isapprox(x, ArbFloat{P}(y)) -isapprox(x::Real, y::ArbFloat{P}) where {P} = isapprox(ArbFloat{P}(x), y) -isapprox(x::ArbReal{P}, y::Real) where {P} = isapprox(x, ArbReal{P}(y)) -isapprox(x::Number, y::ArbReal{P}) where {P} = isapprox(ArbReal{P}(x), y) -isapprox(x::ArbComplex{P}, y::Number) where {P} = isapprox(x, ArbComplex{P}(y)) -isapprox(x::Number, y::ArbComplex{P}) where {P} = isapprox(ArbComplex{P}(x), y) diff --git a/src/values/constructors.jl b/src/values/constructors.jl index 9fb62620..48c65dd0 100644 --- a/src/values/constructors.jl +++ b/src/values/constructors.jl @@ -1,38 +1,3 @@ -@inline function ArbFloat(x) - prec = DEFAULT_PRECISION.x - res = ArbFloat{prec}(x) - return res -end -@inline function ArbReal(x) - prec = DEFAULT_PRECISION.x - res = ArbReal{prec}(x) - return res -end - -@inline function ArbComplex(x) - prec = DEFAULT_PRECISION.x - res = ArbComplex{prec}(x) - return res -end - -@inline function ArbComplex(x, y) - prec = DEFAULT_PRECISION.x - res = ArbComplex{prec}(x, y) - return res -end - -@inline function ArbComplex(x::I, y::T) where {I<:Integer, T<:AbstractFloat} - prec = DEFAULT_PRECISION.x - res = ArbComplex{prec}(T(x), y) - return res -end - -@inline function ArbComplex(x::T, y::I) where {I<:Integer, T<:AbstractFloat} - prec = DEFAULT_PRECISION.x - res = ArbComplex{prec}(x, T(y)) - return res -end - # IEEEFloat # rounds up (widens) @@ -150,14 +115,9 @@ function ArbReal(x::Irrational{S}) where {S} return setball(mid, rad) end -ArbComplex(x::Irrational{S}) where {S} = ArbComplex(ArbReal(x), ArbReal(0)) -ArbComplex{P}(x::Irrational{S}) where {P,S} = ArbComplex{P}(ArbReal{P}(x), ArbReal{P}(0)) -ArbComplex(x::Irrational{S}, y::Real) where {S} = ArbComplex(ArbReal(x), ArbReal(y)) -ArbComplex{P}(x::Irrational{S}, y::Real) where {P,S} = ArbComplex{P}(ArbReal{P}(x), ArbReal{P}(y)) - # fallback -ArbFloat{P}(x::T) where {P,T<:Real} = ArbFloat{P}(BigFloat(x)) +# ArbFloat{P}(x::T) where {P,T<:Real} = ArbFloat{P}(BigFloat(x)) ArbFloat(x::T) where {T<:Real} = ArbFloat{workingprecision(ArbFloat)}(BigFloat(x)) ArbFloat{P}(x::T) where {P,T<:Complex} = ArbFloat{P}(BigFloat(real(x))) ArbFloat(x::T) where {T<:Complex} = ArbFloat{workingprecision(ArbFloat)}(BigFloat(real(x))) diff --git a/src/values/random.jl b/src/values/random.jl index 7b5138fe..0ce3e3fb 100644 --- a/src/values/random.jl +++ b/src/values/random.jl @@ -1,7 +1,7 @@ -function rand(rng::AbstractRNG, ::Random.SamplerTrivial{Random.CloseOpen01{ArbFloat{P}}},n) where P +function rand(rng::AbstractRNG, ::Random.SamplerTrivial{Random.CloseOpen01{ArbFloat{P}}}, dims::Integer...) where P bfprec = precision(BigFloat) setprecision(BigFloat, P) - bfrand = rand(BigFloat,n) + bfrand = rand(rng, BigFloat, dims...) result = ArbFloat{P}.(bfrand) setprecision(BigFloat, bfprec); println("ok") return result @@ -16,10 +16,10 @@ function rand(::Type{ArbFloat{P}}, n::Integer) where P return result end -function rand(rng::AbstractRNG, ::Random.SamplerTrivial{Random.CloseOpen01{ArbReal{P}}},n) where P +function rand(rng::AbstractRNG, ::Random.SamplerTrivial{Random.CloseOpen01{ArbReal{P}}}, dims::Integer...) where P bfprec = precision(BigFloat) setprecision(BigFloat, P) - bfrand = rand(BigFloat,n) + bfrand = rand(rng, BigFloat, dims...) result = ArbReal{P}.(bfrand) setprecision(BigFloat, bfprec); println("ok") return result diff --git a/src/values/rounding.jl b/src/values/rounding.jl index ce9f11b8..bd1d4a17 100644 --- a/src/values/rounding.jl +++ b/src/values/rounding.jl @@ -144,7 +144,16 @@ rounddigits(x::ArbComplex{P}, roundingmode::RoundingMode, digits::Int) where {P} rounddigits(x::ArbComplex{P}, digits::Int) where {P} = roundbits(x, RoundNearest, bits4digits(digits)) +function round(x::ArbFloat{P}, roundingmode::RoundingMode{:NearestTiesUp}; sigdigits::Integer, base::Integer=10) where {P}v + _round(x, roundingmode, sigdigits, base) +end +function round(x::ArbFloat{P}, roundingmode::RoundingMode{:NearestTiesAway}; sigdigits::Integer, base::Integer=10) where {P} + _round(x, roundingmode, sigdigits, base) +end function round(x::ArbFloat{P}, roundingmode::RoundingMode; sigdigits::Integer, base::Integer=10) where {P} + _round(x, roundingmode, sigdigits, base) +end +function _round(x::ArbFloat{P}, roundingmode::RoundingMode, sigdigits::Integer, base::Integer=10) where {P} r = round!(x, roundingmode, sigdigits=sigdigits, base=base) s = string(r) result = ArbFloat{P}(s) @@ -201,7 +210,7 @@ function round!(x::ArbComplex{P}, roundingmode::RoundingMode; sigdigits::Integer end end - + function roundfrac(x::ArbFloat{P}, roundingmode::RoundingMode, sigdigits::Integer, base::Integer) where {P} if base==10 return rounddigits(x, roundingmode, sigdigits) @@ -214,7 +223,7 @@ end -#= +#= function round(x::ArbFloat{P}, bitprecision::Int) where {P} z = ArbFloat{P}() rounding = match_rounding_mode(RoundNearest) @@ -237,18 +246,18 @@ function round(x::ArbComplex{P}, bitprecision::Int) where {P} end =# -#= +#= for A in (:ArbFloat, :ArbReal, :ArbComplex) @eval begin - + round(::Type{T}, x::$A{P}, rounding::RoundingMode=RoundNearest) where {P,T} = T(round(x, digits=P, base=rounding)) - function round(x::$A{P}, rounding::RoundingMode=RoundNearest; + function round(x::$A{P}, rounding::RoundingMode=RoundNearest; sigdigits::Integer = 0, digits::Integer = 0, base = 10) where {P} - sigdigits = max(sigdigits, digits) + sigdigits = max(sigdigits, digits) if base == 10 - round(x, digits=digits4bits(sigdigits), base=rounding) + round(x, digits=digits4bits(sigdigits), base=rounding) elseif base == 2 round(x, digits=sigdigits, base=rounding) else @@ -256,7 +265,7 @@ for A in (:ArbFloat, :ArbReal, :ArbComplex) end end end - + end =# @@ -280,5 +289,3 @@ end integer, with ties (fractional values of 0.5) being rounded to the nearest even integer. Note that round may give incorrect results if the global rounding mode is changed (see rounding). =# - - diff --git a/test/complex.jl b/test/complex.jl index 7d6aaa8e..d8e1816c 100644 --- a/test/complex.jl +++ b/test/complex.jl @@ -1,6 +1,3 @@ - -import ArbNumerics: sign_bit, sign_bits - @testset "ArbComplex" begin @testset "complex" begin @@ -13,9 +10,28 @@ import ArbNumerics: sign_bit, sign_bits @test precision(c) == BITS end - # issue 34 and similar Cint-related stuff + @testset "missing data" begin + @test ismissing(ArbFloat(missing)) + @test ismissing(ArbReal(missing)) + @test ismissing(ArbComplex(missing)) + + @test ismissing(ArbComplex(missing, 1)) + @test ismissing(ArbComplex(1, missing)) + @test ismissing(ArbComplex(missing, missing)) + @test_throws MethodError ismissing(ArbComplex(1+im, missing)) + @test_throws MethodError ismissing(ArbComplex(missing, 1+im)) + + @test ismissing(ArbFloat{99}(missing)) + @test ismissing(ArbReal{99}(missing)) + @test ismissing(ArbComplex{99}(missing)) + + @test ismissing(ArbComplex{99}(missing, 1)) + @test ismissing(ArbComplex{99}(1, missing)) + @test ismissing(ArbComplex{99}(missing, missing)) + @test_throws MethodError ismissing(ArbComplex{99}(1+im, missing)) + @test_throws MethodError ismissing(ArbComplex{99}(missing, 1+im)) + end @testset "ArbComplex constructors" begin - @test ismissing(ArbComplex{100}(missing)) @test_throws ErrorException ArbComplex(5.0, base=3) @test_throws ErrorException ArbComplex(5.0, 1.0, base=3) @test ArbComplex{128}(1) == 1 @@ -28,24 +44,24 @@ import ArbNumerics: sign_bit, sign_bits @test ArbComplex(5, digits=53) == 5 @test ArbComplex(5, 0, digits=53) == 5 @test ArbComplex{53}(1 + 2im) == 1 + 2im - @test sign_bit(ac) == false - @test sign_bits(ac) == (false, false) + @test signbit(ac) == false + @test signbits(ac) == (false, false) @test ArbComplex(fac) == ArbComplex(rac) @test copy(ac) == deepcopy(ac) @test ArbComplex(rac, fac) == Complex(rac, fac) == complex(rac, fac) @test ArbComplex(rac) == Complex(rac) == complex(rac) - @test ArbComplex(rac, Int128(1)) isa ArbComplex{64} - @test ArbComplex(fac, UInt8(1)) isa ArbComplex{128} # inconsistent for fac::ArbFloat - @test ArbComplex(Float32(1.0), fac) isa ArbComplex{64} - @test ArbComplex{25}(1, fac) isa ArbComplex{25} - @test ArbComplex(π) isa ArbComplex{128} + @test typeof(ArbComplex(rac, Int128(1))) == ArbComplex{64} + @test typeof(ArbComplex(fac, UInt8(1))) == ArbComplex{64} + @test typeof(ArbComplex(Float32(1.0), fac)) == ArbComplex{64} + @test typeof(ArbComplex{25}(1, fac)) == ArbComplex{25} + @test typeof(ArbComplex(π)) == ArbComplex{128} @test ArbComplex(Int32(-1)) == ArbComplex(-1.0) @test ArbComplex(Int32(-1), Int32(-2)) == ArbComplex(-1.0, -2.0) @test ArbComplex(typemax(Int32) + 1) == ArbComplex(float(typemax(Int32) + 1)) @test hash(ArbComplex(-2, -2)) isa UInt - @test_broken ArbComplex{256}(1.0 + im) == ArbComplex(1 + im, bits=256) # TODO missing + @test ArbComplex{256}(1.0 + im) == ArbComplex(1 + im, bits=256) @test ArbComplex{256}(1.0 + im) == ArbComplex(1, 1, bits=256) @test midpoint(ac) == ac @@ -82,4 +98,41 @@ import ArbNumerics: sign_bit, sign_bits end + INF = typemax(Int) + sprec(a::ArbNumber) = workingprecision(a) + sprec(::Number) = workingprecision(ArbFloat) + dprec(a::ArbNumber) = sprec(a) + dprec(::Number) = INF + + TEST_TYPES = (Int8, Int, UInt32, Float16, Float64, BigFloat, ArbFloat{200}, ArbReal{200}) + + @testset "ArbComplex(::$T)" for T in (TEST_TYPES..., ArbComplex{200}) + a = T(10) + DEF = sprec(a) + c = ArbComplex(a) + @test c == 10 + @test typeof(c) == ArbComplex{DEF} + + d = ArbComplex{100}(a) + @test typeof(d) == ArbComplex{100} + end + + @testset "ArbComplex(::$TA, ::$TB)" for TA in TEST_TYPES, TB in TEST_TYPES + a = TA(10) + b = TB(20) + DEF = min(dprec(a), dprec(b)) + DEF == INF && (DEF = sprec(a)) + c = ArbComplex(a, b) + @test c == 10 + 20im + @test typeof(c) == ArbComplex{DEF} || a isa ArbNumber || b isa ArbNumber # TODO should be true also for ArbNumbers + + d = ArbComplex{250}(a, b) + @test typeof(d) == ArbComplex{250} + end + + @testset "ArbComplex(::Int, ::ArbComplex)" begin + @test_throws MethodError ArbComplex(-12, ArbComplex(1)) + @test_throws MethodError ArbComplex(ArbComplex(1), 0) + end + end diff --git a/test/runtests.jl b/test/runtests.jl index b3d00955..9679c3be 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,11 +8,14 @@ setprecision(BigFloat, 512) # end @testset "ArbNumerics tests" begin - include("specialvalues.jl") - include("compare.jl") - include("arithmetic.jl") - include("elementaryfunctions.jl") - include("misc.jl") - include("complex.jl") - include("dft.jl") + @testset "ArbNumerics without ambiguous methods" begin + @test isempty(detect_ambiguities(ArbNumerics)) + end + include("specialvalues.jl") + include("compare.jl") + include("arithmetic.jl") + include("elementaryfunctions.jl") + include("misc.jl") + include("complex.jl") + include("dft.jl") end From 6fd75e0fcdf07eb1c3ff2fbe17a72daad9b0f926 Mon Sep 17 00:00:00 2001 From: KlausC Date: Thu, 2 Nov 2023 18:35:18 +0100 Subject: [PATCH 07/15] fix `angle` --- src/libarb/ArbComplex.jl | 9 +-------- test/complex.jl | 2 +- 2 files changed, 2 insertions(+), 9 deletions(-) diff --git a/src/libarb/ArbComplex.jl b/src/libarb/ArbComplex.jl index 7a81c808..d7a75a40 100644 --- a/src/libarb/ArbComplex.jl +++ b/src/libarb/ArbComplex.jl @@ -344,14 +344,7 @@ fmod(repart::Tuple{ArbReal{P1},ArbReal{P1}}, impart::Tuple{ArbReal{P2},ArbReal{P # phase angle function Base.angle(x::ArbComplex{P}) where {P} - rea, ima = reim(x / hypot(reim(x)...,)) - - y = hypot(rea - 1.0, ima) - x = hypot(rea + 1.0, ima) - - a = 2 * atan(y, x) - T = ArbFloat{P} - !(signbit(a) || signbit(T(pi) - a)) ? a : (signbit(a) ? zero(T) : T(pi)) + atan(imag(x), real(x)) end # a type specific hash function helps the type to 'just work' diff --git a/test/complex.jl b/test/complex.jl index d8e1816c..a0d5b3ac 100644 --- a/test/complex.jl +++ b/test/complex.jl @@ -75,7 +75,7 @@ @test ceil(Int16, api) == (4, -3) @test modf(ArbComplex(3.5, -2.5)) == ((0.5, 3.0), (-0.5, -2.0)) @test fmod(modf(api)...) == api - @test_broken angle(api) ≈ angle(Complex(api)) + @test angle(api) ≈ angle(Complex(api)) @test flipsign(api, -1) == -api @test flipsign(api, 0.0) == api From c6130fbe8e24f091808748046f0eb0ac1bd273d9 Mon Sep 17 00:00:00 2001 From: KlausC Date: Thu, 2 Nov 2023 18:39:37 +0100 Subject: [PATCH 08/15] typo --- src/values/rounding.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/values/rounding.jl b/src/values/rounding.jl index bd1d4a17..66a4b909 100644 --- a/src/values/rounding.jl +++ b/src/values/rounding.jl @@ -144,7 +144,7 @@ rounddigits(x::ArbComplex{P}, roundingmode::RoundingMode, digits::Int) where {P} rounddigits(x::ArbComplex{P}, digits::Int) where {P} = roundbits(x, RoundNearest, bits4digits(digits)) -function round(x::ArbFloat{P}, roundingmode::RoundingMode{:NearestTiesUp}; sigdigits::Integer, base::Integer=10) where {P}v +function round(x::ArbFloat{P}, roundingmode::RoundingMode{:NearestTiesUp}; sigdigits::Integer, base::Integer=10) where {P} _round(x, roundingmode, sigdigits, base) end function round(x::ArbFloat{P}, roundingmode::RoundingMode{:NearestTiesAway}; sigdigits::Integer, base::Integer=10) where {P} From 764278a6f97b77a11c8ceb620de592c38b4e2991 Mon Sep 17 00:00:00 2001 From: KlausC Date: Fri, 3 Nov 2023 00:25:03 +0100 Subject: [PATCH 09/15] simplifications and code coverage --- src/ArbNumerics.jl | 2 +- src/libarb/ArbComplex.jl | 148 +++---------------------------------- src/support/complex.jl | 67 ++++------------- src/values/constructors.jl | 22 ------ test/complex.jl | 3 +- 5 files changed, 30 insertions(+), 212 deletions(-) diff --git a/src/ArbNumerics.jl b/src/ArbNumerics.jl index e1e0dde7..859623d7 100644 --- a/src/ArbNumerics.jl +++ b/src/ArbNumerics.jl @@ -85,7 +85,7 @@ import Base: IEEEFloat, Float16, Float32, Float64, float, UInt8, UInt16, UInt32, UInt64, UInt128, Int8, Int16, Int32, Int64, Int128, - BigInt, BigFloat, Rational, Complex, real, imag, complex, + BigInt, BigFloat, Rational, Complex, real, imag, complex, angle, floatmax, floatmin, typemax, typemin, maxintfloat, rationalize, diff --git a/src/libarb/ArbComplex.jl b/src/libarb/ArbComplex.jl index d7a75a40..96188c6c 100644 --- a/src/libarb/ArbComplex.jl +++ b/src/libarb/ArbComplex.jl @@ -15,7 +15,7 @@ The precision of an `ArbComplex` must be >= 24 bits or >= 8 digits. -Internally, an `ArbComplex` represents as a pair of `ArbReals`, where the real part and the imaginary part are each +Internally, an `ArbComplex` represents as a pair of `ArbReal`, where the real part and the imaginary part are each intervals given by midpoint and radius (half-width of the interval). This representation serves as a _ball_ over the complex numbers. Calculations with `ArbComplex` generate results that are assured to contain the true mathematical result. There is no a priori assurance on the width of the interval that results. @@ -64,6 +64,8 @@ const PtrToPtrToArbComplex = Ptr{Ptr{ArbComplex}} # acb_ptr* const Slong = Int # to accomodate windows const ArbFloatReal{P} = Union{ArbFloat{P},ArbReal{P}} const ArbNumber1{P} = Union{ArbFloatReal{P},ArbComplex{P}} +const ArbInts = Union{Int,Int32,Int16,Int8} # Int may be Int32 +const ArbRealTypes = Real # Union{Integer,AbstractFloat,ArbReal} # finalizer: clear_acb(x::ArbComplex{P}) where {P} = ccall(@libarb(acb_clear), Cvoid, (Ref{ArbComplex},), x) @@ -78,7 +80,7 @@ ArbComplex{P}(::Real, ::Missing) where {P} = missing ArbComplex(x::T; kw...) where {T<:Number} = ArbComplex(real(x), imag(x); kw...) -function ArbComplex(x::T, y::Real=zero(x); bits::Int=0, digits::Int=0, base::Int=iszero(bits) ? 10 : 2) where {T<:Real} +function ArbComplex(x::T, y::Real=0; bits::Int=0, digits::Int=0, base::Int=iszero(bits) ? 10 : 2) where {T<:Real} bits > 0 && bits < MINIMUM_PRECISION_BASE2 && throw(DomainError("bit precision $bits < $MINIMUM_PRECISION_BASE2")) digits > 0 && digits < MINIMUM_PRECISION_BASE10 && throw(DomainError("digit precision $digits < $MINIMUM_PRECISION_BASE10")) if base === 10 @@ -88,7 +90,7 @@ function ArbComplex(x::T, y::Real=zero(x); bits::Int=0, digits::Int=0, base::Int else throw(ErrorException("base expects 2 or 10")) end - ArbComplex{bits}(x, y) + iszero(y) ? ArbComplex{bits}(x) : ArbComplex{bits}(x, y) end ArbComplex(x::ArbNumber1{P}) where {P} = ArbComplex{P}(x) @@ -103,31 +105,18 @@ ArbComplex(x::Real, y::ArbFloatReal{P}) where {P} = ArbComplex{P}(x, y) @inline sign_bit(x::ArbComplex{P}, ::Type{ImagPart}) where {P} = isodd(x.imag_mid_size) # ArbComplex{P}(sqrt(2.0)) -function ArbComplex{P}(rea::Float64) where {P} +function ArbComplex{P}(rea::Base.IEEEFloat) where {P} z = ArbComplex{P}() ccall(@libarb(acb_set_d), Cvoid, (Ref{ArbComplex}, Float64), z, rea) return z end -const ArbInts = Union{Int,Int32,Int16,Int8} # Int may be Int32 - function ArbComplex{P}(rea::ArbInts) where {P} z = ArbComplex{P}() ccall(@libarb(acb_set_si), Cvoid, (Ref{ArbComplex}, Slong), z, rea) return z end - -function ArbComplex{P}(rea::ArbFloat{P}) where {P} - x = ArbReal{P}() - ccall(@libarb(arb_set_arf), Cvoid, (Ref{ArbReal}, Ref{ArbFloat}), x, rea) - y = ArbReal{P}() - ccall(@libarb(arb_set_arf), Cvoid, (Ref{ArbReal}, Ref{ArbFloat}), y, x) - z = ArbComplex{P}() - ccall(@libarb(acb_set_arb), Cvoid, (Ref{ArbComplex}, Ref{ArbReal}), z, y) - return z -end - function ArbComplex{P}(rea::ArbReal{P}) where {P} z = ArbComplex{P}() ccall(@libarb(acb_set_arb), Cvoid, (Ref{ArbComplex}, Ref{ArbReal}), z, rea) @@ -135,7 +124,7 @@ function ArbComplex{P}(rea::ArbReal{P}) where {P} end # ArbComplex{P}(sqrt(2.0), inv(17.0)) -function ArbComplex{P}(re::Float64, im::Float64) where {P} +function ArbComplex{P}(re::Base.IEEEFloat, im::Base.IEEEFloat) where {P} z = ArbComplex{P}() ccall(@libarb(acb_set_d_d), Cvoid, (Ref{ArbComplex}, Float64, Float64), z, re, im) return z @@ -180,118 +169,9 @@ function ArbComplex{P}(x::T) where {P,T<:Integer} end -function ArbComplex{P}(x::T1, y::T2) where {P,T1<:Union{Integer,AbstractFloat},T2<:Union{Integer,AbstractFloat}} - x1 = ArbReal{P}(x) - y1 = ArbReal{P}(y) - z = ArbComplex{P}(x1, y1) - return z -end - -function ArbComplex{P}(x::Irrational{S1}, y::Irrational{S2}) where {P,S1,S2} - x1 = ArbReal{P}(x) - y1 = ArbReal{P}(y) - z = ArbComplex{P}(x1, y1) - return z -end - -function ArbComplex{P}(x::Irrational{S}, y::T) where {P,S,T<:Union{Integer,AbstractFloat}} - x1 = ArbReal{P}(x) - y1 = ArbReal{P}(y) - z = ArbComplex{P}(x1, y1) - return z -end - -function ArbComplex{P}(x::T, y::Irrational{S}) where {P,S,T<:Union{Integer,AbstractFloat}} - x1 = ArbReal{P}(x) - y1 = ArbReal{P}(y) - z = ArbComplex{P}(x1, y1) - return z -end - -function ArbComplex{P}(x::T, y::ArbFloat{P}) where {P,T<:Union{Integer,AbstractFloat}} - x1 = ArbReal{P}(x) - y1 = ArbReal{P}(y) - z = ArbComplex{P}(x1, y1) - return z -end - -ArbFloat{P}(x::ArbComplex{P}, bits::Int) where {P} = ArbFloat(real(x), bits=bits) -ArbReal{P}(x::ArbComplex{P}, bits::Int) where {P} = ArbReal(real(x), bits=bits) - -function ArbComplex{P}(x::ArbFloat{P}, y::T) where {P,T<:Union{Integer,AbstractFloat}} - x1 = ArbReal{P}(x) - y1 = ArbReal{P}(y) - z = ArbComplex{P}(x1, y1) - return z -end - -function ArbComplex{P}(x::T, y::ArbReal{P}) where {P,T<:Union{Integer,AbstractFloat}} - x1 = ArbReal{P}(x) - z = ArbComplex{P}(x1, y) - return z -end -function ArbComplex{P}(x::ArbReal{P}, y::T) where {P,T<:Union{Integer,AbstractFloat}} - y1 = ArbReal{P}(y) - z = ArbComplex{P}(x, y1) - return z -end - -function ArbComplex{P}(x::ArbFloat{P}, y::ArbReal{P}) where {P} - x1 = ArbReal{P}(x) - z = ArbComplex{P}(x1, y) - return z -end -function ArbComplex{P}(x::ArbReal{P}, y::ArbFloat{P}) where {P} - y1 = ArbReal{P}(y) - z = ArbComplex{P}(x, y1) - return z -end - -function ArbComplex{P}(x::T, y::ArbFloat{S}) where {P,S,T<:Union{Integer,AbstractFloat}} - x1 = ArbReal{P}(x) - y1 = ArbReal{P}(y) - z = ArbComplex{P}(x1, y1) - return z -end -function ArbComplex{P}(x::ArbFloat{S}, y::T) where {P,S,T<:Union{Integer,AbstractFloat}} - x1 = ArbReal{P}(x) - y1 = ArbReal{P}(y) - z = ArbComplex{P}(x1, y1) - return z -end -function ArbComplex{P}(x::T, y::ArbReal{S}) where {P,S,T<:Union{Integer,AbstractFloat}} - x1 = ArbReal{P}(x) - y1 = ArbReal{P}(y) - z = ArbComplex{P}(x1, y1) - return z -end -function ArbComplex{P}(x::ArbReal{S}, y::T) where {P,S,T<:Union{Integer,AbstractFloat}} - x1 = ArbReal{P}(x) - y1 = ArbReal{P}(y) - z = ArbComplex{P}(x1, y1) - return z -end -function ArbComplex{P}(x::ArbFloat{S1}, y::ArbFloat{S2}) where {P,S1,S2} - x1 = ArbReal{P}(x) - y1 = ArbReal{P}(y) - z = ArbComplex{P}(x1, y1) - return z -end -function ArbComplex{P}(x::ArbReal{S1}, y::ArbReal{S2}) where {P,S1,S2} - x1 = ArbReal{P}(x) - y1 = ArbReal{P}(y) - z = ArbComplex{P}(x1, y1) - return z -end -function ArbComplex{P}(x::ArbFloat{S1}, y::ArbReal{S2}) where {P,S1,S2} - x1 = ArbReal{P}(x) - y1 = ArbReal{P}(y) - z = ArbComplex{P}(x1, y1) - return z -end -function ArbComplex{P}(x::ArbReal{S1}, y::ArbFloat{S2}) where {P,S1,S2} - x1 = ArbReal{P}(x) - y1 = ArbReal{P}(y) +function ArbComplex{P}(x::T1, y::T2) where {P,T1<:ArbRealTypes,T2<:ArbRealTypes} + x1 = convert(ArbReal{P}, x) + y1 = convert(ArbReal{P}, y) z = ArbComplex{P}(x1, y1) return z end @@ -322,8 +202,6 @@ function radius(x::ArbComplex{P}) where {P} return ArbComplex{P}(zreal, zimag) end - - trunc(x::ArbComplex{P}) where {P} = ArbComplex{P}(trunc(real(x)), trunc(imag(x))) trunc(::Type{T}, x::ArbComplex{P}) where {P,T<:Integer} = T(trunc(real(x))), T(trunc(imag(x))) @@ -339,12 +217,8 @@ function modf(x::ArbComplex{P}) where {P} return (repart, impart) end -fmod(repart::Tuple{ArbReal{P1},ArbReal{P1}}, impart::Tuple{ArbReal{P2},ArbReal{P2}}) where {P1,P2} = +function fmod(repart::Tuple{ArbReal{P1},ArbReal{P1}}, impart::Tuple{ArbReal{P2},ArbReal{P2}}) where {P1,P2} ArbComplex(fmod(repart), fmod(impart)) - -# phase angle -function Base.angle(x::ArbComplex{P}) where {P} - atan(imag(x), real(x)) end # a type specific hash function helps the type to 'just work' diff --git a/src/support/complex.jl b/src/support/complex.jl index 5ac0a1c1..ee79118f 100644 --- a/src/support/complex.jl +++ b/src/support/complex.jl @@ -1,76 +1,41 @@ -float(::Type{ArbFloat}) = ArbFloat -float(::Type{ArbReal}) = ArbFloat -float(::Type{ArbComplex}) = ArbFloat +float(::Type{T}) where {T<:ArbNumber1} = ArbFloat +float(::Type{T}) where {P,T<:ArbNumber1{P}} = ArbFloat{P} -float(::Type{ArbFloat{P}}) where {P} = ArbFloat{P} -float(::Type{ArbReal{P}}) where {P} = ArbFloat{P} -float(::Type{ArbComplex{P}}) where {P} = ArbFloat{P} +float(x::ArbNumber1{P}) where {P} = convert(ArbFloat{P}, x) # TODO senseless for Complex? -float(x::ArbFloat{P}) where {P} = x -float(x::ArbReal{P}) where {P} = ArbFloat{P}(x) -float(z::ArbComplex{P}) where {P} = ArbFloat{P}(real(z)) +real(::Type{ArbFloat}) = ArbFloat +real(::Type{T}) where {T<:ArbNumber1} = ArbReal +real(::Type{ArbFloat{P}}) where {P} = ArbFloat{P} +real(::Type{T}) where {P,T<:ArbNumber1{P}} = ArbReal{P} -real(::Type{ArbFloat}) = ArbFloat -real(::Type{ArbReal}) = ArbReal -real(::Type{ArbComplex}) = ArbReal - -real(::Type{ArbFloat{P}}) where {P} = ArbFloat{P} -real(::Type{ArbReal{P}}) where {P} = ArbReal{P} -real(::Type{ArbComplex{P}}) where {P} = ArbReal{P} - -real(x::ArbFloat{P}) where {P} = x -real(x::ArbReal{P}) where {P} = x +real(x::ArbNumber1) = x function real(x::ArbComplex{P}) where {P} z = ArbReal{P}() ccall(@libarb(acb_get_real), Cvoid, (Ref{ArbReal}, Ref{ArbComplex}), z, x) return z end +imag(::Type{T}) where {T<:ArbNumber1} = real(T) -imag(::Type{ArbFloat}) = ArbFloat -imag(::Type{ArbReal}) = ArbReal -imag(::Type{ArbComplex}) = ArbReal - -imag(::Type{ArbFloat{P}}) where {P} = ArbFloat{P} -imag(::Type{ArbReal{P}}) where {P} = ArbReal{P} -imag(::Type{ArbComplex{P}}) where {P} = ArbReal{P} - -imag(x::ArbFloat{P}) where {P} = x -imag(x::ArbReal{P}) where {P} = x +imag(::T) where {T<:ArbNumber1} = zero(T) function imag(x::ArbComplex{P}) where {P} z = ArbReal{P}() ccall(@libarb(acb_get_imag), Cvoid, (Ref{ArbReal}, Ref{ArbComplex}), z, x) return z end - -complex(::Type{ArbFloat}) = ArbComplex -complex(::Type{ArbReal}) = ArbComplex -complex(::Type{ArbComplex}) = ArbComplex - -complex(::Type{ArbFloat{P}}) where {P} = ArbComplex{P} -complex(::Type{ArbReal{P}}) where {P} = ArbComplex{P} -complex(::Type{ArbComplex{P}}) where {P} = ArbComplex{P} - -complex(x::ArbFloat{P}) where {P} = ArbComplex{P}(x) -complex(x::ArbReal{P}) where {P} = ArbComplex{P}(x) -complex(z::ArbComplex{P}) where {P} = z - +complex(::Type{T}) where {T<:ArbNumber1} = ArbComplex +complex(::Type{T}) where {P,T<:ArbNumber1{P}} = ArbComplex{P} +complex(x::T) where {T<:ArbNumber1} = convert(ArbComplex, x) # other parts -angle(x::ArbFloat{P}) where {P} = zero(ArbFloat{P}) -angle(x::ArbReal{P}) where {P} = zero(ArbReal{P}) +angle(x::T) where {T<:ArbNumber1} = signbit(x) ? T(pi) : zero(T) function angle(x::ArbComplex{P}) where {P} z = ArbReal{P}() - ccall(@libarb(acb_arg), Cvoid, (Ref{ArbReal}, Ref{ArbComplex}, Clong), z, x, P) + ccall(@libarb(acb_arg), Cvoid, (Ref{ArbReal}, Ref{ArbComplex}, Slong), z, x, P) return z end -magnitude(x::T) where {T<:Complex} = hypot(reim(x)...) - -magnitude(x::ArbFloat{P}) where {P} = x -magnitude(x::ArbReal{P}) where {P} = x -magnitude(x::ArbComplex{P}) where {P} = hypot(real(x), imag(x)) - +magnitude(x::Number) = abs(x) diff --git a/src/values/constructors.jl b/src/values/constructors.jl index 48c65dd0..21d25d61 100644 --- a/src/values/constructors.jl +++ b/src/values/constructors.jl @@ -127,39 +127,17 @@ ArbReal(x::T) where {T<:Real} = ArbReal{workingprecision(ArbReal)}(BigFloat(real ArbReal{P}(x::T) where {P,T<:Complex} = ArbReal{P}(BigFloat(x)) ArbReal(x::T) where {T<:Complex} = ArbReal{workingprecision(ArbReal)}(BigFloat(real(x))) -ArbComplex{P}(x::BigInt) where {P} = ArbComplex{P}(ArbReal{P}(x)) -ArbComplex{P}(x::BigFloat) where {P} = ArbComplex{P}(ArbReal{P}(x)) -ArbComplex{P}(x::T) where {P,T<:Real} = ArbComplex{P}(BigFloat(x)) -#ArbComplex(x::T) where {T<:Real} = ArbComplex{workingprecision(ArbComplex)}(BigFloat(x)) -ArbComplex{P}(x::T, y::T) where {P,T<:Real} = ArbComplex{P}(BigFloat(x), BigFloat(y)) -#ArbComplex(x::T, y::T) where {T<:Real} = ArbComplex{workingprecision(ArbComplex)}(BigFloat(x), BigFloat(y)) -ArbComplex{P}(x::T) where {P,T<:Complex} = ArbComplex{P}(BigFloat(real(x)), BigFloat(imag(x))) -ArbComplex(x::T) where {T<:Complex} = ArbComplex{workingprecision(ArbComplex)}(BigFloat(real(x)), BigFloat(imag(x))) - # retype ArbFloat(x::ArbReal{P}) where {P} = ArbFloat{P}(x) ArbFloat(x::ArbComplex{P}) where {P} = ArbFloat{P}(real(x)) ArbReal(x::ArbFloat{P}) where {P} = ArbReal{P}(x) ArbReal(x::ArbComplex{P}) where {P} = ArbReal{P}(real(x)) -# ArbComplex(x::ArbFloat{P}) where {P} = ArbComplex{P}(ArbReal{P}(x)) -# ArbComplex(x::ArbReal{P}) where {P} = ArbComplex{P}(x) ArbFloat{Q}(x::ArbReal{P}) where {P,Q} = ArbFloat{Q}(ArbReal{Q}(x)) ArbFloat{Q}(x::ArbComplex{P}) where {P,Q} = ArbFloat{Q}(ArbReal{Q}(real(x))) ArbReal{Q}(x::ArbFloat{P}) where {P,Q} = ArbReal{Q}(ArbFloat{Q}(x)) ArbReal{Q}(x::ArbComplex{P}) where {P,Q} = ArbReal{Q}(real(x)) -ArbComplex{Q}(x::ArbFloat{P}) where {P,Q} = ArbComplex{Q}(ArbReal{Q}(ArbFloat{Q}(x))) -ArbComplex{Q}(x::ArbReal{P}) where {P,Q} = ArbComplex{Q}(ArbReal{Q}(x)) - - -@inline convert(::Type{ArbFloat{Q}}, x::ArbReal{P}) where {P,Q} = ArbFloat{Q}(x) -@inline convert(::Type{ArbFloat{Q}}, x::ArbComplex{P}) where {P,Q} = ArbFloat{Q}(x) -@inline convert(::Type{ArbReal{Q}}, x::ArbFloat{P}) where {P,Q} = ArbReal{Q}(x) -@inline convert(::Type{ArbReal{Q}}, x::ArbComplex{P}) where {P,Q} = ArbReal{Q}(x) -@inline convert(::Type{ArbComplex{Q}}, x::ArbFloat{P}) where {P,Q} = ArbComplex{Q}(x) -@inline convert(::Type{ArbComplex{Q}}, x::ArbReal{P}) where {P,Q} = ArbComplex{Q}(x) - # change precision diff --git a/test/complex.jl b/test/complex.jl index a0d5b3ac..d98bc598 100644 --- a/test/complex.jl +++ b/test/complex.jl @@ -66,7 +66,7 @@ @test midpoint(ac) == ac @test 1e-38 <= abs(radius(ArbComplex(π))) <= 2e-38 - api = ArbComplex(pi, -pi) + api = ArbComplex(pi, -pi) # note `-pi is Float64` @test trunc(api) == ArbComplex(3 - 3im) @test floor(api) == ArbComplex(3 - 4im) @test ceil(api) == ArbComplex(4 - 3im) @@ -76,6 +76,7 @@ @test modf(ArbComplex(3.5, -2.5)) == ((0.5, 3.0), (-0.5, -2.0)) @test fmod(modf(api)...) == api @test angle(api) ≈ angle(Complex(api)) + @test api' ≈ ArbComplex(pi, pi) rtol=1e-16 @test flipsign(api, -1) == -api @test flipsign(api, 0.0) == api From c130f3a9534febfcfbfb2ab88602390f2a698718 Mon Sep 17 00:00:00 2001 From: KlausC Date: Fri, 3 Nov 2023 18:06:58 +0100 Subject: [PATCH 10/15] tested ArbFloat --- src/ArbNumerics.jl | 6 +- src/float/arith.jl | 6 +- src/float/prearith.jl | 5 +- src/intervals/eps_ulp.jl | 8 +-- src/libarb/ArbComplex.jl | 38 +++--------- src/libarb/ArbFloat.jl | 122 +++++++++++++++++++------------------ src/support/complex.jl | 4 -- src/values/constructors.jl | 17 ++---- test/complex.jl | 35 ++++++++++- test/misc.jl | 96 ++++++++++++++++++++++++++++- 10 files changed, 218 insertions(+), 119 deletions(-) diff --git a/src/ArbNumerics.jl b/src/ArbNumerics.jl index 859623d7..ddc27675 100644 --- a/src/ArbNumerics.jl +++ b/src/ArbNumerics.jl @@ -50,7 +50,7 @@ export ArbNumber, addmul, submul, mulsub, square, cube, rsqrt, rcbrt, pow, root, loghypot, risingfactorial, doublefactorial, - tanpi, cotpi, + tanpi, cotpi, swap!, # special functions agm1, agm, @@ -153,6 +153,10 @@ using Libdl using Random using Random: SamplerType, SamplerTrivial, CloseOpen01 +const Slong = Int # to accomodate windows +const USlong = unsigned(Slong) +const ArbInts = Union{Int,Int32,Int16,Int8} # Int may be Int32 + include("support/abstractions.jl") include("support/matrices.jl") diff --git a/src/float/arith.jl b/src/float/arith.jl index 9cd95a85..ea23e8ff 100644 --- a/src/float/arith.jl +++ b/src/float/arith.jl @@ -123,11 +123,11 @@ end # interact with Bool @inline ArbFloat{P}(b::Bool) where {P} = b ? one(ArbFloat{P}) : zero(ArbFloat{P}) -@inline Base.Bool(x::ArbFloat{P}) where {P} = iszero(x) ? false : (isone(x) ? true : throw(InexactError(""))) +@inline Base.Bool(x::ArbFloat{P}) where {P} = iszero(x) ? false : (isone(x) ? true : throw(InexactError(:Bool, Bool, x))) @inline ArbReal{P}(b::Bool) where {P} = b ? one(ArbReal{P}) : zero(ArbReal{P}) -@inline Base.Bool(x::ArbReal{P}) where {P} = iszero(x) ? false : (isone(x) ? true : throw(InexactError(""))) +@inline Base.Bool(x::ArbReal{P}) where {P} = iszero(x) ? false : (isone(x) ? true : throw(InexactError(:Bool, Bool, x))) @inline ArbComplex{P}(b::Bool) where {P} = b ? one(ArbComplex{P}) : zero(ArbComplex{P}) -@inline Base.Bool(x::ArbComplex{P}) where {P} = iszero(x) ? false : (isone(x) ? true : throw(InexactError(""))) +@inline Base.Bool(x::ArbComplex{P}) where {P} = iszero(x) ? false : (isone(x) ? true : throw(InexactError(:Bool, Bool, x))) (+)(x::ArbFloat{P}, b::Bool) where {P} = b ? one(ArbFloat{P}) + x : zero(ArbFloat{P}) (+)(b::Bool, x::ArbFloat{P}) where {P} = b ? x + one(ArbFloat{P}) : zero(ArbFloat{P}) diff --git a/src/float/prearith.jl b/src/float/prearith.jl index 3ef0405a..4075bc9b 100644 --- a/src/float/prearith.jl +++ b/src/float/prearith.jl @@ -4,12 +4,11 @@ signbit(x::ArbFloat{P}) where {P} = isfinite(x) ? sign_bit(x) : isneginf(x) signbit(x::ArbReal{P}) where {P} = isfinite(x) ? sign_bit(x) : isneginf(x) -signbit(x::ArbComplex{P}, ::Type{RealPart}) where {P} = signbit(real(x)) -signbit(x::ArbComplex{P}, ::Type{ImagPart}) where {P} = signbit(imag(x)) +signbit(x::ArbComplex{P}, ::Type{RealPart}) where {P} = sign_bit(x, RealPart) +signbit(x::ArbComplex{P}, ::Type{ImagPart}) where {P} = sign_bit(x, ImagPart) signbit(x::ArbComplex{P}) where {P} = signbit(x, RealPart) signbits(x::ArbComplex{P}) where {P} = signbit(x, RealPart), signbit(x, ImagPart) - function sign(x::ArbFloat{P}) where {P} thesgn = ccall(@libarb(arf_sgn), Cint, (Ref{ArbFloat},), x) return ArbFloat{P}(thesgn) diff --git a/src/intervals/eps_ulp.jl b/src/intervals/eps_ulp.jl index b9f17b9a..977bd4a5 100644 --- a/src/intervals/eps_ulp.jl +++ b/src/intervals/eps_ulp.jl @@ -108,16 +108,16 @@ end void arf_mag_set_ulp(mag_t res, const arf_t x, slong prec) Sets res to the magnitude of the unit in the last place (ulp) of x at precision prec. =# - + function ulp(x::ArbFloat{P}) where {P} if isnan(x) || isinf(x) return ArbFloat{P}(NaN) - end + end if !iszero(x) w = Mag() # Sets z to the magnitude of the unit in the last place (ulp) of x at precision P. - ccall(@libarb(arf_mag_set_ulp), Cvoid, (Ref{Mag}, Ref{ArbFloat}, Clong), w, x, P) + ccall(@libarb(arf_mag_set_ulp), Cvoid, (Ref{Mag}, Ref{ArbFloat}, Slong), w, x, P) z = ArbFloat{P}(w) else z = ArbFloat{P}(2.0)^(-workingprecision(x)) @@ -128,7 +128,7 @@ end function eps(x::ArbFloat{P}) where {P} if isnan(x) || isinf(x) return ArbFloat{P}(NaN) - end + end if !iszero(x) w = Mag() # Sets z to the magnitude of the unit in the last place (ulp) of x at precision P. diff --git a/src/libarb/ArbComplex.jl b/src/libarb/ArbComplex.jl index 96188c6c..871a7917 100644 --- a/src/libarb/ArbComplex.jl +++ b/src/libarb/ArbComplex.jl @@ -61,11 +61,8 @@ ArbComplex{P}(x::T) where {P,T<:Number} = ArbComplex{P}(real(x), imag(x)) const PtrToArbComplex = Ptr{ArbComplex} # acb_ptr const PtrToPtrToArbComplex = Ptr{Ptr{ArbComplex}} # acb_ptr* -const Slong = Int # to accomodate windows const ArbFloatReal{P} = Union{ArbFloat{P},ArbReal{P}} const ArbNumber1{P} = Union{ArbFloatReal{P},ArbComplex{P}} -const ArbInts = Union{Int,Int32,Int16,Int8} # Int may be Int32 -const ArbRealTypes = Real # Union{Integer,AbstractFloat,ArbReal} # finalizer: clear_acb(x::ArbComplex{P}) where {P} = ccall(@libarb(acb_clear), Cvoid, (Ref{ArbComplex},), x) @@ -81,15 +78,7 @@ ArbComplex{P}(::Real, ::Missing) where {P} = missing ArbComplex(x::T; kw...) where {T<:Number} = ArbComplex(real(x), imag(x); kw...) function ArbComplex(x::T, y::Real=0; bits::Int=0, digits::Int=0, base::Int=iszero(bits) ? 10 : 2) where {T<:Real} - bits > 0 && bits < MINIMUM_PRECISION_BASE2 && throw(DomainError("bit precision $bits < $MINIMUM_PRECISION_BASE2")) - digits > 0 && digits < MINIMUM_PRECISION_BASE10 && throw(DomainError("digit precision $digits < $MINIMUM_PRECISION_BASE10")) - if base === 10 - bits = digits > 0 ? bits4digits(digits) + extrabits() : (bits > 0 ? bits + extrabits() : DEFAULT_PRECISION.x) - elseif base === 2 - bits = bits > 0 ? bits + extrabits() : (digits > 0 ? digits + extrabits() : DEFAULT_PRECISION.x) - else - throw(ErrorException("base expects 2 or 10")) - end + bits = get_bits(bits, digits, base) iszero(y) ? ArbComplex{bits}(x) : ArbComplex{bits}(x, y) end @@ -99,8 +88,6 @@ ArbComplex(x::ArbFloatReal{P}, y::ArbFloatReal{S}) where {P,S} = ArbComplex{min( ArbComplex(x::ArbFloatReal{P}, y::Real) where {P} = ArbComplex{P}(x, y) ArbComplex(x::Real, y::ArbFloatReal{P}) where {P} = ArbComplex{P}(x, y) -@inline sign_bit(x::ArbComplex{P}) where {P} = isodd(x.real_mid_size) -@inline sign_bits(x::ArbComplex{P}) where {P} = isodd(x.real_mid_size), isodd(x.imag_mid_size) @inline sign_bit(x::ArbComplex{P}, ::Type{RealPart}) where {P} = isodd(x.real_mid_size) @inline sign_bit(x::ArbComplex{P}, ::Type{ImagPart}) where {P} = isodd(x.imag_mid_size) @@ -136,7 +123,6 @@ function ArbComplex{P}(x::ArbInts, y::ArbInts) where {P} return z end - function copy(x::ArbComplex{P}) where {P} z = ArbComplex{P}() ccall(@libarb(acb_set), Cvoid, (Ref{ArbComplex}, Ref{ArbComplex}), z, x) @@ -145,22 +131,16 @@ end deepcopy(x::ArbComplex{P}) where {P} = copy(x) -Base.Complex(x::ArbFloat{P}) where {P} = ArbComplex(x) -Base.Complex(x::ArbReal{P}) where {P} = ArbComplex(x) -Base.Complex(re::ArbFloat{P}, im::ArbFloat{P}) where {P} = ArbComplex(re, im) -Base.Complex(re::ArbReal{P}, im::ArbReal{P}) where {P} = ArbComplex(re, im) - -Base.complex(re::ArbFloat{P}, im::ArbFloat{P}) where {P} = ArbComplex(re, im) -Base.complex(re::ArbReal{P}, im::ArbReal{P}) where {P} = ArbComplex(re, im) +Complex(x::ArbNumber1) = convert(ArbComplex, x) +Complex(re::ArbFloatReal, im::ArbFloatReal) = ArbComplex(re, im) -Base.Complex(re::ArbFloat{P}, im::ArbReal{P}) where {P} = ArbComplex(ArbReal{P}(re), im) -Base.Complex(re::ArbReal{P}, im::ArbFloat{P}) where {P} = ArbComplex(re, ArbReal{P}(im)) -Base.complex(re::ArbFloat{P}, im::ArbReal{P}) where {P} = ArbComplex(ArbReal{P}(re), im) -Base.complex(re::ArbReal{P}, im::ArbFloat{P}) where {P} = ArbComplex(re, ArbReal{P}(im)) +complex(re::ArbNumber1) = convert(ArbComplex, re) +complex(re::ArbFloatReal, im::ArbFloatReal) = ArbComplex(re, im) -Base.Complex{ArbFloat{P}}(x::ArbComplex{P}) where {P} = x -Base.Complex{ArbReal{P}}(x::ArbComplex{P}) where {P} = x +complex(::Type{T}) where {T<:ArbNumber1} = ArbComplex +complex(::Type{T}) where {P,T<:ArbNumber1{P}} = ArbComplex{P} +Base.Complex{T}(x::ArbComplex{P}) where {P,T<:ArbFloatReal{P}} = x function ArbComplex{P}(x::T) where {P,T<:Integer} y = ArbReal{P}(x) @@ -169,7 +149,7 @@ function ArbComplex{P}(x::T) where {P,T<:Integer} end -function ArbComplex{P}(x::T1, y::T2) where {P,T1<:ArbRealTypes,T2<:ArbRealTypes} +function ArbComplex{P}(x::T1, y::T2) where {P,T1<:Real,T2<:Real} x1 = convert(ArbReal{P}, x) y1 = convert(ArbReal{P}, y) z = ArbComplex{P}(x1, y1) diff --git a/src/libarb/ArbFloat.jl b/src/libarb/ArbFloat.jl index 708a0d3f..7cdba128 100644 --- a/src/libarb/ArbFloat.jl +++ b/src/libarb/ArbFloat.jl @@ -14,7 +14,8 @@ The precision of an ArbFloat must be >= 24 bits or >= 8 digits. Negative zero, unsigned Infinity, NaNs with distinct payloads are not supported. See also: [`ArbReal`](@ref), [`ArbComplex`](@ref), [`workingprecision`](@ref), [`displayprecision`](@ref) -""" ArbFloat +""" +ArbFloat mutable struct ArbFloat{P} <: AbstractFloat # P is the precision in bits exp::Int # fmpz exponent of 2 (2^exp) @@ -23,7 +24,7 @@ mutable struct ArbFloat{P} <: AbstractFloat # P is the precision in bits d2::UInt # (d1, d2) the final part indicating the significand, or 0 function ArbFloat{P}() where {P} - z = new{P}(0,0,0,0) + z = new{P}(0, 0, 0, 0) ccall(@libarb(arf_init), Cvoid, (Ref{ArbFloat},), z) finalizer(arf_clear, z) return z @@ -39,26 +40,37 @@ arf_clear(x::ArbFloat{P}) where {P} = ccall(@libarb(arf_clear), Cvoid, (Ref{ArbF ArbFloat{P}(x::ArbFloat{P}) where {P} = x ArbFloat(x::ArbFloat{P}) where {P} = x -ArbFloat{P}(x::Missing) where {P} = missing -ArbFloat(x::Missing) = missing +ArbFloat{P}(::Missing) where {P} = missing +ArbFloat(::Missing) = missing @inline sign_bit(x::ArbFloat{P}) where {P} = isodd(x.size) - -function ArbFloat(x::T; bits::Int=0, digits::Int=0, base::Int=iszero(bits) ? 10 : 2) where {T<:Number} +function ArbFloat(x::Union{AbstractFloat,AbstractIrrational,Integer}; kw...) + _ArbFloat(x; kw...) +end +function ArbFloat(x::Union{Rational}; kw...) + _ArbFloat(x; kw...) +end +function _ArbFloat(x::Real; bits::Int=0, digits::Int=0, base::Int=iszero(bits) ? 10 : 2) + bits = get_bits(bits, digits, base) + ArbFloat{bits}(x) +end +function get_bits(bits, digits, base) bits > 0 && bits < MINIMUM_PRECISION_BASE2 && throw(DomainError("bit precision $bits < $MINIMUM_PRECISION_BASE2")) digits > 0 && digits < MINIMUM_PRECISION_BASE10 && throw(DomainError("digit precision $digits < $MINIMUM_PRECISION_BASE10")) if base === 10 - bits = digits > 0 ? bits4digits(digits)+extrabits() : (bits > 0 ? bits+extrabits() : DEFAULT_PRECISION.x) + bits = digits > 0 ? bits4digits(digits) + extrabits() : (bits > 0 ? bits + extrabits() : DEFAULT_PRECISION.x) elseif base === 2 - bits = bits > 0 ? bits+extrabits() : (digits > 0 ? digits+extrabits() : DEFAULT_PRECISION.x) + bits = bits > 0 ? bits + extrabits() : (digits > 0 ? digits + extrabits() : DEFAULT_PRECISION.x) else throw(ErrorException("base expects 2 or 10")) end - ArbFloat{bits}(x) + bits end -swap(x::ArbFloat{P}, y::ArbFloat{P}) where {P} = ccall(@libarb(arf_swap), Cvoid, (Ref{ArbFloat}, Ref{ArbFloat}), x, y) +function swap!(x::ArbFloat{P}, y::ArbFloat{P}) where {P} + ccall(@libarb(arf_swap), Cvoid, (Ref{ArbFloat}, Ref{ArbFloat}), x, y) +end function copy(x::ArbFloat{P}) where {P} z = ArbFloat{P}() @@ -77,78 +89,77 @@ copy(x::ArbFloat{P}, roundingmode::RoundingMode) where {P} = copy(x, P, rounding copy(x::ArbFloat{P}, bitprecision::Int) where {P} = copy(x, bitprecision, RoundNearest) -function ArbFloat{P}(x::Int32) where {P} +function ArbFloat{P}(x::ArbInts) where {P} z = ArbFloat{P}() - ccall(@libarb(arf_set_si), Cvoid, (Ref{ArbFloat}, Clong), z, x) + ccall(@libarb(arf_set_si), Cvoid, (Ref{ArbFloat}, Slong), z, x) return z end -ArbFloat{P}(x::T) where {P, T<:Union{Int8, Int16}} = ArbFloat{P}(Int32(x)) -ArbFloat{P}(x::T) where {P, T<:Union{Int64, Int128}} = ArbFloat{P}(BigInt(x)) -function ArbFloat{P}(x::UInt32) where {P} +function ArbFloat{P}(x::USlong) where {P} z = ArbFloat{P}() - ccall(@libarb(arf_set_ui), Cvoid, (Ref{ArbFloat}, Culong), z, x) + ccall(@libarb(arf_set_ui), Cvoid, (Ref{ArbFloat}, USlong), z, x) return z end -ArbFloat{P}(x::T) where {P, T<:Union{UInt8, UInt16}} = ArbFloat{P}(UInt32(x)) -ArbFloat{P}(x::T) where {P, T<:Union{UInt64, UInt128}} = ArbFloat{P}(BigInt(x)) -function ArbFloat{P}(x::Float64) where {P} +ArbFloat{P}(x::T) where {P,T<:Integer} = ArbFloat{P}(BigFloat(x, precision=P + 2)) + +function ArbFloat{P}(x::Base.IEEEFloat) where {P} z = ArbFloat{P}() ccall(@libarb(arf_set_d), Cvoid, (Ref{ArbFloat}, Cdouble), z, x) return z end -ArbFloat{P}(x::T) where {P, T<:Union{Float16, Float32}} = ArbFloat{P}(Float64(x)) function ArbFloat{P}(x::BigFloat) where {P} z = ArbFloat{P}() ccall(@libarb(arf_set_mpfr), Cvoid, (Ref{ArbFloat}, Ref{BigFloat}), z, x) return z end -ArbFloat{P}(x::BigInt) where {P} = ArbFloat{P}(BigFloat(x)) -ArbFloat{P}(x::Rational{T}) where {P, T<:Signed} = ArbFloat{P}(BigFloat(x)) + +ArbFloat{P}(x::Rational{T}) where {P,T<:Signed} = ArbFloat{P}(BigFloat(x, precision=P + 2)) function ArbFloat{P}(x::Irrational{S}) where {P,S} - prec = precision(BigFloat) - newprec = max(prec, P + 32) - setprecision(BigFloat, newprec) - y = BigFloat(x) + y = BigFloat(x, precision=P + 32) z = ArbFloat{P}(y) - setprecision(BigFloat, prec) return z end -function Int64(x::ArbFloat{P}, roundingmode::RoundingMode) where {P} +function (::Type{T})(x::ArbFloat, roundingmode::RoundingMode=RoundNearest) where {T<:ArbInts} rounding = match_rounding_mode(roundingmode) - z = ccall(@libarb(arf_get_si), Clong, (Ref{ArbFloat}, Cint), x, rounding) - return z + if !( typemin(T) <= x <= typemax(T)) + throw(InexactError(nameof(T), T, x)) + # attention: ccall segfaults, if result would become too large for Slong + end + z = ccall(@libarb(arf_get_si), Slong, (Ref{ArbFloat}, Cint), x, rounding) + return convert(T, z) +end + +function (::Type{T})(x::ArbFloat{P}, roundingmode::RoundingMode=RoundNearest) where {P,T<:Integer} + y = round(x, roundingmode) + z = BigFloat(y; precision=P + 2) + return convert(T, z) end -Int32(x::ArbFloat{P}, roundingmode::RoundingMode) where {P} = Int32(Int64(x), roundingmode) -Int16(x::ArbFloat{P}, roundingmode::RoundingMode) where {P} = Int16(Int64(x), roundingmode) """ BigFloat(::ArbFloat; [precision=workingprecision(x), roundingmode=RoundNearest]) Construct a `BigFloat`from an `ArbFloat`. """ -function BigFloat(x::ArbFloat{P}; precision::Int=workingprecision(x), roundingmode::RoundingMode=RoundNearest) where {P} +function BigFloat(x::ArbFloat; precision::Int=workingprecision(x), roundingmode::RoundingMode=RoundNearest) rounding = match_rounding_mode(roundingmode) z = BigFloat(0; precision) roundingdir = ccall(@libarb(arf_get_mpfr), Cint, (Ref{BigFloat}, Ref{ArbFloat}, Cint), z, x, rounding) return z end -BigInt(x::ArbFloat{P}) where {P} = BigInt(trunc(BigFloat(x))) - function Base.Integer(x::ArbFloat{P}) where {P} if isinteger(x) - abs(x) <= typemax(Int64) ? Int64(x) : BigInt(x) + typemin(Int64) <= abs(x) <= typemax(Int64) ? Int64(x, RoundNearest) : BigInt(x) else - throw(InexactError("$x")) + throw(InexactError(:Integer, Integer, x)) end end -for (F,A) in ((:floor, :arf_floor), (:ceil, :arf_ceil)) +for (F, A) in ((:floor, :arf_floor), (:ceil, :arf_ceil)) @eval begin function $F(x::ArbFloat{P}) where {P} z = ArbFloat{P}() @@ -156,16 +167,16 @@ for (F,A) in ((:floor, :arf_floor), (:ceil, :arf_ceil)) return z end $F(::Type{Bool}, x::ArbFloat{P}) where {P} = Bool($F(x)) - $F(::Type{T}, x::ArbFloat{P}) where {P, T<:Integer} = T($F(x)) + $F(::Type{T}, x::ArbFloat{P}) where {P,T<:Integer} = T($F(x)) end end trunc(x::ArbFloat{P}) where {P} = signbit(x) ? ceil(x) : floor(x) trunc(::Type{Bool}, x::ArbFloat{P}) where {P} = Bool(trunc(x)) -trunc(::Type{T}, x::ArbFloat{P}) where {P, T<:Integer} = T(trunc(x)) +trunc(::Type{T}, x::ArbFloat{P}) where {P,T<:Integer} = T(trunc(x)) midpoint(x::ArbFloat{P}) where {P} = x -radius(x::ArbFloat{P}) where {P} = zero(ArbFloat{P}) +radius(::ArbFloat{P}) where {P} = zero(ArbFloat{P}) function modf(x::ArbFloat{P}) where {P} ipart = trunc(x) @@ -173,42 +184,35 @@ function modf(x::ArbFloat{P}) where {P} return (fpart, ipart) end -fmod(fpartipart::Tuple{ArbFloat{P}, ArbFloat{P}}) where {P} = - fpartipart[1] + fpartipart[2] -fmod(fpart::ArbFloat{P}, ipart::ArbFloat{P}) where {P} = - fpart + ipart +fmod(fpartipart::Tuple{ArbFloat{P},ArbFloat{P}}) where {P} = fpartipart[1] + fpartipart[2] +fmod(fpart::ArbFloat{P}, ipart::ArbFloat{P}) where {P} = fpart + ipart -div(x::ArbFloat{P}, y::ArbFloat{P}) where {P} = - trunc(x / y) +div(x::ArbFloat{P}, y::ArbFloat{P}, rm::RoundingMode) where {P} = round(x / y, rm) -rem(x::ArbFloat{P}, y::ArbFloat{P}) where {P} = - x - (div(x,y) * y) +rem(x::ArbFloat{P}, y::ArbFloat{P}) where {P} = x - (div(x, y) * y) function divrem(x::ArbFloat{P}, y::ArbFloat{P}) where {P} - dv = div(x,y) + dv = div(x, y) rm = x - (dv * y) return (dv, rm) end -fld(x::ArbFloat{P}, y::ArbFloat{P}) where {P} = - floor(x / y) +# fld(x::ArbFloat{P}, y::ArbFloat{P}) where {P} = floor(x / y) -mod(x::ArbFloat{P}, y::ArbFloat{P}) where {P} = - x - (fld(x,y) * y) +mod(x::ArbFloat{P}, y::ArbFloat{P}) where {P} = x - (fld(x, y) * y) function fldmod(x::ArbFloat{P}, y::ArbFloat{P}) where {P} - fd = fld(x,y) + fd = fld(x, y) md = x - (fd * y) return (fd, md) end -cld(x::ArbFloat{P}, y::ArbFloat{P}) where {P} = - ceil(x / y) +# cld(x::ArbFloat{P}, y::ArbFloat{P}) where {P} = ceil(x / y) # a type specific hash function helps the type to 'just work' const hash_arbfloat_lo = (UInt === UInt64) ? 0x37e642589da3416a : 0x5d46a6b4 const hash_0_arbfloat_lo = hash(zero(UInt), hash_arbfloat_lo) hash(z::ArbFloat{P}, h::UInt) where {P} = - hash(reinterpret(UInt,z.d1) ⊻ z.exp, - (h ⊻ hash(z.d2 ⊻ (~reinterpret(UInt,P)), hash_arbfloat_lo) ⊻ hash_0_arbfloat_lo)) + hash(reinterpret(UInt, z.d1) ⊻ z.exp, + (h ⊻ hash(z.d2 ⊻ (~reinterpret(UInt, P)), hash_arbfloat_lo) ⊻ hash_0_arbfloat_lo)) diff --git a/src/support/complex.jl b/src/support/complex.jl index ee79118f..6a96865b 100644 --- a/src/support/complex.jl +++ b/src/support/complex.jl @@ -25,10 +25,6 @@ function imag(x::ArbComplex{P}) where {P} return z end -complex(::Type{T}) where {T<:ArbNumber1} = ArbComplex -complex(::Type{T}) where {P,T<:ArbNumber1{P}} = ArbComplex{P} -complex(x::T) where {T<:ArbNumber1} = convert(ArbComplex, x) - # other parts angle(x::T) where {T<:ArbNumber1} = signbit(x) ? T(pi) : zero(T) diff --git a/src/values/constructors.jl b/src/values/constructors.jl index 21d25d61..30984ae3 100644 --- a/src/values/constructors.jl +++ b/src/values/constructors.jl @@ -69,13 +69,6 @@ Complex(x::ArbComplex{P}) where {P} = Complex{Float64}(x, RoundNearest) for I in (:Int8, :Int16, :Int32, :Int64, :Int128) @eval begin - function $I(x::ArbFloat{P}) where {P} - !isinteger(x) && throw(InexactError("$(x) is not an integer")) - bi = BigInt(BigFloat(x)) - !(typemin($I) <= bi <= typemax($I)) && throw(InexactError("$(x)")) - return $I(bi) - end - function $I(x::ArbReal{P}) where {P} (!isexact(x) | !isinteger(x)) && throw(InexactError("$(x) is not an integer")) bi = BigInt(BigFloat(x)) @@ -95,10 +88,10 @@ end # Rational -ArbFloat(x::T) where {S, T<:Rational{S}} = ArbFloat(x.num)/ArbFloat(x.den) -ArbReal(x::T) where {S, T<:Rational{S}} = ArbReal(x.num)/ArbReal(x.den) -ArbComplex(x::T) where {S, T<:Rational{S}} = ArbComplex(ArbReal(x), ArbReal(0)) -ArbComplex(x::T, y::T) where {S, T<:Rational{S}} = ArbComplex(ArbReal(x), ArbReal(y)) +1#ArbFloat(x::T) where {S, T<:Rational{S}} = ArbFloat(x.num)/ArbFloat(x.den) +#ArbReal(x::T) where {S, T<:Rational{S}} = ArbReal(x.num)/ArbReal(x.den) +#ArbComplex(x::T) where {S, T<:Rational{S}} = ArbComplex(ArbReal(x), ArbReal(0)) +#ArbComplex(x::Rational, y::Rational) = ArbComplex(ArbReal(x), ArbReal(y)) # Irrational @@ -118,7 +111,7 @@ end # fallback # ArbFloat{P}(x::T) where {P,T<:Real} = ArbFloat{P}(BigFloat(x)) -ArbFloat(x::T) where {T<:Real} = ArbFloat{workingprecision(ArbFloat)}(BigFloat(x)) +# ArbFloat(x::T) where {T<:Real} = ArbFloat{workingprecision(ArbFloat)}(BigFloat(x)) ArbFloat{P}(x::T) where {P,T<:Complex} = ArbFloat{P}(BigFloat(real(x))) ArbFloat(x::T) where {T<:Complex} = ArbFloat{workingprecision(ArbFloat)}(BigFloat(real(x))) diff --git a/test/complex.jl b/test/complex.jl index d98bc598..e73bf183 100644 --- a/test/complex.jl +++ b/test/complex.jl @@ -35,10 +35,20 @@ @test_throws ErrorException ArbComplex(5.0, base=3) @test_throws ErrorException ArbComplex(5.0, 1.0, base=3) @test ArbComplex{128}(1) == 1 + ac = ArbComplex{64}(1.0 + 2im) rac = real(ac) fac = float(rac) @test abs(ac) ≈ sqrt(abs2(fac) + abs2(imag(ac))) + @test angle(ac) ≈ angle(Complex(ac)) + @test magnitude(ac) ≈ sqrt(5) + @test angle(fac) == 0 + @test angle(-rac) ≈ pi + @test midpoint(ac) == ac + @test radius(ac) ≈ 0 + @test sign(fac) == 1 + @test sign(rac) == 1 + @test ArbComplex{64}(ac) === ac @test ArbComplex{64}(fac) == ArbComplex(fac) @test ArbComplex(5, digits=53) == 5 @@ -46,10 +56,14 @@ @test ArbComplex{53}(1 + 2im) == 1 + 2im @test signbit(ac) == false @test signbits(ac) == (false, false) + @test sign(ac) == 1.0 # TODO that is incompatible with sign(::Complex); maybe rename + @test signs(ac) == (1, 1) + @test_broken sign(ac) ≈ ac / abs(ac) @test ArbComplex(fac) == ArbComplex(rac) @test copy(ac) == deepcopy(ac) @test ArbComplex(rac, fac) == Complex(rac, fac) == complex(rac, fac) @test ArbComplex(rac) == Complex(rac) == complex(rac) + @test Complex{typeof(fac)}(ac) === ac @test typeof(ArbComplex(rac, Int128(1))) == ArbComplex{64} @test typeof(ArbComplex(fac, UInt8(1))) == ArbComplex{64} @@ -64,8 +78,8 @@ @test ArbComplex{256}(1.0 + im) == ArbComplex(1 + im, bits=256) @test ArbComplex{256}(1.0 + im) == ArbComplex(1, 1, bits=256) - @test midpoint(ac) == ac @test 1e-38 <= abs(radius(ArbComplex(π))) <= 2e-38 + api = ArbComplex(pi, -pi) # note `-pi is Float64` @test trunc(api) == ArbComplex(3 - 3im) @test floor(api) == ArbComplex(3 - 4im) @@ -73,9 +87,10 @@ @test trunc(Int, api) == (3, -3) @test floor(Int8, api) == (3, -4) @test ceil(Int16, api) == (4, -3) + @test modf(ArbComplex(3.5, -2.5)) == ((0.5, 3.0), (-0.5, -2.0)) @test fmod(modf(api)...) == api - @test angle(api) ≈ angle(Complex(api)) + @test api' ≈ ArbComplex(pi, pi) rtol=1e-16 @test flipsign(api, -1) == -api @@ -99,6 +114,22 @@ end + TYPE_TESTS = [ + (complex, ArbFloat, ArbComplex), + (complex, ArbReal{42}, ArbComplex{42}), + (complex, ArbComplex{42}, ArbComplex{42}), + (float, ArbReal, ArbFloat), + (float, ArbFloat{42}, ArbFloat{42}), + (real, ArbFloat, ArbFloat), + (real, ArbComplex, ArbReal), + (real, ArbFloat{42}, ArbFloat{42}), + (real, ArbComplex{42}, ArbReal{42}), + (imag, ArbComplex, ArbReal) + ] + @testset "$f(::Type{$T})" for (f, T, R) in TYPE_TESTS + @test f(T) == R + end + INF = typemax(Int) sprec(a::ArbNumber) = workingprecision(a) sprec(::Number) = workingprecision(ArbFloat) diff --git a/test/misc.jl b/test/misc.jl index 9a53dfcc..533919b4 100644 --- a/test/misc.jl +++ b/test/misc.jl @@ -1,3 +1,95 @@ -@testset "rounding mode" begin - @test ArbNumerics.match_rounding_mode(RoundNearest) == 4 +@testset "misc.jl" begin + @testset "rounding mode" begin + @test ArbNumerics.match_rounding_mode(RoundNearest) == 4 + end + + @testset "ArbNumerics.Mag" begin + @test signbit(ArbNumerics.Mag(0)) == false + @test ArbNumerics.Mag(-1e-5) == ArbNumerics.Mag(1e-5) + @test_throws InexactError ArbNumerics.Mag(-1) + end + + @testset "$T($v)" for T in (ArbFloat, ArbFloat{64}, ArbReal, ArbReal{64}), + (v, c, f, t) in [(0.6, 1, 0, 0), (-0.6, 0, -1, 0)] + + x = T(v) + @test ceil(x) == float(c) + 0 <= c <= 1 && @test ceil(Bool, x) == (c == 1) + @test ceil(Int, x) === c + + @test floor(x) == float(f) + if 0 <= f <= 1 + @test floor(Bool, x) == (f == 1) + else + @test_throws InexactError floor(Bool, x) + end + @test floor(Int, x) === f + + @test trunc(x) == float(t) + 0 <= t <= 1 && @test trunc(Bool, x) == (t == 1) + @test trunc(Int, x) === t + end + + @testset "ArbFloat" begin + a = ArbFloat(10) + b = ArbFloat(20) + swap!(a, b) + @test (a, b) == (20, 10) + + a = ArbFloat(pi) + @test ArbFloat(a) === a + @test copy(a) == a + @test copy(a) !== a + @test copy(a, 10, RoundDown) == 804 / 256 + @test copy(a, 10, RoundUp) == 805 / 256 + @test copy(a, RoundUp) > a + @test a > copy(a, RoundDown) + @test copy(a, 10) == copy(a, 10, RoundDown) + + @test ArbFloat(typemax(UInt)) == typemax(UInt) + @test ArbFloat(typemax(Int128)) == typemax(Int128) + @test ArbFloat(1 // 3) ≈ ArbFloat(1) / ArbFloat(3) + end + + @testset "convert ArbFloat to Integer" begin + v = typemax(Int32) + 1.5 + w = round(v) + a = ArbFloat(v) + @test_throws InexactError Int32(a) + @test Int64(a) == w + @test Int128(a) == w + @test BigInt(a) == w + + @test_throws InexactError Integer(a) # TODO inconsistent with Int64, BigInt, ... + v = typemin(Int32) - 1 + a = ArbFloat(v) + x = Integer(a) + @test typeof(x) == Int64 + @test x == v + v = Int128(typemax(Int64)) + 1 + a = ArbFloat(v) + x = Integer(a) + @test typeof(x) == BigInt + @test x == v + end + + @testset "functions with ArbFloat" begin + a = ArbFloat(ℯ) + @test midpoint(a) === a + @test iszero(radius(a)) + + @test modf(a) == (a - 2, 2) + @test fmod(modf(a)) ≈ a + @test fmod(modf(a)...) ≈ a + + b = ArbFloat(1) + @test div(-a, b) == -2 + @test rem(-a, b) ≈ -a + 2 + @test fld(-a, b) == -3 + @test mod(-a, b) ≈ -a + 3 + @test divrem(a, b) == (div(a, b), rem(a, b)) + @test fldmod(a, b) == (fld(a, b), mod(a, b)) + + @test hash(a) isa UInt + end end From fc08471e4bd56598ddd7cf30eefe7bde4210774a Mon Sep 17 00:00:00 2001 From: KlausC Date: Sun, 5 Nov 2023 18:13:14 +0100 Subject: [PATCH 11/15] ARBReal tested - except string functions --- src/intervals/eps_ulp.jl | 2 +- src/libarb/ArbComplex.jl | 1 + src/libarb/ArbFloat.jl | 2 + src/libarb/ArbReal.jl | 64 +++++++++++-------------------- src/libarb/string.jl | 33 ++++++---------- src/support/minprec.jl | 18 +++++---- src/values/constructors.jl | 23 +++++++----- src/values/rounding.jl | 6 +++ test/complex.jl | 2 + test/misc.jl | 77 ++++++++++++++++++++++++++++++++++++-- 10 files changed, 142 insertions(+), 86 deletions(-) diff --git a/src/intervals/eps_ulp.jl b/src/intervals/eps_ulp.jl index 977bd4a5..01ab6baf 100644 --- a/src/intervals/eps_ulp.jl +++ b/src/intervals/eps_ulp.jl @@ -86,7 +86,7 @@ trim_bits(x::Mag) = x function trim_bits(x::ArbFloat{P}, roundingmode::RoundingMode=RoundFromZero) where {P} nbits = significand_bits(x) - nbits == P && return x + (nbits == P || nbits == 0 ) && return x rounding = match_rounding_mode(roundingmode) z = ArbFloat{nbits} res = ccall(@libarb(arf_set_round), Cint, (Ref{ArbFloat}, Ref{ArbFloat}, Clong, Cint), z, x, nbits, rounding) diff --git a/src/libarb/ArbComplex.jl b/src/libarb/ArbComplex.jl index 871a7917..bb1f5672 100644 --- a/src/libarb/ArbComplex.jl +++ b/src/libarb/ArbComplex.jl @@ -48,6 +48,7 @@ mutable struct ArbComplex{P} <: Number # P is the precision in bits function ArbComplex{P}() where {P} + minprec(P, ArbComplex) z = new{P}(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) ccall(@libarb(acb_init), Cvoid, (Ref{ArbComplex},), z) finalizer(clear_acb, z) diff --git a/src/libarb/ArbFloat.jl b/src/libarb/ArbFloat.jl index 7cdba128..b4449433 100644 --- a/src/libarb/ArbFloat.jl +++ b/src/libarb/ArbFloat.jl @@ -24,6 +24,7 @@ mutable struct ArbFloat{P} <: AbstractFloat # P is the precision in bits d2::UInt # (d1, d2) the final part indicating the significand, or 0 function ArbFloat{P}() where {P} + minprec(P, ArbFloat) z = new{P}(0, 0, 0, 0) ccall(@libarb(arf_init), Cvoid, (Ref{ArbFloat},), z) finalizer(arf_clear, z) @@ -79,6 +80,7 @@ function copy(x::ArbFloat{P}) where {P} end function copy(x::ArbFloat{P}, bitprecision::Int, roundingmode::RoundingMode) where {P} + minprec(bitprecision, ArbFloat, copy) z = ArbFloat{P}() rounding = match_rounding_mode(roundingmode) rounddir = ccall(@libarb(arf_set_round), Cint, (Ref{ArbFloat}, Ref{ArbFloat}, Clong, Cint), z, x, bitprecision, rounding) diff --git a/src/libarb/ArbReal.jl b/src/libarb/ArbReal.jl index 802c430d..0906f6a5 100644 --- a/src/libarb/ArbReal.jl +++ b/src/libarb/ArbReal.jl @@ -33,6 +33,7 @@ mutable struct ArbReal{P} <: Real # P is the precision in bits rad_man::UInt # mp_limb_t radius, unsigned by definition function ArbReal{P}() where {P} + minprec(P, ArbReal) z = new{P}(0,0,0,0,0,0) ccall(@libarb(arb_init), Cvoid, (Ref{ArbReal},), z) finalizer(arb_clear, z) @@ -44,31 +45,24 @@ end const PtrToArbReal = Ptr{ArbReal} # arb_ptr const PtrToPtrToArbReal = Ptr{Ptr{ArbReal}} # arb_ptr* - arb_clear(x::ArbReal{P}) where {P} = ccall(@libarb(arb_clear), Cvoid, (Ref{ArbReal},), x) ArbReal{P}(x::ArbReal{P}) where {P} = x -ArbReal(x::ArbReal{P}) where {P} = x +ArbReal(x::ArbReal{P}) where {P} = ArbReal{P}(x) -ArbReal{P}(x::Missing) where {P} = missing -ArbReal(x::Missing) = missing +ArbReal{P}(::Missing) where {P} = missing +ArbReal(::Missing) = missing @inline sign_bit(x::ArbReal{P}) where {P} = isodd(x.mid_size) function ArbReal(x::T; bits::Int=0, digits::Int=0, base::Int=iszero(bits) ? 10 : 2) where {T<:Number} - bits > 0 && bits < MINIMUM_PRECISION_BASE2 && throw(DomainError("bit precision $bits < $MINIMUM_PRECISION_BASE2")) - digits > 0 && digits < MINIMUM_PRECISION_BASE10 && throw(DomainError("digit precision $digits < $MINIMUM_PRECISION_BASE10")) - if base === 10 - bits = digits > 0 ? bits4digits(digits)+extrabits() : (bits > 0 ? bits+extrabits() : DEFAULT_PRECISION.x) - elseif base === 2 - bits = bits > 0 ? bits+extrabits() : (digits > 0 ? digits+extrabits() : DEFAULT_PRECISION.x) - else - throw(ErrorException("base expects 2 or 10")) - end + bits = get_bits(bits, digits, base) ArbReal{bits}(x) end -swap(x::ArbReal{P}, y::ArbReal{P}) where {P} = ccall(@libarb(arb_swap), Cvoid, (Ref{ArbReal}, Ref{ArbReal}), x, y) +function swap!(x::ArbReal{P}, y::ArbReal{P}) where {P} + ccall(@libarb(arb_swap), Cvoid, (Ref{ArbReal}, Ref{ArbReal}), x, y) +end function copy(x::ArbReal{P}) where {P} z = ArbReal{P}() @@ -76,45 +70,41 @@ function copy(x::ArbReal{P}) where {P} return z end -function ArbReal{P}(x::Int32) where {P} +function ArbReal{P}(x::ArbInts) where {P} z = ArbReal{P}() - ccall(@libarb(arb_set_si), Cvoid, (Ref{ArbReal}, Clong), z, x) + ccall(@libarb(arb_set_si), Cvoid, (Ref{ArbReal}, Slong), z, x) return z end -ArbReal{P}(x::T) where {P, T<:Union{Int8, Int16}} = ArbReal{P}(Int32(x)) -ArbReal{P}(x::T) where {P, T<:Union{Int64, Int128}} = ArbReal{P}(BigInt(x)) -function ArbReal{P}(x::UInt32) where {P} +function ArbReal{P}(x::USlong) where {P} z = ArbReal{P}() - ccall(@libarb(arb_set_ui), Cvoid, (Ref{ArbReal}, Culong), z, x) + ccall(@libarb(arb_set_ui), Cvoid, (Ref{ArbReal}, USlong), z, x) return z end -ArbReal{P}(x::T) where {P, T<:Union{UInt8, UInt16}} = ArbReal{P}(UInt32(x)) -ArbReal{P}(x::T) where {P, T<:Union{UInt64, UInt128}} = ArbReal{P}(BigInt(x)) +ArbReal{P}(x::T) where {P, T<:Integer} = ArbReal{P}(BigFloat(x, precision=P + 2)) -function ArbReal{P}(x::Float64) where {P} +function ArbReal{P}(x::Base.IEEEFloat) where {P} z = ArbReal{P}() ccall(@libarb(arb_set_d), Cvoid, (Ref{ArbReal}, Cdouble), z, x) return z end -ArbReal{P}(x::T) where {P, T<:Union{Float16, Float32}} = ArbReal{P}(Float64(x)) function ArbReal{P}(x::BigFloat) where {P} z = ArbReal{P}() ccall(@libarb(arb_set_interval_mpfr), Cvoid, (Ref{ArbReal}, Ref{BigFloat}, Ref{BigFloat}, Clong), z, x, x, P) return z end -ArbReal{P}(x::BigInt) where {P} = ArbReal{P}(BigFloat(x)) -ArbReal{P}(x::Rational{T}) where {P, T<:Signed} = ArbReal{P}(BigFloat(x)) + +ArbReal{P}(x::Rational{T}) where {P, T<:Signed} = ArbReal{P}(BigFloat(x; precision=P + 2)) BigFloat(x::ArbReal{P}) where {P} = BigFloat(ArbFloat{P}(x)) BigInt(x::ArbReal{P}) where {P} = BigInt(trunc(BigFloat(ArbFloat{P}(x)))) -function Base.Integer(x::ArbReal{P}) where {P} +function Base.Integer(x::ArbReal) if isinteger(x) - abs(x) <= typemax(Int64) ? Int64(x) : BigInt(x) + typemin(Int64) <= x <= typemax(Int64) ? Int64(x) : BigInt(x) else - throw(InexactError("$x")) + throw(InexactError(:Integer, Integer, x)) end end @@ -124,9 +114,6 @@ function ArbReal{P}(x::ArbFloat{P}) where {P} return z end -ArbFloat{P}(x::ArbReal{P}, bits::Int) where P = ArbFloat(x, bits=bits) -ArbReal{P}(x::ArbFloat{P}, bits::Int) where P = ArbReal(x, bits=bits) - for (F,A) in ((:floor, :arf_floor), (:ceil, :arf_ceil)) @eval begin function $F(x::ArbReal{P}) where {P} @@ -147,7 +134,6 @@ end trunc(::Type{T}, x::ArbReal{P}) where {P, T<:Integer} = T(trunc(x)) - function modf(x::ArbReal{P}) where {P} ipart = trunc(x) fpart = x - ipart @@ -159,7 +145,6 @@ fmod(fpartipart::Tuple{ArbReal{P}, ArbReal{P}}) where {P} = fmod(fpart::ArbReal{P}, ipart::ArbReal{P}) where {P}= fpart + ipart - div(x::ArbReal{P}, y::ArbReal{P}) where {P} = trunc(x / y) @@ -172,7 +157,6 @@ function divrem(x::ArbReal{P}, y::ArbReal{P}) where {P} return (dv, rm) end - fld(x::ArbReal{P}, y::ArbReal{P}) where {P} = floor(x / y) @@ -185,9 +169,6 @@ function fldmod(x::ArbReal{P}, y::ArbReal{P}) where {P} return (fd, md) end -cld(x::ArbReal{P}, y::ArbReal{P}) where {P} = - ceil(x / y) - function midpoint(x::ArbReal{P}) where {P} z = ArbReal{P}() ccall(@libarb(arb_get_mid_arb), Cvoid, (Ref{ArbReal}, Ref{ArbReal}), z, x) @@ -224,14 +205,13 @@ function radius(x::ArbReal{P}, ::Type{ArbFloat{P}}) where {P} return z end - # a type specific hash function helps the type to 'just work' const hash_arbreal_lo = (UInt === UInt64) ? 0x37e642589da3416a : 0x5d46a6b4 const hash_0_arbreal_lo = hash(zero(UInt), hash_arbreal_lo) # two values of the same precision -# with identical midpoint significands and identical radus_exponentOf2s hash equal +# with identical midpoint significands and identical radius_exponentOf2s hash equal # they are the same value, one is less accurate yet centered about the other Base.hash(z::ArbReal{P}, h::UInt) where {P} = - hash(z.d1 ⊻ z.exp, - (h ⊻ hash(z.d2 ⊻ (~reinterpret(UInt,P)), hash_arbreal_lo) + hash(z.mid_d1 ⊻ z.mid_exp, + (h ⊻ hash(z.mid_d2 ⊻ (~reinterpret(UInt,P)), hash_arbreal_lo) ⊻ hash_0_arbreal_lo)) diff --git a/src/libarb/string.jl b/src/libarb/string.jl index 0ab3708b..d75d8c07 100644 --- a/src/libarb/string.jl +++ b/src/libarb/string.jl @@ -14,7 +14,7 @@ const NO_RADIUS = ARB_STR_NO_RADIUS @inline function digit_precision(bitprecision::Int) bitprecision = evincedbits(bitprecision) - return maximin_digits(bitprecision) + return max(maximin_digits(bitprecision), 1, (bitprecision + 1) ÷ 2) end function trimzeros(str::String) @@ -25,21 +25,21 @@ function trimzeros(str::String) str1b = trimallzeros(str1b) if str1b === "" str1b = "0" - end + end str1 = join((str1a, str1b), '.') end str = join((str1, str2), 'e') - elseif occursin('.', str) + elseif occursin('.', str) str1a, str1b = String.(split(str, '.')) str1b = trimallzeros(str1b) if str1b === "" str1b = "0" - end + end str = join((str1a, str1b), '.') - end + end return str end - + function trimallzeros(str::String) m = n = length(str) (n === 0 || str[n] !== '0') && return str @@ -69,7 +69,7 @@ function string(x::ArbFloat{P}; midpoint::Bool=false) where {P} end stringall(x::ArbFloat{P}) where {P} = return arbstring(x, floor(Int,log10(2)*precision(x))+10; flags=NO_RADIUS) - + function arbstring(x::ArbFloat{P}, maxdigits::Int=digit_precision(P); flags::UInt = NO_RADIUS) where {P} z = ArbReal{P}() ccall(@libarb(arb_set_arf), Cvoid, (Ref{ArbReal}, Ref{ArbFloat}), z, x) @@ -88,38 +88,29 @@ function string(x::ArbReal{P}; midpoint::Bool=false, radius::Bool=false) where { end stringall(x::ArbReal{P}; radius::Bool=false) where {P} = - radius ? arbstring(x, floor(Int,log10(2)*precision(x))+10; flags=ARB_STR_RADIUS) : - arbstring(x, floor(Int,log10(2)*precision(x))+10; flags=ARB_STR_NO_RADIUS) + radius ? arbstring(x, ceil(Int,log10(2)*workingprecision(x)); flags=ARB_STR_RADIUS) : + arbstring(x, ceil(Int,log10(2)*workingprecision(x)); flags=ARB_STR_NO_RADIUS) function arbstring(x::ArbReal{P}, maxdigits::Int=digit_precision(P); flags::UInt = NO_FLAGS) where {P} - unsafestr = ccall(@libarb(arb_get_str), Cstring, - (Ref{ArbReal}, Clong, Culong), x, maxdigits, flags) + unsafestr = ccall(@libarb(arb_get_str), Cstring, (Ref{ArbReal}, Clong, Culong), x, maxdigits, flags) str = deepcopy( unsafe_string(pointer(unsafestr)) ) str = trimzeros(str) ccall(@libflint(flint_free), Cvoid, (Cstring,), unsafestr) return str end - function string(x::ArbComplex{P}; midpoint::Bool=false, radius::Bool=false) where {P} prec = midpoint ? digits4bits(P) : digit_precision(P) flags = radius ? ARB_STR_RADIUS : ARB_STR_NO_RADIUS return arbstring(x, prec, flags=flags) end -stringall(x::ArbComplex{P}; radius::Bool=false) where {P} = - radius ? arbstring(x, floor(Int,log10(2)*precision(x))+10; flags=ARB_STR_RADIUS) : - arbstring(x, floor(Int,log10(2)*precision(x))+10; flags=ARB_STR_NO_RADIUS) - function arbstring(x::ArbComplex{P}, maxdigits::Int=digit_precision(P); flags::UInt = NO_FLAGS) where {P} # rea, ima = real(x), imag(x) - rea = ArbReal{P}() - ima = ArbReal{P}() - ccall(@libarb(acb_get_real), Cvoid, (Ref{ArbReal}, Ref{ArbComplex}), rea, x) - ccall(@libarb(acb_get_imag), Cvoid, (Ref{ArbReal}, Ref{ArbComplex}), ima, x) + rea = ArbReal(x) + ima = ArbReal(x) ima_isneg = signbit(ima) - ima_abs = ima_isneg ? -ima : ima connection = ima_isneg ? " - " : " + " rea_str = arbstring(rea, maxdigits, flags=flags) ima_str = arbstring(abs(ima), maxdigits, flags=flags) diff --git a/src/support/minprec.jl b/src/support/minprec.jl index 0868ea28..f7cc9658 100644 --- a/src/support/minprec.jl +++ b/src/support/minprec.jl @@ -1,11 +1,13 @@ # protect minimum precision -minprecerror() = throw(DomainError("minimum precision is 24 bits")) +@noinline function minprecerror(n, minp, t, f) + text = t === f ? "" : " for $t" + throw(DomainError("$f: minimum precision$text is $minp bits but was $n")) +end -for T in (:Float32, :Float64, :Int16, :Int32, :Int64, :Int128, :BigFloat) - for N in collect(0:23) - @eval ArbFloat{$N}(x::$T) = minprecerror() - @eval ArbReal{$N}(x::$T) = minprecerror() - @eval ArbComplex{$N}(x::$T) = minprecerror() - @eval ArbComplex{$N}(x::$T, y::$T) = minprecerror() - end +function minprec(n, t, f=t) + minp = max(extrabits(ArbFloat), 1) + if n < minp + minprecerror(n, minp, t, f) + end + nothing end diff --git a/src/values/constructors.jl b/src/values/constructors.jl index 30984ae3..827c2fe0 100644 --- a/src/values/constructors.jl +++ b/src/values/constructors.jl @@ -70,16 +70,16 @@ for I in (:Int8, :Int16, :Int32, :Int64, :Int128) @eval begin function $I(x::ArbReal{P}) where {P} - (!isexact(x) | !isinteger(x)) && throw(InexactError("$(x) is not an integer")) + (!isexact(x) || !isinteger(x)) && throw(InexactError(nameof($I), $I, x)) bi = BigInt(BigFloat(x)) - !(typemin($I) <= bi <= typemax($I)) && throw(InexactError("$(x)")) + !(typemin($I) <= bi <= typemax($I)) && throw(InexactError(nameof($I), $I, x)) return $I(bi) end function $I(x::ArbComplex{P}) where {P} - (!isexact(x) | !isinteger(x) | !iszero(imag(x))) && throw(InexactError("$(x) is not an integer")) + (!isexact(x) | !isinteger(x) | !iszero(imag(x))) && throw(InexactError(nameof($I), $I, x)) bi = BigInt(BigFloat(x)) - !(typemin($I) <= bi <= typemax($I)) && throw(InexactError("$(x)")) + !(typemin($I) <= bi <= typemax($I)) && throw(InexactError(nameof($I), $I, x)) return $I(bi) end @@ -113,19 +113,19 @@ end # ArbFloat{P}(x::T) where {P,T<:Real} = ArbFloat{P}(BigFloat(x)) # ArbFloat(x::T) where {T<:Real} = ArbFloat{workingprecision(ArbFloat)}(BigFloat(x)) ArbFloat{P}(x::T) where {P,T<:Complex} = ArbFloat{P}(BigFloat(real(x))) -ArbFloat(x::T) where {T<:Complex} = ArbFloat{workingprecision(ArbFloat)}(BigFloat(real(x))) +ArbFloat(x::T) where {T<:Complex} = ArbFloat{workingprecision(ArbFloat)}(BigFloat(real(x))) # TODO shouldn't that InexactError -ArbReal{P}(x::T) where {P,T<:Real} = ArbReal{P}(BigFloat(x)) -ArbReal(x::T) where {T<:Real} = ArbReal{workingprecision(ArbReal)}(BigFloat(real(x))) -ArbReal{P}(x::T) where {P,T<:Complex} = ArbReal{P}(BigFloat(x)) +#ArbReal{P}(x::T) where {P,T<:Real} = ArbReal{P}(BigFloat(x)) +#ArbReal(x::T) where {T<:Real} = ArbReal{workingprecision(ArbReal)}(BigFloat(real(x))) +ArbReal{P}(x::T) where {P,T<:Complex} = ArbReal{P}(BigFloat(x)) # TODO shouldn't that InexactError ArbReal(x::T) where {T<:Complex} = ArbReal{workingprecision(ArbReal)}(BigFloat(real(x))) # retype ArbFloat(x::ArbReal{P}) where {P} = ArbFloat{P}(x) -ArbFloat(x::ArbComplex{P}) where {P} = ArbFloat{P}(real(x)) +ArbFloat(x::ArbComplex{P}) where {P} = ArbFloat{P}(real(x)) # TODO shouldn't that InexactError ArbReal(x::ArbFloat{P}) where {P} = ArbReal{P}(x) -ArbReal(x::ArbComplex{P}) where {P} = ArbReal{P}(real(x)) +ArbReal(x::ArbComplex{P}) where {P} = ArbReal{P}(real(x)) # TODO shouldn't that InexactError ArbFloat{Q}(x::ArbReal{P}) where {P,Q} = ArbFloat{Q}(ArbReal{Q}(x)) ArbFloat{Q}(x::ArbComplex{P}) where {P,Q} = ArbFloat{Q}(ArbReal{Q}(real(x))) @@ -135,6 +135,7 @@ ArbReal{Q}(x::ArbComplex{P}) where {P,Q} = ArbReal{Q}(real(x)) # change precision function ArbFloat{P}(x::ArbFloat{Q}, roundingmode::RoundingMode) where {P,Q} + minprec(P, ArbFloat) rounding = match_rounding_mode(roundingmode) z = ArbFloat{P}() rnd = ccall(@libarb(arf_set_round), Cint, @@ -144,6 +145,7 @@ end ArbFloat{P}(x::ArbFloat{Q}) where {P,Q} = ArbFloat{P}(x, RoundNearest) function ArbReal{P}(x::ArbReal{Q}) where {P,Q} + minprec(P, ArbReal) z = ArbReal{P}() ccall(@libarb(arb_set_round), Cvoid, (Ref{ArbReal}, Ref{ArbReal}, Clong), z, x, P) @@ -151,6 +153,7 @@ function ArbReal{P}(x::ArbReal{Q}) where {P,Q} end function ArbComplex{P}(x::ArbComplex{Q}) where {P,Q} + minprec(P, ArbComplex) z = ArbComplex{P}() ccall(@libarb(acb_set_round), Cvoid, (Ref{ArbComplex}, Ref{ArbComplex}, Clong), z, x, P) diff --git a/src/values/rounding.jl b/src/values/rounding.jl index 66a4b909..6ccee7dc 100644 --- a/src/values/rounding.jl +++ b/src/values/rounding.jl @@ -73,9 +73,11 @@ end # int arf_set_round(arf_t res, const arf_t x, slong prec, arf_rnd_t rnd) function roundbits(x::ArbFloat{P}, roundingmode::RoundingMode, sigbits::Int) where {P} + minprec(sigbits, ArbFloat{P}, roundbits) sigbits >= P && return ArbFloat{sigbits}(x) z = ArbFloat{P}() rounding = match_rounding_mode(roundingmode) + # sigbits <= 0 evokes SIGSEGFAULT rounddir = ccall(@libarb(arf_set_round), Cint, (Ref{ArbFloat}, Ref{ArbFloat}, Clong, Cint), z, x, sigbits, rounding) return z end @@ -100,8 +102,10 @@ rounddigits!(x::ArbFloat{P}, roundingmode::RoundingMode, digits::Int) where {P} rounddigits!(x::ArbFloat{P}, digits::Int) where {P} = roundbits!(x, RoundNearest, bits4digits(digits)) function roundbits(x::ArbReal{P}, roundingmode::RoundingMode, sigbits::Int) where {P} + minprec(sigbits, ArbReal{P}, roundbits) sigbits >= P && return ArbReal{sigbits}(x) z = ArbReal{P}() + # sigbits <= 0 evokes SIGSEGFAULT ccall(@libarb(arb_set_round), Cvoid, (Ref{ArbReal}, Ref{ArbReal}, Clong), z, x, sigbits) return z end @@ -125,8 +129,10 @@ rounddigits!(x::ArbReal{P}, roundingmode::RoundingMode, digits::Int) where {P} = rounddigits!(x::ArbReal{P}, digits::Int) where {P} = roundbits!(x, RoundNearest, bits4digits(digits)) function roundbits(x::ArbComplex{P}, roundingmode::RoundingMode, sigbits::Int) where {P} + minprec(sigbits, ArbComplex{P}, roundbits) sigbits >= P && return ArbComplex{sigbits}(x) z = ArbComplex{P}() + # sigbits <= 0 evokes SIGSEGFAULT ccall(@libarb(acb_set_round), Cvoid, (Ref{ArbComplex}, Ref{ArbComplex}, Clong), z, x, sigbits) return z end diff --git a/test/complex.jl b/test/complex.jl index e73bf183..871d3a81 100644 --- a/test/complex.jl +++ b/test/complex.jl @@ -32,6 +32,7 @@ @test_throws MethodError ismissing(ArbComplex{99}(missing, 1+im)) end @testset "ArbComplex constructors" begin + @test_throws DomainError ArbComplex{23}(0) @test_throws ErrorException ArbComplex(5.0, base=3) @test_throws ErrorException ArbComplex(5.0, 1.0, base=3) @test ArbComplex{128}(1) == 1 @@ -75,6 +76,7 @@ @test ArbComplex(Int32(-1), Int32(-2)) == ArbComplex(-1.0, -2.0) @test ArbComplex(typemax(Int32) + 1) == ArbComplex(float(typemax(Int32) + 1)) @test hash(ArbComplex(-2, -2)) isa UInt + @test_broken hash(ArbComplex(2, 1)) == hash(Complex(2, 1)) @test ArbComplex{256}(1.0 + im) == ArbComplex(1 + im, bits=256) @test ArbComplex{256}(1.0 + im) == ArbComplex(1, 1, bits=256) diff --git a/test/misc.jl b/test/misc.jl index 533919b4..c6a0b41c 100644 --- a/test/misc.jl +++ b/test/misc.jl @@ -9,7 +9,7 @@ @test_throws InexactError ArbNumerics.Mag(-1) end - @testset "$T($v)" for T in (ArbFloat, ArbFloat{64}, ArbReal, ArbReal{64}), + @testset "$T($v)" for T in (ArbReal, ArbReal{64}, ArbReal, ArbReal{64}), (v, c, f, t) in [(0.6, 1, 0, 0), (-0.6, 0, -1, 0)] x = T(v) @@ -31,6 +31,7 @@ end @testset "ArbFloat" begin + @test_throws DomainError ArbFloat{23}(0) a = ArbFloat(10) b = ArbFloat(20) swap!(a, b) @@ -40,11 +41,12 @@ @test ArbFloat(a) === a @test copy(a) == a @test copy(a) !== a - @test copy(a, 10, RoundDown) == 804 / 256 - @test copy(a, 10, RoundUp) == 805 / 256 + @test copy(a, 24, RoundDown) == 13176794 / 2^22 + @test copy(a, 24, RoundUp) == 13176795 / 2^22 @test copy(a, RoundUp) > a @test a > copy(a, RoundDown) - @test copy(a, 10) == copy(a, 10, RoundDown) + @test copy(a, 24) == copy(a, 24, RoundUp) + @test_throws DomainError copy(a, 10) @test ArbFloat(typemax(UInt)) == typemax(UInt) @test ArbFloat(typemax(Int128)) == typemax(Int128) @@ -91,5 +93,72 @@ @test fldmod(a, b) == (fld(a, b), mod(a, b)) @test hash(a) isa UInt + @test_broken hash(b) == hash(1) # because isequal(b, 1), see ?hash end + + @testset "ArbReal" begin + @test_throws DomainError ArbReal{23}(0) + a = ArbReal(10) + b = ArbReal(20) + swap!(a, b) + @test (a, b) == (20, 10) + + a = ArbReal(pi) + @test ArbReal(a) === a + @test copy(a) == a + @test copy(a) !== a + + @test ArbReal(typemax(UInt)) == typemax(UInt) + @test ArbReal(typemax(Int128)) == typemax(Int128) + @test ArbReal(1 // 3) ≈ ArbReal(1) / ArbReal(3) + end + + @testset "convert ArbReal to Integer" begin + v = typemax(Int32) + 1 + w = round(v) + a = ArbReal(v) + @test_throws InexactError Int32(a) + @test Int64(a) == w + @test Int128(a) == w + @test BigInt(a) == w + + a = ArbReal(1.5) + @test_throws InexactError Int(a) + @test_throws InexactError Integer(a) + v = typemin(Int32) - 1 + a = ArbReal(v) + x = Integer(a) + @test typeof(x) == Int64 + @test x == v + v = Int128(typemax(Int64)) + 1 + a = ArbReal(v) + x = Integer(a) + @test typeof(x) == BigInt + @test x == v + end + + @testset "functions with ArbReal" begin + a = ArbReal(ℯ) + @test midpoint(a) ≈ a + @test !iszero(radius(a)) + @test radius(a, ArbFloat) == radius(a) + @test radius(a, ArbNumerics.Mag) == radius(a) + + @test modf(a) == (a - 2, 2) + @test fmod(modf(a)) ≈ a + @test fmod(modf(a)...) ≈ a + + b = ArbReal(1) + @test iszero(radius(b)) + @test div(-a, b) == -2 + @test rem(-a, b) ≈ -a + 2 + @test fld(-a, b) == -3 + @test mod(-a, b) ≈ -a + 3 + @test divrem(a, b) == (div(a, b), rem(a, b)) + @test fldmod(a, b) == (fld(a, b), mod(a, b)) + + @test hash(a) isa UInt + @test_broken hash(b) == hash(1) # because isequal(b, 1), see ?hash + end + end From 2049e1354d75336fa06e3c26a34a44eddc7c220a Mon Sep 17 00:00:00 2001 From: KlausC Date: Mon, 6 Nov 2023 13:55:33 +0100 Subject: [PATCH 12/15] fixed `radius --- src/libarb/ArbReal.jl | 4 ++-- test/misc.jl | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/libarb/ArbReal.jl b/src/libarb/ArbReal.jl index 0906f6a5..8a21d95d 100644 --- a/src/libarb/ArbReal.jl +++ b/src/libarb/ArbReal.jl @@ -198,8 +198,8 @@ end midpoint(x::ArbReal{P}, ::Type{ArbFloat{P}}) where {P} = midpoint_byref(x) radius(x::ArbReal{P}, ::Type{Mag}) where {P} = radius_byref(x) -function radius(x::ArbReal{P}, ::Type{ArbFloat{P}}) where {P} - z = ArbFloat{P}() +function radius(x::ArbReal, ::Type{T}) where {T<:ArbFloat} + z = ArbFloat{workingprecision(T)}() mag = radius_byref(x) ccall(@libarb(arf_set_mag), Cvoid, (Ref{ArbFloat}, Ref{Mag}), z, mag) return z diff --git a/test/misc.jl b/test/misc.jl index c6a0b41c..6d989fa5 100644 --- a/test/misc.jl +++ b/test/misc.jl @@ -141,7 +141,7 @@ a = ArbReal(ℯ) @test midpoint(a) ≈ a @test !iszero(radius(a)) - @test radius(a, ArbFloat) == radius(a) + @test radius(a, ArbFloat) ≈ radius(a) rtol=1e-5 @test radius(a, ArbNumerics.Mag) == radius(a) @test modf(a) == (a - 2, 2) From 05ffbb67461b108f65e1dc36994f8b5c3cb64d0e Mon Sep 17 00:00:00 2001 From: KlausC Date: Mon, 6 Nov 2023 18:52:09 +0100 Subject: [PATCH 13/15] new ArbTypes.jl, rounding, tests --- src/ArbNumerics.jl | 8 +-- src/ArbTypes.jl | 94 ++++++++++++++++++++++++++ src/libarb/ArbComplex.jl | 53 ++------------- src/libarb/ArbComplexMatrix.jl | 39 ++++++----- src/libarb/ArbFloat.jl | 43 ++++++------ src/libarb/ArbFloatMatrix.jl | 26 ++++---- src/libarb/ArbReal.jl | 27 -------- src/libarb/ArbRealMatrix.jl | 8 +-- src/support/complex.jl | 22 +++--- src/values/constructors.jl | 118 +++++++++------------------------ test/complex.jl | 9 ++- test/misc.jl | 28 +++++++- 12 files changed, 234 insertions(+), 241 deletions(-) create mode 100644 src/ArbTypes.jl diff --git a/src/ArbNumerics.jl b/src/ArbNumerics.jl index ddc27675..d0e00cd1 100644 --- a/src/ArbNumerics.jl +++ b/src/ArbNumerics.jl @@ -85,7 +85,7 @@ import Base: IEEEFloat, Float16, Float32, Float64, float, UInt8, UInt16, UInt32, UInt64, UInt128, Int8, Int16, Int32, Int64, Int128, - BigInt, BigFloat, Rational, Complex, real, imag, complex, angle, + BigInt, BigFloat, Rational, Complex, real, imag, isreal, complex, angle, floatmax, floatmin, typemax, typemin, maxintfloat, rationalize, @@ -153,9 +153,7 @@ using Libdl using Random using Random: SamplerType, SamplerTrivial, CloseOpen01 -const Slong = Int # to accomodate windows -const USlong = unsigned(Slong) -const ArbInts = Union{Int,Int32,Int16,Int8} # Int may be Int32 +include("ArbTypes.jl") include("support/abstractions.jl") include("support/matrices.jl") @@ -173,8 +171,6 @@ include("support/ArblibVector.jl") include("support/random.jl") -const ArbNumber{P} = Union{ArbFloat{P}, ArbReal{P}, ArbComplex{P}} - include("libarb/ArbMatrix.jl") # must preceed ArbRealMatrix include("libarb/ArbRealMatrix.jl") # must preceed ArbFloatMatrix include("libarb/ArbFloatMatrix.jl") # must preceed ArbComplexMatrix diff --git a/src/ArbTypes.jl b/src/ArbTypes.jl new file mode 100644 index 00000000..53467a26 --- /dev/null +++ b/src/ArbTypes.jl @@ -0,0 +1,94 @@ + +mutable struct ArbFloat{P} <: AbstractFloat # P is the precision in bits + exp::Int # fmpz exponent of 2 (2^exp) + size::UInt # mp_size_t nwords and sign (lsb holds sign of significand) + d1::UInt # significand unsigned, immediate value or the initial span + d2::UInt # (d1, d2) the final part indicating the significand, or 0 + + function ArbFloat{P}() where {P} + minprec(P, ArbFloat) + z = new{P}(0, 0, 0, 0) + ccall(@libarb(arf_init), Cvoid, (Ref{ArbFloat},), z) + finalizer(arf_clear, z) + return z + end +end + +# for use within structs, e.g ArbFloatMatrix +const PtrToArbFloat = Ptr{ArbFloat} # arf_ptr +const PtrToPtrToArbFloat = Ptr{Ptr{ArbFloat}} # arf_ptr* + +arf_clear(x::ArbFloat{P}) where {P} = ccall(@libarb(arf_clear), Cvoid, (Ref{ArbFloat},), x) + +# ArbReal structs hold balls given as midpoint, radius +mutable struct ArbReal{P} <: Real # P is the precision in bits + # midpoint + mid_exp::Int # fmpz exponent of 2 (2^exp) + mid_size::UInt # mp_size_t nwords and sign (lsb holds sign of significand) + mid_d1::UInt # significand unsigned, immediate value or the initial span + mid_d2::UInt # (d1, d2) the final part indicating the significand, or 0 + # radius + rad_exp::Int # fmpz exponent of 2 (2^exp) + rad_man::UInt # mp_limb_t radius, unsigned by definition + + function ArbReal{P}() where {P} + minprec(P, ArbReal) + z = new{P}(0, 0, 0, 0, 0, 0) + ccall(@libarb(arb_init), Cvoid, (Ref{ArbReal},), z) + finalizer(arb_clear, z) + return z + end +end + +# for use within structs, e.g. ArbRealMatrix +const PtrToArbReal = Ptr{ArbReal} # arb_ptr +const PtrToPtrToArbReal = Ptr{Ptr{ArbReal}} # arb_ptr* + +arb_clear(x::ArbReal{P}) where {P} = ccall(@libarb(arb_clear), Cvoid, (Ref{ArbReal},), x) + +# ArbComplex structs hold complex balls given as ArbReal pairs + +mutable struct ArbComplex{P} <: Number # P is the precision in bits + # real midpoint + real_mid_exp::Int # fmpz exponent of 2 (2^exp) + real_mid_size::UInt # mp_size_t nwords and sign (lsb holds sign of significand) + real_mid_d1::UInt # significand unsigned, immediate value or the initial span + real_mid_d2::UInt # (d1, d2) the final part indicating the significand, or 0 + # real radius + real_rad_exp::Int # fmpz exponent of 2 (2^exp) + real_rad_man::UInt # mp_limb_t radius, unsigned by definition + # imaginary midpoint + imag_mid_exp::Int # fmpz exponent of 2 (2^exp) + imag_mid_size::UInt # mp_size_t nwords and sign (lsb holds sign of significand) + imag_mid_d1::UInt # significand unsigned, immediate value or the initial span + imag_mid_d2::UInt # (d1, d2) the final part indicating the significand, or 0 + # imaginary radius + imag_rad_exp::Int # fmpz exponent of 2 (2^exp) + imag_rad_man::UInt # mp_limb_t radius, unsigned by definition + + + function ArbComplex{P}() where {P} + minprec(P, ArbComplex) + z = new{P}(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) + ccall(@libarb(acb_init), Cvoid, (Ref{ArbComplex},), z) + finalizer(clear_acb, z) + return z + end +end + +# finalizer: +clear_acb(x::ArbComplex{P}) where {P} = ccall(@libarb(acb_clear), Cvoid, (Ref{ArbComplex},), x) + +ArbComplex{P}(x::T) where {P,T<:Number} = ArbComplex{P}(real(x), imag(x)) + +# for use within a struct, eg. ArbComplexMatrix +const PtrToArbComplex = Ptr{ArbComplex} # acb_ptr +const PtrToPtrToArbComplex = Ptr{Ptr{ArbComplex}} # acb_ptr* + +const ArbFloatReal{P} = Union{ArbFloat{P},ArbReal{P}} +const ArbNumber{P} = Union{ArbFloat{P}, ArbReal{P}, ArbComplex{P}} + +# Types supporting C-API +const Slong = Int # to accomodate windows +const USlong = unsigned(Slong) +const ArbInts = Union{Int,Int32,Int16,Int8} # Int may be Int32 diff --git a/src/libarb/ArbComplex.jl b/src/libarb/ArbComplex.jl index bb1f5672..e099a98a 100644 --- a/src/libarb/ArbComplex.jl +++ b/src/libarb/ArbComplex.jl @@ -26,48 +26,6 @@ See also: [`ArbFloat`](@ref), [`ArbReal`](@ref), [`ball`](@ref), [`setball`](@re """ ArbComplex -# ArbComplex structs hold complex balls given as ArbReal pairs - -mutable struct ArbComplex{P} <: Number # P is the precision in bits - # real midpoint - real_mid_exp::Int # fmpz exponent of 2 (2^exp) - real_mid_size::UInt # mp_size_t nwords and sign (lsb holds sign of significand) - real_mid_d1::UInt # significand unsigned, immediate value or the initial span - real_mid_d2::UInt # (d1, d2) the final part indicating the significand, or 0 - # real radius - real_rad_exp::Int # fmpz exponent of 2 (2^exp) - real_rad_man::UInt # mp_limb_t radius, unsigned by definition - # imaginary midpoint - imag_mid_exp::Int # fmpz exponent of 2 (2^exp) - imag_mid_size::UInt # mp_size_t nwords and sign (lsb holds sign of significand) - imag_mid_d1::UInt # significand unsigned, immediate value or the initial span - imag_mid_d2::UInt # (d1, d2) the final part indicating the significand, or 0 - # imaginary radius - imag_rad_exp::Int # fmpz exponent of 2 (2^exp) - imag_rad_man::UInt # mp_limb_t radius, unsigned by definition - - - function ArbComplex{P}() where {P} - minprec(P, ArbComplex) - z = new{P}(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) - ccall(@libarb(acb_init), Cvoid, (Ref{ArbComplex},), z) - finalizer(clear_acb, z) - return z - end -end - -ArbComplex{P}(x::T) where {P,T<:Number} = ArbComplex{P}(real(x), imag(x)) - -# for use within a struct, eg. ArbComplexMatrix -const PtrToArbComplex = Ptr{ArbComplex} # acb_ptr -const PtrToPtrToArbComplex = Ptr{Ptr{ArbComplex}} # acb_ptr* - -const ArbFloatReal{P} = Union{ArbFloat{P},ArbReal{P}} -const ArbNumber1{P} = Union{ArbFloatReal{P},ArbComplex{P}} - -# finalizer: -clear_acb(x::ArbComplex{P}) where {P} = ccall(@libarb(acb_clear), Cvoid, (Ref{ArbComplex},), x) - ArbComplex{P}(x::ArbComplex{P}) where {P} = x ArbComplex(x::ArbComplex{P}) where {P} = x @@ -83,7 +41,7 @@ function ArbComplex(x::T, y::Real=0; bits::Int=0, digits::Int=0, base::Int=iszer iszero(y) ? ArbComplex{bits}(x) : ArbComplex{bits}(x, y) end -ArbComplex(x::ArbNumber1{P}) where {P} = ArbComplex{P}(x) +ArbComplex(x::ArbNumber{P}) where {P} = ArbComplex{P}(x) ArbComplex(x::ArbFloatReal{P}, y::ArbFloatReal{S}) where {P,S} = ArbComplex{min(P, S)}(x, y) ArbComplex(x::ArbFloatReal{P}, y::Real) where {P} = ArbComplex{P}(x, y) @@ -132,14 +90,14 @@ end deepcopy(x::ArbComplex{P}) where {P} = copy(x) -Complex(x::ArbNumber1) = convert(ArbComplex, x) +Complex(x::ArbNumber) = convert(ArbComplex, x) Complex(re::ArbFloatReal, im::ArbFloatReal) = ArbComplex(re, im) -complex(re::ArbNumber1) = convert(ArbComplex, re) +complex(re::ArbNumber) = convert(ArbComplex, re) complex(re::ArbFloatReal, im::ArbFloatReal) = ArbComplex(re, im) -complex(::Type{T}) where {T<:ArbNumber1} = ArbComplex -complex(::Type{T}) where {P,T<:ArbNumber1{P}} = ArbComplex{P} +complex(::Type{T}) where {T<:ArbNumber} = ArbComplex +complex(::Type{T}) where {P,T<:ArbNumber{P}} = ArbComplex{P} Base.Complex{T}(x::ArbComplex{P}) where {P,T<:ArbFloatReal{P}} = x @@ -149,7 +107,6 @@ function ArbComplex{P}(x::T) where {P,T<:Integer} return z end - function ArbComplex{P}(x::T1, y::T2) where {P,T1<:Real,T2<:Real} x1 = convert(ArbReal{P}, x) y1 = convert(ArbReal{P}, y) diff --git a/src/libarb/ArbComplexMatrix.jl b/src/libarb/ArbComplexMatrix.jl index 71e01743..ad24e9e7 100644 --- a/src/libarb/ArbComplexMatrix.jl +++ b/src/libarb/ArbComplexMatrix.jl @@ -225,7 +225,7 @@ function matmul(x::ArbComplexMatrix{P}, y::ArbComplexMatrix{P}) where {P} z = ArbComplexMatrix{P}(rowcount(x), colcount(y)) ccall(@libarb(acb_mat_mul), Cvoid, (Ref{ArbComplexMatrix}, Ref{ArbComplexMatrix}, Ref{ArbComplexMatrix}, Cint), z, x, y, P) return Matrix(z) -end +end function matmul(x::Array{ArbComplex, 2}, y::Array{ArbComplex, 2}) xx = ArbComplexMatrix(x) @@ -253,7 +253,7 @@ end @inline function Base.:(*)(x::Array{ArbComplex,2}, y::Array{ArbComplex,2}) checkmulable(x, y) - P = workingprecision(ArbComplex) + P = workingprecision(ArbComplex) return matmul(ArbComplexMatrix{P}(x), ArbComplexMatrix{P}(y)) end @@ -457,7 +457,7 @@ function exp(x::Array{ArbComplex, 2}) end function exp(x::Array{ArbComplex{P}, 2}) where {P} - checksquare(x) + checksquare(x) y = ArbComplexMatrix(x) z = ArbComplexMatrix{P}(rowcount(x), colcount(x)) ccall(@libarb(acb_mat_exp), Cvoid, (Ref{ArbComplexMatrix}, Ref{ArbComplexMatrix}, Cint), z, y, P) @@ -465,7 +465,7 @@ function exp(x::Array{ArbComplex{P}, 2}) where {P} end function exp(x::ArbComplexMatrix{P}) where {P} - checksquare(x) + checksquare(x) z = ArbComplexMatrix{P}(rowcount(x), colcount(x)) ccall(@libarb(acb_mat_exp), Cvoid, (Ref{ArbComplexMatrix}, Ref{ArbComplexMatrix}, Cint), z, x, P) return Matrix(z) @@ -490,17 +490,17 @@ int acb_mat_approx_eig_qr(acb_ptr E, acb_mat_t L, acb_mat_t R, const acb_mat_t A acb_mat_approx_eig_qr(E, NULL, R, C, NULL, 0, prec) Computes floating-point approximations of all the n eigenvalues (and optionally eigenvectors) of the given n by n matrix A. -The approximations of the eigenvalues are written to the vector E, in no particular order. -If L is not NULL, approximations of the corresponding left eigenvectors are written to the rows of L. +The approximations of the eigenvalues are written to the vector E, in no particular order. +If L is not NULL, approximations of the corresponding left eigenvectors are written to the rows of L. If R is not NULL, approximations of the corresponding right eigenvectors are written to the columns of R. The parameters tol and maxiter can be used to control the target numerical error and the maximum number of iterations allowed before giving up. Passing NULL and 0 respectively results in default values being used. -Uses the implicitly shifted QR algorithm with reduction to Hessenberg form. -No guarantees are made about the accuracy of the output. -A nonzero return value indicates that the QR iteration converged numerically, -but this is only a heuristic termination test and does not imply any statement whatsoever about error bounds. +Uses the implicitly shifted QR algorithm with reduction to Hessenberg form. +No guarantees are made about the accuracy of the output. +A nonzero return value indicates that the QR iteration converged numerically, +but this is only a heuristic termination test and does not imply any statement whatsoever about error bounds. The output may also be accurate even if this function returns zero. =# @@ -510,11 +510,11 @@ function LinearAlgebra.eigvals(m::ArbComplexMatrix{P}) where {P} checksquare(m) n = rowcount(m) eigvalues = ccall(@libarb(_acb_vec_init), Ptr{ArbComplex}, (Cint,), n) - eigvectors = ArbComplexMatrix(rowcount(m), colcount(m)) - eigvectors2 = copy(m) + eigvectors = ArbComplexMatrix(rowcount(m), colcount(m)) + eigvectors2 = copy(m) tol = Base.C_NULL maxiter = 0 - result = ccall(@libarb(acb_mat_approx_eig_qr), Cint, + result = ccall(@libarb(acb_mat_approx_eig_qr), Cint, (Ptr{ArbComplex}, Ptr{Nothing}, Ref{ArbComplexMatrix}, Ref{ArbComplexMatrix}, Ptr{Nothing}, Clong, Clong), eigvalues, Base.C_NULL, eigvectors, eigvectors2, Base.C_NULL, maxiter, P) eigvalues = sort(eigvalues, lt=complex_lt) @@ -524,30 +524,29 @@ end function LinearAlgebra.eigvecs(m::ArbComplexMatrix{P}) where {P} checksquare(m) eigvalues = zeros(ArbComplex, rowcount(m)) - eigvectors = ArbComplexMatrix(rowcount(m), colcount(m)) + eigvectors = ArbComplexMatrix(rowcount(m), colcount(m)) tol = Base.C_NULL maxiter = 0 - result = ccall(@libarb(acb_mat_approx_eig_qr), Cint, + result = ccall(@libarb(acb_mat_approx_eig_qr), Cint, (Ref{Vector{ArbComplex}}, Ref{ArbComplexMatrix}, Ref{ArbComplexMatrix}, Ref{ArbComplexMatrix}, Ptr{Nothing}, Clong, Clong), eigvalues, Base.C_NULL, eigvectors, m, tol, maxiter, P) - return eigvectors + return eigvectors end function LinearAlgebra.eigen(m::ArbComplexMatrix{P}) where {P} checksquare(m) eigvalues = zeros(ArbComplex, rowcount(m)) - eigvectors = ArbComplexMatrix(rowcount(m), colcount(m)) + eigvectors = ArbComplexMatrix(rowcount(m), colcount(m)) tol = Base.C_NULL maxiter = 0 - result = ccall(@libarb(acb_mat_approx_eig_qr), Cint, + result = ccall(@libarb(acb_mat_approx_eig_qr), Cint, (Ref{Vector{ArbComplex}}, Ref{ArbComplexMatrix}, Ref{ArbComplexMatrix}, Ref{ArbComplexMatrix}, Ptr{Nothing}, Clong, Clong), eigvalues, Base.C_NULL, eigvectors, m, tol, maxiter, P) eigvalues = sort(eigvalues, lt=complex_lt) - return eigvalues, eigvectors + return eigvalues, eigvectors end LinearAlgebra.eigvals(m::Array{ArbComplex{P},2}) where {P} = eigvals(ArbComplexMatrix(m)) LinearAlgebra.eigvecs(m::Array{ArbComplex{P},2}) where {P} = eigvecs(ArbComplexMatrix(m)) LinearAlgebra.eigen(m::Array{ArbComplex{P},2}) where {P} = eigs(ArbComplexMatrix(m)) =# - diff --git a/src/libarb/ArbFloat.jl b/src/libarb/ArbFloat.jl index b4449433..ef6abaab 100644 --- a/src/libarb/ArbFloat.jl +++ b/src/libarb/ArbFloat.jl @@ -17,27 +17,6 @@ See also: [`ArbReal`](@ref), [`ArbComplex`](@ref), [`workingprecision`](@ref), [ """ ArbFloat -mutable struct ArbFloat{P} <: AbstractFloat # P is the precision in bits - exp::Int # fmpz exponent of 2 (2^exp) - size::UInt # mp_size_t nwords and sign (lsb holds sign of significand) - d1::UInt # significand unsigned, immediate value or the initial span - d2::UInt # (d1, d2) the final part indicating the significand, or 0 - - function ArbFloat{P}() where {P} - minprec(P, ArbFloat) - z = new{P}(0, 0, 0, 0) - ccall(@libarb(arf_init), Cvoid, (Ref{ArbFloat},), z) - finalizer(arf_clear, z) - return z - end -end - -# for use within structs, e.g ArbFloatMatrix -const PtrToArbFloat = Ptr{ArbFloat} # arf_ptr -const PtrToPtrToArbFloat = Ptr{Ptr{ArbFloat}} # arf_ptr* - -arf_clear(x::ArbFloat{P}) where {P} = ccall(@libarb(arf_clear), Cvoid, (Ref{ArbFloat},), x) - ArbFloat{P}(x::ArbFloat{P}) where {P} = x ArbFloat(x::ArbFloat{P}) where {P} = x @@ -125,17 +104,33 @@ function ArbFloat{P}(x::Irrational{S}) where {P,S} return z end -function (::Type{T})(x::ArbFloat, roundingmode::RoundingMode=RoundNearest) where {T<:ArbInts} +function (::Type{T})(x::ArbFloat, roundingmode::RoundingMode) where {T<:ArbInts} + if !(typemin(T) <= x <= typemax(T)) + throw(InexactError(nameof(T), T, x)) + # attention: ccall segfaults, if result would become too large for Slong + end rounding = match_rounding_mode(roundingmode) - if !( typemin(T) <= x <= typemax(T)) + z = ccall(@libarb(arf_get_si), Slong, (Ref{ArbFloat}, Cint), x, rounding) + return convert(T, z) +end +function (::Type{T})(x::ArbFloat) where {T<:ArbInts} + if !( isinteger(x) && typemin(T) <= x <= typemax(T)) throw(InexactError(nameof(T), T, x)) # attention: ccall segfaults, if result would become too large for Slong end + rounding = match_rounding_mode(RoundNearest) z = ccall(@libarb(arf_get_si), Slong, (Ref{ArbFloat}, Cint), x, rounding) return convert(T, z) end -function (::Type{T})(x::ArbFloat{P}, roundingmode::RoundingMode=RoundNearest) where {P,T<:Integer} +function (::Type{T})(x::ArbFloat{P}, roundingmode::RoundingMode) where {P,T<:Integer} + y = round(x, roundingmode) + z = BigFloat(y; precision=P + 2) + return convert(T, z) +end +function (::Type{T})(x::ArbFloat{P}) where {P,T<:Integer} + !isinteger(x) && throw(InexactError(nameof(T), T, x)) + roundingmode = RoundNearest y = round(x, roundingmode) z = BigFloat(y; precision=P + 2) return convert(T, z) diff --git a/src/libarb/ArbFloatMatrix.jl b/src/libarb/ArbFloatMatrix.jl index 1c6007c2..bf322795 100644 --- a/src/libarb/ArbFloatMatrix.jl +++ b/src/libarb/ArbFloatMatrix.jl @@ -9,7 +9,7 @@ mutable struct ArbFloatMatrix{P} <: AbstractArbMatrix{P, ArbFloat} arbrealmatrix::ArbRealMatrix{P} - + function ArbFloatMatrix{P}(rowcount::Int, colcount::Int) where {P} z = ArbRealMatrix{P}(rowcount, colcount) return new{P}(z) @@ -166,7 +166,7 @@ function transpose(x::Array{ArbFloat,2}) y = transpose(Array{ArbReal,2}(x)) z = Array{ArbFloat,2}(undef, size(y)) z[:,:] = y[:,:] - return (Transpose{ArbFloat,Array{ArbFloat,2}}).(z) + return (Transpose{ArbFloat,Array{ArbFloat,2}}).(z) end =# @@ -175,14 +175,14 @@ end if Sys.CPU_THREADS <= 2 function mul!(z::ArbFloatMatrix{P}, x::ArbFloatMatrix{P}, y::ArbFloatMatrix{P}) where {P} - ccall(@libarb(arb_mat_approx_mul), Cvoid, (Ref{ArbRealMatrix}, Ref{ArbRealMatrix}, Ref{ArbRealMatrix}, Cint), + ccall(@libarb(arb_mat_approx_mul), Cvoid, (Ref{ArbRealMatrix}, Ref{ArbRealMatrix}, Ref{ArbRealMatrix}, Cint), z.arbrealmatrix, x.arbrealmatrix, y.arbrealmatrix, P) return nothing end function matmul(x::ArbFloatMatrix{P}, y::ArbFloatMatrix{P}) where {P} z = ArbFloatMatrix{P}(rowcount(x), colcount(y)) - ccall(@libarb(arb_mat_approx_mul), Cvoid, (Ref{ArbRealMatrix}, Ref{ArbRealMatrix}, Ref{ArbRealMatrix}, Cint), + ccall(@libarb(arb_mat_approx_mul), Cvoid, (Ref{ArbRealMatrix}, Ref{ArbRealMatrix}, Ref{ArbRealMatrix}, Cint), z.arbrealmatrix, x.arbrealmatrix, y.arbrealmatrix, P) return ArbFloat{P}.(Matrix(z)) end @@ -190,14 +190,14 @@ if Sys.CPU_THREADS <= 2 else function mul!(z::ArbFloatMatrix{P}, x::ArbFloatMatrix{P}, y::ArbFloatMatrix{P}) where {P} - ccall(@libarb(arb_mat_mul_threaded), Cvoid, (Ref{ArbRealMatrix}, Ref{ArbRealMatrix}, Ref{ArbRealMatrix}, Cint), + ccall(@libarb(arb_mat_mul_threaded), Cvoid, (Ref{ArbRealMatrix}, Ref{ArbRealMatrix}, Ref{ArbRealMatrix}, Cint), z.arbrealmatrix, x.arbrealmatrix, y.arbrealmatrix, P) return nothing end function matmul(x::ArbFloatMatrix{P}, y::ArbFloatMatrix{P}) where {P} z = ArbFloatMatrix{P}(rowcount(x), colcount(y)) - ccall(@libarb(arb_mat_mul_threaded), Cvoid, (Ref{ArbRealMatrix}, Ref{ArbRealMatrix}, Ref{ArbRealMatrix}, Cint), + ccall(@libarb(arb_mat_mul_threaded), Cvoid, (Ref{ArbRealMatrix}, Ref{ArbRealMatrix}, Ref{ArbRealMatrix}, Cint), z.arbrealmatrix, x.arbrealmatrix, y.arbrealmatrix, P) return ArbFloat{P}.(Matrix(z)) end @@ -230,7 +230,7 @@ end @inline function Base.:(*)(x::Array{ArbFloat,2}, y::Array{ArbFloat,2}) checkmulable(x, y) - P = workingprecision(ArbFloat) + P = workingprecision(ArbFloat) return matmul(ArbFloatMatrix{P}(x), ArbFloatMatrix{P}(y)) end @@ -269,9 +269,9 @@ end function exp(x::Array{ArbFloat{P},2}) where {P} checksquare(x) - y = ArbFloatMatrix{P}(x) + y = ArbFloatMatrix{P}(x) z = ArbFloatMatrix{P}(rowcount(x), colcount(x)) - ccall(@libarb(arb_mat_exp), Cvoid, (Ref{ArbRealMatrix}, Ref{ArbRealMatrix}, Cint), + ccall(@libarb(arb_mat_exp), Cvoid, (Ref{ArbRealMatrix}, Ref{ArbRealMatrix}, Cint), z.arbrealmatrix, y.arbrealmatrix, P) return ArbFloat{P}.(Matrix(z)) end @@ -279,18 +279,18 @@ end function exp(x::Array{ArbFloat,2}) checksquare(x) P = workingprecision(ArbFloat) - y = ArbFloatMatrix{P}(x) + y = ArbFloatMatrix{P}(x) z = ArbFloatMatrix{P}(rowcount(x), colcount(x)) - ccall(@libarb(arb_mat_exp), Cvoid, (Ref{ArbRealMatrix}, Ref{ArbRealMatrix}, Cint), + ccall(@libarb(arb_mat_exp), Cvoid, (Ref{ArbRealMatrix}, Ref{ArbRealMatrix}, Cint), z.arbrealmatrix, y.arbrealmatrix, P) return ArbFloat{P}.(Matrix(z)) end #= function exp(x::ArbFloatMatrix{P}) where {P} - y = ArbFloatMatrix{P}(x) + y = ArbFloatMatrix{P}(x) z = ArbFloatMatrix{P}(rowcount(x), colcount(x)) - ccall(@libarb(arb_mat_exp), Cvoid, (Ref{ArbRealMatrix}, Ref{ArbRealMatrix}, Cint), + ccall(@libarb(arb_mat_exp), Cvoid, (Ref{ArbRealMatrix}, Ref{ArbRealMatrix}, Cint), z.arbrealmatrix, y.arbrealmatrix, P) return ArbFloat{P}.(Matrix(z)) end diff --git a/src/libarb/ArbReal.jl b/src/libarb/ArbReal.jl index 8a21d95d..8f2e314c 100644 --- a/src/libarb/ArbReal.jl +++ b/src/libarb/ArbReal.jl @@ -20,33 +20,6 @@ of the resultant width required, and to check the radius of results. See also: [`ArbFloat`](@ref), [`ArbComplex`](@ref), [`ball`](@ref), [`setball`](@ref), [`midpoint`](@ref), [`radius`](@ref) """ ArbReal -# ArbReal structs hold balls given as midpoint, radius - -mutable struct ArbReal{P} <: Real # P is the precision in bits - # midpoint - mid_exp::Int # fmpz exponent of 2 (2^exp) - mid_size::UInt # mp_size_t nwords and sign (lsb holds sign of significand) - mid_d1::UInt # significand unsigned, immediate value or the initial span - mid_d2::UInt # (d1, d2) the final part indicating the significand, or 0 - # radius - rad_exp::Int # fmpz exponent of 2 (2^exp) - rad_man::UInt # mp_limb_t radius, unsigned by definition - - function ArbReal{P}() where {P} - minprec(P, ArbReal) - z = new{P}(0,0,0,0,0,0) - ccall(@libarb(arb_init), Cvoid, (Ref{ArbReal},), z) - finalizer(arb_clear, z) - return z - end -end - -# for use within structs, e.g. ArbRealMatrix -const PtrToArbReal = Ptr{ArbReal} # arb_ptr -const PtrToPtrToArbReal = Ptr{Ptr{ArbReal}} # arb_ptr* - -arb_clear(x::ArbReal{P}) where {P} = ccall(@libarb(arb_clear), Cvoid, (Ref{ArbReal},), x) - ArbReal{P}(x::ArbReal{P}) where {P} = x ArbReal(x::ArbReal{P}) where {P} = ArbReal{P}(x) diff --git a/src/libarb/ArbRealMatrix.jl b/src/libarb/ArbRealMatrix.jl index b4f0032c..6581287b 100644 --- a/src/libarb/ArbRealMatrix.jl +++ b/src/libarb/ArbRealMatrix.jl @@ -231,7 +231,7 @@ end # matrix multiply - + function mul!(z::ArbRealMatrix{P}, x::ArbRealMatrix{P}, y::ArbRealMatrix{P}) where {P} ccall(@libarb(arb_mat_mul_threaded), Cvoid, (Ref{ArbRealMatrix}, Ref{ArbRealMatrix}, Ref{ArbRealMatrix}, Cint), z, x, y, P) return nothing @@ -277,7 +277,7 @@ end @inline function Base.:(*)(x::Array{ArbReal,2}, y::Array{ArbReal,2}) checkmulable(x, y) - P = workingprecision(ArbReal) + P = workingprecision(ArbReal) return matmul(ArbRealMatrix{P}(x), ArbRealMatrix{P}(y)) end @@ -481,9 +481,9 @@ function matmul(x::ArbRealMatrix{P}, y::ArbRealMatrix{P}) where {P} return Matrix(z) end =# - + function exp(x::Array{ArbReal, 2}) - checksquare(x) + checksquare(x) P = workingprecision(ArbReal) y = ArbRealMatrix(x) z = ArbRealMatrix{P}(rowcount(x), colcount(x)) diff --git a/src/support/complex.jl b/src/support/complex.jl index 6a96865b..26b6f4fa 100644 --- a/src/support/complex.jl +++ b/src/support/complex.jl @@ -1,33 +1,37 @@ -float(::Type{T}) where {T<:ArbNumber1} = ArbFloat -float(::Type{T}) where {P,T<:ArbNumber1{P}} = ArbFloat{P} +float(::Type{T}) where {T<:ArbNumber} = ArbFloat +float(::Type{T}) where {P,T<:ArbNumber{P}} = ArbFloat{P} -float(x::ArbNumber1{P}) where {P} = convert(ArbFloat{P}, x) # TODO senseless for Complex? +float(x::ArbNumber{P}) where {P} = convert(ArbFloat{P}, x) # TODO senseless for Complex? real(::Type{ArbFloat}) = ArbFloat -real(::Type{T}) where {T<:ArbNumber1} = ArbReal +real(::Type{T}) where {T<:ArbNumber} = ArbReal real(::Type{ArbFloat{P}}) where {P} = ArbFloat{P} -real(::Type{T}) where {P,T<:ArbNumber1{P}} = ArbReal{P} +real(::Type{T}) where {P,T<:ArbNumber{P}} = ArbReal{P} -real(x::ArbNumber1) = x +real(x::ArbNumber) = x function real(x::ArbComplex{P}) where {P} z = ArbReal{P}() ccall(@libarb(acb_get_real), Cvoid, (Ref{ArbReal}, Ref{ArbComplex}), z, x) return z end -imag(::Type{T}) where {T<:ArbNumber1} = real(T) +imag(::Type{T}) where {T<:ArbNumber} = real(T) -imag(::T) where {T<:ArbNumber1} = zero(T) +imag(::T) where {T<:ArbNumber} = zero(T) function imag(x::ArbComplex{P}) where {P} z = ArbReal{P}() ccall(@libarb(acb_get_imag), Cvoid, (Ref{ArbReal}, Ref{ArbComplex}), z, x) return z end +function isreal(x::ArbComplex) + iszero(imag(x)) +end + # other parts -angle(x::T) where {T<:ArbNumber1} = signbit(x) ? T(pi) : zero(T) +angle(x::T) where {T<:ArbNumber} = signbit(x) ? T(pi) : zero(T) function angle(x::ArbComplex{P}) where {P} z = ArbReal{P}() ccall(@libarb(acb_arg), Cvoid, (Ref{ArbReal}, Ref{ArbComplex}, Slong), z, x, P) diff --git a/src/values/constructors.jl b/src/values/constructors.jl index 827c2fe0..a7fbcc8f 100644 --- a/src/values/constructors.jl +++ b/src/values/constructors.jl @@ -1,112 +1,60 @@ # IEEEFloat # rounds up (widens) -function Float64(x::Mag) +function (::Type{T})(x::Mag) where {T<:IEEEFloat} z = ccall(@libarb(mag_get_d), Float64, (Ref{Mag},), x) - return z + return convert(T, z) end -function Float64(x::ArbFloat{P}) where {P} - rounding = match_rounding_mode(RoundNearest) +function (::Type{T})(x::ArbFloat{P}, rm::RoundingMode=RoundNearest) where {P,T<:IEEEFloat} + rounding = match_rounding_mode(rm) z = ccall(@libarb(arf_get_d), Float64, (Ref{ArbFloat}, Cint), x, rounding) - return z + return convert(T, z) end -Float32(x::ArbFloat{P}) where {P} = Float32(Float64(x)) -Float16(x::ArbFloat{P}) where {P} = Float16(Float64(x)) -function Float64(x::ArbReal{P}) where {P} - x = midpoint(x) - y = ArbFloat{P}(x) - Float64(y) -end -Float32(x::ArbReal{P}) where {P} = Float32(Float64(x)) -Float16(x::ArbReal{P}) where {P} = Float16(Float64(x)) - -function Float64(x::ArbComplex{P}) where {P} - x = midpoint(real(x)) - y = ArbFloat{P}(x) - Float64(y) -end -Float32(x::ArbComplex{P}) where {P} = Float32(Float64(x)) -Float16(x::ArbComplex{P}) where {P} = Float16(Float64(x)) - -function Float64(x::ArbFloat{P}, roundingmode::RoundingMode) where {P} - rounding = match_rounding_mode(roundingmode) - z = ccall(@libarb(arf_get_d), Float64, (Ref{ArbFloat}, Cint), x, rounding) - return z +function (::Type{T})(x::ArbReal{P}, rm::RoundingMode=RoundNearest) where {P,T<:IEEEFloat} + T(ArbFloat(midpoint(x)), rm) end -Float32(x::ArbFloat{P}, roundingmode::RoundingMode) where {P} = Float32(Float64(x, roundingmode)) -Float16(x::ArbFloat{P}, roundingmode::RoundingMode) where {P} = Float16(Float64(x, roundingmode)) -function Float64(x::ArbReal{P}, roundingmode::RoundingMode) where {P} - x = midpoint(x) - y = ArbFloat{P}(x) - Float64(y, roundingmode) +function (::Type{T})(x::ArbComplex{P}, rm::RoundingMode=RoundNearest) where {P,T<:IEEEFloat} + isreal(x) || throw(InexactError(nameof(T), T, x)) + T(ArbFloat(midpoint(real(x))), rm) end -Float32(x::ArbReal{P}, roundingmode::RoundingMode) where {P} = Float32(Float64(x, roundingmode)) -Float16(x::ArbReal{P}, roundingmode::RoundingMode) where {P} = Float16(Float64(x, roundingmode)) -function Float64(x::ArbComplex{P}, roundingmode::RoundingMode) where {P} - x = midpoint(real(x)) - y = ArbFloat{P}(x) - Float64(y, roundingmode) +function Complex{T}(x::ArbComplex{P}, rm::RoundingMode=RoundNearest) where {P,T<:Real} + re = T(real(x), rm) + im = T(imag(x), rm) + return Complex{T}(re, im) end -Float32(x::ArbComplex{P}, roundingmode::RoundingMode) where {P} = Float32(Float64(x, roundingmode)) -Float16(x::ArbComplex{P}, roundingmode::RoundingMode) where {P} = Float16(Float64(x, roundingmode)) - - -function Complex{T}(x::ArbComplex{P}, roundingmode::RoundingMode) where {P, T<:IEEEFloat} - re = T(real(x), roundingmode) - im = T(imag(x), roundingmode) - return Complex{T}(re, im) -end - -Complex{T}(x::ArbComplex{P}) where {P, T<:IEEEFloat} = Complex{T}(x, RoundNearest) -Complex(x::ArbComplex{P}) where {P} = Complex{Float64}(x, RoundNearest) # integers -for I in (:Int8, :Int16, :Int32, :Int64, :Int128) - @eval begin - - function $I(x::ArbReal{P}) where {P} - (!isexact(x) || !isinteger(x)) && throw(InexactError(nameof($I), $I, x)) - bi = BigInt(BigFloat(x)) - !(typemin($I) <= bi <= typemax($I)) && throw(InexactError(nameof($I), $I, x)) - return $I(bi) - end - - function $I(x::ArbComplex{P}) where {P} - (!isexact(x) | !isinteger(x) | !iszero(imag(x))) && throw(InexactError(nameof($I), $I, x)) - bi = BigInt(BigFloat(x)) - !(typemin($I) <= bi <= typemax($I)) && throw(InexactError(nameof($I), $I, x)) - return $I(bi) - end +function (::Type{T})(x::ArbReal) where {T<:Integer} + isexact(x) || throw(InexactError(nameof(T), T, x)) + T(ArbFloat(x)) +end - end +function (::Type{T})(x::ArbComplex) where {T<:Integer} + (isexact(x) && iszero(imag(x))) || throw(InexactError(nameof(T), T, x)) + T(ArbFloat(x)) end # Rational -1#ArbFloat(x::T) where {S, T<:Rational{S}} = ArbFloat(x.num)/ArbFloat(x.den) +#ArbFloat(x::T) where {S, T<:Rational{S}} = ArbFloat(x.num)/ArbFloat(x.den) #ArbReal(x::T) where {S, T<:Rational{S}} = ArbReal(x.num)/ArbReal(x.den) #ArbComplex(x::T) where {S, T<:Rational{S}} = ArbComplex(ArbReal(x), ArbReal(0)) #ArbComplex(x::Rational, y::Rational) = ArbComplex(ArbReal(x), ArbReal(y)) # Irrational -function ArbReal{P}(x::Irrational{S}) where {P,S} - mid = ArbFloat{P}(x) - rad = ulp(mid) - return setball(mid, rad) +function ArbReal{P}(x::Irrational) where {P} + mid = ArbFloat{P}(x) + rad = ulp(mid) + return setball(mid, rad) end -function ArbReal(x::Irrational{S}) where {S} - P = workingprecision(ArbReal) - mid = ArbFloat{P}(x) - rad = ulp(mid) - return setball(mid, rad) -end +ArbReal(x::Irrational) = ArbReal{workingprecision(ArbReal)}(x) # fallback @@ -139,7 +87,7 @@ function ArbFloat{P}(x::ArbFloat{Q}, roundingmode::RoundingMode) where {P,Q} rounding = match_rounding_mode(roundingmode) z = ArbFloat{P}() rnd = ccall(@libarb(arf_set_round), Cint, - (Ref{ArbFloat}, Ref{ArbFloat}, Clong, Cint), z, x, P, rounding) + (Ref{ArbFloat}, Ref{ArbFloat}, Clong, Cint), z, x, P, rounding) return z end ArbFloat{P}(x::ArbFloat{Q}) where {P,Q} = ArbFloat{P}(x, RoundNearest) @@ -148,7 +96,7 @@ function ArbReal{P}(x::ArbReal{Q}) where {P,Q} minprec(P, ArbReal) z = ArbReal{P}() ccall(@libarb(arb_set_round), Cvoid, - (Ref{ArbReal}, Ref{ArbReal}, Clong), z, x, P) + (Ref{ArbReal}, Ref{ArbReal}, Clong), z, x, P) return z end @@ -156,12 +104,12 @@ function ArbComplex{P}(x::ArbComplex{Q}) where {P,Q} minprec(P, ArbComplex) z = ArbComplex{P}() ccall(@libarb(acb_set_round), Cvoid, - (Ref{ArbComplex}, Ref{ArbComplex}, Clong), z, x, P) + (Ref{ArbComplex}, Ref{ArbComplex}, Clong), z, x, P) return z end # widen -widen(::Type{ArbFloat{P}}) where {P} = ArbFloat{P+P+4} -widen(::Type{ArbReal{P}}) where {P} = ArbReal{P+P+4} -widen(::Type{ArbComplex{P}}) where {P} = ArbComplex{P+P+4} +widen(::Type{ArbFloat{P}}) where {P} = ArbFloat{P + P + 4} +widen(::Type{ArbReal{P}}) where {P} = ArbReal{P + P + 4} +widen(::Type{ArbComplex{P}}) where {P} = ArbComplex{P + P + 4} diff --git a/test/complex.jl b/test/complex.jl index 871d3a81..aede8b2d 100644 --- a/test/complex.jl +++ b/test/complex.jl @@ -114,6 +114,7 @@ @test abs2(api) ≈ abs2(real(api)) + abs2(imag(api)) + @test widen(ArbComplex{24}) == ArbComplex{52} end TYPE_TESTS = [ @@ -138,7 +139,7 @@ dprec(a::ArbNumber) = sprec(a) dprec(::Number) = INF - TEST_TYPES = (Int8, Int, UInt32, Float16, Float64, BigFloat, ArbFloat{200}, ArbReal{200}) + TEST_TYPES = (Int8, Int, UInt32, Float16, Float32, Float64, BigFloat, ArbFloat{200}, ArbReal{200}) @testset "ArbComplex(::$T)" for T in (TEST_TYPES..., ArbComplex{200}) a = T(10) @@ -158,10 +159,14 @@ DEF == INF && (DEF = sprec(a)) c = ArbComplex(a, b) @test c == 10 + 20im - @test typeof(c) == ArbComplex{DEF} || a isa ArbNumber || b isa ArbNumber # TODO should be true also for ArbNumbers + @test typeof(c) == ArbComplex{DEF} || a isa ArbNumber || b isa ArbNumber # TODO should be true also for ArbNumber d = ArbComplex{250}(a, b) @test typeof(d) == ArbComplex{250} + + @test_throws InexactError Float64(c) + @test isreal(c) == false + @test isreal(ArbComplex(pi)) == true end @testset "ArbComplex(::Int, ::ArbComplex)" begin diff --git a/test/misc.jl b/test/misc.jl index 6d989fa5..303a1c35 100644 --- a/test/misc.jl +++ b/test/misc.jl @@ -9,7 +9,7 @@ @test_throws InexactError ArbNumerics.Mag(-1) end - @testset "$T($v)" for T in (ArbReal, ArbReal{64}, ArbReal, ArbReal{64}), + @testset "$T($v)" for T in (ArbFloat, ArbFloat{64}, ArbReal, ArbReal{64}), (v, c, f, t) in [(0.6, 1, 0, 0), (-0.6, 0, -1, 0)] x = T(v) @@ -54,15 +54,28 @@ end @testset "convert ArbFloat to Integer" begin - v = typemax(Int32) + 1.5 + v = typemax(Int32) + 1.51 w = round(v) a = ArbFloat(v) @test_throws InexactError Int32(a) + @test_throws InexactError Int64(a) == w + @test_throws InexactError Int128(a) == w + @test_throws InexactError BigInt(a) == w + + @test_throws InexactError Int32(a, RoundUp) == w + @test Int64(a, RoundDown) == w - 1 + @test Int128(a, RoundNearest) == w + @test Int128(a, RoundUp) == w + @test BigInt(a, RoundToZero) == w - 1 + + @test_throws InexactError Integer(a) + + a = ArbFloat(w) + @test_throws InexactError Int32(a) == w @test Int64(a) == w @test Int128(a) == w @test BigInt(a) == w - @test_throws InexactError Integer(a) # TODO inconsistent with Int64, BigInt, ... v = typemin(Int32) - 1 a = ArbFloat(v) x = Integer(a) @@ -94,6 +107,8 @@ @test hash(a) isa UInt @test_broken hash(b) == hash(1) # because isequal(b, 1), see ?hash + + @test widen(ArbFloat{24}) == ArbFloat{52} end @testset "ArbReal" begin @@ -108,6 +123,11 @@ @test copy(a) == a @test copy(a) !== a + ma = ArbNumerics.Mag(a) + @test Float64(ma) ≈ Float64(pi) + @test Float32(a) ≈ Float32(pi) + @test Float16(a) ≈ Float16(pi) + @test ArbReal(typemax(UInt)) == typemax(UInt) @test ArbReal(typemax(Int128)) == typemax(Int128) @test ArbReal(1 // 3) ≈ ArbReal(1) / ArbReal(3) @@ -159,6 +179,8 @@ @test hash(a) isa UInt @test_broken hash(b) == hash(1) # because isequal(b, 1), see ?hash + + @test widen(ArbReal{24}) == ArbReal{52} end end From 0ca5d7faf33ba04ec8bd68c4dd1fe978129da75f Mon Sep 17 00:00:00 2001 From: KlausC Date: Wed, 8 Nov 2023 10:40:26 +0100 Subject: [PATCH 14/15] coverage --- src/ArbNumerics.jl | 2 +- src/float/arith.jl | 76 ++----------------------------------- src/float/prearith.jl | 18 ++++++++- src/libarb/ArbFloat.jl | 12 +++--- src/libarb/ArbReal.jl | 1 - src/libarb/roundingmodes.jl | 15 +++++++- src/libarb/string.jl | 4 +- src/values/constructors.jl | 33 ++++++++-------- src/values/intraconvert.jl | 34 +---------------- test/complex.jl | 7 ++-- 10 files changed, 64 insertions(+), 138 deletions(-) diff --git a/src/ArbNumerics.jl b/src/ArbNumerics.jl index d0e00cd1..ada4b81f 100644 --- a/src/ArbNumerics.jl +++ b/src/ArbNumerics.jl @@ -37,7 +37,7 @@ export ArbNumber, inf, posinf, neginf, nan, typemax, typemin, floatmax, floatmin, - magnitude, # complex magnitude, `angle` gives phase + magnitude, csign, # complex magnitude, `angle` gives phase signs, signbits, significand_bits, relerror_bits, relaccuracy_bits, diff --git a/src/float/arith.jl b/src/float/arith.jl index ea23e8ff..9cd55ec0 100644 --- a/src/float/arith.jl +++ b/src/float/arith.jl @@ -46,80 +46,12 @@ for (A,F) in ((:(+), :acb_add), (:(-), :acb_sub), (:(*), :acb_mul), (:(/), :acb_ end end -(+)(x::ArbComplex{P}, y::ArbReal{P}) where {P} = (+)(promote(x, y)...,) -(-)(x::ArbComplex{P}, y::ArbReal{P}) where {P} = (-)(promote(x, y)...,) -(*)(x::ArbComplex{P}, y::ArbReal{P}) where {P} = (*)(promote(x, y)...,) -(/)(x::ArbComplex{P}, y::ArbReal{P}) where {P} = (/)(promote(x, y)...,) -(^)(x::ArbComplex{P}, y::ArbReal{P}) where {P} = (^)(promote(x, y)...,) -(+)(x::ArbReal{P}, y::ArbComplex{P}) where {P} = (+)(promote(x, y)...,) -(-)(x::ArbReal{P}, y::ArbComplex{P}) where {P} = (-)(promote(x, y)...,) -(*)(x::ArbReal{P}, y::ArbComplex{P}) where {P} = (*)(promote(x, y)...,) -(/)(x::ArbReal{P}, y::ArbComplex{P}) where {P} = (/)(promote(x, y)...,) -(^)(x::ArbReal{P}, y::ArbComplex{P}) where {P} = (^)(promote(x, y)...,) - -(+)(x::ArbComplex{P}, y::ArbFloat{P}) where {P} = (+)(promote(x, y)...,) -(-)(x::ArbComplex{P}, y::ArbFloat{P}) where {P} = (-)(promote(x, y)...,) -(*)(x::ArbComplex{P}, y::ArbFloat{P}) where {P} = (*)(promote(x, y)...,) -(/)(x::ArbComplex{P}, y::ArbFloat{P}) where {P} = (/)(promote(x, y)...,) -(^)(x::ArbComplex{P}, y::ArbFloat{P}) where {P} = (^)(promote(x, y)...,) -(+)(x::ArbFloat{P}, y::ArbComplex{P}) where {P} = (+)(promote(x, y)...,) -(-)(x::ArbFloat{P}, y::ArbComplex{P}) where {P} = (-)(promote(x, y)...,) -(*)(x::ArbFloat{P}, y::ArbComplex{P}) where {P} = (*)(promote(x, y)...,) -(/)(x::ArbFloat{P}, y::ArbComplex{P}) where {P} = (/)(promote(x, y)...,) -(^)(x::ArbFloat{P}, y::ArbComplex{P}) where {P} = (^)(promote(x, y)...,) - -(+)(x::ArbReal{P}, y::ArbFloat{P}) where {P} = (+)(promote(x, y)...,) -(-)(x::ArbReal{P}, y::ArbFloat{P}) where {P} = (-)(promote(x, y)...,) -(*)(x::ArbReal{P}, y::ArbFloat{P}) where {P} = (*)(promote(x, y)...,) -(/)(x::ArbReal{P}, y::ArbFloat{P}) where {P} = (/)(promote(x, y)...,) -(^)(x::ArbReal{P}, y::ArbFloat{P}) where {P} = (^)(promote(x, y)...,) -(+)(x::ArbFloat{P}, y::ArbReal{P}) where {P} = (+)(promote(x, y)...,) -(-)(x::ArbFloat{P}, y::ArbReal{P}) where {P} = (-)(promote(x, y)...,) -(*)(x::ArbFloat{P}, y::ArbReal{P}) where {P} = (*)(promote(x, y)...,) -(/)(x::ArbFloat{P}, y::ArbReal{P}) where {P} = (/)(promote(x, y)...,) -(^)(x::ArbFloat{P}, y::ArbReal{P}) where {P} = (^)(promote(x, y)...,) - -(+)(x::Mag, y::ArbFloat{P}) where {P} = (+)(promote(x, y)...,) -(-)(x::Mag, y::ArbFloat{P}) where {P} = (-)(promote(x, y)...,) -(*)(x::Mag, y::ArbFloat{P}) where {P} = (*)(promote(x, y)...,) -(/)(x::Mag, y::ArbFloat{P}) where {P} = (/)(promote(x, y)...,) -(^)(x::Mag, y::ArbFloat{P}) where {P} = (^)(promote(x, y)...,) -(+)(x::ArbFloat{P}, y::Mag) where {P} = (+)(promote(x, y)...,) -(-)(x::ArbFloat{P}, y::Mag) where {P} = (-)(promote(x, y)...,) -(*)(x::ArbFloat{P}, y::Mag) where {P} = (*)(promote(x, y)...,) -(/)(x::ArbFloat{P}, y::Mag) where {P} = (/)(promote(x, y)...,) -(^)(x::ArbFloat{P}, y::Mag) where {P} = (^)(promote(x, y)...,) - -(+)(x::Mag, y::ArbReal{P}) where {P} = (+)(promote(x, y)...,) -(-)(x::Mag, y::ArbReal{P}) where {P} = (-)(promote(x, y)...,) -(*)(x::Mag, y::ArbReal{P}) where {P} = (*)(promote(x, y)...,) -(/)(x::Mag, y::ArbReal{P}) where {P} = (/)(promote(x, y)...,) -(^)(x::Mag, y::ArbReal{P}) where {P} = (^)(promote(x, y)...,) -(+)(x::ArbReal{P}, y::Mag) where {P} = (+)(promote(x, y)...,) -(-)(x::ArbReal{P}, y::Mag) where {P} = (-)(promote(x, y)...,) -(*)(x::ArbReal{P}, y::Mag) where {P} = (*)(promote(x, y)...,) -(/)(x::ArbReal{P}, y::Mag) where {P} = (/)(promote(x, y)...,) -(^)(x::ArbReal{P}, y::Mag) where {P} = (^)(promote(x, y)...,) - -for F in (:(+), :(-), :(*), :(/), :(^)) - @eval begin - ($F)(x::ArbFloat{P}, y::Integer) where {P} = ($F)(promote(x, y)...,) - ($F)(x::ArbReal{P}, y::Integer) where {P} = ($F)(promote(x, y)...,) - ($F)(x::ArbComplex{P}, y::Integer) where {P} = ($F)(promote(x, y)...,) - ($F)(x::ArbFloat{P}, y::AbstractFloat) where {P} = ($F)(promote(x, y)...,) - ($F)(x::ArbReal{P}, y::AbstractFloat) where {P} = ($F)(promote(x, y)...,) - ($F)(x::ArbComplex{P}, y::AbstractFloat) where {P} = ($F)(promote(x, y)...,) - ($F)(x::ArbComplex{P}, y::Complex) where {P} = ($F)(promote(x, y)...,) - ($F)(x::Integer, y::ArbFloat{P}) where {P} = ($F)(promote(x, y)...,) - ($F)(x::Integer, y::ArbReal{P}) where {P} = ($F)(promote(x, y)...,) - ($F)(x::Integer, y::ArbComplex{P}) where {P} = ($F)(promote(x, y)...,) - ($F)(x::IEEEFloat, y::ArbFloat{P}) where {P} = ($F)(promote(x, y)...,) - ($F)(x::IEEEFloat, y::ArbReal{P}) where {P} = ($F)(promote(x, y)...,) - ($F)(x::IEEEFloat, y::ArbComplex{P}) where {P} = ($F)(promote(x, y)...,) - ($F)(x::Complex, y::ArbComplex{P}) where {P} = ($F)(promote(x, y)...,) - end +function (^)(a::ArbNumber, b::Integer) + x = Base.power_by_squaring(a, abs(b)) + b < 0 ? inv(x) : x end + # interact with Bool @inline ArbFloat{P}(b::Bool) where {P} = b ? one(ArbFloat{P}) : zero(ArbFloat{P}) diff --git a/src/float/prearith.jl b/src/float/prearith.jl index 4075bc9b..1ebbe96b 100644 --- a/src/float/prearith.jl +++ b/src/float/prearith.jl @@ -40,16 +40,32 @@ end and the sign of the imaginary part when z is on the imaginary axis. =# function sign(x::ArbComplex{P}) where {P} + z = ArbComplex{P}() + ccall(@libarb(acb_sgn), Cvoid, (Ref{ArbComplex}, Ref{ArbComplex}, Slong), z, x, P) + return z +end +""" + csign(::ArbComplex) + +Return the extension of the real sign function taking the value 1 +for z strictly in the right half plane, -1 for z strictly in the left half plane, +and the sign of the imaginary part when z is on the imaginary axis. +""" +function csign(x::ArbComplex{P}) where {P} z = ArbReal{P}() ccall(@libarb(acb_csgn), Cvoid, (Ref{ArbReal}, Ref{ArbComplex}), z, x) return z end +""" + signs(z::ArbComplex) + +Return (sign(real(z)), sign(imag(z))) +""" function signs(x::ArbComplex{P}) where {P} return sign(real(x)), sign(imag(x)) end - function (-)(x::ArbFloat{P}) where {P} z = ArbFloat{P}() ccall(@libarb(arf_neg), Cvoid, (Ref{ArbFloat}, Ref{ArbFloat}), z, x) diff --git a/src/libarb/ArbFloat.jl b/src/libarb/ArbFloat.jl index ef6abaab..eaf84896 100644 --- a/src/libarb/ArbFloat.jl +++ b/src/libarb/ArbFloat.jl @@ -114,7 +114,7 @@ function (::Type{T})(x::ArbFloat, roundingmode::RoundingMode) where {T<:ArbInts} return convert(T, z) end function (::Type{T})(x::ArbFloat) where {T<:ArbInts} - if !( isinteger(x) && typemin(T) <= x <= typemax(T)) + if !(isinteger(x) && typemin(T) <= x <= typemax(T)) throw(InexactError(nameof(T), T, x)) # attention: ccall segfaults, if result would become too large for Slong end @@ -130,19 +130,17 @@ function (::Type{T})(x::ArbFloat{P}, roundingmode::RoundingMode) where {P,T<:Int end function (::Type{T})(x::ArbFloat{P}) where {P,T<:Integer} !isinteger(x) && throw(InexactError(nameof(T), T, x)) - roundingmode = RoundNearest - y = round(x, roundingmode) - z = BigFloat(y; precision=P + 2) + z = BigFloat(x, precision=P + 2) return convert(T, z) end """ - BigFloat(::ArbFloat; [precision=workingprecision(x), roundingmode=RoundNearest]) + BigFloat(::ArbFloat, [roundingmode=RoundNearest]; [precision=workingprecision(x)) Construct a `BigFloat`from an `ArbFloat`. """ -function BigFloat(x::ArbFloat; precision::Int=workingprecision(x), roundingmode::RoundingMode=RoundNearest) - rounding = match_rounding_mode(roundingmode) +function BigFloat(x::ArbFloat, roundingmode::RoundingMode=RoundNearest; precision::Int=workingprecision(x)) + rounding = mpfr_rounding_mode(roundingmode) z = BigFloat(0; precision) roundingdir = ccall(@libarb(arf_get_mpfr), Cint, (Ref{BigFloat}, Ref{ArbFloat}, Cint), z, x, rounding) return z diff --git a/src/libarb/ArbReal.jl b/src/libarb/ArbReal.jl index 8f2e314c..97a4b617 100644 --- a/src/libarb/ArbReal.jl +++ b/src/libarb/ArbReal.jl @@ -70,7 +70,6 @@ end ArbReal{P}(x::Rational{T}) where {P, T<:Signed} = ArbReal{P}(BigFloat(x; precision=P + 2)) -BigFloat(x::ArbReal{P}) where {P} = BigFloat(ArbFloat{P}(x)) BigInt(x::ArbReal{P}) where {P} = BigInt(trunc(BigFloat(ArbFloat{P}(x)))) function Base.Integer(x::ArbReal) diff --git a/src/libarb/roundingmodes.jl b/src/libarb/roundingmodes.jl index a4a38e66..263624c9 100644 --- a/src/libarb/roundingmodes.jl +++ b/src/libarb/roundingmodes.jl @@ -27,23 +27,34 @@ const MpfrRoundFromZero = Cint(4) # RoundingMode{:FromZero}() # matched with Julia's terminology +#= match_rounding_mode(::Type{Val{ArbRoundNearest}}) = RoundNearest match_rounding_mode(::Type{Val{ArbRoundDown}}) = RoundDown match_rounding_mode(::Type{Val{ArbRoundUp}}) = RoundUp match_rounding_mode(::Type{Val{ArbRoundToZero}}) = RoundToZero match_rounding_mode(::Type{Val{ArbRoundFromZero}}) = RoundFromZero +=# match_rounding_mode(::Type{Val{RoundNearest}}) = ArbRoundNearest match_rounding_mode(::Type{Val{RoundDown}}) = ArbRoundDown match_rounding_mode(::Type{Val{RoundUp}}) = ArbRoundUp match_rounding_mode(::Type{Val{RoundToZero}}) = ArbRoundToZero match_rounding_mode(::Type{Val{RoundFromZero}}) = ArbRoundFromZero -match_rounding_mode(::Type{Val{RoundNearestTiesAway}}) = ArbRoundNearest -match_rounding_mode(::Type{Val{RoundNearestTiesUp}}) = throw(ErrorException("RoundNearestTiesUp is not supported")) +match_rounding_mode(::Type{Val{RoundNearestTiesAway}}) = ArbRoundFromZero +match_rounding_mode(::Type{Val{RoundNearestTiesUp}}) = ArbRoundUp +mpfr_rounding_mode(::Type{Val{RoundNearest}}) = MpfrRoundNearest +mpfr_rounding_mode(::Type{Val{RoundDown}}) = MpfrRoundDown +mpfr_rounding_mode(::Type{Val{RoundUp}}) = MpfrRoundUp +mpfr_rounding_mode(::Type{Val{RoundToZero}}) = MpfrRoundToZero +mpfr_rounding_mode(::Type{Val{RoundFromZero}}) = MpfrRoundFromZero +mpfr_rounding_mode(::Type{Val{RoundNearestTiesAway}}) = MpfrRoundFromZero +mpfr_rounding_mode(::Type{Val{RoundNearestTiesUp}}) = MpfrRoundUp @inline match_rounding_mode(mode::Cint) = match_rounding_mode(Val{mode}) @inline match_rounding_mode(mode::RoundingMode{S}) where {S} = match_rounding_mode(Val{mode}) + +@inline mpfr_rounding_mode(mode::RoundingMode) = mpfr_rounding_mode(Val{mode}) diff --git a/src/libarb/string.jl b/src/libarb/string.jl index d75d8c07..4c76391c 100644 --- a/src/libarb/string.jl +++ b/src/libarb/string.jl @@ -107,8 +107,8 @@ end function arbstring(x::ArbComplex{P}, maxdigits::Int=digit_precision(P); flags::UInt = NO_FLAGS) where {P} # rea, ima = real(x), imag(x) - rea = ArbReal(x) - ima = ArbReal(x) + rea = real(x) + ima = imag(x) ima_isneg = signbit(ima) connection = ima_isneg ? " - " : " + " diff --git a/src/values/constructors.jl b/src/values/constructors.jl index a7fbcc8f..25a6635d 100644 --- a/src/values/constructors.jl +++ b/src/values/constructors.jl @@ -15,6 +15,9 @@ end function (::Type{T})(x::ArbReal{P}, rm::RoundingMode=RoundNearest) where {P,T<:IEEEFloat} T(ArbFloat(midpoint(x)), rm) end +function BigFloat(x::ArbReal{P}, rm::RoundingMode=RoundNearest) where {P} + BigFloat(convert(ArbFloat, midpoint(x)), rm) +end function (::Type{T})(x::ArbComplex{P}, rm::RoundingMode=RoundNearest) where {P,T<:IEEEFloat} isreal(x) || throw(InexactError(nameof(T), T, x)) @@ -56,29 +59,27 @@ end ArbReal(x::Irrational) = ArbReal{workingprecision(ArbReal)}(x) -# fallback - -# ArbFloat{P}(x::T) where {P,T<:Real} = ArbFloat{P}(BigFloat(x)) -# ArbFloat(x::T) where {T<:Real} = ArbFloat{workingprecision(ArbFloat)}(BigFloat(x)) -ArbFloat{P}(x::T) where {P,T<:Complex} = ArbFloat{P}(BigFloat(real(x))) -ArbFloat(x::T) where {T<:Complex} = ArbFloat{workingprecision(ArbFloat)}(BigFloat(real(x))) # TODO shouldn't that InexactError - -#ArbReal{P}(x::T) where {P,T<:Real} = ArbReal{P}(BigFloat(x)) -#ArbReal(x::T) where {T<:Real} = ArbReal{workingprecision(ArbReal)}(BigFloat(real(x))) -ArbReal{P}(x::T) where {P,T<:Complex} = ArbReal{P}(BigFloat(x)) # TODO shouldn't that InexactError -ArbReal(x::T) where {T<:Complex} = ArbReal{workingprecision(ArbReal)}(BigFloat(real(x))) - # retype ArbFloat(x::ArbReal{P}) where {P} = ArbFloat{P}(x) -ArbFloat(x::ArbComplex{P}) where {P} = ArbFloat{P}(real(x)) # TODO shouldn't that InexactError ArbReal(x::ArbFloat{P}) where {P} = ArbReal{P}(x) -ArbReal(x::ArbComplex{P}) where {P} = ArbReal{P}(real(x)) # TODO shouldn't that InexactError +ArbReal(x::ArbComplex{P}) where {P} = ArbReal{P}(x) +ArbFloat(x::ArbComplex{P}) where {P} = ArbFloat{P}(x) +ArbReal(x::Complex) = ArbReal{workingprecision(ArbReal)}(x) +ArbFloat(x::Complex) = ArbFloat{workingprecision(ArbFloat)}(x) + +function _ArbFloatReal(::Type{T}, x::ArbComplex{P}) where {P,T<:ArbFloatReal} + if isreal(x) + convert(T, real(x)) + else + throw(InexactError(:ArbReal, ArbReal, x)) + end +end ArbFloat{Q}(x::ArbReal{P}) where {P,Q} = ArbFloat{Q}(ArbReal{Q}(x)) -ArbFloat{Q}(x::ArbComplex{P}) where {P,Q} = ArbFloat{Q}(ArbReal{Q}(real(x))) ArbReal{Q}(x::ArbFloat{P}) where {P,Q} = ArbReal{Q}(ArbFloat{Q}(x)) -ArbReal{Q}(x::ArbComplex{P}) where {P,Q} = ArbReal{Q}(real(x)) +ArbReal{Q}(x::ArbComplex) where {Q} = _ArbFloatReal(ArbReal{Q}, x) +ArbFloat{Q}(x::ArbComplex) where {Q} = _ArbFloatReal(ArbFloat{Q}, x) # change precision diff --git a/src/values/intraconvert.jl b/src/values/intraconvert.jl index f1f0d3a4..6cc0a029 100644 --- a/src/values/intraconvert.jl +++ b/src/values/intraconvert.jl @@ -90,40 +90,8 @@ function ArbFloat{P}(x::ArbReal{P}, ::Type{LowerBound}) where {P} (Ref{ArbFloat}, Ref{ArbReal}, Clong), z, x, P) return z end -convert(::Type{ArbFloat{P}}, x::ArbReal{P}, ::Type{LowerBound}) where {P} = ArbFloat{P}(x) - -ArbFloat{P}(x::ArbComplex{P}) where {P} = ArbFloat{P}(x, UpperBound) -function ArbFloat{P}(x::ArbComplex{P}, ::Type{UpperBound}, ) where {P} - z = ArbFloat{P}() - ccall(@libarb(arb_get_ubound_arf), Cvoid, - (Ref{ArbFloat}, Ref{ArbComplex}, Clong), z, x, P) - return z -end - -function ArbFloat{P}(x::ArbComplex{P}, ::Type{LowerBound}) where {P} - z = ArbFloat{P}() - ccall(@libarb(arb_get_lbound_arf), Cvoid, - (Ref{ArbFloat}, Ref{ArbComplex}, Clong), z, x, P) - return z -end - -ArbFloat{P}(x::ArbComplex{P}, ::Type{Radius}) where {P} = ArbFloat{P}(x, UpperBound, Radius) -convert(::Type{ArbFloat{P}}, x::ArbComplex{P}, ::Type{Radius}) where {P} = ArbFloat{P}(x, Radius) - -function ArbFloat{P}(x::ArbComplex{P}, ::Type{UpperBound}, ::Type{Radius}) where {P} - z = ArbFloat{P}() - ccall(@libarb(arb_get_rad_ubound_arf), Cvoid, - (Ref{ArbFloat}, Ref{ArbComplex}, Clong), z, x, P) - return z -end - -function ArbFloat{P}(x::ArbComplex{P}, ::Type{LowerBound}, ::Type{Radius}) where {P} - z = ArbFloat{P}() - ccall(@libarb(arb_get_rad_lbound_arf), Cvoid, - (Ref{ArbFloat}, Ref{ArbComplex}, Clong), z, x, P) - return z -end +convert(::Type{ArbFloat{P}}, x::ArbReal{P}, ::Type{LowerBound}) where {P} = ArbFloat{P}(x) Mag(x::ArbComplex{P}) where {P} = Mag(x, UpperBound) diff --git a/test/complex.jl b/test/complex.jl index aede8b2d..0f081636 100644 --- a/test/complex.jl +++ b/test/complex.jl @@ -57,9 +57,9 @@ @test ArbComplex{53}(1 + 2im) == 1 + 2im @test signbit(ac) == false @test signbits(ac) == (false, false) - @test sign(ac) == 1.0 # TODO that is incompatible with sign(::Complex); maybe rename + @test csign(ac) == 1.0 # the name `acb_csgn` is used in Arb.jll @test signs(ac) == (1, 1) - @test_broken sign(ac) ≈ ac / abs(ac) + @test sign(ac) ≈ ac / abs(ac) @test ArbComplex(fac) == ArbComplex(rac) @test copy(ac) == deepcopy(ac) @test ArbComplex(rac, fac) == Complex(rac, fac) == complex(rac, fac) @@ -105,9 +105,10 @@ @test copysign(-api, ArbFloat(-1)) == -api @test copysign(-api, ArbReal(0.0)) == api - @test sign(api) == 1 + @test csign(api) == 1 @test signs(api) == (1, -1) @test signs(ArbComplex(0)) == (0, 0) + @test sign(api) ≈ (1 - im) * sqrt(0.5) rtol=1e-16 @test signbit(api) == false @test signbits(api) == (false, true) From 538e6e64be6fee58f6a42fa9035b94be96cd6efb Mon Sep 17 00:00:00 2001 From: KlausC Date: Thu, 9 Nov 2023 20:58:14 +0100 Subject: [PATCH 15/15] coverage `constructors.jl` --- .gitignore | 2 ++ Project.toml | 2 +- src/libarb/ArbComplex.jl | 10 ++-------- src/libarb/ArbFloat.jl | 41 +++++++++++++++++++++++++++++--------- src/libarb/ArbReal.jl | 8 ++++++-- src/values/constructors.jl | 28 ++++++++++++++++---------- test/complex.jl | 14 +++++++++++++ test/misc.jl | 22 ++++++++++++++++++-- 8 files changed, 94 insertions(+), 33 deletions(-) diff --git a/.gitignore b/.gitignore index 0d4b7fea..33484143 100644 --- a/.gitignore +++ b/.gitignore @@ -7,3 +7,5 @@ lcov.info Manifest.toml settings.json deps/deps.jl +gen/ +src/clang/arblib_api.* diff --git a/Project.toml b/Project.toml index 9e37ac60..9c32f750 100644 --- a/Project.toml +++ b/Project.toml @@ -4,7 +4,7 @@ repo = "https://github.com/JeffreySarnoff/ArbNumerics.jl.git" keywords = ["math", "floating-point", "extended-precision", "accuracy", "precision"] author = "Jeffrey Sarnoff" license = "MIT" -version = "v1.3.5" +version = "v1.4" [deps] Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" diff --git a/src/libarb/ArbComplex.jl b/src/libarb/ArbComplex.jl index e099a98a..f25e20da 100644 --- a/src/libarb/ArbComplex.jl +++ b/src/libarb/ArbComplex.jl @@ -36,17 +36,11 @@ ArbComplex{P}(::Real, ::Missing) where {P} = missing ArbComplex(x::T; kw...) where {T<:Number} = ArbComplex(real(x), imag(x); kw...) -function ArbComplex(x::T, y::Real=0; bits::Int=0, digits::Int=0, base::Int=iszero(bits) ? 10 : 2) where {T<:Real} - bits = get_bits(bits, digits, base) +function ArbComplex(x::T, y::Real; bits::Int=0, digits::Int=0, base::Int=iszero(bits) ? 10 : 2) where {T<:Real} + bits = get_bits(bits, digits, base, x, y) iszero(y) ? ArbComplex{bits}(x) : ArbComplex{bits}(x, y) end -ArbComplex(x::ArbNumber{P}) where {P} = ArbComplex{P}(x) - -ArbComplex(x::ArbFloatReal{P}, y::ArbFloatReal{S}) where {P,S} = ArbComplex{min(P, S)}(x, y) -ArbComplex(x::ArbFloatReal{P}, y::Real) where {P} = ArbComplex{P}(x, y) -ArbComplex(x::Real, y::ArbFloatReal{P}) where {P} = ArbComplex{P}(x, y) - @inline sign_bit(x::ArbComplex{P}, ::Type{RealPart}) where {P} = isodd(x.real_mid_size) @inline sign_bit(x::ArbComplex{P}, ::Type{ImagPart}) where {P} = isodd(x.imag_mid_size) diff --git a/src/libarb/ArbFloat.jl b/src/libarb/ArbFloat.jl index eaf84896..83c16277 100644 --- a/src/libarb/ArbFloat.jl +++ b/src/libarb/ArbFloat.jl @@ -21,31 +21,54 @@ ArbFloat{P}(x::ArbFloat{P}) where {P} = x ArbFloat(x::ArbFloat{P}) where {P} = x ArbFloat{P}(::Missing) where {P} = missing -ArbFloat(::Missing) = missing +ArbFloat(::Missing; kw...) = missing @inline sign_bit(x::ArbFloat{P}) where {P} = isodd(x.size) -function ArbFloat(x::Union{AbstractFloat,AbstractIrrational,Integer}; kw...) +function ArbFloat(x::Union{AbstractFloat,AbstractIrrational,Integer,ArbNumber}; kw...) _ArbFloat(x; kw...) end -function ArbFloat(x::Union{Rational}; kw...) +function ArbFloat(x::Rational; kw...) _ArbFloat(x; kw...) end -function _ArbFloat(x::Real; bits::Int=0, digits::Int=0, base::Int=iszero(bits) ? 10 : 2) - bits = get_bits(bits, digits, base) +function _ArbFloat(x::Number; bits::Int=0, digits::Int=0, base::Int=iszero(bits) ? 10 : 2) + bits = get_bits(bits, digits, base, x) ArbFloat{bits}(x) end -function get_bits(bits, digits, base) + +function get_bits(bits, digits, base, x::ArbNumber{P}, y::ArbNumber{Q}) where {P,Q} + if bits == digits == 0 + max(P, Q) + else + get_bits(bits, digits, base) + end +end +function get_bits(bits, digits, base, x::ArbNumber{P}, y::Any=0) where {P} + if bits == digits == 0 + P + else + get_bits(bits, digits, base) + end +end +function get_bits(bits, digits, base, y::Any, x::ArbNumber{P}) where {P} + if bits == digits == 0 + P + else + get_bits(bits, digits, base) + end +end + +function get_bits(bits, digits, base, ::Any...) bits > 0 && bits < MINIMUM_PRECISION_BASE2 && throw(DomainError("bit precision $bits < $MINIMUM_PRECISION_BASE2")) digits > 0 && digits < MINIMUM_PRECISION_BASE10 && throw(DomainError("digit precision $digits < $MINIMUM_PRECISION_BASE10")) if base === 10 - bits = digits > 0 ? bits4digits(digits) + extrabits() : (bits > 0 ? bits + extrabits() : DEFAULT_PRECISION.x) + bi = digits > 0 ? bits4digits(digits) + extrabits() : (bits > 0 ? bits + extrabits() : DEFAULT_PRECISION.x) elseif base === 2 - bits = bits > 0 ? bits + extrabits() : (digits > 0 ? digits + extrabits() : DEFAULT_PRECISION.x) + bi = bits > 0 ? bits + extrabits() : (digits > 0 ? digits + extrabits() : DEFAULT_PRECISION.x) else throw(ErrorException("base expects 2 or 10")) end - bits + bi end function swap!(x::ArbFloat{P}, y::ArbFloat{P}) where {P} diff --git a/src/libarb/ArbReal.jl b/src/libarb/ArbReal.jl index 97a4b617..cbb6e2ec 100644 --- a/src/libarb/ArbReal.jl +++ b/src/libarb/ArbReal.jl @@ -28,8 +28,12 @@ ArbReal(::Missing) = missing @inline sign_bit(x::ArbReal{P}) where {P} = isodd(x.mid_size) -function ArbReal(x::T; bits::Int=0, digits::Int=0, base::Int=iszero(bits) ? 10 : 2) where {T<:Number} - bits = get_bits(bits, digits, base) +function ArbReal(x::Real; bits::Int=0, digits::Int=0, base::Int=iszero(bits) ? 10 : 2) + bits = get_bits(bits, digits, base, x) + ArbReal{bits}(x) +end +function ArbReal(x::Complex; bits::Int=0, digits::Int=0, base::Int=iszero(bits) ? 10 : 2) + bits = get_bits(bits, digits, base, x) ArbReal{bits}(x) end diff --git a/src/values/constructors.jl b/src/values/constructors.jl index 25a6635d..d835e2e7 100644 --- a/src/values/constructors.jl +++ b/src/values/constructors.jl @@ -21,7 +21,11 @@ end function (::Type{T})(x::ArbComplex{P}, rm::RoundingMode=RoundNearest) where {P,T<:IEEEFloat} isreal(x) || throw(InexactError(nameof(T), T, x)) - T(ArbFloat(midpoint(real(x))), rm) + T(convert(ArbFloat, midpoint(real(x))), rm) +end +function BigFloat(x::ArbComplex{P}, rm::RoundingMode=RoundNearest) where {P} + isreal(x) || throw(InexactError(nameof(T), T, x)) + BigFloat(convert(ArbFloat, midpoint(real(x))), rm) end function Complex{T}(x::ArbComplex{P}, rm::RoundingMode=RoundNearest) where {P,T<:Real} @@ -61,25 +65,27 @@ ArbReal(x::Irrational) = ArbReal{workingprecision(ArbReal)}(x) # retype -ArbFloat(x::ArbReal{P}) where {P} = ArbFloat{P}(x) -ArbReal(x::ArbFloat{P}) where {P} = ArbReal{P}(x) -ArbReal(x::ArbComplex{P}) where {P} = ArbReal{P}(x) -ArbFloat(x::ArbComplex{P}) where {P} = ArbFloat{P}(x) -ArbReal(x::Complex) = ArbReal{workingprecision(ArbReal)}(x) -ArbFloat(x::Complex) = ArbFloat{workingprecision(ArbFloat)}(x) +#ArbFloat(x::ArbNumber{P}) where {P} = ArbFloat{P}(x) +#ArbReal(x::ArbNumber{P}) where {P} = ArbReal{P}(x) +#ArbFloat(x::ArbReal{P}) where {P} = ArbFloat{P}(x) +#ArbReal(x::ArbFloat{P}) where {P} = ArbReal{P}(x) +#ArbReal(x::ArbComplex{P}) where {P} = ArbReal{P}(x) +#ArbFloat(x::ArbComplex{P}) where {P} = ArbFloat{P}(x) +#ArbReal(x::Complex) = ArbReal{workingprecision(ArbReal)}(x) +#ArbFloat(x::Complex) = ArbFloat{workingprecision(ArbFloat)}(x) -function _ArbFloatReal(::Type{T}, x::ArbComplex{P}) where {P,T<:ArbFloatReal} +function convert_to_real(::Type{T}, x::ArbComplex{P}) where {P,T<:Real} if isreal(x) convert(T, real(x)) else - throw(InexactError(:ArbReal, ArbReal, x)) + throw(InexactError(nameof(T), T, x)) end end ArbFloat{Q}(x::ArbReal{P}) where {P,Q} = ArbFloat{Q}(ArbReal{Q}(x)) ArbReal{Q}(x::ArbFloat{P}) where {P,Q} = ArbReal{Q}(ArbFloat{Q}(x)) -ArbReal{Q}(x::ArbComplex) where {Q} = _ArbFloatReal(ArbReal{Q}, x) -ArbFloat{Q}(x::ArbComplex) where {Q} = _ArbFloatReal(ArbFloat{Q}, x) +ArbReal{Q}(x::ArbComplex) where {Q} = convert_to_real(ArbReal{Q}, x) +ArbFloat{Q}(x::ArbComplex) where {Q} = convert_to_real(ArbFloat{Q}, x) # change precision diff --git a/test/complex.jl b/test/complex.jl index 0f081636..0c267365 100644 --- a/test/complex.jl +++ b/test/complex.jl @@ -80,6 +80,8 @@ @test ArbComplex{256}(1.0 + im) == ArbComplex(1 + im, bits=256) @test ArbComplex{256}(1.0 + im) == ArbComplex(1, 1, bits=256) + @test typeof(ArbComplex(1, ArbReal(2), bits=25)) == ArbComplex{49} + @test 1e-38 <= abs(radius(ArbComplex(π))) <= 2e-38 api = ArbComplex(pi, -pi) # note `-pi is Float64` @@ -175,4 +177,16 @@ @test_throws MethodError ArbComplex(ArbComplex(1), 0) end + @testset "(...)(::ArbComplex)" begin + ac = ArbComplex{64}(1 + pi*im) + @test Complex{BigFloat}(ac) ≈ ac + @test_throws InexactError Float64(ac) + @test_throws InexactError Int(ac) + + ac = ArbComplex{64}(120000) + @test BigFloat(ac) == 120000 + @test Int64(ac) == 120000 + @test_throws InexactError Int16(ac) + end + end diff --git a/test/misc.jl b/test/misc.jl index 303a1c35..f6ed7da4 100644 --- a/test/misc.jl +++ b/test/misc.jl @@ -48,9 +48,19 @@ @test copy(a, 24) == copy(a, 24, RoundUp) @test_throws DomainError copy(a, 10) + @test Float64(a) ≈ Float64(pi) + @test Float32(a) ≈ Float32(pi) + @test Float16(a) ≈ Float16(pi) + @test BigFloat(a) ≈ BigFloat(pi, precision=128) + @test ArbFloat(typemax(UInt)) == typemax(UInt) @test ArbFloat(typemax(Int128)) == typemax(Int128) @test ArbFloat(1 // 3) ≈ ArbFloat(1) / ArbFloat(3) + @test ArbFloat(a, bits=25) isa ArbFloat{49} + + @test ArbFloat{25}(ArbReal{26}(42)) == 42 + @test ArbFloat{25}(ArbComplex{26}(42)) == 42 + @test_throws InexactError ArbFloat{25}(ArbComplex{26}(1+im)) end @testset "convert ArbFloat to Integer" begin @@ -123,14 +133,22 @@ @test copy(a) == a @test copy(a) !== a - ma = ArbNumerics.Mag(a) - @test Float64(ma) ≈ Float64(pi) + @test Float64(a) ≈ Float64(pi) @test Float32(a) ≈ Float32(pi) @test Float16(a) ≈ Float16(pi) + @test BigFloat(a) ≈ BigFloat(pi, precision=128) + + ma = ArbNumerics.Mag(a) + @test Float64(ma) ≈ Float64(pi) + @test Float32(ma) ≈ Float32(pi) + @test Float16(ma) ≈ Float16(pi) + @test_throws MethodError BigFloat(ma) ≈ BigFloat(pi, precision=128) @test ArbReal(typemax(UInt)) == typemax(UInt) @test ArbReal(typemax(Int128)) == typemax(Int128) @test ArbReal(1 // 3) ≈ ArbReal(1) / ArbReal(3) + @test ArbReal(1 + 0im) == 1 + @test ArbReal{25}(ArbComplex{26}(42)) == 42 end @testset "convert ArbReal to Integer" begin