{-# LINE 1 "src/Data/Number/Flint/Arb/Calc/FFI.hsc" #-} {-| module : Data.Number.Flint.Arb.Calc.FFI copyright : (c) 2022 Hartmut Monien license : GNU GPL, version 2 or above (see LICENSE) maintainer : hmonien@uni-bonn.de -} module Data.Number.Flint.Arb.Calc.FFI ( -- * Calculus with real-valued functions ArfInterval (..) , CArfInterval (..) , newArfInterval , withArfInterval , withNewArfInterval , CArbCalcFunc , arf_interval_init , arf_interval_clear , _arf_interval_vec_init , _arf_interval_vec_clear , arf_interval_set , arf_interval_swap , arf_interval_get_arb , arf_interval_printd , arf_interval_fprintd , arb_calc_isolate_roots , arb_calc_refine_root_bisect -- * Newton-based root finding , arb_calc_newton_conv_factor , arb_calc_newton_step , arb_calc_refine_root_newton -- * Return types , arb_calc_success , arb_calc_imprecise_input , arb_calc_no_convergence ) where -- Calculus with real-valued functions ----------------------------------------- import Control.Monad import Foreign.C.String import Foreign.C.Types import Foreign.ForeignPtr import Foreign.Ptr ( Ptr, FunPtr, plusPtr, nullPtr, castPtr ) import Foreign.Storable import Foreign.Marshal ( free ) import Data.Int ( Int64 ) import Data.Functor ((<&>)) import Data.Number.Flint.Flint import Data.Number.Flint.Arb import Data.Number.Flint.Arb.Types -- arf_interval_t -------------------------------------------------------------- data ArfInterval = ArfInterval {-# UNPACK #-} !(ForeignPtr CArfInterval) data CArfInterval = CArfInterval CArf CArf instance Storable CArfInterval where {-# INLINE sizeOf #-} sizeOf :: CArfInterval -> Int sizeOf CArfInterval _ = (Int 64) {-# LINE 66 "src/Data/Number/Flint/Arb/Calc/FFI.hsc" #-} {-# INLINE alignment #-} alignment :: CArfInterval -> Int alignment CArfInterval _ = Int 8 {-# LINE 68 "src/Data/Number/Flint/Arb/Calc/FFI.hsc" #-} peek = error "CArfInterval.peek: Not defined" poke :: Ptr CArfInterval -> CArfInterval -> IO () poke = [Char] -> Ptr CArfInterval -> CArfInterval -> IO () forall a. HasCallStack => [Char] -> a error [Char] "CArfInterval.poke: Not defined" -- | Create a new `ArfInterval` structure. newArfInterval :: IO ArfInterval newArfInterval = do ForeignPtr CArfInterval x <- IO (ForeignPtr CArfInterval) forall a. Storable a => IO (ForeignPtr a) mallocForeignPtr ForeignPtr CArfInterval -> (Ptr CArfInterval -> IO ()) -> IO () forall a b. ForeignPtr a -> (Ptr a -> IO b) -> IO b withForeignPtr ForeignPtr CArfInterval x Ptr CArfInterval -> IO () arf_interval_init FinalizerPtr CArfInterval -> ForeignPtr CArfInterval -> IO () forall a. FinalizerPtr a -> ForeignPtr a -> IO () addForeignPtrFinalizer FinalizerPtr CArfInterval p_arf_interval_clear ForeignPtr CArfInterval x ArfInterval -> IO ArfInterval forall a. a -> IO a forall (m :: * -> *) a. Monad m => a -> m a return (ArfInterval -> IO ArfInterval) -> ArfInterval -> IO ArfInterval forall a b. (a -> b) -> a -> b $ ForeignPtr CArfInterval -> ArfInterval ArfInterval ForeignPtr CArfInterval x -- | Use `ArfInterval` structure. {-# INLINE withArfInterval #-} withArfInterval :: ArfInterval -> (Ptr CArfInterval -> IO a) -> IO (ArfInterval, a) withArfInterval (ArfInterval ForeignPtr CArfInterval x) Ptr CArfInterval -> IO a f = ForeignPtr CArfInterval -> (Ptr CArfInterval -> IO (ArfInterval, a)) -> IO (ArfInterval, a) forall a b. ForeignPtr a -> (Ptr a -> IO b) -> IO b withForeignPtr ForeignPtr CArfInterval x ((Ptr CArfInterval -> IO (ArfInterval, a)) -> IO (ArfInterval, a)) -> (Ptr CArfInterval -> IO (ArfInterval, a)) -> IO (ArfInterval, a) forall a b. (a -> b) -> a -> b $ \Ptr CArfInterval xp -> Ptr CArfInterval -> IO a f Ptr CArfInterval xp IO a -> (a -> (ArfInterval, a)) -> IO (ArfInterval, a) forall (f :: * -> *) a b. Functor f => f a -> (a -> b) -> f b <&> (ForeignPtr CArfInterval -> ArfInterval ArfInterval ForeignPtr CArfInterval x,) -- | Use new `ArfInterval` structure. {-# INLINE withNewArfInterval #-} withNewArfInterval :: (Ptr CArfInterval -> IO a) -> IO (ArfInterval, a) withNewArfInterval Ptr CArfInterval -> IO a f = IO ArfInterval newArfInterval IO ArfInterval -> (ArfInterval -> IO (ArfInterval, a)) -> IO (ArfInterval, a) forall a b. IO a -> (a -> IO b) -> IO b forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b >>= (ArfInterval -> (Ptr CArfInterval -> IO a) -> IO (ArfInterval, a)) -> (Ptr CArfInterval -> IO a) -> ArfInterval -> IO (ArfInterval, a) forall a b c. (a -> b -> c) -> b -> a -> c flip ArfInterval -> (Ptr CArfInterval -> IO a) -> IO (ArfInterval, a) forall {a}. ArfInterval -> (Ptr CArfInterval -> IO a) -> IO (ArfInterval, a) withArfInterval Ptr CArfInterval -> IO a f -- arb_calc_func_t ------------------------------------------------------------- type CArbCalcFunc = Ptr CArb -> Ptr CArb -> Ptr () -> CLong -> CLong -> IO CInt -- arb_calc_returns ------------------------------------------------------------ type ArbCalcReturn = CInt arb_calc_success, arb_calc_imprecise_input, arb_calc_no_convergence :: ArbCalcReturn arb_calc_success :: ArbCalcReturn arb_calc_success = ArbCalcReturn 0 {-# LINE 98 "src/Data/Number/Flint/Arb/Calc/FFI.hsc" #-} arb_calc_imprecise_input = 1 {-# LINE 99 "src/Data/Number/Flint/Arb/Calc/FFI.hsc" #-} arb_calc_no_convergence = 2 {-# LINE 100 "src/Data/Number/Flint/Arb/Calc/FFI.hsc" #-} -- Subdivision-based root finding ---------------------------------------------- -- | /arf_interval_init/ /v/ foreign import ccall "arb_calc.h arf_interval_init_" arf_interval_init :: Ptr CArfInterval -> IO () -- | /arf_interval_clear/ /v/ foreign import ccall "arb_calc.h arf_interval_clear_" arf_interval_clear :: Ptr CArfInterval -> IO () foreign import ccall "arb_calc.h &arf_interval_clear_" p_arf_interval_clear :: FunPtr (Ptr CArfInterval -> IO ()) -- | /_arf_interval_vec_init/ /n/ foreign import ccall "arb_calc.h _arf_interval_vec_init_" _arf_interval_vec_init :: CLong -> IO (Ptr CArfInterval) -- | /_arf_interval_vec_clear/ /v/ /n/ foreign import ccall "arb_calc.h _arf_interval_vec_clear_" _arf_interval_vec_clear :: Ptr CArfInterval -> CLong -> IO () -- | /arf_interval_set/ /v/ /u/ foreign import ccall "arb_calc.h arf_interval_set_" arf_interval_set :: Ptr CArfInterval -> Ptr CArfInterval -> IO () -- | /arf_interval_swap/ /v/ /u/ foreign import ccall "arb_calc.h arf_interval_swap_" arf_interval_swap :: Ptr CArfInterval -> Ptr CArfInterval -> IO () -- | /arf_interval_get_arb/ /x/ /v/ /prec/ foreign import ccall "arb_calc.h arf_interval_get_arb_" arf_interval_get_arb :: Ptr CArb -> Ptr CArfInterval -> CLong -> IO () foreign import ccall "arb_calc.h arf_interval_get_strd" arf_interval_get_strd :: Ptr CArfInterval -> CLong -> IO CString -- | /arf_interval_printd/ /v/ /n/ -- -- Helper functions for endpoint-based intervals. arf_interval_printd :: Ptr CArfInterval -> CLong -> IO () arf_interval_printd :: Ptr CArfInterval -> CLong -> IO () arf_interval_printd Ptr CArfInterval x CLong digits = do (Ptr CArfInterval -> IO CString) -> Ptr CArfInterval -> IO ArbCalcReturn forall a. (Ptr a -> IO CString) -> Ptr a -> IO ArbCalcReturn printCStr (Ptr CArfInterval -> CLong -> IO CString `arf_interval_get_strd` CLong digits) Ptr CArfInterval x () -> IO () forall a. a -> IO a forall (m :: * -> *) a. Monad m => a -> m a return () -- | /arf_interval_fprintd/ /file/ /v/ /n/ -- -- Helper functions for endpoint-based intervals. foreign import ccall "arb_calc.h arf_interval_fprintd" arf_interval_fprintd :: Ptr CFile -> Ptr CArfInterval -> CLong -> IO () -- | /arb_calc_isolate_roots/ /found/ /flags/ /func/ /param/ /interval/ /maxdepth/ /maxeval/ /maxfound/ /prec/ -- -- Rigorously isolates single roots of a real analytic function on the -- interior of an interval. -- -- This routine writes an array of /n/ interesting subintervals of -- /interval/ to /found/ and corresponding flags to /flags/, returning the -- integer /n/. The output has the following properties: -- -- - The function has no roots on /interval/ outside of the output -- subintervals. -- - Subintervals are sorted in increasing order (with no overlap except -- possibly starting and ending with the same point). -- - Subintervals with a flag of 1 contain exactly one (single) root. -- - Subintervals with any other flag may or may not contain roots. -- -- If no flags other than 1 occur, all roots of the function on /interval/ -- have been isolated. If there are output subintervals on which the -- existence or nonexistence of roots could not be determined, the user may -- attempt further searches on those subintervals (possibly with increased -- precision and\/or increased bounds for the breaking criteria). Note that -- roots of multiplicity higher than one and roots located exactly at -- endpoints cannot be isolated by the algorithm. -- -- The following breaking criteria are implemented: -- -- - At most /maxdepth/ recursive subdivisions are attempted. The -- smallest details that can be distinguished are therefore about -- \(2^{-\text{maxdepth}}\) times the width of /interval/. A typical, -- reasonable value might be between 20 and 50. -- - If the total number of tested subintervals exceeds /maxeval/, the -- algorithm is terminated and any untested subintervals are added to -- the output. The total number of calls to /func/ is thereby -- restricted to a small multiple of /maxeval/ (the actual count can be -- slightly higher depending on implementation details). A typical, -- reasonable value might be between 100 and 100000. -- - The algorithm terminates if /maxfound/ roots have been isolated. In -- particular, setting /maxfound/ to 1 can be used to locate just one -- root of the function even if there are numerous roots. To try to -- find all roots, /LONG_MAX/ may be passed. -- -- The argument /prec/ denotes the precision used to evaluate the function. -- It is possibly also used for some other arithmetic operations performed -- internally by the algorithm. Note that it probably does not make sense -- for /maxdepth/ to exceed /prec/. -- -- Warning: it is assumed that subdivision points of /interval/ can be -- represented exactly as floating-point numbers in memory. Do not pass -- \(1 \pm 2^{-10^{100}}\) as input. foreign import ccall "arb_calc.h arb_calc_isolate_roots" arb_calc_isolate_roots :: Ptr (Ptr CArfInterval) -> Ptr (Ptr CInt)-> FunPtr CArbCalcFunc -> Ptr () -> Ptr CArfInterval -> CLong -> CLong -> CLong -> CLong -> IO CLong -- | /arb_calc_refine_root_bisect/ /r/ /func/ /param/ /start/ /iter/ /prec/ -- -- Given an interval /start/ known to contain a single root of /func/, -- refines it using /iter/ bisection steps. The algorithm can return a -- failure code if the sign of the function at an evaluation point is -- ambiguous. The output /r/ is set to a valid isolating interval (possibly -- just /start/) even if the algorithm fails. foreign import ccall "arb_calc.h arb_calc_refine_root_bisect" arb_calc_refine_root_bisect :: Ptr CArfInterval -> FunPtr CArbCalcFunc -> Ptr () -> Ptr CArfInterval -> CLong -> CLong -> IO CInt -- Newton-based root finding --------------------------------------------------- -- | /arb_calc_newton_conv_factor/ /conv_factor/ /func/ /param/ /conv_region/ /prec/ -- -- Given an interval \(I\) specified by /conv_region/, evaluates a bound -- for \(C = \sup_{t,u \in I} \frac{1}{2} |f''(t)| / |f'(u)|\), where \(f\) -- is the function specified by /func/ and /param/. The bound is obtained -- by evaluating \(f'(I)\) and \(f''(I)\) directly. If \(f\) is -- ill-conditioned, \(I\) may need to be extremely precise in order to get -- an effective, finite bound for /C/. foreign import ccall "arb_calc.h arb_calc_newton_conv_factor" arb_calc_newton_conv_factor :: Ptr CArf -> FunPtr CArbCalcFunc -> Ptr () -> Ptr CArb -> CLong -> IO () -- | /arb_calc_newton_step/ /xnew/ /func/ /param/ /x/ /conv_region/ /conv_factor/ /prec/ -- -- Performs a single step with an interval version of Newton\'s method. The -- input consists of the function \(f\) specified by /func/ and /param/, a -- ball \(x = [m-r, m+r]\) known to contain a single root of \(f\), a ball -- \(I\) (/conv_region/) containing \(x\) with an associated bound -- (/conv_factor/) for -- \(C = \sup_{t,u \in I} \frac{1}{2} |f''(t)| / |f'(u)|\), and a working -- precision /prec/. -- -- The Newton update consists of setting \(x' = [m'-r', m'+r']\) where -- \(m' = m - f(m) / f'(m)\) and \(r' = C r^2\). The expression -- \(m - f(m) / f'(m)\) is evaluated using ball arithmetic at a working -- precision of /prec/ bits, and the rounding error during this evaluation -- is accounted for in the output. We now check that \(x' \in I\) and -- \(r' < r\). If both conditions are satisfied, we set /xnew/ to \(x'\) -- and return /ARB_CALC_SUCCESS/. If either condition fails, we set /xnew/ -- to \(x\) and return /ARB_CALC_NO_CONVERGENCE/, indicating that no -- progress is made. foreign import ccall "arb_calc.h arb_calc_newton_step" arb_calc_newton_step :: Ptr CArb -> FunPtr CArbCalcFunc -> Ptr () -> Ptr CArb -> Ptr CArb -> Ptr CArf -> CLong -> IO CInt -- | /arb_calc_refine_root_newton/ /r/ /func/ /param/ /start/ /conv_region/ /conv_factor/ /eval_extra_prec/ /prec/ -- -- Refines a precise estimate of a single root of a function to high -- precision by performing several Newton steps, using nearly optimally -- chosen doubling precision steps. -- -- The inputs are defined as for /arb_calc_newton_step/, except for the -- precision parameters: /prec/ is the target accuracy and -- /eval_extra_prec/ is the estimated number of guard bits that need to be -- added to evaluate the function accurately close to the root (for -- example, if the function is a polynomial with large coefficients of -- alternating signs and Horner\'s rule is used to evaluate it, the extra -- precision should typically be approximately the bit size of the -- coefficients). -- -- This function returns /ARB_CALC_SUCCESS/ if all attempted Newton steps -- are successful (note that this does not guarantee that the computed root -- is accurate to /prec/ bits, which has to be verified by the user), only -- that it is more accurate than the starting ball. -- -- On failure, /ARB_CALC_IMPRECISE_INPUT/ or /ARB_CALC_NO_CONVERGENCE/ may -- be returned. In this case, /r/ is set to a ball for the root which is -- valid but likely does have full accuracy (it can possibly just be equal -- to the starting ball). foreign import ccall "arb_calc.h arb_calc_refine_root_newton" arb_calc_refine_root_newton :: Ptr CArb -> FunPtr CArbCalcFunc -> Ptr () -> Ptr CArb -> Ptr CArb -> Ptr CArf -> CLong -> CLong -> IO CInt