diff --git a/Numeric/SpecFunctions/Internal.hs b/Numeric/SpecFunctions/Internal.hs --- a/Numeric/SpecFunctions/Internal.hs +++ b/Numeric/SpecFunctions/Internal.hs @@ -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 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. diff --git a/bench/bench.hs b/bench/bench.hs new file mode 100644 --- /dev/null +++ b/bench/bench.hs @@ -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 + ] + ] diff --git a/changelog.md b/changelog.md --- a/changelog.md +++ b/changelog.md @@ -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 diff --git a/math-functions.cabal b/math-functions.cabal --- a/math-functions.cabal +++ b/math-functions.cabal @@ -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 diff --git a/tests/Tests/SpecFunctions.hs b/tests/Tests/SpecFunctions.hs --- a/tests/Tests/SpecFunctions.hs +++ b/tests/Tests/SpecFunctions.hs @@ -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 + ] diff --git a/tests/tables/generate.py b/tests/tables/generate.py --- a/tests/tables/generate.py +++ b/tests/tables/generate.py @@ -1,4 +1,4 @@ -#!/usr/bin/python +#!/usr/bin/env python3 """ """ import itertools diff --git a/tests/tables/igamma.dat b/tests/tables/igamma.dat --- a/tests/tables/igamma.dat +++ b/tests/tables/igamma.dat @@ -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 diff --git a/tests/tables/inputs/igamma.dat b/tests/tables/inputs/igamma.dat --- a/tests/tables/inputs/igamma.dat +++ b/tests/tables/inputs/igamma.dat @@ -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