{-# LINE 1 "src/Data/Number/Flint/Hypgeom/FFI.hsc" #-}
{-|
module      :  Data.Number.Flint.Hypgeom.FFI
copyright   :  (c) 2022 Hartmut Monien
license     :  GNU GPL, version 2 or above (see LICENSE)
maintainer  :  hmonien@uni-bonn.de
-}
module Data.Number.Flint.Hypgeom.FFI (
  -- * Support for hypergeometric series
  -- * Types
    Hypgeom (..)
  , CHypgeom (..)
  , newHypgeom
  , withHypgeom
  , withNewHypgeom
  -- * Memory management
  , hypgeom_init
  , hypgeom_clear
  -- * Error bounding
  , hypgeom_estimate_terms
  , hypgeom_bound
  , hypgeom_precompute
  -- * Summation
  , arb_hypgeom_sum
  , arb_hypgeom_infsum
) where

-- Support for hypergeometric series -------------------------------------------

import Foreign.Ptr
import Foreign.ForeignPtr
import Foreign.Storable
import Foreign.C.Types

import Data.Number.Flint.Flint
import Data.Number.Flint.Arb.Types



-- hypgeom_t -------------------------------------------------------------------

data Hypgeom = Hypgeom {-# UNPACK #-} !(ForeignPtr CHypgeom)
type CHypgeom = CFlint Hypgeom

instance Storable CHypgeom where
  sizeOf :: CHypgeom -> Int
sizeOf    CHypgeom
_ = (Int
152)
{-# LINE 46 "src/Data/Number/Flint/Hypgeom/FFI.hsc" #-}
  alignment _ = 8
{-# LINE 47 "src/Data/Number/Flint/Hypgeom/FFI.hsc" #-}
  peek = undefined
  poke :: Ptr CHypgeom -> CHypgeom -> IO ()
poke = Ptr CHypgeom -> CHypgeom -> IO ()
forall a. HasCallStack => a
undefined

newHypgeom :: IO Hypgeom
newHypgeom = do
  ForeignPtr CHypgeom
x <- IO (ForeignPtr CHypgeom)
forall a. Storable a => IO (ForeignPtr a)
mallocForeignPtr
  ForeignPtr CHypgeom -> (Ptr CHypgeom -> IO ()) -> IO ()
forall a b. ForeignPtr a -> (Ptr a -> IO b) -> IO b
withForeignPtr ForeignPtr CHypgeom
x Ptr CHypgeom -> IO ()
hypgeom_init
  FinalizerPtr CHypgeom -> ForeignPtr CHypgeom -> IO ()
forall a. FinalizerPtr a -> ForeignPtr a -> IO ()
addForeignPtrFinalizer FinalizerPtr CHypgeom
p_hypgeom_clear ForeignPtr CHypgeom
x
  Hypgeom -> IO Hypgeom
forall a. a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return (Hypgeom -> IO Hypgeom) -> Hypgeom -> IO Hypgeom
forall a b. (a -> b) -> a -> b
$ ForeignPtr CHypgeom -> Hypgeom
Hypgeom ForeignPtr CHypgeom
x

withHypgeom :: Hypgeom -> (Ptr CHypgeom -> IO a) -> IO (Hypgeom, a)
withHypgeom (Hypgeom ForeignPtr CHypgeom
p) Ptr CHypgeom -> IO a
f = do
  ForeignPtr CHypgeom
-> (Ptr CHypgeom -> IO (Hypgeom, a)) -> IO (Hypgeom, a)
forall a b. ForeignPtr a -> (Ptr a -> IO b) -> IO b
withForeignPtr ForeignPtr CHypgeom
p ((Ptr CHypgeom -> IO (Hypgeom, a)) -> IO (Hypgeom, a))
-> (Ptr CHypgeom -> IO (Hypgeom, a)) -> IO (Hypgeom, a)
forall a b. (a -> b) -> a -> b
$ \Ptr CHypgeom
fp -> (ForeignPtr CHypgeom -> Hypgeom
Hypgeom ForeignPtr CHypgeom
p,) (a -> (Hypgeom, a)) -> IO a -> IO (Hypgeom, a)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Ptr CHypgeom -> IO a
f Ptr CHypgeom
fp

withNewHypgeom :: (Ptr CHypgeom -> IO a) -> IO (Hypgeom, a)
withNewHypgeom Ptr CHypgeom -> IO a
f = do
  Hypgeom
x <- IO Hypgeom
newHypgeom
  Hypgeom -> (Ptr CHypgeom -> IO a) -> IO (Hypgeom, a)
forall {a}. Hypgeom -> (Ptr CHypgeom -> IO a) -> IO (Hypgeom, a)
withHypgeom Hypgeom
x Ptr CHypgeom -> IO a
f
  
-- Memory management -----------------------------------------------------------

-- | /hypgeom_init/ /hyp/ 
--
foreign import ccall "hypgeom.h hypgeom_init"
  hypgeom_init :: Ptr CHypgeom -> IO ()

-- | /hypgeom_clear/ /hyp/ 
--
foreign import ccall "hypgeom.h hypgeom_clear"
  hypgeom_clear :: Ptr CHypgeom -> IO ()

foreign import ccall "hypgeom.h &hypgeom_clear"
  p_hypgeom_clear :: FunPtr (Ptr CHypgeom -> IO ())

-- Error bounding --------------------------------------------------------------

-- | /hypgeom_estimate_terms/ /z/ /r/ /d/ 
--
-- Computes an approximation of the largest \(n\) such that
-- \(|z|^n/(n!)^r = 2^{-d}\), giving a first-order estimate of the number
-- of terms needed to approximate the sum of a hypergeometric series of
-- weight \(r \ge 0\) and argument \(z\) to an absolute precision of
-- \(d \ge 0\) bits. If \(r = 0\), the direct solution of the equation is
-- given by \(n = (\log(1-z) - d \log 2) / \log z\). If \(r > 0\), using
-- \(\log n! \approx n \log n - n\) gives an equation that can be solved in
-- terms of the Lambert /W/-function as \(n = (d \log 2) / (r\,W\!(t))\)
-- where \(t = (d \log 2) / (e r z^{1/r})\).
-- 
-- The evaluation is done using double precision arithmetic. The function
-- aborts if the computed value of \(n\) is greater than or equal to
-- LONG_MAX \/ 2.
foreign import ccall "hypgeom.h hypgeom_estimate_terms"
  hypgeom_estimate_terms :: Ptr CMag -> CInt -> CLong -> IO CLong

-- | /hypgeom_bound/ /error/ /r/ /C/ /D/ /K/ /TK/ /z/ /prec/ 
--
-- Computes a truncation parameter sufficient to achieve /prec/ bits of
-- absolute accuracy, according to the strategy described above. The input
-- consists of \(r\), \(C\), \(D\), \(K\), precomputed bound for \(T(K)\),
-- and \(\tilde z = z (a_p / b_q)\), such that for \(k > K\), the
-- hypergeometric term ratio is bounded by
-- 
-- \[`\]
-- \[\frac{\tilde z}{k^r} \frac{k(k-D)}{(k-C)(k-2D)}.\]
-- 
-- Given this information, we compute a \(\varepsilon\) and an integer
-- \(n\) such that
-- \(\left| \sum_{k=n}^{\infty} T(k) \right| \le \varepsilon \le 2^{-\mathrm{prec}}\).
-- The output variable /error/ is set to the value of \(\varepsilon\), and
-- \(n\) is returned.
foreign import ccall "hypgeom.h hypgeom_bound"
  hypgeom_bound :: Ptr CMag -> CInt -> CLong -> CLong -> CLong -> Ptr CMag -> Ptr CMag -> CLong -> IO CLong

-- | /hypgeom_precompute/ /hyp/ 
--
-- Precomputes the bounds data \(C\), \(D\), \(K\) and an upper bound for
-- \(T(K)\).
foreign import ccall "hypgeom.h hypgeom_precompute"
  hypgeom_precompute :: Ptr CHypgeom -> IO ()

-- Summation -------------------------------------------------------------------

-- | /arb_hypgeom_sum/ /P/ /Q/ /hyp/ /n/ /prec/ 
--
-- Computes \(P, Q\) such that \(P / Q = \sum_{k=0}^{n-1} T(k)\) where
-- \(T(k)\) is defined by /hyp/, using binary splitting and a working
-- precision of /prec/ bits.
foreign import ccall "hypgeom.h arb_hypgeom_sum"
  arb_hypgeom_sum :: Ptr CArb -> Ptr CArb -> Ptr CHypgeom -> CLong -> CLong -> IO ()

-- | /arb_hypgeom_infsum/ /P/ /Q/ /hyp/ /tol/ /prec/ 
--
-- Computes \(P, Q\) such that \(P / Q = \sum_{k=0}^{\infty} T(k)\) where
-- \(T(k)\) is defined by /hyp/, using binary splitting and working
-- precision of /prec/ bits. The number of terms is chosen automatically to
-- bound the truncation error by at most \(2^{-\mathrm{tol}}\). The bound
-- for the truncation error is included in the output as part of /P/.
foreign import ccall "hypgeom.h arb_hypgeom_infsum"
  arb_hypgeom_infsum :: Ptr CArb -> Ptr CArb -> Ptr CHypgeom -> CLong -> CLong -> IO ()