Factorisation including units

[b-mehta]

Is it possible to have something like

factorise' :: UniqueFactorisation a => a -> ([(Prime a, Word)], a)

for which the second argument is the remaining unit, with the invariant

unit * product (map (\(p,k) -> unPrime p ^ k) powers) == n
  where (powers, unit) = factorise' n

and snd (factorise' n) is always a unit?

The invariant for factorise is significantly weaker, since 3+4i and 5 have the same absolute value in the Gaussian integers, but ought to have different factorisations, so just checking that

abs n == abs (product (map (\(p, k) -> unPrime p ^ k) (factorise n)))

seems disappointing for domains with possibly non-real parts.

In addition, it would be nice to have that factorise' has a right inverse, and one could reconstruct a number from its factorisation.

[Bodigrim]

It is true that norm 5 == norm (3 + 4i), but abs 5 /= norm (3 + 4i).

[b-mehta]

Ah fair point, my mistake. I still feel like this may be a valuable function in any case.

[Bodigrim]

Do you have any particular use case in mind? In my personal experience such function was rarely needed.

[rockbmb]

If the unit is necessary, can't one just do

let (powers, unit) = (factorise n, signum n)
in ...

since factorise = factorise . abs?

[b-mehta]

That approach sounds reasonable - seems like I didn't think this through well enough. Thanks both.

[Bodigrim]

Even if n == abs n, it is not guaranteed that product (map (\(p, k) -> unPrime p ^ k) (factorise n) == n. For example, take n = 2 :: GaussianInteger. We have factorise n = [(Prime 1+ι,2)], but (1+ι)^2 = 2*ι /= 2.

[rockbmb]

@b-mehta my idea's not entirely correct, but I found the correct way to express what I meant. Let g be some GaussianInteger.

let g' = product . map (\(p, k) -> unPrime p ^ k) . factorise $ g
    sig' = signum g'
    sig = signum g
    -- This is the sign you want
    s = div sig sig'

Now we have that

s * g' == g

You can try running it with the test-suite by adding this to GaussianIntegersTests:

factoriseProperty1 :: GaussianInteger -> Bool
factoriseProperty1 g = g == 0 || s * g' == g
  where
    sig = signum g
    sig' = signum g'
    s = div sig sig'
    factors = factorise g
    g' = product $ map (\(p, k) -> unPrime p ^ k) factors

[b-mehta]

@rockbmb Right - my initial concern had been that this approach requires actually computing the product again which could be more costly than necessary if factorise throws away the unit information anyway.

@Bodigrim I don't have a use case in mind, it just struck me as odd that the units were thrown away, given that the statement of unique factorisation I'm used to gives something like n = u p1^k1 ... pi^ki.

[Bodigrim]

For Int/Word/Integer/Natural one can just take signum of an argument, as pointed above by @rockbmb. It is true that things get more complicated for quadratic fields, but here tracking of units is not actually free of charge, it makes factorisation algorithm more involved. I believe that for most of applications one can spend time reconstructing the product back.

[rockbmb]

@b-mehta if the cost of multiplicating the prime factors is what you're worried about, you can instead map (signum . fst) over factorise n and then multiply the units produced. Much cheaper, and the unit you get ~is the same as what my previous code would have given~ can then be divided by the original one as I described above.