<commit_before>{-# LANGUAGE ExistentialQuantification #-}
-- |
-- Code for working with series
module Numeric.Series (
    -- * Data type for series
    Series(..)
  , next
    -- * Constructors
  , enumSeriesFrom
  , enumSeriesFromStep
  , scanSeries
    -- * Evaluation of series
  , sumSeries
  , sumPowerSeries
  , seriesToList
    -- * Evaluation of continued fractions
  , evalModLentz
  ) where

import Control.Applicative
import Data.List (unfoldr)

import Numeric.MathFunctions.Constants (m_epsilon)


----------------------------------------------------------------

-- | Representation of infinite series
data Series a = forall s. Series s (s -> (a,s))

instance Functor Series where
  fmap f (Series s0 step) = Series s0 (\s -> let (a,s') = step s in (f a, s'))
  {-# INLINE fmap #-}

instance Applicative Series where
  pure a = Series () (\_ -> (a,()))
  Series sA fA <*> Series sB fB = Series (sA,sB) $ \(sa,sb) ->
    let (a,sa') = fA sa
        (b,sb') = fB sb
    in (a b, (sa',sb'))
  {-# INLINE pure  #-}
  {-# INLINE (<*>) #-}

next :: Series a -> (a,Series a)
{-# INLINE next #-}
next (Series s f) = case f s of
  (a,s') -> (a,Series s' f)


----------------------------------------------------------------
-- Constructors
----------------------------------------------------------------0

enumSeriesFrom :: Num a => a -> Series a
enumSeriesFrom i = Series i (\n -> (n,n+1))
{-# INLINE enumSeriesFrom #-}

enumSeriesFromStep :: Num a => a -> a -> Series a
enumSeriesFromStep n d = Series n (\i -> (i,i+d))
{-# INLINE enumSeriesFromStep #-}

scanSeries :: (b -> a -> b) -> b -> Series a -> Series b
{-# INLINE scanSeries #-}
scanSeries f b0 (Series s0 step) = Series (b0,s0) $ \(b,s) ->
  let (a,s') = step s
      b'     = f b a
  in (b,(b',s'))


----------------------------------------------------------------
-- Evaluation of series
----------------------------------------------------------------

sumSeries :: Series Double -> Double
{-# INLINE sumSeries #-}
sumSeries ser0
  = go x0 ser
  where 
    (x0,ser) = next ser0
    go x s | abs (d/x) < m_epsilon = x'
           | otherwise             = go x' s'
      where (d,s') = next s
            x' = x + d

sumPowerSeries :: Double -> Series Double -> Double
sumPowerSeries x ser = sumSeries $ liftA2 (*) (scanSeries (*) 1 (pure x)) ser
{-# INLINE sumPowerSeries #-}
            
seriesToList :: Series a -> [a]
seriesToList (Series s f) = unfoldr (Just . f) s



----------------------------------------------------------------
-- Evaluation of continued fractions
----------------------------------------------------------------

evalModLentz :: Series (Double,Double) -> Double
{-# INLINE evalModLentz #-}
evalModLentz ser0
  = let ((_,b0), ser) = next ser0
        f0            = maskZero b0
    in  go f0 f0 0 ser
  where
    tiny = 1e-60
    maskZero 0 = tiny
    maskZero x = x
    
    go f c d ser
      | abs (delta - 1) < m_epsilon = f'
      | otherwise                   = go f' c' d' ser'
      where     
          ((a,b),ser') = next ser
          d'    = recip $ maskZero $ b + a*d
          c'    = maskZero $ b + a/c 
          delta = c'*d'
          f'    = f*delta
<commit_msg>Optimize functions for evaluation of series<commit_after>{-# LANGUAGE BangPatterns              #-}
{-# LANGUAGE ExistentialQuantification #-}
-- |
-- Code for working with series
module Numeric.Series (
    -- * Data type for series
    Series(..)
  , next
    -- * Constructors
  , enumSeriesFrom
  , enumSeriesFromStep
  , scanSeries
    -- * Evaluation of series
  , sumSeries
  , sumPowerSeries
  , seriesToList
    -- * Evaluation of continued fractions
  , evalModLentz
  ) where

import Control.Applicative
import Data.List (unfoldr)

import Numeric.MathFunctions.Constants (m_epsilon)


----------------------------------------------------------------

-- | Representation of infinite series
data Series a = forall s. Series s (s -> (a,s))

instance Functor Series where
  fmap f (Series s0 step) = Series s0 (\s -> let (a,s') = step s in (f a, s'))
  {-# INLINE fmap #-}

instance Applicative Series where
  pure a = Series () (\() -> (a,()))
  Series sA fA <*> Series sB fB = Series (sA,sB) $ \(!sa,!sb) ->
    let (a,sa') = fA sa
        (b,sb') = fB sb
    in (a b, (sa',sb'))
  {-# INLINE pure  #-}
  {-# INLINE (<*>) #-}

next :: Series a -> (a,Series a)
{-# INLINE next #-}
next (Series s f) = case f s of
  (a,s') -> (a,Series s' f)


----------------------------------------------------------------
-- Constructors
----------------------------------------------------------------

enumSeriesFrom :: Num a => a -> Series a
enumSeriesFrom i = Series i (\n -> (n,n+1))
{-# INLINE enumSeriesFrom #-}

enumSeriesFromStep :: Num a => a -> a -> Series a
enumSeriesFromStep n d = Series n (\i -> (i,i+d))
{-# INLINE enumSeriesFromStep #-}

scanSeries :: (b -> a -> b) -> b -> Series a -> Series b
{-# INLINE scanSeries #-}
scanSeries f b0 (Series s0 step) = Series (b0,s0) $ \(b,s) ->
  let (a,s') = step s
      b'     = f b a
  in (b,(b',s'))


----------------------------------------------------------------
-- Evaluation of series
----------------------------------------------------------------

sumSeries :: Series Double -> Double
{-# INLINE sumSeries #-}
sumSeries (Series sInit step)
  = go x0 s0
  where 
    (x0,s0) = step sInit
    go x s | abs (d/x) < m_epsilon = x'
           | otherwise             = go x' s'
      where (d,s') = step s
            x'     = x + d

sumPowerSeries :: Double -> Series Double -> Double
sumPowerSeries x ser = sumSeries $ liftA2 (*) (scanSeries (*) 1 (pure x)) ser
{-# INLINE sumPowerSeries #-}
            
seriesToList :: Series a -> [a]
seriesToList (Series s f) = unfoldr (Just . f) s



----------------------------------------------------------------
-- Evaluation of continued fractions
----------------------------------------------------------------

evalModLentz :: Series (Double,Double) -> Double
{-# INLINE evalModLentz #-}
evalModLentz (Series sInit step)
  = let ((_,b0),s0) = step sInit
        f0          = maskZero b0
    in  go f0 f0 0 s0
  where
    tiny = 1e-60
    maskZero 0 = tiny
    maskZero x = x
    
    go f c d s
      | abs (delta - 1) < m_epsilon = f'
      | otherwise                   = go f' c' d' s'
      where     
          ((a,b),s') = step s
          d'    = recip $ maskZero $ b + a*d
          c'    = maskZero $ b + a/c 
          delta = c'*d'
          f'    = f*delta

{'fraction_non_alphanumeric': 0.26695571384329514, 'fraction_numerical': 0.008206875193558377, 'mean_word_length': 1.430022573363431, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 55, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3900772', 'n_tokens_mistral': 2062, 'n_tokens_neox': 2008, 'n_words': 649}