Enh/acb integrate

deps/build.jl

@@ -6,7 +6,7 @@ oldwdir = pwd()
 @show MPFR_VERSION = "4.0.0"
 @show ANTIC_VERSION = "9fb5b8d5ccfad13d1ec5b59d4fd13a9fde94c78e"
 @show FLINT_VERSION = "84ee095d0cdb66ca228cf3ff67f0f3a206184971"
-@show ARB_VERSION = "6035ee2420b7a3fa0259c92dcfa5de4bc76a4b95"
+@show ARB_VERSION = "232135f35eeebb74afcea5dd7be436142bee9227"
 
 pkgdir = dirname(dirname(@__FILE__))
 wdir = joinpath(pkgdir, "deps")

src/Fields.jl

@@ -18,10 +18,12 @@ include("arb/arb.jl")
 
 include("arb/acb.jl")
 
+include("arb/acb_calc.jl")
+
 //(a::T, b::T) where {T <: FieldElem} = divexact(a, b)
 
 //(x::T, y::Union{Integer, Rational}) where {T <: RingElem} = x//parent(x)(y)
- 
+
 //(x::Union{Integer, Rational}, y::T) where {T <: RingElem} = parent(y)(x)//y
 
 function gcd(x::T, y::T) where {T <: FieldElem}
@@ -30,5 +32,3 @@ function gcd(x::T, y::T) where {T <: FieldElem}
 end
 
 characteristic(R::Field) = 0
-
-

src/arb/ArbTypes.jl

@@ -26,32 +26,10 @@ const ARB_RND_NEAR = Cint(4) # to nearest
 
 ################################################################################
 #
-# Types and memory management for ArbField
+# Structs for shallow operations
 #
 ################################################################################
 
-mutable struct ArbField <: Field
- prec::Int
-
- function ArbField(p::Int = 256, cached::Bool = true)
- arb_check_prec(p)
- if haskey(ArbFieldID, p)
- return ArbFieldID[p]
- else
- z = new(p)
- if cached
- ArbFieldID[p] = z
- end
- return z
- end
- end
-end
-
-const ArbFieldID = Dict{Int, ArbField}()
-
-prec(x::ArbField) = x.prec
-
-# these may be used for shallow operations
 mutable struct arf_struct
  exp::Int # fmpz
  size::UInt # mp_size_t
@@ -62,6 +40,17 @@ end
 mutable struct mag_struct
  exp::Int # fmpz
  man::UInt # mp_limb_t
+
+ function mag_struct()
+ m = new()
+ ccall((:mag_init, :libarb), Void, (Ref{mag_struct},), m)
+ finalizer(m, _mag_struct_fn)
+ return m
+ end
+end
+
+function _mag_struct_fn(m::mag_struct)
+ ccall((:mag_clear, :libarb), Void, (Ref{mag_struct}, ), m)
 end
 
 mutable struct arb_struct
@@ -86,8 +75,46 @@ mutable struct acb_struct
  imag_mid_d2::Int
  imag_rad_exp::Int # fmpz
  imag_rad_man::UInt
+
+ function acb_struct()
+ z = new()
+ ccall((:acb_init, :libarb), Void, (Ref{Nemo.acb_struct},), z)
+ finalizer(z, _acb_clear_fn)
+ return z
+ end
+end
+
+function _acb_clear_fn(x::acb_struct)
+ ccall((:acb_clear, :libarb), Void, (Ref{acb_struct}, ), x)
 end
 
+################################################################################
+#
+# Types and memory management for ArbField
+#
+################################################################################
+
+mutable struct ArbField <: Field
+ prec::Int
+
+ function ArbField(p::Int = 256, cached::Bool = true)
+ arb_check_prec(p)
+ if haskey(ArbFieldID, p)
+ return ArbFieldID[p]
+ else
+ z = new(p)
+ if cached
+ ArbFieldID[p] = z
+ end
+ return z
+ end
+ end
+end
+
+const ArbFieldID = Dict{Int, ArbField}()
+
+prec(x::ArbField) = x.prec
+
 mutable struct arb <: FieldElem
  mid_exp::Int # fmpz
  mid_size::UInt # mp_size_t
@@ -230,6 +257,25 @@ function _acb_clear_fn(x::acb)
  ccall((:acb_clear, :libarb), Void, (Ref{acb}, ), x)
 end
 
