Faster primitive root

[b-mehta]

(As mentioned in #119)

Unfinished so far, but coprimality tests for Int and Word fail seemingly due to overflow: Is it okay to just remove these tests? I think the logic in the new isPrimitiveRoot' is correct, but sanity checks would be appreciated also; the idea is that to check if g is a primitive root mod p^2 it suffices to check that g is a primitive root mod p and that g^(p-1) mod p^2 is not 1. In addition, these same conditions are exactly the conditions for g to be a primitive root mod p^k for any k. I've kept the method for testing primitive roots mod p as wikipedia describes. Similarly to check primitive roots mod 2p^k it suffices to check g is a primitive root mod p^k and that g mod 2p^k is odd, so just that g is odd.

I'd also like to add a test counting the number of primitive roots and possibly laziness tests. Are there any other appropriate tests?

[Bodigrim]

AFAIU neither tests nor arbitrary instances were changed, so they run through the same range of values as before. And if they succeeded before and fail now, then something went wrong.

A wild guess is to add

{-# SPECIALIZE isPrimitiveRoot' :: CyclicGroup Int -> Int -> Bool #-}
{-# SPECIALIZE isPrimitiveRoot' :: CyclicGroup Word -> Word -> Bool #-}

[b-mehta]

My suspicion is that overflow occurred in both the implementation and the tests - it's notable that the Integer and Natural versions pass. In addition, the test cases that fail always seem to be those which have overflow, and the new implementation has no chance of internal overflow since prefValue is never used. 

[b-mehta]

Confirming my suspicion, the test which fails (isPrimitiveRoot'Property1) tests gcd n (prefValue (cyclicGroupToModulo cg)) == 1 and the function isPrimitiveRoot' asserts gcd n (prefValue m) == 1 where m = cyclicGroupToModulo cg (line 121). The overflow occurs in calculating prefValue m, since this is a prime to a (probably) large power, but the new implementation does not require the large power to be calculated.

[Bodigrim]

Yes, you are right: previously both implementation and test were incorrect. Thanks for noticing this. I propose to amend the test, avoiding overflow in arguments of gcd. I often use Int/Word instead of Integer/Natural for the sake of performance, so it is worth to check that primitive root is correct.

import qualified Math.NumberTheory.GCD as GCD

isPrimitiveRoot'Property1
  :: (Eq a, Integral a, UniqueFactorisation a)
  => AnySign a -> CyclicGroup a -> Bool
isPrimitiveRoot'Property1 (AnySign n) cg
  = gcd (toInteger n) (prefValue (castPrefactored (cyclicGroupToModulo cg))) == 1
  || not (isPrimitiveRoot' cg n)

castPrefactored :: Integral a => Prefactored a -> Prefactored Integer
castPrefactored = fromFactors . GCD.splitIntoCoprimes . map (first toInteger) . GCD.toList . prefFactors

castPrefactored is very ugly, of course.

---

I'd also like to add a test counting the number of primitive roots and possibly laziness tests.

Sounds good. It is also worth to add a benchmark, showing the achieved performance improvement.

[b-mehta]

Should be all completed now.

Benchmarking: I've added benchmarks for cyclicGroupFromModulo and isPrimitiveRoot', as these are the two time-intensive functions in this module. I haven't benchmarked isPrimitiveRoot since it is essentially just a combination of the above two functions.

The new version of isPrimitiveRoot' speeds up primitive root checking for primes to a large power by a huge amount: checking a primitive root modulo 3^20000 went from around 12s to around 1.4μs and similarly modulo 2*3^20000 went from around 11s to around 1.5μs. For large primes, there was essentially no change, and there was a small change for large prime squares: from around 1.3ms to around 850μs.

[Bodigrim]

Let us rather derive NFData instance for CyclicGroup, similar to benchmarks of Gaussian integers: one and two.

[b-mehta]

Fixed.
EDIT: maybe not. The problem is that this requires NFData PrimeNat (or something analogous). This can't be done as a standalone/orphan instance since PrimeNat is in Math.NumberTheory.Primes.Types, which is a hidden module. I'm hesitant to include the instance in the hidden module since this involves including deepseq in the library dependencies.

[Bodigrim]

They are not dubious anymore, right?

[b-mehta]

Fixed.

[Bodigrim]

If you need to promote m to type level, there is a shorter approach:

case someNatVal m of
  SomeNat (_ :: Proxy t) -> ...

[b-mehta]

Fixed.

[Bodigrim]

Let us use Test.QuickCheck.Small newtype here. Positive does not mean much for Natural and I am worried that if the Arbitrary instance of Natural gets changed, this test may take a really long time.

[b-mehta]

I'm not entirely sure what your proposed fix here is, could you elaborate? Positive Natural seems to be used fairly often across the test suites here, including in isPrimitiveRootProperty2.

[Bodigrim]

My proposal was to write

isPrimitiveRootProperty5 :: Small Natural -> Bool
isPrimitiveRootProperty5 (Smalll m) = ...

But now I realised that there is no Small newtype in smallcheck and my approach actually would not work. Nevermind.

[Bodigrim]

checking a primitive root modulo 3^20000 went from around 12s to around 1.4μs and similarly modulo 2*3^20000 went from around 11s to around 1.5μs.

Wow! This is impressive.

[Bodigrim]

Overall looks good to me, but Travis build fails with

benchmark/Math/NumberTheory/PrimitiveRootsBench.hs:15:1: error:
    Could not load module ‘Math.NumberTheory.Primes.Types’
    it is a hidden module in the package ‘arithmoi-0.7.0.0’
    Use -v to see a list of the files searched for.
   |
15 | import Math.NumberTheory.Primes.Types
   | ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
cabal: Failed to build bench:criterion from arithmoi-0.7.0.0.

I do not want this module to be publicly exposed.
Is it possible to write instance NFData (Prime Integer)? If no, then let us move NDFata instances from benchmarks to the library, under relevant types.

[b-mehta]

Everything should be working now!

[Bodigrim]

Thanks for your contribution! Merged.