FFI.hsc

{-|
module      :  Data.Number.Flint.Bernoulli.FFI
copyright   :  (c) 2022 Hartmut Monien
license     :  GNU GPL, version 2 or above (see LICENSE)
maintainer  :  hmonien@uni-bonn.de
-}
module Data.Number.Flint.Bernoulli.FFI (
  -- * Support for Bernoulli numbers
  -- * Generation of Bernoulli numbers
    BernoulliRev (..)
  , CBernoulliRev (..)
  , newBernoulliRev
  , withBernoulliRev
  , withNewBernoulliRev
  , bernoulli_rev_init
  , bernoulli_rev_next
  , bernoulli_rev_clear
  , bernoulli_fmpq_vec_no_cache
  -- * Caching
  , bernoulli_cache_compute
  -- * Bounding
  , bernoulli_bound_2exp_si
  -- * Isolated Bernoulli numbers
  , bernoulli_mod_p_harvey
  , _bernoulli_fmpq_ui_zeta
  , _bernoulli_fmpq_ui
) where

-- Support for Bernoulli numbers -----------------------------------------------

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

import Data.Number.Flint.Flint
import Data.Number.Flint.Fmpz
import Data.Number.Flint.Fmpq

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

#include <flint/bernoulli.h>

-- bernoulli_rev_t -------------------------------------------------------------