+mutable struct acb_calc_integrate_opts
+ deg_limit::Int # <=0: default of 0.5*min(prec, rel_goal) + 10
+ eval_limit::Int # <=0: default of 1000*prec*prec^2
+ depth_limit::Int # <=0: default of 2*prec
+ use_heap::Int32 # 0 append to the top of a stack; 1 binary heap
+ verbose::Int32 # 1 less verbose; 2 more verbose
+
+ function acb_calc_integrate_opts(deg_limit::Int, eval_limit::Int,
+ depth_limit::Int, use_heap::Int32, verbose::Int32)
+ return new(deg_limit, eval_limit, depth_limit, use_heap, verbose)
+ end
+
+ function acb_calc_integrate_opts()
+ opts = new()
+ ccall((:acb_calc_integrate_opt_init, :libarb),
+ Void, (Ref{acb_calc_integrate_opts}, ), opts)
+ return opts
+ end
+end
 
 ################################################################################
 #

src/arb/acb.jl

@@ -1582,6 +1582,16 @@ for (typeofx, passtoc) in ((acb, Ref{acb}), (Ptr{acb}, Ptr{acb}))
  end
 end
 
+function _acb_set(x::acb_struct, y::acb_struct)
+ ccall((:acb_set, :libarb), Void, (Ref{acb_struct}, Ref{acb_struct}), x, y)
+end
+
+function mag_set_ui_2exp_si(m::mag_struct, x::Int, y::Int)
+ ccall((:mag_set_ui_2exp_si, :libarb), Void,
+ (Ref{mag_struct}, UInt, Int), m, abs(x), y)
+ return m
+end
+
 ################################################################################
 #
 # Parent object overload
@@ -1616,6 +1626,13 @@ end
 (r::AcbField)(x::Rational{S}, y::Rational{T}) where {S <: Integer, T <: Integer} =
  r(fmpq(x), fmpq(y))
 
+function (T::Type{acb_struct})(z::acb)
+ res = acb_struct()
+ ccall((:acb_set, :libarb), Void,
+ (Ref{acb_struct}, Ref{acb}), res, z)
+ return res
+end
+
 ################################################################################
 #
 # AcbField constructor

src/arb/acb_calc.jl

@@ -0,0 +1,44 @@
+function acb_calc_func_wrap(res::acb_struct, x::acb, param::Ptr{Void}, order::Int, prec::Int)
+
+ acbF = unsafe_pointer_to_objref(param)::Function
+ FF = AcbField(prec)
+ z = acbF(FF(x))
+
+ _acb_set(res, acb_struct(z))
+ return Cint(0)::Cint
+end
+
+acb_calc_func_wrap_c() = cfunction(acb_calc_func_wrap, Cint,
+ (Ref{acb_struct}, Ref{acb}, Ptr{Void}, Int, Int))
+
+const ARB_CALC_SUCCESS = UInt(0)
+const ARB_CALC_NO_CONVERGENCE = UInt(2)
+
+function integrate(F::Function, a::acb, b::acb;
+ rel_goal::Int=prec(parent(a)),
+ abs_tol::mag_struct=mag_set_ui_2exp_si(mag_struct(), 1, -prec(parent(a))),
+ opts::acb_calc_integrate_opts=acb_calc_integrate_opts(),
+ precision::Int=prec(parent(a)))
+
+ res = AcbField(precision)()
+ ptrF = pointer_from_objref(F)
+
+ status = ccall((:acb_calc_integrate, :libarb), UInt,
+ (Ref{acb}, #res
+ Ptr{Void}, #func
+ Ptr{Void}, #params
+ Ref{acb}, #a
+ Ref{acb}, #b
+ Int, #rel_goal
+ Ref{mag_struct}, #abs_tol
+ Ref{acb_calc_integrate_opts}, #opts
+ Int),
+ res, acb_calc_func_wrap_c(), ptrF, a, b, rel_goal, abs_tol, opts, precision)
+
+ if status == ARB_CALC_SUCCESS
+ nothing
+ elseif status == ARB_CALC_NO_CONVERGENCE
+ warn("Integration did not obtained convergence")
+ end
+ return res
+end

test/arb/acb-test.jl

@@ -417,6 +417,27 @@ function test_acb_functions()
  println("PASS")
 end
 
+function test_acb_integration()
+ print("acb.integration...")
+
+ res = Nemo.integrate(x->x, CC(-1), CC(1))
+ @test contains(res, CC(0))
+ @test imag(res) == CC(0)
+ @test radius(real(res)) < 3e-19
+
+ res = Nemo.integrate(x->x^2,CC(-1), CC(1))
+ @test contains(res, CC(2//3))
+ @test imag(res) == CC(0)
+ @test radius(real(res)) < 7e-19
+
+ res = Nemo.integrate(sin, CC(0), const_pi(CC))
+ @test overlaps(res, CC(2))
+ @test imag(res) == CC(0)
+ @test radius(real(res)) < 4e-18
+
+ println("PASS")
+end
+
 function test_acb()
  test_acb_constructors()
  test_acb_printing()
@@ -429,6 +450,7 @@ function test_acb()
  test_acb_unsafe_ops()
  test_acb_constants()
  test_acb_functions()
+ test_acb_integration()
 
  println("")
 end