Prefactored binomial coefficients

[Bodigrim]

There are applications, when it is desirable to represent binomial coefficients binomial via their factorisation. Binomials easily get huge, so it is vastly inefficient to compute them as numbers first and factorise afterwards.

Define

factorialFactors :: (UniqueFactorisation a, Integral a) => a -> [(Prime a, Word)]
binomialFactors :: (UniqueFactorisation a, Integral a) => a -> a -> [(Prime a, Word)]

Earlier approach, rejected now

There is a Prefactored data type, designed to keep information about factors under the hood by using a clever Num instance. For example,

factorial :: Integer -> Prefactored Integer 
factorial n = product $ map fromValue [1..10]

> factorial 10
Prefactored {prefValue = 3628800, prefFactors = Coprimes {unCoprimes = [(7,1),(5,2),(3,4),(2,8)]}}

Unfortunately, computation of binomial coefficients involves not only multiplication, but also division (without remainder). And Prefactored is not an instance of Integral (or Euclidean).

Is is possible to write a clever Integral instance for Prefactored, such that could preserve information about factors, when remainder appears to be zero?

[gilgamec]

Since Coprimes handles negative exponents exactly as one would hope, it seems that it'd be pretty simple. I'd point out that we can create a RationalPrefactored, allowing both positive and negative exponents, which is a full member of Fractional:

import Data.Ratio

import Math.NumberTheory.Euclidean
import Math.NumberTheory.Euclidean.Coprimes

mkCoprimes :: (Euclidean a, Eq b, Num b) => [(a,b)] -> Coprimes a b
mkCoprimes = foldr (uncurry insert) mempty

data RationalPrefactored a = RationalPrefactored
  { prefNumerator, prefDenominator :: a
  , prefFactors :: Coprimes a Int
  } deriving Show

fromValue :: (Num a, Eq a) => a -> RationalPrefactored a
fromValue a = RationalPrefactored a 1 (singleton a 1)

fromFactors :: Num a => Coprimes a Int -> RationalPrefactored a
fromFactors as = RationalPrefactored
  { prefNumerator   = product [ k^e | (k,e) <- unCoprimes as, e > 0 ]
  , prefDenominator = product [ k^(-e) | (k,e) <- unCoprimes as, e < 0 ]
  , prefFactors = as }

denProduct :: Euclidean a => Coprimes a Int -> Coprimes a Int -> Coprimes a Int
denProduct fs1 fs2 = mkCoprimes [ (k,-e) | (k,e) <- unCoprimes fs1, e < 0 ]
                  <> mkCoprimes [ (k,-e) | (k,e) <- unCoprimes fs2, e < 0 ]

instance Euclidean a => Num (RationalPrefactored a) where
  fromInteger = fromValue . fromInteger
  negate (RationalPrefactored n d f) = RationalPrefactored (negate n) d f
  abs (RationalPrefactored n d f) = RationalPrefactored (abs n) (abs d) f
  signum _ = undefined
  (RationalPrefactored n1 d1 fs1) + (RationalPrefactored n2 d2 fs2)
    = fromValue (n1*d2 + n2*d1) / fromFactors (denProduct fs1 fs2)
  (RationalPrefactored n1 d1 fs1) - (RationalPrefactored n2 d2 fs2)
    = fromValue (n1*d2 - n2*d1) / fromFactors (denProduct fs1 fs2)
  (RationalPrefactored _ _ fs1) * (RationalPrefactored _ _ fs2)
    = fromFactors (fs1 <> fs2)

instance Euclidean a => Fractional (RationalPrefactored a) where
  fromRational r = (fromValue . fromInteger $ numerator r) /
                   (fromValue . fromInteger $ denominator r)
  recip pref = RationalPrefactored
    { prefNumerator = prefDenominator pref
    , prefDenominator = prefNumerator pref
    , prefFactors = mkCoprimes [ (k,-e) | (k,e) <- unCoprimes (prefFactors pref) ]
    }

factorial :: Integer -> RationalPrefactored Integer
factorial n = product $ map fromValue [1..n]

choose :: Integer -> Integer -> RationalPrefactored Integer
choose m n = factorial m / (factorial n * factorial (m - n))

It could undoubtedly be made more efficient (in particular, by having insert remove factors with zero exponent), but the gist is there.

[Bodigrim]

Yep, RationalPrefactored looks cool.

I have a similar draft of Prefactored representation of binomial, but unfortunately it is too slow. Basically, Semigroup Coprimes is not terribly fast (for a good reason) and (<>) may take roughly up to O(n^2) time. This could be improved by employing the fact that prefFactors of factorial consists of primes only. It would be interesting to experiment with

data Prefactored a = Prefactored 
  { prefValue :: a 
  , prefPrimeFactors :: Map (Prime a) Word 
  , prefPossiblyCompositeFactors :: Coprimes a Word
  }

[Bodigrim]

After some deliberarion, now I think that the most profitable approach is just to define

factorialFactors :: UniqueFactorisation a => Int -> [(Prime a, Word)]
binomialFactors :: UniqueFactorisation a => Int -> Int -> [(Prime a, Word)]