This is your Haskel (not) only for factorials and fit

When it comes to favorite languages, I usually say that, other things being equal, I prefer C ++ for number of processors and haskel for everything else. It is useful to periodically check whether such a division is justified, and here one idle and very simple question arose recently: how will the sum of all the divisors of a number behave with the growth of this very number, say, for the first billion numbers. This task can simply be programmed (it’s already embarrassing to call the resulting number crusher), so it looks like a great option for such a check.

In addition, I still do not have the skill of accurately predicting the performance of a Haskel code, so it is useful to try obviously bad approaches to see how performance will degrade.

Well, in addition, you can easily show off a more efficient algorithm than the frontal search for divisors for each number from $$$1$ before $$$n$ .

Algorithm

So let's start with the algorithm.

How to find the sum of all divisors of a number $$$n$ ? You can go for all $$$k_1 \ in \ {1 \ dots \ lfloor \ sqrt n \ rfloor \}$ and for each such $$$k_1$ check the remainder of the division $$$n$ on $$$k_1$ . If the remainder is $$$0$ then add to the battery $$$k_1 + k_2$ where $$$k_2 = \ frac {n} {k_1}$ , if a $$$k_1 \ neq k_2$ , and just $$$k_1$ otherwise.

Is it possible to apply this algorithm $$$n$ times for each number from $$$1$ before $$$n$ ? Of course it is possible. What will be the complexity? Easy to see that order $$$O (n ^ \ frac {3} {2})$ divisions - for each number we do exactly the root-of-it divisions, and we have numbers $$$n$ . Can we do better? It turns out that yes.

One of the problems with this method is that we waste too much energy. Too many divisions do not lead us to success, giving a non-zero remainder. It is natural to try to be a little more lazy and approach the task from the other side: let's just generate all sorts of candidates for dividers and look at what numbers do they satisfy?

So let now we need in one fell swoop for each number from $$$1$ before $$$n$ calculate the sum of all its divisors. To do this, go through all $$$k_1 \ in \ {1 \ dots \ lfloor \ sqrt n \ rfloor \}$ and for each such $$$k_1$ go through all $$$k_2 \ in \ {k_1 \ dots \ lfloor \ frac {n} {k} \ rfloor \}$ . For each pair $$$(k_1, k_2)$ add to the cell with the index $$$k_1 \ cdot k_2$ value $$$k_1 + k_2$ , if a $$$k_1 \ neq k_2$ and $$$k_1$ otherwise.

This algorithm makes exactly $$$n ^ \ frac {1} {2}$ divisions, and each multiplication (which is cheaper than divisions) leads us to success: at each iteration we increase something. It is much more effective than the frontal approach.

In addition, with this very frontal approach, you can compare both implementations and make sure that they give the same results for fairly small numbers, which should add a bit of confidence.

First implementation

And, by the way, this is almost the pseudo-code of the initial implementation on Haskel:

module Divisors.Multi(divisorSums) where

import Data.IntMap.Strict as IM

divisorSums :: Int -> Int
divisorSums n = IM.fromListWith (+) premap IM.! n
  where premap = [ (k1 * k2, if k1 /= k2 then k1 + k2 else k1)
                 | k1 <- [ 1 .. floor $ sqrt $ fromIntegral n ]
                 , k2 <- [ k1 .. n `quot` k1 ]
                 ]

Main- , .

, $$$n$ . , — , ( ), , - .

? i7 3930k 100'000 0.4 . 0.15 0.25 — GC. 8 , , — 8 , 800 .

( ). , , ́? 1'000'000 7.5 , 4.5 GC, 80 ( 10 , ). Senior Java Software Developer' GC, . . , , : 64 , 80, .

,

C++

, , .

, , :

#include <vector>
#include <string>
#include <cmath>
#include <iostream>

int main(int argc, char **argv)
{
    if (argc != 2)
    {
        std::cerr << "Usage: " << argv[0] << " maxN" << std::endl;
        return 1;
    }
    int64_t n = std::stoi(argv[1]);

    std::vector<int64_t> arr;
    arr.resize(n + 1);

    for (int64_t k1 = 1; k1 <= static_cast<int64_t>(std::sqrt(n)); ++k1)
    {
        for (int64_t k2 = k1; k2 <= n / k1; ++k2)
        {
            auto val = k1 != k2 ? k1 + k2 : k1;
            arr[k1 * k2] += val;
        }
    }

    std::cout << arr.back() << std::endl;
}

-O3 -march=native, clang 8, 0.024 , 8 . — 155 , 8 , . . . . ! ?

, IntMap, , , — , (, , ). , C++?

:

module Divisors.Multi(divisorSums) where

import qualified Data.Array.IArray as A
import qualified Data.Array.Unboxed as A

divisorSums :: Int -> Int
divisorSums n = arr A.! n
  where arr = A.accumArray (+) 0 (1, n) premap :: A.UArray Int Int
        premap = [ (k1 * k2, if k1 /= k2 then k1 + k2 else k1)
                 | k1 <- [ 1 .. floor bound ]
                 , k2 <- [ k1 .. n `quot` k1 ]
                 ]
        bound = sqrt $ fromIntegral n :: Double

unboxed- , Int , . Boxed- arr, . , bound, , LICM, , defaulting' floor.

0.045 ( !). 8 , GC (!). — , C++, . ! ?

, . accumArray , — . accumArray unsafeAccumArray:

module Divisors.Multi(divisorSums) where

import qualified Data.Array.Base as A
import qualified Data.Array.IArray as A
import qualified Data.Array.Unboxed as A

divisorSums :: Int -> Int
divisorSums n = arr A.! (n - 1)
  where arr = A.unsafeAccumArray (+) 0 (0, n - 1) premap :: A.UArray Int Int
        premap = [ (k1 * k2 - 1, if k1 /= k2 then k1 + k2 else k1)
                 | k1 <- [ 1 .. floor bound ]
                 , k2 <- [ k1 .. n `quot` k1 ]
                 ]
        bound = sqrt $ fromIntegral n :: Double

, , (, , API , ). ?

— 0.021 (, , , !). , 8 , 0 GC.

— 152 (, !). 8 . 0 GC. - . , , .

-, , accumArray unsafe- . 10-20 ( , operator[] at() ), !

-, , , .

-, , , . , , . , , ( ) . LLVM JIT . , , .

-, : . unsafe , , k_1 * k_2 <= n k_1, k_2, . , . , , , , ( ), array .

-: , , . , . , - Go, Rust, Julia, D, Java, Malbolge , , C++- — , , .

P.S.: . .