import GHC.Conc import System.IO.Unsafe import System.CPUTime import Options.Applicative import Control.Monad import Control.Applicative ((<|>)) import Foreign.C.Types import Foreign.C.String import Foreign.Ptr import Foreign.Marshal.Alloc (free) import Foreign.Marshal.Utils (with) import Foreign.Storable import Text.ParserCombinators.ReadP hiding (option) import Text.Read (readMaybe) import Text.Printf import Data.Bits import Data.Char import Data.List (intercalate) import Data.Number.Flint import Integrands main = run =<< execParser opts where hDesc = "Calculate integrals using acb_calculate." desc = "Calculate integrals from list in range 0:" ++ show (length description) ++ ". " ++ "For a list of implemented integrals use --list." opts = info (parameters <**> helper) ( fullDesc <> progDesc desc <> header hDesc) run params = do if list params then do putStrLn "List of implemented integrals:\n" mapM_ (uncurry (printf " %2d %s\n")) $ zip [0 :: Int ..] description else do calc params calc params@(Parameters list (Range (start, end)) prec opt_goal tol twice heap verbose deg eval depth num_threads) = do numProcessors <- getNumProcessors if num_threads <= fromIntegral numProcessors then do flint_set_num_threads num_threads else do let msg = "number of threads (=" ++ show num_threads ++ ") greater than " ++ "number of processors (=" ++ show numProcessors ++ ").\n\n" ++ "choose [1.." ++ show numProcessors ++ "] threads." error msg let use_heap = if heap then 1 else 0 goal = if opt_goal == 0 then prec else opt_goal opts <- newAcbCalcIntegrateOpt_ deg eval depth use_heap verbose print opts withAcbCalcIntegrateOpt opts $ \opts -> do withMag tol $ \tol -> do flag <- mag_is_zero tol when (flag == 1) $ do mag_set_ui_2exp_si tol 1 (-prec) withNewAcb $ \a -> do withNewAcb $ \b -> do withNewAcb $ \s -> do withNewAcb $ \t -> do forM_ [start .. end] $ \j -> do let desc = description !! j startTime <- getCPUTime replicateM_ (if twice then 2 else 1) $ do case j of 0 -> do f <- makeFunPtr f_sin acb_set_si a 0 acb_set_si b 100 acb_calc_integrate s f nullPtr a b goal tol opts prec return () 1 -> do f <- makeFunPtr f_atanderiv acb_set_si a 0 acb_set_si b 1 acb_calc_integrate s f nullPtr a b goal tol opts prec acb_mul_2exp_si s s 2 return () 2 -> do f <- makeFunPtr f_atanderiv acb_set_si a 0 acb_set_si b 1 acb_mul_2exp_si b b goal acb_calc_integrate s f nullPtr a b goal tol opts prec arb_add_error_2exp_si (acb_realref s) (-goal) acb_mul_2exp_si s s 1 return () 3 -> do f <- makeFunPtr f_circle acb_set_si a 0 acb_set_si b 1 acb_calc_integrate s f nullPtr a b goal tol opts prec acb_mul_2exp_si s s 2 return () 4 -> do f <- makeFunPtr f_rump acb_set_si a 0 acb_set_si b 8 acb_calc_integrate s f nullPtr a b goal tol opts prec return () 5 -> do f <- makeFunPtr f_floor acb_set_si a 1 acb_set_si b 101 acb_calc_integrate s f nullPtr a b goal tol opts prec return () 6 -> do f <- makeFunPtr f_helfgott acb_set_si a 0 acb_set_si b 1 acb_calc_integrate s f nullPtr a b goal tol opts prec return () 7 -> do f <- makeFunPtr f_zeta let points = [(-1, -1), (2, -1), (2, 1), (-1, 1)] path = zip points $ tail $ cycle points acb_zero s forM_ path $ \((x1, y1), (x2, y2)) -> do acb_set_si_si a x1 y1 acb_set_si_si b x2 y2 acb_calc_integrate t f nullPtr a b goal tol opts prec acb_add s s t prec acb_const_pi t prec acb_div s s t prec acb_mul_2exp_si s s (-1) acb_div_onei s s return() 8 -> do f <- makeFunPtr f_essing acb_set_si a 0 acb_set_si b 1 acb_calc_integrate s f nullPtr a b goal tol opts prec return () 9 -> do f <- makeFunPtr f_essing2 acb_set_si a 0 acb_set_si b 1 acb_calc_integrate s f nullPtr a b goal tol opts prec return () 10 -> do f <- makeFunPtr f_factorial1000 acb_set_si a 0 acb_set_si b 10000 acb_calc_integrate s f nullPtr a b goal tol opts prec return () 11 -> do f <- makeFunPtr f_gamma acb_set_si_si a 1 0 acb_set_si_si b 1 1000 acb_calc_integrate s f nullPtr a b goal tol opts prec return () 12 -> do f <- makeFunPtr f_sin_plus_small acb_set_si a (-10) acb_set_si b 10 acb_calc_integrate s f nullPtr a b goal tol opts prec return () 13 -> do f <- makeFunPtr f_exp acb_set_si a (-1020) acb_set_si b (-1010) acb_calc_integrate s f nullPtr a b goal tol opts prec return () 14 -> do f <- makeFunPtr f_gaussian acb_zero a acb_set_si b $ fromIntegral $ ceiling $ sqrt (fromIntegral goal * log 2) + 1 acb_calc_integrate s f nullPtr a b goal tol opts prec acb_mul b b b prec acb_neg b b acb_exp b b prec arb_add_error (acb_realref s) (acb_realref b) return () 15 -> do f <- makeFunPtr f_spike acb_zero a acb_one b acb_calc_integrate s f nullPtr a b goal tol opts prec return () 16 -> do f <- makeFunPtr f_monster acb_zero a acb_set_si b 8 acb_calc_integrate s f nullPtr a b goal tol opts prec return () 17 -> do f <- makeFunPtr f_sech acb_set_d a 0 acb_set_d b $ fromIntegral $ ceiling $ fromIntegral goal * log 2 + 1 acb_calc_integrate s f nullPtr a b goal tol opts prec acb_neg b b acb_exp b b prec acb_mul_2exp_si b b 1 arb_add_error (acb_realref s) (acb_realref b) return () 18 -> do f <- makeFunPtr f_sech3 acb_zero a acb_set_d b $ fromIntegral $ ceiling $ fromIntegral goal * log 2 / 3 + 2 acb_calc_integrate s f nullPtr a b goal tol opts prec acb_neg b b acb_mul_ui b b 3 prec acb_exp b b prec acb_mul_2exp_si b b 3 acb_div_ui b b 3 prec arb_add_error (acb_realref s) (acb_realref b) return () 19 -> do when (goal < 0) $ do error "goal < 0" -- error bound 2^-N (1+N) when truncated at 2^-N let bitCount x = finiteBitSize x - countLeadingZeros x n = goal + fromIntegral (bitCount goal) acb_one a acb_mul_2exp_si a a (-n) acb_one b f <- makeFunPtr f_log_div1p acb_calc_integrate s f nullPtr a b goal tol opts prec acb_set_si b (n+1) acb_mul_2exp_si b b (-n) arb_add_error (acb_realref s) (acb_realref b) return () 20 -> do f <- makeFunPtr f_log_div1p_transformed when (goal < 0) $ do error "goal < 0" -- error bound (N+1) exp(-N) when truncated at N let bitCount x = finiteBitSize x - countLeadingZeros x n = goal + fromIntegral (bitCount goal) acb_zero a acb_set_si b n acb_calc_integrate s f nullPtr a b goal tol opts prec acb_neg b b acb_exp b b prec acb_mul_si b b (n+1) prec arb_add_error (acb_realref s) (acb_realref b) return () 21 -> do f <- makeFunPtr f_elliptic_p_laurent_n let n = 10 ::CLong points = [(1, -1), (1, 1), (-1, 1), (-1, -1)] path = zip points $ tail $ cycle points acb_zero s with n $ \np -> do let p = castPtr np :: Ptr () forM_ path $ \((x1, y1), (x2, y2)) -> do acb_set_si_si a x1 y1 acb_set_si_si b x2 y2 acb_div_ui a a 2 prec acb_div_ui b b 2 prec acb_calc_integrate t f p a b goal tol opts prec acb_add s s t prec acb_const_pi t prec acb_div s s t prec acb_mul_2exp_si s s (-1) acb_div_onei s s return () 22 -> do f <- makeFunPtr f_zeta_frac let n = 1000 :: CLong acb_zero s acb_set_si_si a 100 0 acb_set_si_si b 100 n acb_calc_integrate t f nullPtr a b goal tol opts prec acb_add s s t prec acb_set_si_si a 100 n acb_set_si_si b 1 (2*n) acb_div_ui b b 2 prec acb_calc_integrate t f nullPtr a b goal tol opts prec acb_add s s t prec acb_div_onei s s arb_zero (acb_imagref s) acb_set_si t n acb_dirichlet_hardy_theta t t nullPtr nullPtr 1 prec acb_add s s t prec acb_const_pi t prec acb_div s s t prec acb_add_ui s s 1 prec return () 23 -> do acb_set_si a 0 acb_set_si b 1000 f <- makeFunPtr f_lambertw acb_calc_integrate s f nullPtr a b goal tol opts prec return () 24 -> do acb_zero a acb_const_pi b prec f <- makeFunPtr f_max_sin_cos acb_calc_integrate s f nullPtr a b goal tol opts prec return () 25 -> do acb_set_si a (-1) acb_set_si b 1 f <- makeFunPtr f_erf_bent acb_calc_integrate s f nullPtr a b goal tol opts prec return () 26 -> do acb_set_si a (-10) acb_set_si b 10 f <- makeFunPtr f_airy_ai acb_calc_integrate s f nullPtr a b goal tol opts prec return () 27 -> do acb_set_si a 0 acb_set_si b 10 f <- makeFunPtr f_horror acb_calc_integrate s f nullPtr a b goal tol opts prec return () 28 -> do acb_set_si_si a (-1) (-1) acb_set_si_si b (-1) 1 f <- makeFunPtr f_sqrt acb_calc_integrate s f nullPtr a b goal tol opts prec return () 29 -> do f <- makeFunPtr f_gaussian_twist acb_zero a acb_set_si b $ ceiling $ sqrt (fromIntegral goal * log 2) + 1 acb_calc_integrate s f nullPtr a b goal tol opts prec acb_mul b b b prec acb_neg b b acb_exp b b prec arb_add_error (acb_realref s) (acb_realref b) arb_add_error (acb_imagref s) (acb_realref b) return () 30 -> do f <- makeFunPtr f_exp_airy acb_zero a acb_set_si b $ ceiling $ fromIntegral goal * log 2 + 1 acb_calc_integrate s f nullPtr a b goal tol opts prec acb_mul b b b prec acb_neg b b acb_exp b b prec acb_mul_2exp_si b b 1 arb_add_error (acb_realref s) (acb_realref b) return () 31 -> do f <- makeFunPtr f_sin_cos_frac acb_zero a acb_const_pi b prec acb_calc_integrate s f nullPtr a b goal tol opts prec return () 32 -> do f <- makeFunPtr f_sin_near_essing acb_zero a acb_set_ui b 3 acb_calc_integrate s f nullPtr a b goal tol opts prec return () 33 -> do f <- makeFunPtr f_scaled_bessel acb_zero a acb_zero b let k = 3 scaled_bessel_select_N (acb_realref b) k prec with k $ \kp -> do let p = castPtr kp :: Ptr () acb_calc_integrate s f p a b goal tol opts prec scaled_bessel_tail_bound (acb_realref a) k (acb_realref b) prec arb_add_error (acb_realref s) (acb_realref a) return () 34 -> do f <- makeFunPtr f_scaled_bessel acb_zero a acb_zero b let k = 15 scaled_bessel_select_N (acb_realref b) k prec with k $ \kp -> do let p = castPtr kp :: Ptr () acb_calc_integrate s f p a b goal tol opts prec scaled_bessel_tail_bound (acb_realref a) k (acb_realref b) prec arb_add_error (acb_realref s) (acb_realref a) return () 35 -> do f <- makeFunPtr f_rsqrt acb_set_si_si a (-1) (-1) acb_set_si_si b (-1) 1 acb_calc_integrate s f nullPtr a b goal tol opts prec return () 36 -> do f <- makeFunPtr f_rgamma when (goal < 0) $ do error "goal < 0." acb_zero a acb_set_si b (4 + (goal + 1) `div` 2) acb_calc_integrate s f nullPtr a b goal tol opts prec return () _ -> do error "case does not exists." return () endTime <- getCPUTime putStrLn $ "I" ++ show j ++ " = " ++ desc let diff = fromIntegral (endTime - startTime) / 10^12 :: Double printf "computation time: %0.3f sec\n" diff let digits = round (0.333 * fromIntegral prec) :: CLong acb_printn s digits arb_str_none putStr "\n\n" return () -- auxiliary routines ---------------------------------------------------------- scaled_bessel_tail_bound b k n prec = do arb_const_pi b prec arb_mul b b n prec arb_ui_div b k b prec arb_sqrt b b prec arb_pow_ui b b k prec arb_mul b b n prec arb_mul_ui b b 2 prec arb_div_ui b b (k-2) prec scaled_bessel_select_N n k prec = do let f = log (fromIntegral k / pi) / log 2 :: CDouble c e = fromIntegral e * fromIntegral k / 2 - f - fromIntegral e < fromIntegral prec + 5 e = last $ takeWhile c [1, 2 ..] arb_one n arb_mul_2exp_si n n e -- Parser ---------------------------------------------------------------------- data Parameters = Parameters { list :: Bool , range :: Range , prec :: CLong , goal :: CLong , tol :: Mag , twice :: Bool , heap :: Bool , verbosity :: CInt , deg :: CLong , eval :: CLong , depth :: CLong , num_threads :: CInt } deriving Show parameters :: Parser Parameters parameters = Parameters <$> switch ( help "show list of implemented integrals." <> long "list" <> short 'l') <*> option rng ( help "range (for example --range 2:6 or --range 7)" <> long "range" <> short 'r' <> value (Range (0, length description - 1)) <> metavar "range") <*> option pos ( help "precision in bits (default p = 64)" <> long "prec" <> short 'p' <> value 64 <> metavar "p") <*> option pos ( help "approximate relative accuracy goal (p)" <> long "goal" <> value 0 <> metavar "goal") <*> option mag ( help "approximate absolute accuracy goal (default 2^-p)" <> long "tol" <> value (read "0") <> metavar "abstol") <*> switch ( help "run twice (to see overhead of computing nodes)" <> long "twice") <*> switch ( help "use heap for subinterval queue" <> long "heap") <*> option pos ( help "verbosity level" <> short 'v' <> long "verbosity" <> value 0 <> metavar "verbosity") <*> option pos ( help "use quadrature degree up to n" <> long "deg" <> value 0 <> metavar "degree") <*> option pos ( help "limit number of function evaluations to n" <> long "eval" <> value 0 <> metavar "eval") <*> option pos ( help "limit subinterval queue size to n" <> long "depth" <> value 0 <> metavar "depth") <*> option auto ( help "number of threads" <> long "threads" <> value 1 <> metavar "threads") -- ReadM parsers --------------------------------------------------------------- rng :: ReadM Range rng = eitherReader $ \s -> do let result@(Range (a, b)) = read s :: Range if 0 <= a && a <= b && b < length description then Right result else Left "could not parse range" pos :: (Read a, Integral a) => ReadM a pos = eitherReader $ \s -> do let result = read s if result >= 0 then Right result else Left "expected positive number" between a b = eitherReader $ \s -> do let result = read s if a <= result && result <= b then Right result else Left $ "expected number in range [" ++ show a ++ " .. " ++ show b ++ "]." mag = eitherReader $ \s -> do case readMaybe s of Just result -> Right result _ -> Left $ "parsing " ++ show s ++ " failed." -- instances Mag --------------------------------------------------------------- instance Read Mag where readsPrec _ = readP_to_S (mkMag 10 <$> parseArb) -- instances Range ------------------------------------------------------------- newtype Range = Range (Int, Int) instance Read Range where readsPrec _ = readP_to_S (parseRange <|> parseIndex) instance Show Range where show (Range (a, b)) = "[" ++ show a ++ ":" ++ show b ++ "]" -- parsers -------------------------------------------------------------------- parseArb = do (res, _) <- gather $ choice [ char '[' *> parseFloat *> pm *> parseFloat <* char ']' , char '[' *> pm *> parseFloat <* char ']' , char '[' *> parseFloat <* char ']' , parseFloat ] return res where pm = skipSpaces *> string "+/-" <* many1 (char ' ') mkMag prec s = unsafePerformIO $ do (result, flag) <- withNewMag $ \m -> do withNewArb $ \x -> do withCString s $ \s -> do arb_set_str x s prec arb_get_mag m x return result parseFloat = do munch (\x -> x == '+' || x == '-') choice [nan, inf, num *> e, num] where nan = string "nan" inf = string "inf" num = munch1 isNumber *> munch (== '.') *> munch isNumber e = do char 'e' munch (\x -> x == '+' || x == '-') munch1 isNumber parseRange :: ReadP Range parseRange = do a <- read <$> munch1 isNumber char ':' b <- read <$> munch1 isNumber return $ Range (a, b) parseIndex :: ReadP Range parseIndex = do a <- read <$> munch1 isNumber return $ Range (a, a) -------------------------------------------------------------------------------- instance Show AcbCalcIntegrateOpt where show x = unsafePerformIO $ do (_, result) <- withAcbCalcIntegrateOpt x $ \x -> do CAcbCalcIntegrateOpt deg eval depth use_heap verbosity <- peek x return $ "options:" ++ " deg=" ++ show deg ++ " eval=" ++ show eval ++ " depth=" ++ show depth ++ " heap=" ++ show use_heap ++ " verbosity=" ++ show verbosity return result