Step by Step solution about How to resolve the algorithm Chowla numbers step by step in the Haskell programming language

The Haskell code you provided implements the computation of the Chowla function for a range of numbers, up to 35 million.

The chowla function takes a positive integer n as input and returns the value of the Chowla function for that number. The Chowla function is defined as follows:

chowla(n) = n - product(sumFactor(p^e) for (p, e) in factorisation(n))

where factorisation(n) is the prime factorisation of n, sumFactor(p^e) is the sum of p^k for k in [1, e], and product is the product of the values obtained by applying sumFactor to each prime factor of n.

The chowlas function takes a list of positive integers as input and returns a list of tuples, where each tuple consists of a number and the value of the Chowla function for that number. The chowlas function uses parallel processing to compute the Chowla function for each number in the input list.

The chowlaPrimes function takes a list of tuples, where each tuple consists of a number and the value of the Chowla function for that number, and a range of numbers as input. The chowlaPrimes function returns the number of primes in the given range that have a Chowla function value that is equal to the second element of the given tuple.

The chowlaPerfects function takes a list of tuples, where each tuple consists of a number and the value of the Chowla function for that number, as input. The chowlaPerfects function returns a list of the numbers in the input list that are perfect numbers. A perfect number is a number that is equal to the sum of its proper divisors.

The commas function takes a positive integer as input and returns a string representation of that number with commas separating the thousands, millions, etc.

The main function is the entry point of the program. The main function first sets the number of threads to use to the number of cores available on the computer. Then, the main function computes the Chowla function for the first 37 numbers, the number of primes in several ranges of numbers that have a Chowla function value that is equal to the Chowla function value of the perfect number closest to the upper bound of the range, and the list of perfect numbers less than 35 million. Finally, the main function prints the results to the console.

import Control.Concurrent       (setNumCapabilities)
import Control.Monad.Par        (runPar, get, spawnP)
import Control.Monad            (join, (>=>))
import Data.List.Split          (chunksOf)
import Data.List                (intercalate, mapAccumL, genericTake, genericDrop)
import Data.Bifunctor           (bimap)
import GHC.Conc                 (getNumProcessors)
import Math.NumberTheory.Primes (factorise, unPrime)
import Text.Printf              (printf)

chowla :: Word -> Word
chowla 1 = 0
chowla n = f n
  where
    f = (-) =<< pred . product . fmap sumFactor . factorise
    sumFactor (n, e) = foldr (\p s -> s + unPrime n^p) 1 [1..e]

chowlas :: [Word] -> [(Word, Word)]
chowlas [] = []
chowlas xs = runPar $ join <$>
  (mapM (spawnP . fmap ((,) <*> chowla)) >=> mapM get) (chunksOf (10^6) xs)

chowlaPrimes :: [(Word, Word)] -> (Word, Word) -> (Word, Word)
chowlaPrimes chowlas range = (count chowlas, snd range)
  where
    isPrime (1, n) = False
    isPrime (_, n) = n == 0
    count = fromIntegral . length . filter isPrime . between range
    between (min, max) = genericTake (max - pred min) . genericDrop (pred min)

chowlaPerfects :: [(Word, Word)] -> [Word]
chowlaPerfects = fmap fst . filter isPerfect
  where 
    isPerfect (1, _) = False
    isPerfect (n, c) = c == pred n

commas :: (Show a, Integral a) => a -> String
commas = reverse . intercalate "," . chunksOf 3 . reverse . show

main :: IO ()
main = do
  cores <- getNumProcessors
  setNumCapabilities cores
  printf "Using %d cores\n" cores

  mapM_ (uncurry (printf "chowla(%2d) = %d\n")) $ take 37 allChowlas
  mapM_ (uncurry (printf "There are %8s primes < %10s\n"))
    (chowlaP
       [ (1, 10^2)
       , (succ $ 10^2, 10^3)
       , (succ $ 10^3, 10^4)
       , (succ $ 10^4, 10^5)
       , (succ $ 10^5, 10^6)
       , (succ $ 10^6, 10^7) ])

  mapM_ (printf "%10s is a perfect number.\n" . commas) perfects
  printf "There are %2d perfect numbers < 35,000,000\n" $ length perfects
  where
    chowlaP = fmap (bimap commas commas) . snd
      . mapAccumL (\total (count, max) -> (total + count, (total + count, max))) 0
      . fmap (chowlaPrimes $ take (10^7) allChowlas)
    perfects = chowlaPerfects allChowlas
    allChowlas = chowlas [1..35*10^6]