numerics-aez-notes

firstDerivative :: Floating a => a -> (a -> a) -> a -> a
firstDerivative h f x =
  (-f (x + 2 * h) + 8 * f (x + h) - 8 * f (x - h) + f (x - 2 * h)) / (12 * h)

secondDerivative :: Floating a => a -> (a -> a) -> a -> a
secondDerivative h f x =
  (-f (x + 2 * h) + 16 * f (x + h) - 30 * f x + 16 * f (x - h) - f (x - 2 * h)) / (12 * h ** 2)

logSumExp :: (Floating a, Ord a) => [a] -> a
logSumExp xs = x' + log (sum [exp (x - x') | x <- xs])
               where x' = maximum xs

----------------------------------------------------------------
-- Combinatorics
----------------------------------------------------------------

-- |
-- Quickly compute the natural logarithm of /n/ @`choose`@ /k/, with
-- no checking.
--
-- Less numerically stable:
--
-- > exp $ lg (n+1) - lg (k+1) - lg (n-k+1)
-- >   where lg = logGamma . fromIntegral
logChooseFast :: Double -> Double -> Double
logChooseFast n k = -log (n + 1) - logBeta (n - k + 1) (k + 1)

-- | Calculate binomial coefficient using exact formula
chooseExact :: Int -> Int -> Double
n `chooseExact` k
  = U.foldl' go 1 $ U.enumFromTo 1 k
  where
    go a i      = a * (nk + j) / j
        where j = fromIntegral i :: Double
    nk = fromIntegral (n - k)

-- | Compute logarithm of the binomial coefficient.
logChoose :: Int -> Int -> Double
n `logChoose` k
    | k  > n    = (-1) / 0
      -- For very large N exact algorithm overflows double so we
      -- switch to beta-function based one
    | k' < 50 && (n < 20000000) = log $ chooseExact n k'
    | otherwise                 = logChooseFast (fromIntegral n) (fromIntegral k)
  where
    k' = min k (n-k)

-- | Compute the binomial coefficient /n/ @\``choose`\`@ /k/. For
-- values of /k/ > 50, this uses an approximation for performance
-- reasons.  The approximation is accurate to 12 decimal places in the
-- worst case
--
-- Example:
--
-- > 7 `choose` 3 == 35
choose :: Int -> Int -> Double
n `choose` k
    | k  > n         = 0
    | k' < 50        = chooseExact n k'
    | approx < max64 = fromIntegral . round64 $ approx
    | otherwise      = approx
  where
    k'             = min k (n-k)
    approx         = exp $ logChooseFast (fromIntegral n) (fromIntegral k')
    max64          = fromIntegral (maxBound :: Int64)
    round64 x      = round x :: Int64

df <- data.frame(t = c(38,106,255,847), n = c(180,190,200,210))

ggplot(data = df, mapping = aes(x = n, y = t / 60)) +
  geom_point() +
  scale_y_log10(n.breaks = 10) +
  geom_smooth(method = "lm", fullrange = TRUE) +
  xlim(c(170, 220))

public class AAA {
    public static void main(String[] args) {
        double size = 100000;
        double answer = 1051287.7089736538;
        int repeats = Integer.parseInt(args[0]);
        for (int m = 0; m < repeats; m++) {
            double x = 0;
            for (double n = 1.0; n < size; n+=1.0) {
                x += Math.log(n);
            }
            if (x != answer) {
                throw new RuntimeException("incorrect value");
            }
        }
    }
}

public class BBB {
    public static void main(String[] args) {
        int size = 100000;
        double answer = 1051287.7089736538;
        int repeats = Integer.parseInt(args[0]);
        double[] ls;
        ls = new double[size+1];
        for (int ix = 0; ix <= size; ix++) {
            ls[ix] = Math.log(ix);
        }

        for (int m = 0; m < repeats; m++) {
            double x = 0;
            for (int n = 1; n < size; n++) {
                x += ls[n];
            }
            if (x != answer) {
                throw new RuntimeException("incorrect value");
            }
        }
    }
}

method,repeats,seconds
AAA,100,0.286
AAA,200,0.474
AAA,300,0.612
AAA,400,0.781
AAA,500,0.945
AAA,600,1.154
AAA,700,1.312
AAA,800,1.521
AAA,900,1.712
AAA,1000,1.764
AAA,1100,2.084
AAA,1200,2.177
AAA,1300,2.397
AAA,1400,2.619
AAA,1500,2.661
BBB,100,0.133
BBB,200,0.128
BBB,300,0.156
BBB,400,0.154
BBB,500,0.169
BBB,600,0.174
BBB,700,0.177
BBB,800,0.186
BBB,900,0.178
BBB,1000,0.240
BBB,1100,0.205
BBB,1200,0.229
BBB,1300,0.248
BBB,1400,0.254
BBB,1500,0.272