math-functions 0.3.4.0 → 0.3.4.1
raw patch · 8 files changed
+374/−30 lines, 8 filesdep +gaugedep +randomdep ~data-default-classdep ~vectorPVP ok
version bump matches the API change (PVP)
Dependencies added: gauge, random
Dependency ranges changed: data-default-class, vector
API changes (from Hackage documentation)
Files
- Numeric/SpecFunctions/Internal.hs +63/−4
- bench/bench.hs +125/−0
- changelog.md +6/−0
- math-functions.cabal +28/−3
- tests/Tests/SpecFunctions.hs +42/−22
- tests/tables/generate.py +1/−1
- tests/tables/igamma.dat +99/−0
- tests/tables/inputs/igamma.dat +10/−0
Numeric/SpecFunctions/Internal.hs view
@@ -409,12 +409,21 @@ || (abs mu < 0.4) -- Gautschi's algorithm. --- -- Evaluate series for P(a,x). See [Temme1994] Eq. 5.5- --- -- FIXME: Term `exp (log x * z - x - logGamma (z+1))` doesn't give full precision+ -- Evaluate series for P(a,x). See [Temme1994] Eq. 5.5 and [NOTE:+ -- incompleteGamma.taylorP]+ factorP+ | a < 10 = x ** a+ / (exp x * exp (logGamma (a + 1)))+ | a < 1182.5 = (x * exp 1 / a) ** a+ / exp x+ / sqrt (2*pi*a)+ / exp (logGammaCorrection a)+ | otherwise = (x * exp 1 / a * exp (-x/a)) ** a+ / sqrt (2*pi*a)+ / exp (logGammaCorrection a) taylorSeriesP = sumPowerSeries x (scanSequence (/) 1 $ enumSequenceFrom (a+1))- * exp (log x * a - x - logGamma (a+1))+ * factorP -- Series for 1-Q(a,x). See [Temme1994] Eq. 5.5 taylorSeriesComplQ = sumPowerSeries (-x) (scanSequence (/) 1 (enumSequenceFrom 1) / enumSequenceFrom a)@@ -1322,3 +1331,53 @@ , 4.269068009004705e304 , 7.257415615307998e306 ]+++-- [NOTE: incompleteGamma.taylorP]+--+-- Incompltete gamma uses several algorithms for different parts of+-- parameter space. Most troublesome is P(a,x) Taylor series+-- [Temme1994,Eq.5.5] which requires to evaluate rather nasty+-- expression:+--+-- x^a x^a+-- ------------- = -------------+-- exp(x)·Γ(a+1) exp(x)·a·Γ(a)+--+-- Conditions:+-- | 0.5<x<1.1 = x < 4/3*a+-- | otherwise = x < a+--+-- For small `a` computation could be performed directly. However for+-- largish values of `a` it's possible some of factor in the+-- expression overflow. Values below take into account ranges for+-- Taylor P approximation:+--+-- · a > 155 - x^a could overflow+-- · a > 1182.5 - exp(x) could overflow+--+-- Usual way to avoid overflow problem is to perform calculations in+-- the log domain. It however doesn't work very well in this case+-- since we encounter catastrophic cancellations and could easily lose+-- up to 6(!) digits for large `a`.+--+-- So we take another approach and use Stirling approximation with+-- correction (logGammaCorrection).+--+-- x^a / x·e \^a 1+-- ≈ ------------------------- = | --- | · ----------------+-- exp(x)·sqrt(2πa)·(a/e)^a) \ a / exp(x)·sqrt(2πa)+--+-- We're using this approach as soon as logGammaCorrection starts+-- working (a>10) because we don't have implementation for gamma+-- function and exp(logGamma z) results in errors for large a.+--+-- Once we get into region when exp(x) could overflow we rewrite+-- expression above once more:+--+-- / x·e \^a 1+-- | --- · e^(-x/a) | · ---------+-- \ a / sqrt(2πa)+--+-- This approach doesn't work very well but it's still big improvement+-- over calculations in the log domain.
+ bench/bench.hs view
@@ -0,0 +1,125 @@+{-# LANGUAGE NumDecimals #-}+import Gauge.Main+import Data.Default.Class+import qualified Data.Vector.Unboxed as U+import Text.Printf+import System.Random (randomIO)++import qualified Numeric.Sum as Sum+import Numeric.SpecFunctions+import Numeric.Polynomial+import Numeric.RootFinding++++-- Uniformly sample logGamma performance between 10^-6 to 10^6+benchmarkLogGamma logG =+ [ bench (printf "%.3g" x) $ nf logG x+ | x <- [ m * 10**n | n <- [ -8 .. 8 ]+ , m <- [ 10**(i / tics) | i <- [0 .. tics-1] ]+ ]+ ]+ where tics = 3+{-# INLINE benchmarkLogGamma #-}+++-- Power of polynomial to be evaluated (In other words length of coefficients vector)+coef_size :: [Int]+coef_size = [ 1,2,3,4,5,6,7,8,9+ , 10, 30+ , 100, 300+ , 1000, 3000+ , 10000, 30000+ ]+{-# INLINE coef_size #-}++-- Precalculated coefficients+coef_list :: [U.Vector Double]+coef_list = [ U.replicate n 1.2 | n <- coef_size]+{-# NOINLINE coef_list #-}++++main :: IO ()+main = do+ v <- U.replicateM 1e6 randomIO :: IO (U.Vector Double)+ defaultMain+ [ bgroup "logGamma" $+ benchmarkLogGamma logGamma+ , bgroup "incompleteGamma" $+ [ bench (show p) $ nf (incompleteGamma p) p+ | p <- [ 0.1+ , 1, 3+ , 10, 30+ , 100, 300+ , 999, 1000+ ]+ ]+ , bgroup "factorial"+ [ bench (show n) $ nf factorial n+ | n <- [ 0, 1, 3, 6, 9, 11, 15+ , 20, 30, 40, 50, 60, 70, 80, 90, 100+ ]+ ]+ , bgroup "incompleteBeta"+ [ bench (show (p,q,x)) $ nf (incompleteBeta p q) x+ | (p,q,x) <- [ (10, 10, 0.5)+ , (101, 101, 0.5)+ , (1010, 1010, 0.5)+ , (10100, 10100, 0.5)+ , (100100, 100100, 0.5)+ , (1001000, 1001000, 0.5)+ , (10010000,10010000,0.5)+ ]+ ]+ , bgroup "log1p"+ [ bench (show x) $ nf log1p x+ | x <- [ -0.9+ , -0.5+ , -0.1+ , 0.1+ , 0.5+ , 1+ , 10+ , 100+ ] :: [Double]+ ]+ , bgroup "sinc" $+ bench "sin" (nf sin (0.55 :: Double))+ : [ bench (show x) $ nf sinc x+ | x <- [0, 1e-6, 1e-3, 0.5]+ ]+ , bgroup "erf & Co"+ [ bgroup "erf"+ [ bench (show x) $ nf erf x+ | x <- [0, 1.1, 100, 1000]+ ]+ , bgroup "erfc"+ [ bench (show x) $ nf erfc x+ | x <- [0, 1.1, 100, 1000]+ ]+ , bgroup "invErfc"+ [ bench (show x) $ nf erfc x+ | x <- [1e-9, 1e-6, 1e-3, 0.1, 1]+ ]+ ]+ , bgroup "expm1"+ [ bench (show x) $ nf expm1 (x :: Double)+ | x <- [-0.1, 0, 1, 19]+ ]+ , bgroup "poly"+ $ [ bench ("vector_"++show (U.length coefs)) $ nf (\x -> evaluatePolynomial x coefs) (1 :: Double)+ | coefs <- coef_list ]+ ++ [ bench ("unpacked_"++show n) $ nf (\x -> evaluatePolynomialL x (map fromIntegral [1..n])) (1 :: Double)+ | n <- coef_size ]+ , bgroup "RootFinding"+ [ bench "ridders sin" $ nf (ridders def (0,pi/2)) (\x -> sin x - 0.525)+ , bench "newton sin" $ nf (newtonRaphson def (0,1.2,pi/2)) (\x -> (sin x - 0.525,cos x))+ ]+ , bgroup "Sum"+ [ bench "naive" $ whnf U.sum v+ , bench "kahan" $ whnf (Sum.sumVector Sum.kahan) v+ , bench "kbn" $ whnf (Sum.sumVector Sum.kbn) v+ , bench "kb2" $ whnf (Sum.sumVector Sum.kb2) v+ ]+ ]
changelog.md view
@@ -1,7 +1,13 @@+## Changes in 0.3.4.1++ * Precision of `incompleteGamma` improved.++ ## Changes in 0.3.4.0 * Dependency on `vector-th-unbox` is dropped. All instances are written by hand now.+ ## Changes in 0.3.3.0
math-functions.cabal view
@@ -1,5 +1,5 @@ name: math-functions-version: 0.3.4.0+version: 0.3.4.1 cabal-version: >= 1.10 license: BSD2 license-file: LICENSE@@ -74,7 +74,7 @@ build-depends: base >= 4.5 && < 5 , deepseq , data-default-class >= 0.1.2.0- , vector >= 0.7+ , vector >= 0.11 , primitive if flag(system-expm1) && !os(windows) cpp-options: -DUSE_SYSTEM_EXPM1@@ -94,7 +94,7 @@ other-modules: Numeric.SpecFunctions.Compat -test-suite tests+test-suite math-function-tests default-language: Haskell2010 other-extensions: ViewPatterns @@ -127,6 +127,31 @@ , tasty >= 1.2 , tasty-hunit >= 0.10 , tasty-quickcheck >= 0.10++benchmark math-functions-bench+ type: exitcode-stdio-1.0+ if impl(ghc <= 7.10 ) || impl(ghcjs)+ buildable: False+ default-language: Haskell2010+ other-extensions:+ BangPatterns+ CPP+ DeriveDataTypeable+ FlexibleContexts+ MultiParamTypeClasses+ ScopedTypeVariables+ TemplateHaskell+ TypeFamilies+ DeriveGeneric+ ghc-options: -Wall -O2+ hs-source-dirs: bench+ Main-is: bench.hs+ build-depends: base >= 4.5 && < 5+ , math-functions+ , data-default-class+ , vector+ , random+ , gauge >=0.2.5 source-repository head type: git
tests/Tests/SpecFunctions.hs view
@@ -8,6 +8,7 @@ import Control.Monad import Data.List+import Data.Maybe import qualified Data.Vector as V import Data.Vector ((!)) import qualified Data.Vector.Unboxed as U@@ -35,6 +36,13 @@ erfcLargeTol = 64 #endif +isGHCJS :: Bool+#if defined(__GHCJS__)+isGHCJS = True+#else+isGHCJS = False+#endif+ tests :: TestTree tests = testGroup "Special functions" [ testGroup "erf"@@ -42,14 +50,14 @@ -- large arguments testCase "erfc table" $ forTable "tests/tables/erfc.dat" $ \[x, exact] ->- checkTabular erfcTol (show x) exact (erfc x)+ checkTabularPure erfcTol (show x) exact (erfc x) , testCase "erfc table [large]" $ forTable "tests/tables/erfc-large.dat" $ \[x, exact] ->- checkTabular erfcLargeTol (show x) exact (erfc x)+ checkTabularPure erfcLargeTol (show x) exact (erfc x) -- , testCase "erf table" $ forTable "tests/tables/erf.dat" $ \[x, exact] -> do- checkTabular erfTol (show x) exact (erf x)+ checkTabularPure erfTol (show x) exact (erf x) , testProperty "id = erfc . invErfc" invErfcIsInverse , testProperty "id = invErfc . erfc" invErfcIsInverse2 , testProperty "invErf = erf^-1" invErfIsInverse@@ -58,16 +66,16 @@ , testGroup "log1p & Co" [ testCase "expm1 table" $ forTable "tests/tables/expm1.dat" $ \[x, exact] ->- checkTabular 2 (show x) exact (expm1 x)+ checkTabularPure 2 (show x) exact (expm1 x) , testCase "log1p table" $ forTable "tests/tables/log1p.dat" $ \[x, exact] ->- checkTabular 1 (show x) exact (log1p x)+ checkTabularPure 1 (show x) exact (log1p x) ] ---------------- , testGroup "gamma function" [ testCase "logGamma table [fractional points" $ forTable "tests/tables/loggamma.dat" $ \[x, exact] -> do- checkTabular 2 (show x) exact (logGamma x)+ checkTabularPure 2 (show x) exact (logGamma x) , testProperty "Gamma(x+1) = x*Gamma(x)" $ gammaReccurence , testCase "logGamma is expected to be precise at 1e-15 level" $ forM_ [3..10000::Int] $ \n -> do@@ -79,7 +87,11 @@ , testGroup "incomplete gamma" [ testCase "incompleteGamma table" $ forTable "tests/tables/igamma.dat" $ \[a,x,exact] -> do- checkTabular 16 (show (a,x)) exact (incompleteGamma a x)+ let err | a < 10 = 16+ | a <= 101 = if isGHCJS then 64 else 32+ | a == 201 = 200+ | otherwise = 32+ checkTabularPure err (show (a,x)) exact (incompleteGamma a x) , testProperty "incomplete gamma - increases" $ \(abs -> s) (abs -> x) (abs -> y) -> s > 0 ==> monotonicallyIncreases (incompleteGamma s) x y , testProperty "0 <= gamma <= 1" incompleteGammaInRange@@ -89,7 +101,7 @@ ---------------- , testGroup "beta function" [ testCase "logBeta table" $- forTable "tests/tables/logbeta.dat" $ \[p,q,exact] -> do+ forTable "tests/tables/logbeta.dat" $ \[p,q,exact] -> let errEst -- For Stirling approx. errors are very good | b > 10 = 2@@ -109,9 +121,7 @@ est = ceiling $ abs (logGamma a) + abs (logGamma b) + abs (logGamma (a + b)) / abs (logBeta a b)--- checkTabular errEst (show (p,q)) exact (logBeta p q)+ in checkTabularPure errEst (show (p,q)) exact (logBeta p q) , testCase "logBeta factorial" betaFactorial , testProperty "beta(1,p) = 1/p" beta1p -- , testProperty "beta recurrence" betaRecurrence@@ -137,7 +147,7 @@ -- Relative precision is lost when digamma(x) ≈ 0 , testCase "digamma is expected to be precise at 1e-12" $ forTable "tests/tables/digamma.dat" $ \[x, exact] ->- checkTabular 2048+ checkTabularPure 2048 (show x) (digamma x) exact ] ----------------@@ -159,7 +169,7 @@ $ U.length factorialTable == 171 , testCase "Log factorial table" $ forTable "tests/tables/factorial.dat" $ \[i,exact] ->- checkTabular 3+ checkTabularPure 3 (show i) (logFactorial (round i :: Int)) exact ] ----------------@@ -435,16 +445,26 @@ = fmap (fmap (fmap read . words) . lines) . readFile -forTable :: FilePath -> ([Double] -> IO ()) -> IO ()+forTable :: FilePath -> ([Double] -> Maybe String) -> IO () forTable path fun = do- mapM_ fun =<< readTable path+ rows <- readTable path+ case mapMaybe fun rows of+ [] -> return ()+ errs -> assertFailure $ intercalate "---\n" errs checkTabular :: Int -> String -> Double -> Double -> IO () checkTabular prec x exact val =- assertBool (unlines [ " x = " ++ x- , " expected = " ++ show exact- , " got = " ++ show val- , " ulps diff = " ++ show (ulpDistance exact val)- , " err.est. = " ++ show prec- ])- (within prec exact val)+ case checkTabularPure prec x exact val of+ Nothing -> return ()+ Just s -> assertFailure s++checkTabularPure :: Int -> String -> Double -> Double -> Maybe String+checkTabularPure prec x exact val+ | within prec exact val = Nothing+ | otherwise = Just $ unlines+ [ " x = " ++ x+ , " expected = " ++ show exact+ , " got = " ++ show val+ , " ulps diff = " ++ show (ulpDistance exact val)+ , " err.est. = " ++ show prec+ ]
tests/tables/generate.py view
@@ -1,4 +1,4 @@-#!/usr/bin/python+#!/usr/bin/env python3 """ """ import itertools
tests/tables/igamma.dat view
@@ -7,6 +7,7 @@ 0.000100000000000000005 10 0.999999999584179852012789127241 0.000100000000000000005 100 1.0 0.000100000000000000005 1000 1.0+0.000100000000000000005 3301 1.0 0.00100000000000000002 9.99999999999999955e-07 0.986848133694076665339510700432 0.00100000000000000002 1.00000000000000008e-05 0.989123044695782668850940465364 0.00100000000000000002 0.00100000000000000002 0.993687646708860290096219637037@@ -16,6 +17,7 @@ 0.00100000000000000002 10 0.999999995830692182809738396874 0.00100000000000000002 100 1.0 0.00100000000000000002 1000 1.0+0.00100000000000000002 3301 1.0 0.0100000000000000002 9.99999999999999955e-07 0.875933759832353305070970748339 0.0100000000000000002 1.00000000000000008e-05 0.896336798267197200971193210546 0.0100000000000000002 0.00100000000000000002 0.938570652526128985382985766694@@ -25,6 +27,7 @@ 0.0100000000000000002 10 0.99999995718295022590061122615 0.0100000000000000002 100 1.0 0.0100000000000000002 1000 1.0+0.0100000000000000002 3301 1.0 0.100000000000000006 9.99999999999999955e-07 0.264033654327922324857187514752 0.100000000000000006 1.00000000000000008e-05 0.332398405040503295401903974218 0.100000000000000006 0.00100000000000000002 0.526768568392445111817869842601@@ -34,6 +37,7 @@ 0.100000000000000006 10 0.999999445201428209809392042838 0.100000000000000006 100 1.0 0.100000000000000006 1000 1.0+0.100000000000000006 3301 1.0 0.200000000000000011 9.99999999999999955e-07 0.0687190937987684780583487585964 0.200000000000000011 1.00000000000000008e-05 0.108912260585591831357806184556 0.200000000000000011 0.00100000000000000002 0.273530102033034019134013607911@@ -43,6 +47,7 @@ 0.200000000000000011 10 0.999998540143010797930830697127 0.200000000000000011 100 1.0 0.200000000000000011 1000 1.0+0.200000000000000011 3301 1.0 0.299999999999999989 9.99999999999999955e-07 0.0176595495901936665671778397068 0.299999999999999989 1.00000000000000008e-05 0.0352353606155625769158079986437 0.299999999999999989 0.00100000000000000002 0.140242458924867370902963425856@@ -52,6 +57,7 @@ 0.299999999999999989 10 0.99999715515533278951644372428 0.299999999999999989 100 1.0 0.299999999999999989 1000 1.0+0.299999999999999989 3301 1.0 0.400000000000000022 9.99999999999999955e-07 0.00448690737698480048018769908146 0.400000000000000022 1.00000000000000008e-05 0.0112705727782256817005130809278 0.400000000000000022 0.00100000000000000002 0.0710923978953330760218610475958@@ -61,6 +67,7 @@ 0.400000000000000022 10 0.999995127544578487259167603836 0.400000000000000022 100 1.0 0.400000000000000022 1000 1.0+0.400000000000000022 3301 1.0 0.5 9.99999999999999955e-07 0.00112837879096923635441785924383 0.5 1.00000000000000008e-05 0.00356823633818045042058077366094 0.5 0.00100000000000000002 0.0356705917296798854171108747554@@ -70,6 +77,7 @@ 0.5 10 0.999992255783568955916362323619 0.5 100 1.0 0.5 1000 1.0+0.5 3301 1.0 0.599999999999999978 9.99999999999999955e-07 0.000281123932739927969642852394278 0.599999999999999978 1.00000000000000008e-05 0.00111917075717695837092315556803 0.599999999999999978 0.00100000000000000002 0.0177310780570833845378940622156@@ -79,3 +87,94 @@ 0.599999999999999978 10 0.999988293084421631090774556099 0.599999999999999978 100 1.0 0.599999999999999978 1000 1.0+0.599999999999999978 3301 1.0+2 9.99999999999999955e-07 4.9999966666679162138149041803e-13+2 1.00000000000000008e-05 4.99996666679166715135638665241e-11+2 0.00100000000000000002 0.000000499666791633340297383350611252+2 0.0100000000000000002 0.000049667913340265892415918274953+2 0.100000000000000006 0.00467884016044447002161170213187+2 1 0.264241117657115356808952459677+2 10 0.999500600772612666633108493329+2 100 1.0+2 1000 1.0+2 3301 1.0+3 9.99999999999999955e-07 1.6666654166671664402685631963e-19+3 1.00000000000000008e-05 1.66665416671666693678925483422e-16+3 0.00100000000000000002 1.66541716652780763845435502936e-10+3 0.0100000000000000002 0.000000165421652807487686572525438025+3 0.100000000000000006 0.000154653070264671678619072687835+3 1 0.0803013970713941960111905745964+3 10 0.997230604284488424056328917551+3 100 1.0+3 1000 1.0+3 3301 1.0+6 9.99999999999999955e-07 1.38888769841321886877873406008e-39+6 1.00000000000000008e-05 1.38887698417906798768211894792e-33+6 0.00100000000000000002 1.38769893337745993330072800382e-21+6 0.0100000000000000002 1.37703605634306470721227911913e-15+6 0.100000000000000006 0.00000000127489869222979188498133341591+6 1 0.000594184817581692998827091061365+6 10 0.932914037120968217714240937173+6 100 1.0+6 1000 1.0+6 3301 1.0+12 9.99999999999999955e-07 2.08767377170244306355144979982e-81+12 1.00000000000000008e-05 2.08765642802367929860815380041e-69+12 0.00100000000000000002 2.08574950796625691581696824468e-45+12 0.0100000000000000002 2.06849404029269947208812650771e-33+12 0.100000000000000006 1.90364240064062767847395955076e-21+12 1 8.31610742688233390952391737168e-10+12 10 0.30322385369689331180582927404+12 100 0.99999999999999999999999999999+12 1000 1.0+12 3301 1.0+18 9.99999999999999955e-07 1.56191921714497984569652172566e-124+18 1.00000000000000008e-05 1.56190589978546723049068336087e-106+18 0.00100000000000000002 1.56044168515546614690893456571e-70+18 0.0100000000000000002 1.5471936172459858931960924774e-52+18 0.100000000000000006 1.42075999849733529123506295102e-34+18 1 6.06428067721557331946176715382e-17+18 10 0.0142776135970496129089830965836+18 100 0.999999999999999999999998742916+18 1000 1.0+18 3301 1.0+101 9.99999999999999955e-07 1.06090022487198225461869933883e-766+101 1.00000000000000008e-05 1.06089077042094832542309156458e-665+101 0.00100000000000000002 1.05985129506822713235416044578e-463+101 0.0100000000000000002 1.05044811631749169413217529395e-362+101 0.100000000000000006 9.60885206142996276573449906708e-262+101 1 3.94147589063752014615062291541e-161+101 10 5.33940546071971052337450543217e-64+101 100 0.473437801470001529623393607105+101 1000 1.0+101 3301 1.0+201 9.99999999999999955e-07 6.30833677503352627325382205745e-1584+201 1.00000000000000008e-05 6.30828028132019299365575041933e-1383+201 0.00100000000000000002 6.30206906057329746426353105998e-981+201 0.0100000000000000002 6.24588319207081699412164871293e-780+201 0.100000000000000006 5.71085198693835010992991011704e-579+201 1 2.33225525188187904433135775933e-378+201 10 3.01310889066541622340100430602e-181+201 100 4.62617947019577288096442044516e-19+201 1000 1.0+201 3301 1.0+1000 9.99999999999999955e-07 2.4851656605824546988202041697e-8568+1000 1.00000000000000008e-05 2.48514331653642003877116895294e-7568+1000 0.00100000000000000002 2.48268669750003876617472467308e-5568+1000 0.0100000000000000002 2.46046488714861784228553201292e-4568+1000 0.100000000000000006 2.24889779123046964018322233315e-3568+1000 1 9.15156509116491575720545785542e-2569+1000 10 1.13964958763725134069338684487e-1572+1000 100 1.0270971815582277877665615585e-611+1000 1000 0.504205244180215508503777843602+1000 3301 1.0+3003 9.99999999999999955e-07 8.90812855529069811196158732723e-27160+3003 1.00000000000000008e-05 8.90804840918643793907789460332e-24157+3003 0.00100000000000000002 8.89923673804629728415411068459e-18151+3003 0.0100000000000000002 8.81952937102869844156580901528e-15148+3003 0.100000000000000006 8.0606844309353722733298755906e-12145+3003 1 3.27821191297260784144201433956e-9142+3003 10 4.05779611158376897690039719814e-6143+3003 100 3.42800829838332693574242322961e-3179+3003 1000 6.77752709164469153719156665409e-567+3003 3301 0.999999933219976432069280011406
tests/tables/inputs/igamma.dat view
@@ -8,6 +8,15 @@ 0.4 0.5 0.6+2+3+6+12+18+101+201+1000+3003 x = 1e-6@@ -19,3 +28,4 @@ 10 100 1000+3301