data BernoulliRev = BernoulliRev {-# UNPACK #-} !(ForeignPtr CBernoulliRev)
type CBernoulliRev = CFlint BernoulliRev

instance Storable CBernoulliRev where
  sizeOf    _ = #{size      bernoulli_rev_t}
  alignment _ = #{alignment bernoulli_rev_t}
  peek = error "CBernoulliRev.peek: not implemented."
  poke = error "CBernoulliRev.poke: not implemented."
  
newBernoulliRev n = do
  x <- mallocForeignPtr
  withForeignPtr x $ \x -> do
    bernoulli_rev_init x n
  addForeignPtrFinalizer p_bernoulli_rev_clear x
  return $ BernoulliRev x

withBernoulliRev (BernoulliRev x) f = do
  withForeignPtr x $ \xp -> (BernoulliRev x,) <$> f xp

withNewBernoulliRev n f = do
  x <- newBernoulliRev n
  withBernoulliRev x f
  
-- Generation of Bernoulli numbers ---------------------------------------------

-- | /bernoulli_rev_init/ /iter/ /n/ 
-- 
-- Initializes the iterator /iter/. The first Bernoulli number to be
-- generated by calling @bernoulli_rev_next@ is \(B_n\). It is assumed that
-- \(n\) is even.
foreign import ccall "bernoulli.h bernoulli_rev_init"
  bernoulli_rev_init :: Ptr CBernoulliRev -> CULong -> IO ()

-- | /bernoulli_rev_next/ /numer/ /denom/ /iter/ 
-- 
-- Sets /numer/ and /denom/ to the exact, reduced numerator and denominator
-- of the Bernoulli number \(B_k\) and advances the state of /iter/ so that
-- the next invocation generates \(B_{k-2}\).
foreign import ccall "bernoulli.h bernoulli_rev_next"
  bernoulli_rev_next :: Ptr CFmpz -> Ptr CFmpz -> Ptr CBernoulliRev -> IO ()

-- | /bernoulli_rev_clear/ /iter/ 
-- 
-- Frees all memory allocated internally by /iter/.
foreign import ccall "bernoulli.h bernoulli_rev_clear"
  bernoulli_rev_clear :: Ptr CBernoulliRev -> IO ()

foreign import ccall "bernoulli.h &bernoulli_rev_clear"
  p_bernoulli_rev_clear :: FunPtr (Ptr CBernoulliRev -> IO ())

-- | /bernoulli_fmpq_vec_no_cache/ /res/ /a/ /num/ 
-- 
-- Writes /num/ consecutive Bernoulli numbers to /res/ starting with
-- \(B_a\). This function is not currently optimized for a small count
-- /num/. The entries are not read from or written to the Bernoulli number
-- cache; if retrieving a vector of Bernoulli numbers is needed more than
-- once, use @bernoulli_cache_compute@ followed by @bernoulli_fmpq_ui@
-- instead.
-- 
-- This function is a wrapper for the /rev/ iterators. It can use multiple
-- threads internally.
foreign import ccall "bernoulli.h bernoulli_fmpq_vec_no_cache"
  bernoulli_fmpq_vec_no_cache :: Ptr CFmpq -> CULong -> CLong -> IO ()

-- Caching ---------------------------------------------------------------------

-- | /bernoulli_cache_compute/ /n/ 
-- 
-- Makes sure that the Bernoulli numbers up to at least \(B_{n-1}\) are
-- cached. Calling @flint_cleanup()@ frees the cache.
-- 
-- The cache is extended by calling @bernoulli_fmpq_vec_no_cache@
-- internally.
foreign import ccall "bernoulli.h bernoulli_cache_compute"
  bernoulli_cache_compute :: CLong -> IO ()

-- Bounding --------------------------------------------------------------------

-- | /bernoulli_bound_2exp_si/ /n/ 
-- 
-- Returns an integer \(b\) such that \(|B_n| \le 2^b\). Uses a lookup
-- table for small \(n\), and for larger \(n\) uses the inequality
-- \(|B_n| < 4 n! / (2 \pi)^n < 4 (n+1)^{n+1} e^{-n} / (2 \pi)^n\). Uses
-- integer arithmetic throughout, with the bound for the logarithm being
-- looked up from a table. If \(|B_n| = 0\), returns /LONG_MIN/. Otherwise,
-- the returned exponent \(b\) is never more than one percent larger than
-- the true magnitude.
-- 
-- This function is intended for use when \(n\) small enough that one might
-- comfortably compute \(B_n\) exactly. It aborts if \(n\) is so large that
-- internal overflow occurs.
foreign import ccall "bernoulli.h bernoulli_bound_2exp_si"
  bernoulli_bound_2exp_si :: CULong -> IO CLong

-- Isolated Bernoulli numbers --------------------------------------------------

-- | /bernoulli_mod_p_harvey/ /n/ /p/ 
-- 
-- Returns the \(B_n\) modulo the prime number /p/, computed using
-- Harvey\'s algorithm < [Har2010]>. The running time is linear in /p/. If
-- /p/ divides the numerator of \(B_n\), /UWORD_MAX/ is returned as an
-- error code.
foreign import ccall "bernoulli.h bernoulli_mod_p_harvey"
  bernoulli_mod_p_harvey :: CULong -> CULong -> IO CULong

-- | /_bernoulli_fmpq_ui_zeta/ /num/ /den/ /n/ 
-- 
-- Sets /num/ and /den/ to the reduced numerator and denominator of the
-- Bernoulli number \(B_n\).
-- 
-- The /zeta/ version computes the denominator \(d\) using the von
-- Staudt-Clausen theorem, numerically approximates \(B_n\) using
-- @arb_bernoulli_ui_zeta@, and then rounds \(d B_n\) to the correct
-- numerator.
-- 
-- The /multi_mod/ version reconstructs \(B_n\) by computing the high bits
-- via the Riemann zeta function and the low bits via Harvey\'s
-- multimodular algorithm. The tuning parameter /alpha/ should be a
-- fraction between 0 and 1 controlling the number of bits to compute by
-- the multimodular algorithm. If set to a negative number, a default value
-- will be used.
foreign import ccall "bernoulli.h _bernoulli_fmpq_ui_zeta"
  _bernoulli_fmpq_ui_zeta :: Ptr CFmpz -> Ptr CFmpz -> CULong -> IO ()

-- | /_bernoulli_fmpq_ui/ /num/ /den/ /n/ 
-- 
-- Computes the Bernoulli number \(B_n\) as an exact fraction, for an
-- isolated integer \(n\). This function reads \(B_n\) from the global
-- cache if the number is already cached, but does not automatically extend
-- the cache by itself.
foreign import ccall "bernoulli.h _bernoulli_fmpq_ui"
  _bernoulli_fmpq_ui :: Ptr CFmpz -> Ptr CFmpz -> CULong -> IO ()