Count numbers, coprime with a given one

[Bodigrim]

It would be desirable to have an efficient implementation of the following function:

countCoprimes :: (Enum a, Euclidean a) => a -> a -> Int
countCoprimes n m = length $ filter (coprime n) [0..m]

or even better this one:

countCoprimes :: (Enum a, Euclidean a) => [Prime a] -> a -> Int 
countCoprimes ps m = length $ filter (\i -> all (coprime i . unPrime) ps) [0..m]

emphasing the fact that only prime factors of n matter.

[Bodigrim]

Here is a draft implementation plan for anyone willing to take a stab:

1. Implement in an efficient way:

   foo :: Int -> [a] -> [[a]]
   foo k = filter ((== k) . length) . Data.List.subsequences 

   Bonus points if it is productive on an infinite input list.
   Extra bonus points if it satisfies a prefix property:

   foo k (take n xs) == take (binomial !! n !! k) (foo k xs) 

2. Define

   bar :: Num a => Int -> [a] -> [a]
   bar k = map product . foo k  

   For example, bar 2 [2, 3, 5] = [6, 10, 15].

3. Implement countCoprimes n by finding prime factors of n and using inclusion-exclusion principle.

[CarlEdman]

Well, here is a trivial, but I think probably reasonably efficient, implementation of foo:

foo1 = go
  where
    go 0 _      = [[]]
    go _ []     = []
    go k (l:ls) = map (l:) (go (k-1) (ls)) ++ go k ls

While this is correct and fulfills the bonus property, it is a little dull on infinite input lists. The following implementation is only slightly more complex and has the additional property that any given subsequence of an infinite list will appear in finite time:

foo2 = go
  where
    go 0 _      = [[]]
    go _ []     = []
    go k (l:ls) = map (l:) (go (k-1) (ls)) `alternate` go k ls

    alternate (a:as) b = a:alternate b as
    alternate _ b      = b

Still thinking about how to meet the extra bonus property. It feels like the sort of thing which has a blindingly simple and efficient implementation, but every one which comes to mind fails to meet one or two of those criteria.

[rockbmb]

@Bodigrim isn't this Euler's totient function?

Count numbers, coprime with a given one

Isn't this an infinite quantity unless we bound it by e.g. asking only for numbers up to the given one, thereby arriving at the definition of Euler's totient function?
I think I'm missing something.

[b-mehta]

@rockbmb countCoprimes n m should be the numbers between 0 and m which are coprime to n, so countCoprimes n n == totient n. That said, countCoprimes n m should be a simple function of countCoprimes n (mod n m) and totient n.

[CarlEdman]

I happened to need something like this for another problem, so I constructed a foo which also fulfills the extra bonus property

import Data.List (scanl')

foo3 :: [a] -> [[[a]]]
foo3 = go 0 [[]]
  where
    go k p as
      | null as   = [p]
      | otherwise = p:(go (k+1) [q:qs|(q,b)<-zip as bs,qs<-take b p] (tail as))
      where
        bs = scanl' (\b i -> quot (b*(k+i)) i) 1 [1..]

Note that foo3 takes a single argument, a potentially infinite list. Its result is a list of list of subsets of its argument. In other words, rather than foo k xs, you do (foo xs)!!k.

Note also that the elements in the subsets appear in reverse order from the original list, but that ensures maximum sharing while at the same time still allowing infinite input lists. If this is a problem and you are dealing with finite lists, the best way is to reverse the input list before passing it to foo.

Note finally that each of the list of subsets is already sorted in order of increasing value if by value of a subset you define the integer value of the bit pattern corresponding to the subset (interpreted little-endian). That is useful for my task.

I knew there was a semi-elegant way of getting this!

[Bodigrim]

Sounds amazing! Sorry, I am a bit overloaded at the moment, will take a look next week.

[CarlEdman]

On of my OCDs is that once I get set on a simple Haskell problem which ought to be encapsulated in a reusable function, my mind just gnaws at it until it has achieved its most pure form. In that spirit, I offer foo4 which should be identical to foo3 in output, but is a bit more concise:

foo4 :: [a] -> [[[a]]]
foo4 = scanl' go [[]] . zip pascalTriangleTranspose . init . tails
  where
    go p (bs,as) = [q:qs|(q,b)<-zip as bs,qs<-take b p]

For pascalTriangleTranspose a generally useful function, you can use one of these:

[scanl' (\b i -> quot (b*(k+i)) i) 1 [1..]|k<-[0..]]

or

transpose pascalTriangle
pascalTriangle = iterate (\l -> zipWith (+) (l++[0]) (0:l)) [1]

or

pascalTriangleTranspose = iterate (tail . scanl' (+) 0) (repeat 1)

[rockbmb]

@CarlEdman for pascalTriangleTranspose, the third alternative seems clearer, although I cannot comment on its performance.

[CarlEdman]

@rockbmb Yes, the third is my favorite too. It is clean and simple and performs just a single integer addition to generate each element, while the first one looks messier and does a multiply and a divide for every element.

And yet, when I benchmarked them against each other, the first one was marginally, but consistently faster. Otherwise I would not even have included it.

The difference is too small to decide between the versions, but was a useful, repeat lessons for old hackers like me that the hardware our code runs on these days does really fast integer and floating-point multiplies and divides and hence the habit of trying to avoid them has become obsolete. Meanwhile, harder-to-predict factors, like cache hits and locality of reference, are huge for performance.

[rockbmb]

@CarlEdman a nice alternative would be to use the first, but add a comment above it that says: Equivalent to iterate (tail . scanl' (+) 0) (repeat 1). One gets simplicity in the explanation, and increased performance in practice.

[Bodigrim]

This is really beautiful!

The choice between first and third implementations of pascalTriangleTranspose is the same as for binomial (https://github.com/cartazio/arithmoi/blob/master/Math/NumberTheory/Recurrences/Bilinear.hs#L66). The third one looks nicer, but you need to retain the whole structure in memory, while in the first one sublists are completely independent.

[CarlEdman]

Very nice! And I should have used binomial. My near-identical homebrew pascalTriangle is a leftover from before arithmoi had the function. I'm slowly replacing the functions in my own library with usually superior arithmoi functions, but had missed that.

Two points:

1. Could I suggest adding a binomialTranspose function (from one of the above) to Recurrences.Bilinear? It comes in useful sometimes, as here, and the most efficient solution takes a little while to figure out.

2. One point in favor of the simpler/"elegant" implementations of both: They can be implemented without div and hence with just a Num, rather than Integral constraint. That in turn means that they would work fine with Mods as implemented in Math.NumberTheory.Moduli.Class, which are Nums, but not Integrals. Alternatively, one could add a Integral instance for Mods. That would be convenient, but I'm torn about whether it would be right.

[Bodigrim]

1. Yeah, it seems reasonable. I am thinking about having (note the difference in constraints)

binomial :: Num a => [[a]]
binomialTransposed :: Num a => [[a]]
binomialLine :: Integral a => a -> [a]
binomialColumn :: Integral a => a -> [a]

+ #152 to compute individual binomial coefficients efficiently.

2. I am not aware of any reasonable definition for Integral (Mod m).