atrophy 0.1.0.0 → 0.2.0.0

raw patch · 12 files changed

+1715/−522 lines, 12 filesdep +primitivedep −HUnitdep −contiguousdep −quickcheck-classesdep ~QuickCheckdep ~basedep ~deepseqPVP ok

version bump matches the API change (PVP)

Dependencies added: primitive

Dependencies removed: HUnit, contiguous, quickcheck-classes

Dependency ranges changed: QuickCheck, base, deepseq, random, tasty-bench, wide-word

API changes (from Hackage documentation)

Files

• CHANGELOG.md +41/−1
• README.md +110/−0
• atrophy.cabal +55/−35
• bench/Main.hs +159/−71
• src/Atrophy.hs +40/−9
• src/Atrophy/Internal.hs +424/−141
• src/Atrophy/Internal/LongDivision.hs +0/−84
• src/Atrophy/Internal/Prim.hs +238/−0
• src/Atrophy/Known.hs +261/−0
• src/Atrophy/LongDivision.hs +116/−37
• src/Atrophy/LongMultiplication.hs +68/−82
• tests/Main.hs +203/−62

@@ -1,5 +1,45 @@ # Revision history for atrophy -## 0.1.0.0 -- YYYY-mm-dd+## 0.2.0.0 -- 2026-09-15++A rewrite for performance. The API has changed.++* Division is now built on GHC's double-word primops (`timesWord2#`,+  `quotRemWord2#`) instead of boxed `Word128` arithmetic, and never allocates.+* `Word64` uses Granlund-Montgomery: one multiplication, no branches, no+  special cases for 1 or powers of two.+* `Word32`, `Word16` and `Word8` use a 64-bit reciprocal (Lemire et al.): one+  multiplication.+* `Word128` strength reduction is now implemented. `new` for `Word128` uses+  hardware 128/64 divisions and a 3-by-2 normalized division rather than+  `Integer`.+* New `StrengthReduce` class with an injective `StrengthReduced` type family+  replaces the `HasField`-based functions. `new`, `divRem`, `div'` and `rem'`+  work at every width; `new64`, `div64`, `rem64` and `divRem64` are gone.+  The `StrengthReducedW*` types are now abstract.+* Numerators known at compile time: `divRemConst` and friends (value level, for+  literals) and `divRemN @n` (type level). Zero, one, powers of two and+  `maxBound` avoid the multiplication entirely; `Word128` numerators below+  2^64 need two multiplications instead of eight.+* `divRemNonZero` and friends: hardware division without the zero check.+  `divRemNonZeroN @n` uses a 32-bit division when the numerator fits.+* Divisors known at compile time: `divRemK @d`, `divK @d`, `remK @d`. GHC's+  native code generator does not strength-reduce division by constants; these+  compute libdivide's magic numbers during type checking.+* `Atrophy.LongDivision` now uses Möller-Granlund 2-by-1 division with a+  precomputed reciprocal (`Divisor2By1`), takes little-endian limbs, returns the+  remainder, and works for every divisor. The previous implementation processed+  limbs most significant first and was incorrect for divisors of 2^32 or more.+  `divide128By64Preshifted` and `divide128MaxBy64` have been removed.+* `Atrophy.LongMultiplication` works on `primitive`'s `PrimArray` and+  `MutablePrimArray` directly. `multiply256By128UpperBits` is a straight chain of+  64-bit multiplications.+* Dropped the `contiguous` dependency. Now depends on `primitive` and+  `deepseq`; `NFData` instances are provided.+* `NonZero` gains `Eq`, `Ord` and `NFData` instances, `nonZero` and+  `getNonZero`.+* Requires GHC 9.4 or later.++## 0.1.0.0 -- 2022  * First version. Released on an unsuspecting world.

@@ -0,0 +1,110 @@+# Atrophy+Fast div/mod via arithmetic strength reduction.++Precompute a divisor once, then divide by it with a multiplication and a few+shifts instead of a hardware division.++```haskell+import Atrophy++let d = new (NonZero 7) :: StrengthReduced Word64+div' 100 d  -- 14+rem' 100 d  -- 2+```++Works for `Word8`, `Word16`, `Word32`, `Word64` and `Word128`, all through the+`StrengthReduce` class. Everything is `INLINE`, branchless where it matters, and+never allocates.++## Constants++GHC's native code generator does **not** strength-reduce division by constants:+``x `quot` 7`` compiles to a `div` instruction. `Atrophy.Known` does it for you,+computing the magic numbers during type checking:++```haskell+divK @7 x       -- one multiplication, no division+remK @1000 x+```++It also handles numerators known at compile time, with a runtime divisor:++```haskell+divN @(2 ^ 63) d          -- d :: StrengthReduced Word64; no multiplication at all+divNonZeroN @1000000 d    -- d :: NonZero Word64; a 32-bit hardware division+```++Zero, one, powers of two and `maxBound` skip the multiplication entirely, and+`Word128` numerators below 2^64 need two 64-bit multiplications instead of eight.+`divConst` and friends do the same for literal numerators passed as values; GHC+folds the branches away.++## Algorithms++* `Word64`: Granlund & Montgomery, "Division by Invariant Integers using+  Multiplication". One `mul`, no branches, no special cases for 1 or powers of+  two.+* `Word32` and smaller: Lemire, Kaser & Kurz, "Faster Remainder by Direct+  Computation". One `mul`.+* `Word128`: Granlund & Montgomery on 64-bit limbs. `new` uses hardware 128/64+  divisions and a normalized 3-by-2 division.+* Compile-time divisors: libdivide's unsigned algorithm, with the choice between+  shift, multiply-shift and multiply-add-shift made at compile time.+* `Atrophy.LongDivision`: Möller & Granlund, "Improved division by invariant+  integers", for dividing little-endian multi-limb numbers by a 64-bit divisor.++## Benchmarks++Nanoseconds per operation, averaged over 10000 uniformly random dividends.+Divisors have a uniformly random bit length. GHC 9.14.1, native code generator,+AMD Ryzen 7 7840U. "atrophy 0.1" is the previous release, run on the same data.++| `Word64`                           | GHC `quot` | atrophy 0.1 | atrophy |+|------------------------------------|-----------:|------------:|--------:|+| `new`                              |            |        5.56 |    2.92 |+| one divisor, many dividends        |       1.52 |        17.6 |    0.92 |+| one divisor, remainder             |            |             |    1.13 |+| unique divisors (`new` + `div'`)   |       1.61 |        33.5 |    2.54 |+| constant divisor 7 (`divK`)        |       1.91 |             |    0.85 |+| constant divisor 10^9+7 (`divK`)   |       1.50 |             |    0.88 |+| constant numerator 10^6, `NonZero` divisor (`divNonZeroN`) | 1.50 | | 0.52 |+| constant numerator 2^63, `NonZero` divisor (`divNonZeroN`) | 3.05 | | 1.61 |++| `Word32`                           | GHC `quot` | atrophy 0.1 | atrophy |+|------------------------------------|-----------:|------------:|--------:|+| `new`                              |            |        1.56 |    1.88 |+| one divisor, many dividends        |       1.29 |        1.05 |    0.70 |+| unique divisors (`new` + `div'`)   |       1.29 |        1.61 |    1.72 |+| constant divisor 7 (`divK`)        |       1.29 |             |    0.69 |++| `Word128`                          | wide-word `quot` | atrophy |+|------------------------------------|-----------------:|--------:|+| `new`                              |                  |    10.8 |+| one divisor, many dividends        |              121 |    4.89 |+| unchecked hardware division (`divNonZero`) |          |    4.12 |+| unique divisors (`new` + `div'`)   |             39.2 |    16.3 |+| constant divisor 10^19 (`divK`)    |             4.68 |    3.91 |+| constant numerator 10^18, `NonZero` divisor (`divNonZeroN`) | 5.13 | 0.92 |++Constant numerators with a precomputed divisor, cycling through 64 divisors:++|                                    | `div'` | `divN` |+|------------------------------------|-------:|-------:|+| `Word64`, numerator 2^63           |   2.26 |   1.78 |+| `Word64`, numerator 2^64 - 1       |   2.23 |   1.78 |+| `Word128`, numerator 10^18         |   5.40 |   2.80 |++Dividing a 64-limb number by a 64-bit divisor with `longDivision` takes 281 ns,+including allocating the quotient; GMP's hand-written assembly, via `Integer`,+takes 166 ns.++Performance is *heavily* platform dependent. Zen 4 has an unusually fast+hardware divider, so these numbers understate the gains on most other CPUs,+where a 64-bit `div` costs 35 to 90 cycles rather than 10 to 20. On this machine+a hardware division is still faster than `new` followed by a single `div'`, so+strength reduction pays off when a divisor is reused.++Run them yourself with `cabal bench`.++## Special thanks+Originally based on https://github.com/ejmahler/strength_reduce

@@ -1,68 +1,88 @@-cabal-version: 2.4+cabal-version: 3.0 name: atrophy-version: 0.1.0.0-synopsis: Faster integer division and modulus operations +version: 0.2.0.0+synopsis: Faster integer division and modulus operations description:   Fast div/mod via arithmetic strength reduction.--  Good compilers already perform this optimization for divisors that are known at compile time; this library enables this optimization for divisors that are only known at runtime.+  .+  Precompute a divisor that is only known at runtime, and divide by it with a+  multiplication and a few shifts instead of a hardware division. Also provides+  strength reduction for divisors and numerators known at compile time, which+  GHC does not do on its own, and multi-limb division and multiplication. license: MIT license-file: LICENSE author: Zachary Churchill <zacharyachurchill@gmail.com> maintainer: Zachary Churchill <zacharyachurchill@gmail.com> copyright: 2022 Zachary Churchill category: Math-extra-source-files: CHANGELOG.md+build-type: Simple+extra-doc-files:+  CHANGELOG.md+  README.md+tested-with: GHC == 9.4.8 || == 9.10.3 || == 9.14.1  source-repository head   type: git   location: https://github.com/goolord/atrophy/ +common warnings+  ghc-options:+    -Wall+    -Wcompat+    -Widentities+    -Wincomplete-record-updates+    -Wincomplete-uni-patterns+    -Wredundant-constraints+    -Wpartial-fields+ library+  import: warnings   exposed-modules:-      Atrophy-    , Atrophy.LongDivision-    , Atrophy.LongMultiplication+    Atrophy+    Atrophy.Known+    Atrophy.LongDivision+    Atrophy.LongMultiplication   other-modules:-      Atrophy.Internal-      Atrophy.Internal.LongDivision-  -- other-extensions:-  ghc-options: -Wall -O2 -Wredundant-constraints+    Atrophy.Internal+    Atrophy.Internal.Prim   build-depends:-      base >=4.10.0.0 && <5-    , wide-word-    , contiguous >= 0.6.0.0+    , base >=4.17 && <5+    , deepseq >=1.4 && <1.7+    , primitive >=0.7 && <0.10+    , wide-word >=0.1 && <0.2   hs-source-dirs: src-  default-language: Haskell2010+  default-language: GHC2021+  ghc-options: -O2  test-suite atrophy-test-  default-language: Haskell2010+  import: warnings   type: exitcode-stdio-1.0   hs-source-dirs: tests   main-is: Main.hs+  default-language: GHC2021+  ghc-options: -O2 -threaded -rtsopts -with-rtsopts=-N   build-depends:-      base >=4.10.0.0 && <5-    , wide-word-    , contiguous >= 0.6.0.0-    , HUnit-    , tasty-    , QuickCheck-    , quickcheck-classes+    , atrophy+    , base+    , primitive+    , QuickCheck >=2.14     , tasty     , tasty-hunit     , tasty-quickcheck-    , atrophy+    , wide-word  benchmark atrophy-bench+  import: warnings   type: exitcode-stdio-1.0-  build-depends:-      base >=4.10.0.0 && <5-    , tasty-    , tasty-bench-    , deepseq-    , atrophy-    , random-  ghc-options: -Wall -O2-  default-language: Haskell2010   hs-source-dirs: bench   main-is: Main.hs+  default-language: GHC2021+  ghc-options: -O2 -rtsopts+  build-depends:+    , atrophy+    , base+    , deepseq+    , primitive+    , random >=1.2+    , tasty-bench >=0.3+    , wide-word

@@ -1,87 +1,175 @@-{-# LANGUAGE-    BangPatterns-  , GADTs-  , DeriveGeneric-  , StandaloneDeriving-  , MagicHash-  , DataKinds-  , GeneralizedNewtypeDeriving-  , TypeApplications-  , ScopedTypeVariables-  , NumericUnderscores-#-}--{-# OPTIONS_GHC-    -fno-warn-orphans-#-}-{-# LANGUAGE DerivingStrategies #-}+{-# LANGUAGE AllowAmbiguousTypes #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE ExplicitNamespaces #-} -module Main where+module Main (main) where -import Test.Tasty.Bench (bench, bgroup, defaultMain, nf, Benchmark, envWithCleanup, nfIO)-import Control.DeepSeq (NFData, force) import Atrophy-import GHC.Generics-import Data.Word (Word64, Word32)-import Control.Exception (evaluate)-import Test.Tasty (withResource)-import System.Mem (performMajorGC)+import Control.DeepSeq (NFData (..))+import Control.Monad (replicateM)+import Control.Monad.ST (runST)+import Data.Bits+import Data.Primitive.PrimArray+import Data.Primitive.SmallArray+import Data.WideWord.Word128+import Data.Primitive.Types (Prim)+import Data.Word+import GHC.TypeNats (type (-), type (^)) import System.Random.Stateful-import Data.Proxy (Proxy (Proxy))+import Test.Tasty.Bench -deriving instance Generic StrengthReducedW64-instance NFData StrengthReducedW64+size :: Int+size = 10000 -deriving instance Generic StrengthReducedW32-instance NFData StrengthReducedW32+-- | Dividends, divisors (never zero) and one fixed divisor.+data Env a = Env !(PrimArray a) !(PrimArray a) !a -instance (Bounded a, Num a, UniformRange a) => Uniform (NonZero a) where-  uniformM g = NonZero <$> uniformRM (1, maxBound) g+instance NFData (Env a) where+  rnf !_ = () -deriving newtype instance NFData a => NFData (NonZero a)+-- | Precomputed divisors, for the known-numerator benchmarks.+newtype Reduced a = Reduced (SmallArray (StrengthReduced a)) -manyRandom :: forall a. (Uniform a, NFData a) => IO [a]-manyRandom = uniformListM 10_000 globalStdGen+instance NFData (Reduced a) where+  rnf (Reduced !_) = () -randomEnv :: NFData b => IO b-  -> (b -> Benchmark)-  -> Benchmark-randomEnv a = envWithCleanup (a >>= evaluate . force) (const performMajorGC)+class (Prim a, Integral a, FiniteBits a, Bounded a, StrengthReduce a) => BenchWord a where+  uniformW :: StdGen -> (a, StdGen) -randomEnv' :: NFData b => IO b-  -> (IO b -> Benchmark)-  -> Benchmark-randomEnv' a = withResource (a >>= evaluate . force) (const performMajorGC)+instance BenchWord Word32 where uniformW = uniform+instance BenchWord Word64 where uniformW = uniform+instance BenchWord Word128 where+  uniformW g0 = let (h, g1) = uniform g0; (l, g2) = uniform g1 in (Word128 h l, g2) -main :: IO ()-main = do-  defaultMain $-    [ bgroup "atrophy"-        [ bgroup "Word64" $ atrophyBench div64 new64-        , bgroup "Word32" $ atrophyBench div' (new StrengthReducedW32)-        ]-    , bgroup "ghc"-        [ bgroup "Word64" $ ghcBench (Proxy @Word64)-        , bgroup "Word32" $ ghcBench (Proxy @Word32)-        ]-    ]+-- | Dividends are uniform. Divisors have a uniform bit length, then are uniform+-- within it: this exercises both small and large quotients, which matters for+-- hardware division latency.+mkEnv :: forall a. BenchWord a => IO (Env a)+mkEnv = do+  let g0 = mkStdGen 2026+      (ns, g1) = go size g0 []+      (ds, g2) = goD size g1 []+      (d, _) = goD 1 g2 []+  pure $ Env (primArrayFromList ns) (primArrayFromList ds) (sum d)+  where+  go 0 g acc = (acc, g)+  go k g acc = let (x, g') = uniformW g in go (k - 1 :: Int) g' (x : acc)+  goD 0 g acc = (acc, g)+  goD k g acc =+    let (bitsLen, g') = uniformR (1, finiteBitSize (0 :: a)) g+        (x, g'') = uniformW g'+    in goD (k - 1 :: Int) g'' (max 1 (x `shiftR` (finiteBitSize (0 :: a) - bitsLen)) : acc) -ghcBench :: forall a. (NFData a, Uniform a, Integral a, Random a, Bounded a) => Proxy a -> [Benchmark]-ghcBench _ =-  [ randomEnv (manyRandom @(a, a)) $ \somePairs ->-      bench "div 10000 uniques" $ nf (\xs -> fmap (\(x, y) -> x `div` y) xs) somePairs-  , randomEnv' ((,) <$> randomRIO (1, maxBound) <*> manyRandom @a) $ \x -> bench "div 10000, 1 divisor" $ nfIO $ do -      (divisor', dividends) <- x-      pure $ fmap (\dividend' -> dividend' `div` divisor') dividends+-- | A small table, so that the benchmark measures division rather than cache+-- misses on 10000 boxed values.+tableSize :: Int+tableSize = 64++mkReduced :: BenchWord a => Env a -> Reduced a+mkReduced (Env _ ds _) = Reduced $ runST $ do+  arr <- newSmallArray tableSize undefined+  let go i | i == tableSize = pure ()+           | otherwise = (writeSmallArray arr i $! new (NonZero (indexPrimArray ds i))) >> go (i + 1)+  go 0+  unsafeFreezeSmallArray arr++{-# INLINE sumMap #-}+sumMap :: (Prim a, Num a) => (a -> a) -> PrimArray a -> a+sumMap f = foldlPrimArray' (\acc x -> acc + f x) 0++{-# INLINE sumIx #-}+sumIx :: Num a => (Int -> a) -> a+sumIx f = go 0 0+  where+  go !i !acc+    | i == size = acc+    | otherwise = go (i + 1) (acc + f i)++{-# INLINE sumReduced #-}+sumReduced :: Num a => (StrengthReduced a -> a) -> Reduced a -> a+sumReduced f (Reduced arr) = sumIx (\i -> f (indexSmallArray arr (i .&. (tableSize - 1))))++{-# NOINLINE runtime #-}+runtime :: a -> a+runtime x = x++main :: IO ()+main = defaultMain+  [ bgroup "Word64" $ runtimeBenches @Word64 +++      [ env (mkEnv @Word64) $ \ ~(Env ns _ _) -> bgroup "known divisor 7"+          [ bench "ghc quot" $ whnf (sumMap (`quot` 7)) ns+          , bench "atrophy divK" $ whnf (sumMap (divK @7)) ns+          , bench "atrophy div' (runtime)" $ whnf (\sr -> sumMap (`div'` sr) ns) (new (NonZero (runtime 7)))+          ]+      , env (mkEnv @Word64) $ \ ~(Env ns _ _) -> bgroup "known divisor 1000000007"+          [ bench "ghc quot" $ whnf (sumMap (`quot` 1000000007)) ns+          , bench "atrophy divK" $ whnf (sumMap (divK @1000000007)) ns+          ]+      , env (mkEnv @Word64) $ \ ~(e@(Env _ ds _)) -> env (pure (mkReduced e)) $ \rs -> bgroup "known numerator 2^63"+          [ bench "ghc quot" $ whnf (sumMap (quot (2 ^ (63 :: Int)))) ds+          , bench "atrophy divNonZeroN" $ whnf (sumMap (divNonZeroN @(2 ^ 63) . NonZero)) ds+          , bench "atrophy div' (runtime)" $ whnf (\n -> sumReduced (div' n) rs) (runtime (2 ^ (63 :: Int)))+          , bench "atrophy divN" $ whnf (sumReduced (divN @(2 ^ 63))) rs+          ]+      , env (mkEnv @Word64) $ \ ~(e@(Env _ ds _)) -> env (pure (mkReduced e)) $ \rs -> bgroup "known numerator 1000000"+          [ bench "ghc quot" $ whnf (sumMap (quot 1000000)) ds+          , bench "atrophy divNonZeroN" $ whnf (sumMap (divNonZeroN @1000000 . NonZero)) ds+          , bench "atrophy div' (runtime)" $ whnf (\n -> sumReduced (div' n) rs) (runtime 1000000)+          , bench "atrophy divN" $ whnf (sumReduced (divN @1000000)) rs+          ]+      , env (mkEnv @Word64) $ \ ~(e@(Env _ ds _)) -> env (pure (mkReduced e)) $ \rs -> bgroup "known numerator maxBound"+          [ bench "ghc quot" $ whnf (sumMap (quot maxBound)) ds+          , bench "atrophy div' (runtime)" $ whnf (\n -> sumReduced (div' n) rs) (runtime maxBound)+          , bench "atrophy divN" $ whnf (sumReduced (divN @(2 ^ 64 - 1))) rs+          ]+      , env mkLimbs $ \ ~(ls, d) -> bgroup "long division 64 limbs"+          [ bench "Integer quotRem" $ whnf (\n -> n `quotRem` fromIntegral d) (limbsToInteger ls)+          , bench "atrophy longDivision" $ whnf (\dv -> runST $ do+              q <- newPrimArray 64+              longDivision dv ls q) (newDivisor2By1 (NonZero d))+          ]+      ]+  , bgroup "Word32" $ runtimeBenches @Word32 +++      [ env (mkEnv @Word32) $ \ ~(Env ns _ _) -> bgroup "known divisor 7"+          [ bench "ghc quot" $ whnf (sumMap (`quot` 7)) ns+          , bench "atrophy divK" $ whnf (sumMap (divK @7)) ns+          ]+      ]+  , bgroup "Word128" $ runtimeBenches @Word128 +++      [ env (mkEnv @Word128) $ \ ~(Env ns _ _) -> bgroup "known divisor 10^19"+          [ bench "wide-word quot" $ whnf (sumMap (`quot` 10000000000000000000)) ns+          , bench "atrophy divK" $ whnf (sumMap (divK @10000000000000000000)) ns+          ]+      , env (mkEnv @Word128) $ \ ~(e@(Env _ ds _)) -> env (pure (mkReduced e)) $ \rs -> bgroup "known numerator 10^18"+          [ bench "wide-word quot" $ whnf (sumMap (quot 1000000000000000000)) ds+          , bench "atrophy divNonZeroN" $ whnf (sumMap (divNonZeroN @1000000000000000000 . NonZero)) ds+          , bench "atrophy div' (runtime)" $ whnf (\n -> sumReduced (div' n) rs) (runtime 1000000000000000000)+          , bench "atrophy divN" $ whnf (sumReduced (divN @1000000000000000000)) rs+          ]+      ]   ] -atrophyBench :: forall base sr out. (NFData base, Uniform base, UniformRange base, NFData sr, NFData out, Random base, Num base, Bounded base) => (base -> sr -> out) -> (NonZero base -> sr) -> [Benchmark]-atrophyBench divF newF =-  [ randomEnv (uniformM globalStdGen) $ \divisor' ->-      bench "new" $ nf newF divisor'-  , randomEnv (manyRandom @(base, base)) $ \somePairs ->-      bench "div 10000 uniques" $ nf (fmap (\(x, y) -> x `divF` newF (NonZero y))) somePairs-  , randomEnv' ((,) <$> (newF <$> (NonZero <$> randomRIO (1, maxBound))) <*> manyRandom @base) $ \x -> bench "div 10000, 1 divisor" $ nfIO $ do -      (divisor', dividends) <- x-      pure $ fmap (\dividend' -> dividend' `divF` divisor') dividends+runtimeBenches :: forall a. BenchWord a => [Benchmark]+runtimeBenches =+  [ env (mkEnv @a) $ \ ~(Env _ _ d) -> bench "new" $ whnf new (NonZero d)+  , env (mkEnv @a) $ \ ~(Env ns _ d) -> bgroup "fixed divisor"+      [ bench "baseline quot" $ whnf (\dv -> sumMap (`quot` dv) ns) d+      , bench "atrophy divNonZero" $ whnf (\dv -> sumMap (`divNonZero` dv) ns) (NonZero d)+      , bench "atrophy div'" $ whnf (\sr -> sumMap (`div'` sr) ns) (new (NonZero d))+      , bench "atrophy rem'" $ whnf (\sr -> sumMap (`rem'` sr) ns) (new (NonZero d))+      ]+  , env (mkEnv @a) $ \ ~(Env ns ds _) -> bgroup "unique divisors"+      [ bench "baseline quot" $ whnf (\xs -> sumIx (\i -> indexPrimArray xs i `quot` indexPrimArray ds i)) ns+      , bench "atrophy div' . new" $ whnf (\xs -> sumIx (\i -> indexPrimArray xs i `div'` new (NonZero (indexPrimArray ds i)))) ns+      ]   ]++mkLimbs :: IO (PrimArray Word64, Word64)+mkLimbs = do+  g <- newIOGenM (mkStdGen 7)+  ls <- replicateM 64 (uniformM g)+  d <- uniformRM (1, maxBound) g+  pure (primArrayFromList ls, d)++limbsToInteger :: PrimArray Word64 -> Integer+limbsToInteger = foldrPrimArray (\l acc -> acc * 2 ^ (64 :: Int) + toInteger l) 0

@@ -1,12 +1,43 @@+-- | Fast division by divisors that are only known at runtime, via arithmetic+-- strength reduction.+--+-- > let d = new (NonZero 7) :: StrengthReduced Word64+-- > div' 100 d == 14+--+-- Precompute a divisor once with 'new', then each 'divRem' is a multiplication+-- and a few shifts instead of a hardware division.+--+-- For constants, see "Atrophy.Known". module Atrophy-  ( module REXPORT-  , StrengthReducedW128 (..)-  , StrengthReducedW64 (..)-  , StrengthReducedW32 (..)-  , StrengthReducedW16 (..)-  , StrengthReducedW8 (..)+  ( -- * Nonzero values+    NonZero (..)+  , nonZero+  , getNonZero+    -- * Strength reduction+  , StrengthReduce (..)+  , div'+  , rem'+    -- ** Compile-time numerators+  , divConst+  , remConst+    -- ** Unchecked hardware division+  , divNonZero+  , remNonZero+  , divNonZeroConst+  , remNonZeroConst+    -- ** Precomputed divisors+  , StrengthReducedW128+  , StrengthReducedW64+  , StrengthReducedW32+  , StrengthReducedW16+  , StrengthReducedW8+    -- * Re-exports+  , module Atrophy.Known+  , module Atrophy.LongDivision+  , module Atrophy.LongMultiplication   ) where -import Atrophy.LongDivision as REXPORT-import Atrophy.LongMultiplication as REXPORT-import Atrophy.Internal as REXPORT+import Atrophy.Internal+import Atrophy.Known+import Atrophy.LongDivision+import Atrophy.LongMultiplication

@@ -1,173 +1,456 @@-{-# LANGUAGE-    CPP-  , TypeApplications-  , DataKinds-  , FlexibleContexts-  , DuplicateRecordFields-  , TypeFamilies-  , BangPatterns-  , NumericUnderscores-  , ScopedTypeVariables-  , DerivingStrategies-  , GeneralizedNewtypeDeriving-#-}+{-# LANGUAGE MagicHash #-}+{-# LANGUAGE UnboxedTuples #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeFamilyDependencies #-}+{-# LANGUAGE DerivingStrategies #-}+{-# OPTIONS_HADDOCK not-home #-}  module Atrophy.Internal where -import Data.WideWord.Word128+import Atrophy.Internal.Prim+import Control.DeepSeq (NFData (..)) import Data.Bits-import Atrophy.Internal.LongDivision-import GHC.Records+import Data.Kind (Type)+import Data.WideWord.Word128 (Word128 (..)) import Data.Word +-- | A value that is known not to be zero. newtype NonZero a = NonZero a-  deriving newtype (Num, Show)+  deriving newtype (Eq, Ord, Show, Num, NFData)  instance (Bounded a, Num a) => Bounded (NonZero a) where-  minBound = 1+  minBound = NonZero 1   maxBound = NonZero maxBound -{-# INLINE isPowerOf2 #-}-isPowerOf2 :: (Bits a, Num a) => a -> Bool-isPowerOf2 x = (x .&. (x - 1)) == 0+{-# INLINE getNonZero #-}+getNonZero :: NonZero a -> a+getNonZero (NonZero a) = a -{-# INLINE new64 #-}-new64 :: NonZero Word64 -> StrengthReducedW64-new64 (NonZero divi) =-  if isPowerOf2 divi-  then StrengthReducedW64 0 divi-  else-    let quotient = divide128MaxBy64 $ fromIntegral divi-    in StrengthReducedW64 (quotient + 1) divi+-- | 'Nothing' for zero.+{-# INLINE nonZero #-}+nonZero :: (Eq a, Num a) => a -> Maybe (NonZero a)+nonZero 0 = Nothing+nonZero a = Just (NonZero a) -{-# INLINE divRem64 #-}-divRem64 ::-  ( HasField "divisor" strRed a-  , HasField "multiplier" strRed Word128-  , Integral a-  , FiniteBits a-  ) => a -> strRed -> (a, a)-divRem64 dividend divis =-  case getField @"multiplier" divis of-    0 ->-      let-        quotient = dividend `unsafeShiftR` (countTrailingZeros $ getField @"divisor" divis)-        remainder = dividend .&. (getField @"divisor" divis - 1)-      in (quotient, remainder)-    multiplier' ->-      let-        numerator128 = fromIntegral @_ @Word128 dividend-        multipliedHi = numerator128 * (upper128 multiplier')-        multipliedLo = upper128 (numerator128 * (lower128 multiplier'))+-- | Unsigned integers that can be divided by a precomputed divisor.+--+-- Every method is @INLINE@; calls at a known type compile down to a handful of+-- instructions with no dictionary passing and no allocation.+class StrengthReduce a where+  -- | The precomputed form of a divisor.+  type StrengthReduced a = (r :: Type) | r -> a -        quotient = fromIntegral (upper128 (multipliedHi + multipliedLo))-        remainder = dividend - quotient * getField @"divisor" divis-      in (quotient, remainder)+  -- | Precompute a divisor. This costs about as much as one hardware division.+  new :: NonZero a -> StrengthReduced a -{-# INLINE divRem #-}-{-# SPECIALIZE divRem :: Word32 -> StrengthReducedW32 -> (Word32, Word32) #-}-divRem :: forall strRed a b.-  ( HasField "divisor" strRed a-  , HasField "multiplier" strRed b-  , Integral a-  , FiniteBits a, Integral b, FiniteBits (Half b), Bits b) => a -> strRed -> (a, a)-divRem dividend divis =-  case getField @"multiplier" divis of-    0 ->-      let-        quotient = dividend `unsafeShiftR` (countTrailingZeros $ getField @"divisor" divis)-        remainder = dividend .&. (getField @"divisor" divis - 1)-      in (quotient, remainder)-    multiplier' ->-      let-        numerator64 = fromIntegral @_ @b dividend-        multipliedHi = numerator64 * (upperHalf multiplier')-        multipliedLo = upperHalf (numerator64 * (lowerHalf multiplier'))+  -- | The original divisor.+  divisor :: StrengthReduced a -> a -        quotient = fromIntegral (upperHalf (multipliedHi + multipliedLo))-        remainder = dividend - quotient * getField @"divisor" divis-      in (quotient, remainder)+  -- | Quotient and remainder, branchless.+  divRem :: a -> StrengthReduced a -> (a, a) -{-# INLINE new #-}-{-# SPECIALIZE new :: (Word64 -> Word32 -> StrengthReducedW32) -> NonZero Word32 -> StrengthReducedW32 #-}-new :: (Bits t, Integral t, Bounded (Multiplier t), Integral (Multiplier t)) =>((Multiplier t) -> t -> a) -> (NonZero t) -> a-new con (NonZero divi) =-  if isPowerOf2 divi-  then con 0 divi-  else-    let quotient = maxBound `div` fromIntegral divi-    in con (quotient + 1) divi+  -- | 'divRem', specialised for numerators that are compile-time constants.+  --+  -- The result is identical to 'divRem', but the implementation branches on+  -- properties of the numerator (zero, a power of two, all ones, fits in half+  -- the width). When the numerator is a literal, GHC folds those branches away+  -- and the multiplication often turns into a shift or disappears entirely.+  -- When it is not a literal you pay for the branches, so use 'divRem'.+  --+  -- See also 'Atrophy.Known.divRemN', which takes the numerator at the type+  -- level and so is guaranteed to be constant.+  divRemConst :: a -> StrengthReduced a -> (a, a)+  divRemConst = divRem+  {-# INLINE divRemConst #-} -{-# INLINE div64 #-}-{-# SPECIALIZE div64 :: Word64 -> StrengthReducedW64 -> Word64 #-}-div64 :: (HasField "divisor" r b, HasField "multiplier" r Word128,- Integral b, FiniteBits b) =>-  b -> r -> b-div64 a rhs = fst $ divRem64 a rhs+  -- | Hardware division without the zero check performed by 'quotRem'.+  divRemNonZero :: a -> NonZero a -> (a, a) -{-# INLINE rem64 #-}-{-# SPECIALIZE rem64 :: Word64 -> StrengthReducedW64 -> Word64 #-}-rem64 :: (HasField "divisor" r b, HasField "multiplier" r Word128,- Integral b, FiniteBits b) =>-  b -> r -> b-rem64 a rhs = snd $ divRem64 a rhs+  -- | 'divRemNonZero', specialised for numerators that are compile-time+  -- constants, e.g. using a narrower hardware division when the numerator is+  -- small. The same caveats as 'divRemConst' apply.+  divRemNonZeroConst :: a -> NonZero a -> (a, a)+  divRemNonZeroConst = divRemNonZero+  {-# INLINE divRemNonZeroConst #-} +-- | Quotient by a precomputed divisor. The remainder is never computed. {-# INLINE div' #-}-{-# SPECIALIZE div' :: Word32 -> StrengthReducedW32 -> Word32 #-}-div' ::-  ( HasField "divisor" strRed b-  , HasField "multiplier" strRed w-  , Integral b, FiniteBits b,  Integral w, FiniteBits (Half w), Bits w) => b -> strRed -> b-div' a rhs = fst $ divRem a rhs+div' :: StrengthReduce a => a -> StrengthReduced a -> a+div' n d = case divRem n d of (q, _) -> q +-- | Remainder by a precomputed divisor. {-# INLINE rem' #-}-{-# SPECIALIZE rem' :: Word32 -> StrengthReducedW32 -> Word32 #-}-rem' ::-  ( HasField "divisor" strRed b-  , HasField "multiplier" strRed w-  , Integral b, FiniteBits b,  Integral w, FiniteBits (Half w), Bits w-  ) => b -> strRed -> b-rem' a rhs = snd $ divRem a rhs+rem' :: StrengthReduce a => a -> StrengthReduced a -> a+rem' n d = case divRem n d of (_, r) -> r -{-# INLINE lower128 #-}-lower128 :: Word128 -> Word128-lower128 (Word128 _hi low) = Word128 0 low+-- | Quotient of a compile-time constant numerator. See 'divRemConst'.+{-# INLINE divConst #-}+divConst :: StrengthReduce a => a -> StrengthReduced a -> a+divConst n d = case divRemConst n d of (q, _) -> q -{-# INLINE upper128 #-}-upper128 :: Word128 -> Word128-upper128 (Word128 hi _low) = Word128 0 hi+-- | Remainder of a compile-time constant numerator. See 'divRemConst'.+{-# INLINE remConst #-}+remConst :: StrengthReduce a => a -> StrengthReduced a -> a+remConst n d = case divRemConst n d of (_, r) -> r -{-# INLINE lowerHalf #-}-lowerHalf :: forall w. ( FiniteBits (Half w), Bits w) =>w -> w-lowerHalf w = (w `unsafeShiftL` halfSize) `unsafeShiftR` halfSize-  where-  halfSize = finiteBitSize @(Half w) zeroBits+-- | Unchecked hardware quotient. See 'divRemNonZero'.+{-# INLINE divNonZero #-}+divNonZero :: StrengthReduce a => a -> NonZero a -> a+divNonZero n d = case divRemNonZero n d of (q, _) -> q -{-# INLINE upperHalf #-}-upperHalf :: forall w. ( Bits w, FiniteBits (Half w)) =>w -> w-upperHalf w = w `unsafeShiftR` halfSize-  where-  halfSize = finiteBitSize @(Half w) zeroBits+-- | Unchecked hardware remainder. See 'divRemNonZero'.+{-# INLINE remNonZero #-}+remNonZero :: StrengthReduce a => a -> NonZero a -> a+remNonZero n d = case divRemNonZero n d of (_, r) -> r -type family Multiplier a where-  Multiplier Word64 = Word128-  Multiplier Word32 = Word64-  Multiplier Word16 = Word32-  Multiplier Word8  = Word16+-- | Unchecked quotient of a compile-time constant numerator. See+-- 'divRemNonZeroConst'.+{-# INLINE divNonZeroConst #-}+divNonZeroConst :: StrengthReduce a => a -> NonZero a -> a+divNonZeroConst n d = case divRemNonZeroConst n d of (q, _) -> q -type family Half a where-  Half Word128 = Word64-  Half Word64 = Word32-  Half Word32 = Word16-  Half Word16 = Word8+-- | Unchecked remainder of a compile-time constant numerator. See+-- 'divRemNonZeroConst'.+{-# INLINE remNonZeroConst #-}+remNonZeroConst :: StrengthReduce a => a -> NonZero a -> a+remNonZeroConst n d = case divRemNonZeroConst n d of (_, r) -> r -data StrengthReducedW64 = StrengthReducedW64 { multiplier :: {-# UNPACK #-} !Word128, divisor :: {-# UNPACK #-} !Word64 }-data StrengthReducedW32 = StrengthReducedW32 { multiplier :: {-# UNPACK #-} !Word64,  divisor :: {-# UNPACK #-} !Word32 }-data StrengthReducedW16 = StrengthReducedW16 { multiplier :: {-# UNPACK #-} !Word32,  divisor :: {-# UNPACK #-} !Word16 }-data StrengthReducedW8  = StrengthReducedW7  { multiplier :: {-# UNPACK #-} !Word16,  divisor :: {-# UNPACK #-} !Word8  }+--------------------------------------------------------------------------------+-- Word64: Granlund-Montgomery, "Division by Invariant Integers using+-- Multiplication", figure 4.1. One multiplication, fully branchless, and+-- handles 1 and powers of two without special cases.+--+--   l  = ceil (log2 d)+--   m  = floor (2^64 * (2^l - d) / d) + 1+--   t  = mulhi m n+--   q  = (t + ((n - t) >> min l 1)) >> max (l - 1) 0+-------------------------------------------------------------------------------- +-- | A precomputed 'Word64' divisor. Four words.+data StrengthReducedW64 = StrengthReducedW64+  {-# UNPACK #-} !Word64 -- multiplier+  {-# UNPACK #-} !Int    -- pre-shift, 0 or 1+  {-# UNPACK #-} !Int    -- post-shift+  {-# UNPACK #-} !Word64 -- divisor+  deriving stock (Eq, Show)++instance NFData StrengthReducedW64 where+  rnf !_ = ()++instance StrengthReduce Word64 where+  type StrengthReduced Word64 = StrengthReducedW64++  {-# INLINE new #-}+  new (NonZero d) =+    let !l = 64 - countLeadingZeros (d - 1)+        !h = if l == 64 then negate d else unsafeShiftL 1 l - d+        !m = if h == 0 then 1 else case quotRem128By64 h 0 d of (# q, _ #) -> q + 1+    in StrengthReducedW64 m (min l 1) (max (l - 1) 0) d++  {-# INLINE divisor #-}+  divisor (StrengthReducedW64 _ _ _ d) = d++  {-# INLINE divRem #-}+  divRem n (StrengthReducedW64 m s1 s2 d) = gm64 n (mulHi64 m n) s1 s2 d++  {-# INLINE divRemConst #-}+  divRemConst n (StrengthReducedW64 m s1 s2 d)+    | n == 0 = (0, 0)+    | n == 1 = gm64 n 0 s1 s2 d+    -- mulhi m (2^64 - 1) == m - 1, since 1 <= m < 2^64+    | n == maxBound = gm64 n (m - 1) s1 s2 d+    -- mulhi m (2^k) == m >> (64 - k)+    | popCount n == 1 = gm64 n (unsafeShiftR m (64 - countTrailingZeros n)) s1 s2 d+    | otherwise = gm64 n (mulHi64 m n) s1 s2 d++  {-# INLINE divRemNonZero #-}+  divRemNonZero n (NonZero d) = case quotRemWord64 n d of (# q, r #) -> (q, r)++  {-# INLINE divRemNonZeroConst #-}+  divRemNonZeroConst n (NonZero d)+    | n == 0 = (0, 0)+    -- a 32-bit hardware division is considerably cheaper on many CPUs+    | n <= 0xffffffff =+        if d > n+          then (0, n)+          else+            let !q = fromIntegral (quotWord32 (fromIntegral n) (fromIntegral d))+            in (q, n - q * d)+    | otherwise = case quotRemWord64 n d of (# q, r #) -> (q, r)++-- | The part of Granlund-Montgomery after the multiplication.+{-# INLINE gm64 #-}+gm64 :: Word64 -> Word64 -> Int -> Int -> Word64 -> (Word64, Word64)+gm64 n t s1 s2 d =+  let !q = (t + unsafeShiftR (n - t) s1) `unsafeShiftR` s2+  in (q, n - q * d)++--------------------------------------------------------------------------------+-- Word32 and smaller: Lemire, Kaser & Kurz, "Faster Remainder by Direct+-- Computation". With a 64-bit reciprocal c = ceil (2^64 / d) the quotient of+-- any 32-bit n is exactly mulhi c n. The reciprocal wraps to 0 for d == 1, so+-- a mask adds n back in that case.+--------------------------------------------------------------------------------++{-# INLINE newSmall #-}+newSmall :: Word64 -> (# Word64, Word64 #)+newSmall d = case quotRemWord64 maxBound d of+  (# q, _ #) -> (# q + 1, negate (ltW d 2) #)++{-# INLINE quotSmall #-}+quotSmall :: Word64 -> Word64 -> Word64 -> Word64+quotSmall c mask n = mulHi64 c n + (n .&. mask)++{-# INLINE quotSmallConst #-}+quotSmallConst :: Word64 -> Word64 -> Word64 -> Word64+quotSmallConst c mask n+  | n == 0 = 0+  | n == 1 = mask .&. 1+  -- mulhi c (2^k) == c >> (64 - k)+  | popCount n == 1 = unsafeShiftR c (64 - countTrailingZeros n) + (n .&. mask)+  | otherwise = quotSmall c mask n++-- | A precomputed 'Word32' divisor.+data StrengthReducedW32 = StrengthReducedW32+  {-# UNPACK #-} !Word64 -- reciprocal+  {-# UNPACK #-} !Word64 -- all ones when the divisor is 1+  {-# UNPACK #-} !Word32 -- divisor+  deriving stock (Eq, Show)++instance NFData StrengthReducedW32 where+  rnf !_ = ()++instance StrengthReduce Word32 where+  type StrengthReduced Word32 = StrengthReducedW32++  {-# INLINE new #-}+  new (NonZero d) = case newSmall (fromIntegral d) of+    (# c, mask #) -> StrengthReducedW32 c mask d++  {-# INLINE divisor #-}+  divisor (StrengthReducedW32 _ _ d) = d++  {-# INLINE divRem #-}+  divRem n (StrengthReducedW32 c mask d) =+    let !q = fromIntegral (quotSmall c mask (fromIntegral n))+    in (q, n - q * d)++  {-# INLINE divRemConst #-}+  divRemConst n (StrengthReducedW32 c mask d) =+    let !q = fromIntegral (quotSmallConst c mask (fromIntegral n))+    in (q, n - q * d)++  {-# INLINE divRemNonZero #-}+  divRemNonZero n (NonZero d) =+    let !q = quotWord32 n d+    in (q, n - q * d)++-- | A precomputed 'Word16' divisor.+data StrengthReducedW16 = StrengthReducedW16+  {-# UNPACK #-} !Word64 -- reciprocal+  {-# UNPACK #-} !Word64 -- all ones when the divisor is 1+  {-# UNPACK #-} !Word16 -- divisor+  deriving stock (Eq, Show)++instance NFData StrengthReducedW16 where+  rnf !_ = ()++instance StrengthReduce Word16 where+  type StrengthReduced Word16 = StrengthReducedW16++  {-# INLINE new #-}+  new (NonZero d) = case newSmall (fromIntegral d) of+    (# c, mask #) -> StrengthReducedW16 c mask d++  {-# INLINE divisor #-}+  divisor (StrengthReducedW16 _ _ d) = d++  {-# INLINE divRem #-}+  divRem n (StrengthReducedW16 c mask d) =+    let !q = fromIntegral (quotSmall c mask (fromIntegral n))+    in (q, n - q * d)++  {-# INLINE divRemConst #-}+  divRemConst n (StrengthReducedW16 c mask d) =+    let !q = fromIntegral (quotSmallConst c mask (fromIntegral n))+    in (q, n - q * d)++  {-# INLINE divRemNonZero #-}+  divRemNonZero n (NonZero d) =+    let !q = fromIntegral (quotWord32 (fromIntegral n) (fromIntegral d))+    in (q, n - q * d)++-- | A precomputed 'Word8' divisor.+data StrengthReducedW8 = StrengthReducedW8+  {-# UNPACK #-} !Word64 -- reciprocal+  {-# UNPACK #-} !Word64 -- all ones when the divisor is 1+  {-# UNPACK #-} !Word8  -- divisor+  deriving stock (Eq, Show)++instance NFData StrengthReducedW8 where+  rnf !_ = ()++instance StrengthReduce Word8 where+  type StrengthReduced Word8 = StrengthReducedW8++  {-# INLINE new #-}+  new (NonZero d) = case newSmall (fromIntegral d) of+    (# c, mask #) -> StrengthReducedW8 c mask d++  {-# INLINE divisor #-}+  divisor (StrengthReducedW8 _ _ d) = d++  {-# INLINE divRem #-}+  divRem n (StrengthReducedW8 c mask d) =+    let !q = fromIntegral (quotSmall c mask (fromIntegral n))+    in (q, n - q * d)++  {-# INLINE divRemConst #-}+  divRemConst n (StrengthReducedW8 c mask d) =+    let !q = fromIntegral (quotSmallConst c mask (fromIntegral n))+    in (q, n - q * d)++  {-# INLINE divRemNonZero #-}+  divRemNonZero n (NonZero d) =+    let !q = fromIntegral (quotWord32 (fromIntegral n) (fromIntegral d))+    in (q, n - q * d)++--------------------------------------------------------------------------------+-- Word128: Granlund-Montgomery again, on 64-bit limbs.+--------------------------------------------------------------------------------++-- | A precomputed 'Word128' divisor. data StrengthReducedW128 = StrengthReducedW128-  { multiplierHi :: {-#UNPACK #-} !Word128-  , multiplierLo :: {-#UNPACK #-} !Word128-  , divisor      :: {-#UNPACK #-} !Word128-  }+  {-# UNPACK #-} !Word64 -- multiplier hi+  {-# UNPACK #-} !Word64 -- multiplier lo+  {-# UNPACK #-} !Int    -- pre-shift, 0 or 1+  {-# UNPACK #-} !Int    -- post-shift+  {-# UNPACK #-} !Word64 -- divisor hi+  {-# UNPACK #-} !Word64 -- divisor lo+  deriving stock (Eq, Show)++instance NFData StrengthReducedW128 where+  rnf !_ = ()++instance StrengthReduce Word128 where+  type StrengthReduced Word128 = StrengthReducedW128++  {-# INLINE new #-}+  new = new128++  {-# INLINE divisor #-}+  divisor (StrengthReducedW128 _ _ _ _ d1 d0) = Word128 d1 d0++  {-# INLINE divRem #-}+  divRem (Word128 n1 n0) (StrengthReducedW128 m1 m0 s1 s2 d1 d0) =+    case mulHi128 m1 m0 n1 n0 of (# t1, t0 #) -> gm128 n1 n0 t1 t0 s1 s2 d1 d0++  {-# INLINE divRemConst #-}+  divRemConst (Word128 n1 n0) (StrengthReducedW128 m1 m0 s1 s2 d1 d0)+    | n1 == 0 && n0 == 0 = (Word128 0 0, Word128 0 0)+    | n1 == maxBound && n0 == maxBound =+        case sub128 m1 m0 0 1 of (# t1, t0 #) -> gm128 n1 n0 t1 t0 s1 s2 d1 d0+    | popCount n1 + popCount n0 == 1 =+        let !k = if n1 == 0 then countTrailingZeros n0 else 64 + countTrailingZeros n1+        in if k == 0+          then gm128 n1 n0 0 0 s1 s2 d1 d0+          else case shr128 m1 m0 (128 - k) of (# t1, t0 #) -> gm128 n1 n0 t1 t0 s1 s2 d1 d0+    -- n < 2^64: the multiplication is only 128x64, the quotient and remainder+    -- fit in 64 bits+    | n1 == 0 =+        -- floor (m * n / 2^128) <= n < 2^64, so it is just the top limb+        case mulHi128By64 m1 m0 n0 of+          (# t, _ #) ->+            let !x = t + unsafeShiftR (n0 - t) s1+                !q = if s2 < 64 then unsafeShiftR x s2 else 0+            in (Word128 0 q, Word128 0 (n0 - q * d0))+    | otherwise =+        case mulHi128 m1 m0 n1 n0 of (# t1, t0 #) -> gm128 n1 n0 t1 t0 s1 s2 d1 d0++  {-# INLINE divRemNonZero #-}+  divRemNonZero = divRemNonZero128++  {-# INLINE divRemNonZeroConst #-}+  divRemNonZeroConst n@(Word128 n1 n0) d@(NonZero (Word128 d1 d0))+    | n1 == 0 && n0 == 0 = (Word128 0 0, Word128 0 0)+    | n1 == 0 =+        if d1 /= 0+          then (Word128 0 0, n)+          else case divRemNonZeroConst n0 (NonZero d0) of+            (q, r) -> (Word128 0 q, Word128 0 r)+    | otherwise = divRemNonZero128 n d++{-# INLINE gm128 #-}+gm128 :: Word64 -> Word64 -> Word64 -> Word64 -> Int -> Int -> Word64 -> Word64 -> (Word128, Word128)+gm128 n1 n0 t1 t0 s1 s2 d1 d0 =+  case sub128 n1 n0 t1 t0 of { (# x1, x0 #) ->+  case shr128Small x1 x0 s1 of { (# y1, y0 #) ->+  case add128 y1 y0 t1 t0 of { (# z1, z0 #) ->+  case shr128 z1 z0 s2 of { (# q1, q0 #) ->+  case mulLo128 q1 q0 d1 d0 of { (# p1, p0 #) ->+  case sub128 n1 n0 p1 p0 of { (# r1, r0 #) ->+  (Word128 q1 q0, Word128 r1 r0) }}}}}}++{-# INLINE new128 #-}+new128 :: NonZero Word128 -> StrengthReducedW128+new128 (NonZero (Word128 d1 d0)) =+  case sub128 d1 d0 0 1 of { (# e1, e0 #) ->+  let !l = 128 - clz128 e1 e0+      !s1 = min l 1+      !s2 = max (l - 1) 0+      sr m1 m0 = StrengthReducedW128 m1 m0 s1 s2 d1 d0+  in+  case pow2Minus l of { (# h1, h0 #) ->+  if h1 == 0 && h0 == 0+    then sr 0 1+  else if d1 == 0+    -- h < d < 2^64: two 128/64 hardware divisions+    then+      case quotRem128By64 h0 0 d0 of { (# q1, r1 #) ->+      case quotRem128By64 r1 0 d0 of { (# q0, _ #) ->+      case add128 q1 q0 0 1 of { (# m1, m0 #) ->+      sr m1 m0 }}}+    -- normalize and do two 3-by-2 divisions of (h << s) * 2^128+    else+      let !s = countLeadingZeros d1 in+      case shl128Small d1 d0 s of { (# dn1, dn0 #) ->+      case shl128Small h1 h0 s of { (# u1, u0 #) ->+      case div3By2 u1 u0 0 dn1 dn0 of { (# q1, r1, r0 #) ->+      case div3By2 r1 r0 0 dn1 dn0 of { (# q0, _, _ #) ->+      case add128 q1 q0 0 1 of { (# m1, m0 #) ->+      sr m1 m0 }}}}}+  }}+  where+  -- 2^l - d, modulo 2^128+  pow2Minus :: Int -> (# Word64, Word64 #)+  pow2Minus l+    | l == 128 = sub128 0 0 d1 d0+    | l >= 64 = sub128 (unsafeShiftL 1 (l - 64)) 0 d1 d0+    | otherwise = sub128 0 (unsafeShiftL 1 l) d1 d0++{-# INLINE divRemNonZero128 #-}+divRemNonZero128 :: Word128 -> NonZero Word128 -> (Word128, Word128)+divRemNonZero128 (Word128 n1 n0) (NonZero (Word128 d1 d0))+  | d1 == 0 =+      if n1 < d0+        then case quotRem128By64 n1 n0 d0 of+          (# q, r #) -> (Word128 0 q, Word128 0 r)+        else+          -- NB. not quotRemWord64: when only the remainder is demanded GHC+          -- turns that into a remWord# expression and sinks it into the high+          -- argument of the quotRemWord2# below, and the x86 NCG (9.14 at+          -- least) then clobbers rax while computing it. quotRemWord2# is+          -- never sunk.+          case quotRem128By64 0 n1 d0 of { (# q1, r1 #) ->+          case quotRem128By64 r1 n0 d0 of { (# q0, r #) ->+          (Word128 q1 q0, Word128 0 r) }}+  | otherwise =+      -- the quotient fits in 64 bits: normalize, one 3-by-2 division+      let !s = countLeadingZeros d1+          !u2 = unsafeShiftR (unsafeShiftR n1 1) (63 - s)+      in+      case shl128Small d1 d0 s of { (# dn1, dn0 #) ->+      case shl128Small n1 n0 s of { (# u1, u0 #) ->+      case div3By2 u2 u1 u0 dn1 dn0 of { (# q, r1, r0 #) ->+      case shr128Small r1 r0 s of { (# rr1, rr0 #) ->+      (Word128 0 q, Word128 rr1 rr0) }}}}

@@ -1,84 +0,0 @@-{-# LANGUAGE-    TypeApplications-  , ScopedTypeVariables-  , LambdaCase-  , NumericUnderscores-#-}--module Atrophy.Internal.LongDivision where--import Data.WideWord.Word128-import Data.Word-import Data.Bits---- divides a 128-bit number by a 64-bit divisor, returning the quotient as a 64-bit number--- assumes that the divisor and numerator have both already been bit-shifted so that countLeadingZeros divisor == 0-{-# INLINE divide128By64Preshifted #-}-divide128By64Preshifted :: Word64 -> Word64 -> Word64 -> Word64-divide128By64Preshifted numeratorHi numeratorLo' divisor =-  let-    numeratorMid = fromIntegral @Word64 @Word128 (numeratorLo' `unsafeShiftR` 32)-    numeratorLo = fromIntegral @Word32 @Word128 (fromIntegral @Word64 @Word32 numeratorLo')-    divisorFull128 = fromIntegral @Word64 @Word128 divisor-    divisorHi = divisor `unsafeShiftR` 32--    -- To get the upper 32 bits of the quotient, we want to divide 'fullUpperNumerator' by 'divisor'-    -- but the problem is, fullUpperNumerator is a 96-bit number, meaning we would need to use u128 to do the division all at once, and the whole point of this is that we don't want to do 128 bit divison because it's slow-    -- so instead, we'll shift both the numerator and divisor right by 32, giving us a 64 bit / 32 bit division. This won't give us the exact quotient -- but it will be close.-    fullUpperNumerator = (Word128 0 numeratorHi `unsafeShiftL` 32) .|. numeratorMid-    quotientHi :: Word64-    quotientHi = min (numeratorHi `div` divisorHi) (fromIntegral $ maxBound @Word32)-    productHi = (Word128 0 quotientHi) * divisorFull128--  -- quotientHi contains our guess at what the quotient is! the problem is that we got this by ignoring the lower 32 bits of the divisor. when we account for that, the quotient might be slightly lower-  -- we will know our quotient is too high if quotient * divisor > numerator. if it is, decrement until it's in range-    (productHi', quotientHi') = clampToFull productHi quotientHi divisorFull128 fullUpperNumerator--    remainderHi = fullUpperNumerator - productHi'--  -- repeat the process using the lower half of the numerator-    fullLowerNumerator = (remainderHi `unsafeShiftL` 32) .|. numeratorLo--    quotientLo = min ((fromIntegral @_ @Word64 remainderHi) `div` divisorHi) (fromIntegral $ maxBound @Word32)-    productLo = (Word128 0 quotientLo) * divisorFull128--  -- again, quotientLo is just a guess at this point, it might be slightly too large-    (_, quotientLo') = clampToFull productLo quotientLo divisorFull128 fullLowerNumerator--  -- We now have our separate quotients, now we just have to add them together-  in (quotientHi' `unsafeShiftL` 32) .|. quotientLo'--divide128MaxBy64 :: Word64 -> Word128-divide128MaxBy64 divisor =-  let-    quotientHi = maxBound @Word64 `div` divisor;-    remainderHi = maxBound @Word64 - quotientHi * divisor;--    leadingZeros = countLeadingZeros divisor-    quotientLo = if leadingZeros >= 32-      then-        let-          numeratorMid = (remainderHi `unsafeShiftL` 32) .|. (fromIntegral (maxBound @Word32))-          quotientMid = numeratorMid `div` divisor;-          remainderMid = numeratorMid - quotientMid * divisor;--          numeratorLo = (remainderMid `unsafeShiftL` 32) .|. (fromIntegral (maxBound @Word32))-          quotientLo' = numeratorLo `div` divisor--        in (quotientMid `unsafeShiftL` 32) .|. quotientLo'-      else-        let-          numeratorHi = if leadingZeros > 0-            then (remainderHi `unsafeShiftL` leadingZeros) .|. (maxBound @Word64 `unsafeShiftR` (64 - leadingZeros))-            else remainderHi-          numeratorLo = maxBound @Word64 `unsafeShiftL` leadingZeros;-        in divide128By64Preshifted numeratorHi numeratorLo (divisor `unsafeShiftL` leadingZeros)-  in ((fromIntegral quotientHi) `unsafeShiftL` 64) .|. (fromIntegral quotientLo)--clampToFull :: Word128 -> Word64 -> Word128 -> Word128 -> (Word128, Word64)-clampToFull product' quotient' divisorFull128 fullUpperNumerator = go product' quotient'-  where-  go prod quotient =-    if prod > fullUpperNumerator-    then go (prod - divisorFull128) (quotient - 1)-    else (prod, quotient)

@@ -0,0 +1,238 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE MagicHash #-}+{-# LANGUAGE UnboxedTuples #-}+{-# LANGUAGE AllowAmbiguousTypes #-}+{-# OPTIONS_HADDOCK hide #-}++#include "MachDeps.h"++-- | Double-word and multi-limb primitives. Everything in here is @INLINE@ and+-- is expected to compile down to a handful of instructions. Limbs are passed+-- most significant first, results are returned as unboxed tuples so that+-- nothing is ever allocated.+module Atrophy.Internal.Prim+  ( -- * 64-bit+    mulHi64+  , mulFull64+  , quotRem128By64+  , quotRemWord64+  , quotWord32+  , ltW+  , addCarry64+  , subBorrow64+  , adc64+    -- * 128-bit, as pairs of limbs+  , add128+  , sub128+  , shr128+  , shr128Small+  , shl128Small+  , clz128+  , mulLo128+  , mulHi128+  , mulHi128By64+  , div3By2+    -- * Type level+  , natWord64+  , natInt+  ) where++import Data.Bits+import GHC.Exts+import GHC.TypeNats (KnownNat, natVal')+import GHC.Word++#if WORD_SIZE_IN_BITS == 64++-- | High 64 bits of the 128-bit product. A single @mul@.+{-# INLINE mulHi64 #-}+mulHi64 :: Word64 -> Word64 -> Word64+mulHi64 (W64# a) (W64# b) = case timesWord2# (word64ToWord# a) (word64ToWord# b) of+  (# h, _ #) -> W64# (wordToWord64# h)++-- | @(# hi, lo #)@ of the 128-bit product. A single @mul@.+{-# INLINE mulFull64 #-}+mulFull64 :: Word64 -> Word64 -> (# Word64, Word64 #)+mulFull64 (W64# a) (W64# b) = case timesWord2# (word64ToWord# a) (word64ToWord# b) of+  (# h, l #) -> (# W64# (wordToWord64# h), W64# (wordToWord64# l) #)++-- | @quotRem128By64 hi lo d@ divides @hi * 2^64 + lo@ by @d@. Requires @hi < d@.+-- A single @div@.+{-# INLINE quotRem128By64 #-}+quotRem128By64 :: Word64 -> Word64 -> Word64 -> (# Word64, Word64 #)+quotRem128By64 (W64# h) (W64# l) (W64# d) =+  case quotRemWord2# (word64ToWord# h) (word64ToWord# l) (word64ToWord# d) of+    (# q, r #) -> (# W64# (wordToWord64# q), W64# (wordToWord64# r) #)++-- | Unchecked division, the divisor must not be zero.+{-# INLINE quotRemWord64 #-}+quotRemWord64 :: Word64 -> Word64 -> (# Word64, Word64 #)+quotRemWord64 (W64# n) (W64# d) = case quotRemWord# (word64ToWord# n) (word64ToWord# d) of+  (# q, r #) -> (# W64# (wordToWord64# q), W64# (wordToWord64# r) #)++-- | 1 if @a < b@, 0 otherwise. Branchless.+{-# INLINE ltW #-}+ltW :: Word64 -> Word64 -> Word64+ltW (W64# a) (W64# b) = W64# (wordToWord64# (int2Word# (ltWord# (word64ToWord# a) (word64ToWord# b))))++#else++{-# INLINE mulHi64 #-}+mulHi64 :: Word64 -> Word64 -> Word64+mulHi64 a b = case mulFull64 a b of (# h, _ #) -> h++{-# INLINE mulFull64 #-}+mulFull64 :: Word64 -> Word64 -> (# Word64, Word64 #)+mulFull64 a b =+  let !aL = a .&. 0xffffffff+      !aH = a `unsafeShiftR` 32+      !bL = b .&. 0xffffffff+      !bH = b `unsafeShiftR` 32+      !ll = aL * bL+      !lh = aL * bH+      !hl = aH * bL+      !hh = aH * bH+      !mid = (ll `unsafeShiftR` 32) + (lh .&. 0xffffffff) + (hl .&. 0xffffffff)+  in (# hh + (lh `unsafeShiftR` 32) + (hl `unsafeShiftR` 32) + (mid `unsafeShiftR` 32), a * b #)++{-# INLINE quotRem128By64 #-}+quotRem128By64 :: Word64 -> Word64 -> Word64 -> (# Word64, Word64 #)+quotRem128By64 h l d =+  case (toInteger h `unsafeShiftL` 64 .|. toInteger l) `quotRem` toInteger d of+    (q, r) -> (# fromInteger q, fromInteger r #)++{-# INLINE quotRemWord64 #-}+quotRemWord64 :: Word64 -> Word64 -> (# Word64, Word64 #)+quotRemWord64 n d = case quotRem n d of (q, r) -> (# q, r #)++{-# INLINE ltW #-}+ltW :: Word64 -> Word64 -> Word64+ltW a b = if a < b then 1 else 0++#endif++-- | Unchecked division, the divisor must not be zero.+{-# INLINE quotWord32 #-}+quotWord32 :: Word32 -> Word32 -> Word32+quotWord32 (W32# n) (W32# d) = W32# (quotWord32# n d)++-- | @(# sum, carry #)@+{-# INLINE addCarry64 #-}+addCarry64 :: Word64 -> Word64 -> (# Word64, Word64 #)+addCarry64 a b = let !s = a + b in (# s, ltW s a #)++-- | @(# difference, borrow #)@+{-# INLINE subBorrow64 #-}+subBorrow64 :: Word64 -> Word64 -> (# Word64, Word64 #)+subBorrow64 a b = (# a - b, ltW a b #)++-- | @a + b + carry@, returning @(# sum, carry #)@. @carry@ must be 0 or 1.+{-# INLINE adc64 #-}+adc64 :: Word64 -> Word64 -> Word64 -> (# Word64, Word64 #)+adc64 a b c =+  let !s1 = a + b+      !s = s1 + c+  in (# s, ltW s1 a + ltW s s1 #)++{-# INLINE add128 #-}+add128 :: Word64 -> Word64 -> Word64 -> Word64 -> (# Word64, Word64 #)+add128 a1 a0 b1 b0 = case addCarry64 a0 b0 of (# s, c #) -> (# a1 + b1 + c, s #)++{-# INLINE sub128 #-}+sub128 :: Word64 -> Word64 -> Word64 -> Word64 -> (# Word64, Word64 #)+sub128 a1 a0 b1 b0 = case subBorrow64 a0 b0 of (# s, c #) -> (# a1 - b1 - c, s #)++-- | Shift right by @0 <= s < 64@, branchless.+{-# INLINE shr128Small #-}+shr128Small :: Word64 -> Word64 -> Int -> (# Word64, Word64 #)+shr128Small x1 x0 s =+  (# x1 `unsafeShiftR` s+  , (x0 `unsafeShiftR` s) .|. ((x1 `unsafeShiftL` 1) `unsafeShiftL` (63 - s))+  #)++-- | Shift left by @0 <= s < 64@, branchless.+{-# INLINE shl128Small #-}+shl128Small :: Word64 -> Word64 -> Int -> (# Word64, Word64 #)+shl128Small x1 x0 s =+  (# (x1 `unsafeShiftL` s) .|. ((x0 `unsafeShiftR` 1) `unsafeShiftR` (63 - s))+  , x0 `unsafeShiftL` s+  #)++-- | Shift right by @0 <= s < 128@.+{-# INLINE shr128 #-}+shr128 :: Word64 -> Word64 -> Int -> (# Word64, Word64 #)+shr128 x1 x0 s+  | s < 64 = shr128Small x1 x0 s+  | otherwise = (# 0, x1 `unsafeShiftR` (s - 64) #)++{-# INLINE clz128 #-}+clz128 :: Word64 -> Word64 -> Int+clz128 x1 x0 = if x1 == 0 then 64 + countLeadingZeros x0 else countLeadingZeros x1++-- | Low 128 bits of the product.+{-# INLINE mulLo128 #-}+mulLo128 :: Word64 -> Word64 -> Word64 -> Word64 -> (# Word64, Word64 #)+mulLo128 a1 a0 b1 b0 = case mulFull64 a0 b0 of+  (# h, l #) -> (# h + a0 * b1 + a1 * b0, l #)++-- | High 128 bits of the 256-bit product.+{-# INLINE mulHi128 #-}+mulHi128 :: Word64 -> Word64 -> Word64 -> Word64 -> (# Word64, Word64 #)+mulHi128 a1 a0 b1 b0 =+  case mulFull64 a0 b0 of { (# h00, _ #) ->+  case mulFull64 a0 b1 of { (# h01, l01 #) ->+  case mulFull64 a1 b0 of { (# h10, l10 #) ->+  case mulFull64 a1 b1 of { (# h11, l11 #) ->+  case addCarry64 h00 l01 of { (# s1, c1a #) ->+  case addCarry64 s1 l10 of { (# _, c1b #) ->+  case addCarry64 h01 h10 of { (# s2a, c2a #) ->+  case addCarry64 s2a l11 of { (# s2b, c2b #) ->+  case addCarry64 s2b (c1a + c1b) of { (# s2, c2c #) ->+  (# h11 + c2a + c2b + c2c, s2 #) }}}}}}}}}++-- | High 128 bits of the 192-bit product of a 128-bit and a 64-bit number.+{-# INLINE mulHi128By64 #-}+mulHi128By64 :: Word64 -> Word64 -> Word64 -> (# Word64, Word64 #)+mulHi128By64 a1 a0 b =+  case mulFull64 a1 b of { (# h1, l1 #) ->+  case mulHi64 a0 b of { h0 ->+  case addCarry64 l1 h0 of { (# s, c #) ->+  (# h1 + c, s #) }}}++-- | Divide the 192-bit @(u2, u1, u0)@ by the normalized (top bit set) 128-bit+-- @(d1, d0)@, given @(u2, u1) < (d1, d0)@. The quotient fits in 64 bits.+-- Returns @(# quotient, remainder hi, remainder lo #)@.+--+-- Knuth's algorithm D, one step: estimate from the top limbs with a hardware+-- division, then correct at most twice.+{-# INLINE div3By2 #-}+div3By2 :: Word64 -> Word64 -> Word64 -> Word64 -> Word64 -> (# Word64, Word64, Word64 #)+div3By2 u2 u1 u0 d1 d0 =+  let !qhat = if u2 >= d1 then maxBound else case quotRem128By64 u2 u1 d1 of (# q, _ #) -> q+  in+  case mulFull64 qhat d0 of { (# a1, a0 #) ->+  case mulFull64 qhat d1 of { (# b1, b0 #) ->+  case addCarry64 b0 a1 of { (# p1, c #) ->+  correct qhat (b1 + c) p1 a0 }}}+  where+  correct :: Word64 -> Word64 -> Word64 -> Word64 -> (# Word64, Word64, Word64 #)+  correct !q !p2 !p1 !p0+    | p2 > u2 || (p2 == u2 && (p1 > u1 || (p1 == u1 && p0 > u0))) =+        case subBorrow64 p0 d0 of { (# p0', b0 #) ->+        case subBorrow64 p1 d1 of { (# t1, b1a #) ->+        case subBorrow64 t1 b0 of { (# p1', b1b #) ->+        correct (q - 1) (p2 - b1a - b1b) p1' p0' }}}+    | otherwise =+        -- the remainder is below d, so its top limb is zero+        case subBorrow64 u0 p0 of { (# r0, b0 #) ->+        (# q, u1 - p1 - b0, r0 #) }++-- | A type level natural as a 'Word64' literal. Folds at compile time.+{-# INLINE natWord64 #-}+natWord64 :: forall n. KnownNat n => Word64+natWord64 = fromIntegral (natVal' (proxy# :: Proxy# n))++-- | A type level natural as an 'Int' literal. Folds at compile time.+{-# INLINE natInt #-}+natInt :: forall n. KnownNat n => Int+natInt = fromIntegral (natVal' (proxy# :: Proxy# n))

@@ -0,0 +1,261 @@+{-# LANGUAGE AllowAmbiguousTypes #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE MagicHash #-}+{-# LANGUAGE NoStarIsType #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE UnboxedTuples #-}+{-# LANGUAGE UndecidableInstances #-}++-- | Division where one side is known at compile time, supplied at the type+-- level so that it is guaranteed to be a constant.+--+-- GHC's native code generator does /not/ strength-reduce division by+-- constants: @x \`quot\` 7@ compiles to a @div@ instruction. 'divK' computes+-- the magic numbers during type checking instead:+--+-- >>> divK @7 (100 :: Word64)+-- 14+--+-- 'divRemN' goes the other way: the divisor is a runtime 'StrengthReduced'+-- value, but the numerator is fixed, which lets e.g. @divN \@(2 ^ 63)@ skip+-- the multiplication entirely.+module Atrophy.Known+  ( -- * Divisors known at compile time+    KnownDivisor (..)+  , divK+  , remK+    -- * Numerators known at compile time+  , KnownNumerator (..)+  , divRemN+  , divN+  , remN+  , divRemNonZeroN+  , divNonZeroN+  , remNonZeroN+  ) where++import Atrophy.Internal+import Atrophy.Internal.Prim+import Data.Bits+import Data.Kind (Constraint, Type)+import Data.Type.Bool (If)+import Data.Type.Equality (type (==))+import Data.WideWord.Word128 (Word128 (..))+import Data.Word+import GHC.TypeLits (ErrorMessage (..), TypeError)+import GHC.TypeNats++type Assert :: Bool -> Constraint -> Constraint+type family Assert b err where+  Assert 'True _ = ()+  Assert _ err = err++type CheckDivisor :: Nat -> Nat -> Constraint+type CheckDivisor bits d =+  ( Assert (1 <=? d) (TypeError ('Text "atrophy: division by zero"))+  , Assert (d <=? 2 ^ bits - 1)+      (TypeError ('Text "atrophy: the divisor " ':<>: 'ShowType d ':<>: 'Text " does not fit in " ':<>: 'ShowType bits ':<>: 'Text " bits"))+  )++type CheckNumerator :: Nat -> Nat -> Constraint+type CheckNumerator bits n =+  Assert (n <=? 2 ^ bits - 1)+    (TypeError ('Text "atrophy: the numerator " ':<>: 'ShowType n ':<>: 'Text " does not fit in " ':<>: 'ShowType bits ':<>: 'Text " bits"))++--------------------------------------------------------------------------------+-- Magic numbers, following libdivide's unsigned branching algorithm. Since the+-- choice of algorithm is made at compile time there is no branch at runtime.+--+-- 0: power of two, shift only+-- 1: q = mulhi m n >> s+-- 2: q = (((n - t) >> 1) + t) >> s  where t = mulhi m n+--------------------------------------------------------------------------------++type IsPow2 d = 2 ^ Log2 d == d+type ProposedM bits d = Div (2 ^ (bits + Log2 d)) d+type ProposedR bits d = Mod (2 ^ (bits + Log2 d)) d+type NoAdd bits d = CmpNat (d - ProposedR bits d) (2 ^ Log2 d) == 'LT++type Algo bits d = If (IsPow2 d) 0 (If (NoAdd bits d) 1 2)++type Magic bits d =+  If (IsPow2 d) 0+    (If (NoAdd bits d)+      (ProposedM bits d + 1)+      (Mod (2 * ProposedM bits d + If (d <=? 2 * ProposedR bits d) 1 0 + 1) (2 ^ bits)))++-- | For 32 bits and below: q = mulhi c n, with c = ceil (2^64 / d).+type Lemire d = Div (2 ^ 64) d + 1++-- | Types that can be divided by the type-level constant @d@.+--+-- @divRemK \@d n == n \`quotRem\` d@, but the division is replaced by at most+-- one multiplication and some shifts, with no branches.+type KnownDivisor :: Nat -> Type -> Constraint+class KnownDivisor d a where+  divRemK :: a -> (a, a)++-- | Quotient by the type-level divisor @d@.+{-# INLINE divK #-}+divK :: forall d a. KnownDivisor d a => a -> a+divK n = case divRemK @d n of (q, _) -> q++-- | Remainder by the type-level divisor @d@.+{-# INLINE remK #-}+remK :: forall d a. KnownDivisor d a => a -> a+remK n = case divRemK @d n of (_, r) -> r++instance+  ( CheckDivisor 64 d+  , KnownNat d+  , KnownNat (Log2 d)+  , KnownNat (Algo 64 d)+  , KnownNat (Magic 64 d)+  ) => KnownDivisor d Word64 where+  {-# INLINE divRemK #-}+  divRemK n = case natWord64 @(Algo 64 d) of+    0 -> (unsafeShiftR n sh, n .&. (dv - 1))+    1 -> let !q = unsafeShiftR (mulHi64 mg n) sh in (q, n - q * dv)+    _ -> let !t = mulHi64 mg n+             !q = unsafeShiftR (unsafeShiftR (n - t) 1 + t) sh+         in (q, n - q * dv)+    where+    sh = natInt @(Log2 d)+    mg = natWord64 @(Magic 64 d)+    dv = natWord64 @d++{-# INLINE smallK #-}+smallK :: Word64 -> Int -> Word64 -> Word64 -> Word64 -> (Word64, Word64)+smallK algo sh c dv n = case algo of+  0 -> (unsafeShiftR n sh, n .&. (dv - 1))+  _ -> let !q = mulHi64 c n in (q, n - q * dv)++instance+  ( CheckDivisor 32 d+  , KnownNat d+  , KnownNat (Log2 d)+  , KnownNat (If (IsPow2 d) 0 1)+  , KnownNat (Lemire d)+  ) => KnownDivisor d Word32 where+  {-# INLINE divRemK #-}+  divRemK n =+    case smallK (natWord64 @(If (IsPow2 d) 0 1)) (natInt @(Log2 d)) (natWord64 @(Lemire d)) (natWord64 @d) (fromIntegral n) of+      (q, r) -> (fromIntegral q, fromIntegral r)++instance+  ( CheckDivisor 16 d+  , KnownNat d+  , KnownNat (Log2 d)+  , KnownNat (If (IsPow2 d) 0 1)+  , KnownNat (Lemire d)+  ) => KnownDivisor d Word16 where+  {-# INLINE divRemK #-}+  divRemK n =+    case smallK (natWord64 @(If (IsPow2 d) 0 1)) (natInt @(Log2 d)) (natWord64 @(Lemire d)) (natWord64 @d) (fromIntegral n) of+      (q, r) -> (fromIntegral q, fromIntegral r)++instance+  ( CheckDivisor 8 d+  , KnownNat d+  , KnownNat (Log2 d)+  , KnownNat (If (IsPow2 d) 0 1)+  , KnownNat (Lemire d)+  ) => KnownDivisor d Word8 where+  {-# INLINE divRemK #-}+  divRemK n =+    case smallK (natWord64 @(If (IsPow2 d) 0 1)) (natInt @(Log2 d)) (natWord64 @(Lemire d)) (natWord64 @d) (fromIntegral n) of+      (q, r) -> (fromIntegral q, fromIntegral r)++instance+  ( CheckDivisor 128 d+  , KnownNat (Log2 d)+  , KnownNat (Algo 128 d)+  , KnownNat (Div (Magic 128 d) (2 ^ 64))+  , KnownNat (Mod (Magic 128 d) (2 ^ 64))+  , KnownNat (Div d (2 ^ 64))+  , KnownNat (Mod d (2 ^ 64))+  ) => KnownDivisor d Word128 where+  {-# INLINE divRemK #-}+  divRemK (Word128 n1 n0) = case natWord64 @(Algo 128 d) of+    0 ->+      case shr128 n1 n0 sh of { (# q1, q0 #) ->+      case sub128 d1 d0 0 1 of { (# e1, e0 #) ->+      (Word128 q1 q0, Word128 (n1 .&. e1) (n0 .&. e0)) }}+    1 ->+      case mulHi128 m1 m0 n1 n0 of { (# t1, t0 #) ->+      case shr128 t1 t0 sh of { (# q1, q0 #) ->+      finish q1 q0 }}+    _ ->+      case mulHi128 m1 m0 n1 n0 of { (# t1, t0 #) ->+      case sub128 n1 n0 t1 t0 of { (# x1, x0 #) ->+      case shr128Small x1 x0 1 of { (# y1, y0 #) ->+      case add128 y1 y0 t1 t0 of { (# z1, z0 #) ->+      case shr128 z1 z0 sh of { (# q1, q0 #) ->+      finish q1 q0 }}}}}+    where+    sh = natInt @(Log2 d)+    m1 = natWord64 @(Div (Magic 128 d) (2 ^ 64))+    m0 = natWord64 @(Mod (Magic 128 d) (2 ^ 64))+    d1 = natWord64 @(Div d (2 ^ 64))+    d0 = natWord64 @(Mod d (2 ^ 64))+    finish q1 q0 =+      case mulLo128 q1 q0 d1 d0 of { (# p1, p0 #) ->+      case sub128 n1 n0 p1 p0 of { (# r1, r0 #) ->+      (Word128 q1 q0, Word128 r1 r0) }}++-- | Type-level numerators.+type KnownNumerator :: Nat -> Type -> Constraint+class KnownNumerator n a where+  -- | The numerator as a value.+  numeratorVal :: a++instance (CheckNumerator 128 n, KnownNat (Div n (2 ^ 64)), KnownNat (Mod n (2 ^ 64))) => KnownNumerator n Word128 where+  {-# INLINE numeratorVal #-}+  numeratorVal = Word128 (natWord64 @(Div n (2 ^ 64))) (natWord64 @(Mod n (2 ^ 64)))++instance (CheckNumerator 64 n, KnownNat n) => KnownNumerator n Word64 where+  {-# INLINE numeratorVal #-}+  numeratorVal = natWord64 @n++instance (CheckNumerator 32 n, KnownNat n) => KnownNumerator n Word32 where+  {-# INLINE numeratorVal #-}+  numeratorVal = fromIntegral (natWord64 @n)++instance (CheckNumerator 16 n, KnownNat n) => KnownNumerator n Word16 where+  {-# INLINE numeratorVal #-}+  numeratorVal = fromIntegral (natWord64 @n)++instance (CheckNumerator 8 n, KnownNat n) => KnownNumerator n Word8 where+  {-# INLINE numeratorVal #-}+  numeratorVal = fromIntegral (natWord64 @n)++-- | @divRemN \@n d == divRem n d@ for a type-level numerator. See 'divRemConst'.+{-# INLINE divRemN #-}+divRemN :: forall n a. (KnownNumerator n a, StrengthReduce a) => StrengthReduced a -> (a, a)+divRemN = divRemConst (numeratorVal @n)++-- | Quotient of the type-level numerator @n@.+{-# INLINE divN #-}+divN :: forall n a. (KnownNumerator n a, StrengthReduce a) => StrengthReduced a -> a+divN d = case divRemN @n d of (q, _) -> q++-- | Remainder of the type-level numerator @n@.+{-# INLINE remN #-}+remN :: forall n a. (KnownNumerator n a, StrengthReduce a) => StrengthReduced a -> a+remN d = case divRemN @n d of (_, r) -> r++-- | @divRemNonZeroN \@n d == divRemNonZero n d@ for a type-level numerator.+-- See 'divRemNonZeroConst'.+{-# INLINE divRemNonZeroN #-}+divRemNonZeroN :: forall n a. (KnownNumerator n a, StrengthReduce a) => NonZero a -> (a, a)+divRemNonZeroN = divRemNonZeroConst (numeratorVal @n)++-- | Unchecked hardware quotient of the type-level numerator @n@.+{-# INLINE divNonZeroN #-}+divNonZeroN :: forall n a. (KnownNumerator n a, StrengthReduce a) => NonZero a -> a+divNonZeroN d = case divRemNonZeroN @n d of (q, _) -> q++-- | Unchecked hardware remainder of the type-level numerator @n@.+{-# INLINE remNonZeroN #-}+remNonZeroN :: forall n a. (KnownNumerator n a, StrengthReduce a) => NonZero a -> a+remNonZeroN d = case divRemNonZeroN @n d of (_, r) -> r

@@ -1,46 +1,125 @@-{-# LANGUAGE-    TypeApplications-  , ScopedTypeVariables-  , LambdaCase-  , NumericUnderscores-#-}+{-# LANGUAGE MagicHash #-}+{-# LANGUAGE UnboxedTuples #-}+{-# LANGUAGE DerivingStrategies #-} +-- | Dividing multi-limb numbers by a single 64-bit limb.+--+-- Uses Möller & Granlund, "Improved division by invariant integers": a 128/64+-- division costs one multiplication-and-a-half and a predictable branch, given+-- a reciprocal precomputed with a single hardware division.+--+-- Limbs are little-endian: index 0 is the least significant. module Atrophy.LongDivision-  ( module X-  , module Atrophy.LongDivision-  )-  where+  ( Divisor2By1+  , newDivisor2By1+  , divisor2By1+  , divRem2By1+  , longDivision+  , longDivisionInPlace+  ) where -import Data.Word-import Atrophy.Internal.LongDivision as X-import Atrophy.Internal-import qualified Data.Primitive.Contiguous as Contiguous-import Data.Primitive.Contiguous (PrimArray, Mutable, Sliced)-import Control.Monad.ST.Strict (ST)-import Data.STRef.Strict (newSTRef, readSTRef, writeSTRef)+import Atrophy.Internal (NonZero (..))+import Atrophy.Internal.Prim+import Control.DeepSeq (NFData (..))+import Control.Monad.ST (ST) import Data.Bits+import Data.Primitive.PrimArray+import Data.Word -{-# NOINLINE longDivision #-}-longDivision :: forall s. Sliced PrimArray Word64 -> StrengthReducedW64 -> Mutable PrimArray s Word64 -> ST s ()-longDivision numeratorSlice reducedDivisor quotient = do-  remainder <- newSTRef 0-  (flip Contiguous.itraverse_) numeratorSlice $ \i numerator -> do-    readSTRef remainder >>= \case-      0 -> do-        -- The remainder is zero, which means we can take a shortcut and only do a single division!-        let (digitQuotient, digitRemainder) = divRem numerator reducedDivisor+-- | A precomputed 64-bit divisor for 128/64 division.+data Divisor2By1 = Divisor2By1+  {-# UNPACK #-} !Word64 -- normalized divisor, top bit set+  {-# UNPACK #-} !Word64 -- reciprocal: floor ((2^128 - 1) / dn) - 2^64+  {-# UNPACK #-} !Int    -- normalization shift+  deriving stock (Eq, Show) -        Contiguous.write quotient i digitQuotient-        writeSTRef remainder digitRemainder+instance NFData Divisor2By1 where+  rnf !_ = () -      remainder' -> do-        -- Do one division that includes the running remainder and the upper half of this numerator element,-        -- then a second division for the first division's remainder combinedwith the lower half-        let upperNumerator = (remainder' `unsafeShiftL` 32) .|. (numerator `unsafeShiftR` 32)-        let (upperQuotient, upperRemainder) = divRem upperNumerator reducedDivisor+{-# INLINE newDivisor2By1 #-}+newDivisor2By1 :: NonZero Word64 -> Divisor2By1+newDivisor2By1 (NonZero d) =+  let !s = countLeadingZeros d+      !dn = unsafeShiftL d s+  in case quotRem128By64 (complement dn) maxBound dn of+    (# v, _ #) -> Divisor2By1 dn v s -        let lowerNumerator = (upperRemainder `unsafeShiftL` 32) .|. (0x00000000_ffffffff .&. numerator)-        let (lowerQuotient, lowerRemainder) = divRem lowerNumerator reducedDivisor+-- | The original divisor.+{-# INLINE divisor2By1 #-}+divisor2By1 :: Divisor2By1 -> Word64+divisor2By1 (Divisor2By1 dn _ s) = unsafeShiftR dn s -        Contiguous.write quotient i $ (upperQuotient `unsafeShiftL` 32) .|. lowerQuotient-        writeSTRef remainder lowerRemainder+-- | Möller-Granlund algorithm 4, on a normalized divisor. Requires @u1 < dn@.+{-# INLINE udivrem2By1 #-}+udivrem2By1 :: Word64 -> Word64 -> Word64 -> Word64 -> (# Word64, Word64 #)+udivrem2By1 u1 u0 dn v =+  case mulFull64 v u1 of { (# p1, p0 #) ->+  case addCarry64 p0 u0 of { (# q0, c #) ->+  let !q1 = p1 + u1 + c + 1+      !r = u0 - q1 * dn+      -- if r > q0 then (q1 - 1, r + dn), branchless since it is unpredictable+      !mask = negate (ltW q0 r)+      !q1' = q1 + mask+      !r' = r + (mask .&. dn)+  in if r' >= dn+    then (# q1' + 1, r' - dn #)+    else (# q1', r' #) }}++-- | @divRem2By1 d hi lo@ divides @hi * 2^64 + lo@ by @d@. Requires @hi < d@.+{-# INLINE divRem2By1 #-}+divRem2By1 :: Divisor2By1 -> Word64 -> Word64 -> (Word64, Word64)+divRem2By1 (Divisor2By1 dn v s) hi lo =+  case shl128Small hi lo s of+    (# u1, u0 #) -> case udivrem2By1 u1 u0 dn v of+      (# q, r #) -> (q, unsafeShiftR r s)++-- | Divide the little-endian number in the first argument, writing the+-- quotient limbs to the mutable array, which must be at least as long.+-- Returns the remainder.+longDivision :: Divisor2By1 -> PrimArray Word64 -> MutablePrimArray s Word64 -> ST s Word64+longDivision d numerator quotient = do+  let !len = sizeofPrimArray numerator+  qlen <- getSizeofMutablePrimArray quotient+  if qlen < len+    then error "Atrophy.LongDivision.longDivision: quotient array is too small"+    else do+      copyPrimArray quotient 0 numerator 0 len+      divInPlace d quotient len++-- | Divide the little-endian number in place, replacing it with the quotient.+-- Returns the remainder.+longDivisionInPlace :: Divisor2By1 -> MutablePrimArray s Word64 -> ST s Word64+longDivisionInPlace d arr = getSizeofMutablePrimArray arr >>= divInPlace d arr++{-# INLINE divInPlace #-}+divInPlace :: Divisor2By1 -> MutablePrimArray s Word64 -> Int -> ST s Word64+divInPlace (Divisor2By1 dn v s) arr len+  | len == 0 = pure 0+  | otherwise = do+      -- Normalize in a separate pass rather than on the fly: two tight loops+      -- beat one loop juggling shifts, since GHC's register allocator spills.+      r0 <- if s == 0 then pure 0 else shiftLeftInPlace arr len s+      let go !i !r+            | i < 0 = pure (unsafeShiftR r s)+            | otherwise = do+                u0 <- readPrimArray arr i+                case udivrem2By1 r u0 dn v of+                  (# q, r' #) -> do+                    writePrimArray arr i q+                    go (i - 1) r'+      go (len - 1) r0++-- | Shift the little-endian number left by @0 < s < 64@, returning the bits+-- shifted out of the top.+shiftLeftInPlace :: MutablePrimArray s Word64 -> Int -> Int -> ST s Word64+shiftLeftInPlace arr len s = do+  top <- readPrimArray arr (len - 1)+  let !rs = 64 - s+      go !i !cur+        | i == 0 = writePrimArray arr 0 (unsafeShiftL cur s)+        | otherwise = do+            next <- readPrimArray arr (i - 1)+            writePrimArray arr i (unsafeShiftL cur s .|. unsafeShiftR next rs)+            go (i - 1) next+  go (len - 1) top+  pure (unsafeShiftR top rs)

@@ -1,90 +1,76 @@-{-# LANGUAGE-    TypeApplications-  , ScopedTypeVariables-  , LambdaCase-#-}+{-# LANGUAGE MagicHash #-}+{-# LANGUAGE UnboxedTuples #-} -module Atrophy.LongMultiplication where+-- | Multi-limb multiplication. Limbs are little-endian: index 0 is the least+-- significant.+module Atrophy.LongMultiplication+  ( multiply256By128UpperBits+  , longMultiply+  ) where -import Data.WideWord.Word128+import Atrophy.Internal.Prim+import Control.Monad.ST (ST)+import Data.Primitive.PrimArray+import Data.WideWord.Word128 (Word128 (..)) import Data.Word-import qualified Data.Primitive.Contiguous as Contiguous-import Data.Primitive.Contiguous (PrimArray, MutableSliced, Mutable)-import Control.Monad.ST.Strict (ST)-import Data.STRef.Strict (newSTRef, modifySTRef, readSTRef)-import Data.Bits-import Data.Foldable (for_) +-- | @multiply256By128UpperBits aHi aLo b@ is bits 256 to 383 of+-- @(aHi * 2^128 + aLo) * b@. {-# INLINE multiply256By128UpperBits #-} multiply256By128UpperBits :: Word128 -> Word128 -> Word128 -> Word128-multiply256By128UpperBits aHi aLo b =-  let-    -- Break a and b into little-endian 64-bit chunks-    aChunks :: PrimArray Word64-    aChunks = Contiguous.quadrupleton-      (word128Lo64 aLo)-      (word128Hi64 aLo)-      (word128Lo64 aHi)-      (word128Hi64 aHi)-    bChunks :: PrimArray Word64-    bChunks = Contiguous.doubleton-      (word128Lo64 b)-      (word128Hi64 b)--    -- Multiply b by a, one chunk of b at a time-    prod :: PrimArray Word64-    prod = Contiguous.create $ do-      prod' <- Contiguous.replicateMut 6 0-      flip Contiguous.itraverse_ bChunks $ \bIndex bDigit -> do-        pSize <- Contiguous.sizeMut prod'-        multiply256By64Helper-            (Contiguous.sliceMut prod' bIndex (pSize - bIndex))-            aChunks-            bDigit-      pure prod'--  in Word128-    { word128Hi64 = Contiguous.index prod 5-    , word128Lo64 = Contiguous.index prod 4-    }--{-# INLINE multiply256By64Helper #-}-multiply256By64Helper :: forall s. MutableSliced PrimArray s Word64 -> PrimArray Word64 -> Word64 -> ST s ()-multiply256By64Helper _ _ 0 = pure ()-multiply256By64Helper prod a b = do-  carry <- newSTRef 0-  productSize <- Contiguous.sizeMut prod-  let-    aSize = Contiguous.size a-    productLo :: MutableSliced PrimArray s Word64-    productLo = Contiguous.sliceMut prod 0 aSize-    productHi :: MutableSliced PrimArray s Word64-    productHi = Contiguous.sliceMut prod aSize (productSize - aSize)-  -- Multiply each of the digits in a by b, adding them into the 'prod' value.-  -- We don't zero out prod, because we this will be called multiple times, so it probably contains a previous iteration's partial prod, and we're adding + carrying on top of it-  for_ [0..aSize - 1] $ \i -> do-    p <- Contiguous.read productLo i-    let aDigit = Contiguous.index a i-    modifySTRef carry $ \x -> x-      + Word128 0 p-      + (Word128 0 aDigit * Word128 0 b)-    Contiguous.write prod i . word128Lo64 =<< readSTRef carry-    modifySTRef carry (`unsafeShiftR` 64)--  let productHiSize = productSize - aSize-  for_ [0..productHiSize - 1] $ \i -> do-    p <- Contiguous.read productHi i-    modifySTRef carry (+ Word128 0 p)-    Contiguous.write productHi i . word128Lo64 =<< readSTRef carry-    modifySTRef carry (`unsafeShiftR` 64)--  readSTRef carry >>= \case-    0 -> pure ()-    _ -> error "carry overflow during multiplication!"+multiply256By128UpperBits (Word128 a3 a2) (Word128 a1 a0) (Word128 b1 b0) =+  -- a * b0, limbs 0..4+  case mulFull64 a0 b0 of { (# h0, _ #) ->+  case mulFull64 a1 b0 of { (# h1, l1 #) ->+  case mulFull64 a2 b0 of { (# h2, l2 #) ->+  case mulFull64 a3 b0 of { (# h3, l3 #) ->+  case addCarry64 l1 h0 of { (# x1, c1 #) ->+  case adc64 l2 h1 c1 of { (# x2, c2 #) ->+  case adc64 l3 h2 c2 of { (# x3, c3 #) ->+  let !x4 = h3 + c3 in+  -- a * b1, limbs 1..5+  case mulFull64 a0 b1 of { (# g0, y1 #) ->+  case mulFull64 a1 b1 of { (# g1, k1 #) ->+  case mulFull64 a2 b1 of { (# g2, k2 #) ->+  case mulFull64 a3 b1 of { (# g3, k3 #) ->+  case addCarry64 k1 g0 of { (# y2, e2 #) ->+  case adc64 k2 g1 e2 of { (# y3, e3 #) ->+  case adc64 k3 g2 e3 of { (# y4, e4 #) ->+  let !y5 = g3 + e4 in+  -- sum+  case addCarry64 x1 y1 of { (# _, s1 #) ->+  case adc64 x2 y2 s1 of { (# _, s2 #) ->+  case adc64 x3 y3 s2 of { (# _, s3 #) ->+  case adc64 x4 y4 s3 of { (# r4, s4 #) ->+  Word128 (y5 + s4) r4 }}}}}}}}}}}}}}}}}} --- compute prod += a * b-{-# INLINE longMultiply #-}-longMultiply :: forall s. PrimArray Word64 -> Word64 -> Mutable PrimArray s Word64 -> ST s ()+-- | @longMultiply a b prod@ computes @prod += a * b@. The product array must be+-- at least as long as @a@; carries propagate through the rest of it, and+-- running off the end is an error.+longMultiply :: PrimArray Word64 -> Word64 -> MutablePrimArray s Word64 -> ST s ()+longMultiply _ 0 _ = pure () longMultiply a b prod = do-  prod' <- Contiguous.toSliceMut prod-  multiply256By64Helper prod' a b+  plen <- getSizeofMutablePrimArray prod+  let !alen = sizeofPrimArray a+      mulLoop !i !carry+        | i == alen = carryLoop i carry+        | otherwise = do+            p <- readPrimArray prod i+            case mulFull64 (indexPrimArray a i) b of+              (# h, l #) -> case addCarry64 l p of+                (# s1, c1 #) -> case addCarry64 s1 carry of+                  (# s, c2 #) -> do+                    writePrimArray prod i s+                    mulLoop (i + 1) (h + c1 + c2)+      carryLoop !i !carry+        | carry == 0 = pure ()+        | i == plen = error "Atrophy.LongMultiplication.longMultiply: carry overflow"+        | otherwise = do+            p <- readPrimArray prod i+            case addCarry64 p carry of+              (# s, c #) -> do+                writePrimArray prod i s+                carryLoop (i + 1) c+  if plen < alen+    then error "Atrophy.LongMultiplication.longMultiply: product array is too small"+    else mulLoop 0 0

@@ -1,78 +1,219 @@-{-# LANGUAGE TypeOperators #-}--{-# OPTIONS_GHC -Wall #-}-{-# OPTIONS_GHC -fno-warn-orphans #-}-{-# OPTIONS_GHC -Wno-unused-top-binds #-}-{-# LANGUAGE StandaloneDeriving #-}-{-# LANGUAGE GeneralizedNewtypeDeriving #-}-{-# LANGUAGE DerivingVia #-}+{-# LANGUAGE AllowAmbiguousTypes #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE ExplicitNamespaces #-}+{-# OPTIONS_GHC -Wno-orphans #-}  module Main (main) where -import Test.QuickCheck hiding (NonZero)-import Test.Tasty-import Test.Tasty.QuickCheck hiding (NonZero)-import qualified Test.Tasty.QuickCheck as QC (NonZero(..))-import Data.WideWord.Word128-import Data.WideWord.Word256 import Atrophy+import Control.Exception (ErrorCall (..), evaluate, try)+import Control.Monad (forM_, unless)+import Control.Monad.ST (runST)+import Data.Bits+import Data.Primitive.PrimArray+import Data.Proxy (Proxy (..))+import Data.WideWord.Word128 import Data.Word+import GHC.TypeNats (KnownNat, Nat, natVal, type (+), type (-), type (^))+import Test.QuickCheck hiding (NonZero (..))+import Test.Tasty+import Test.Tasty.HUnit+import Test.Tasty.QuickCheck (QuickCheckTests (..), testProperty)  main :: IO ()-main = do-  defaultMain tests+main = defaultMain $ localOption (QuickCheckTests 20000) $ testGroup "atrophy"+  [ testGroup "runtime divisors"+      [ widthTests @Word8 "Word8"+      , widthTests @Word16 "Word16"+      , widthTests @Word32 "Word32"+      , widthTests @Word64 "Word64"+      , widthTests @Word128 "Word128"+      , testCase "Word8 exhaustive" exhaustiveWord8+      , testCase "Word16 all divisors, edge numerators" (allDivisors @Word16)+      , testCase "Word32 edge divisors and numerators" (edgeCases @Word32)+      , testCase "Word64 edge divisors and numerators" (edgeCases @Word64)+      , testCase "Word128 edge divisors and numerators" (edgeCases @Word128)+      ]+  , testGroup "known divisors"+      [ testGroup "Word8"+          [ kd @1 @Word8, kd @2 @Word8, kd @3 @Word8, kd @5 @Word8, kd @7 @Word8, kd @10 @Word8+          , kd @64 @Word8, kd @127 @Word8, kd @128 @Word8, kd @129 @Word8, kd @191 @Word8, kd @255 @Word8+          ]+      , testGroup "Word16"+          [ kd @1 @Word16, kd @3 @Word16, kd @7 @Word16, kd @10 @Word16, kd @255 @Word16, kd @256 @Word16+          , kd @257 @Word16, kd @641 @Word16, kd @32767 @Word16, kd @32768 @Word16, kd @32769 @Word16, kd @65535 @Word16+          ]+      , testGroup "Word32"+          [ kd @1 @Word32, kd @3 @Word32, kd @7 @Word32, kd @10 @Word32, kd @641 @Word32, kd @65537 @Word32+          , kd @(2 ^ 31 - 1) @Word32, kd @(2 ^ 31) @Word32, kd @(2 ^ 31 + 1) @Word32, kd @1000000007 @Word32+          , kd @(2 ^ 32 - 1) @Word32+          ]+      , testGroup "Word64"+          [ kd @1 @Word64, kd @2 @Word64, kd @3 @Word64, kd @5 @Word64, kd @6 @Word64, kd @7 @Word64+          , kd @10 @Word64, kd @11 @Word64, kd @641 @Word64, kd @6700417 @Word64, kd @1000000007 @Word64+          , kd @(2 ^ 32 - 1) @Word64, kd @(2 ^ 32) @Word64, kd @(2 ^ 32 + 1) @Word64+          , kd @1000000000000000000 @Word64, kd @14757395258967641293 @Word64+          , kd @(2 ^ 63 - 1) @Word64, kd @(2 ^ 63) @Word64, kd @(2 ^ 63 + 1) @Word64, kd @(2 ^ 64 - 1) @Word64+          ]+      , testGroup "Word128"+          [ kd @1 @Word128, kd @3 @Word128, kd @7 @Word128, kd @10 @Word128, kd @641 @Word128+          , kd @(2 ^ 64 - 1) @Word128, kd @(2 ^ 64) @Word128, kd @(2 ^ 64 + 1) @Word128+          , kd @10000000000000000000 @Word128, kd @(10 ^ 30) @Word128+          , kd @(2 ^ 127 - 1) @Word128, kd @(2 ^ 127) @Word128, kd @(2 ^ 127 + 1) @Word128+          , kd @(2 ^ 128 - 159) @Word128, kd @(2 ^ 128 - 1) @Word128+          ]+      ]+  , testGroup "known numerators"+      [ testGroup "Word8"+          [ kn @0 @Word8, kn @1 @Word8, kn @2 @Word8, kn @3 @Word8, kn @128 @Word8, kn @200 @Word8, kn @255 @Word8 ]+      , testGroup "Word16"+          [ kn @0 @Word16, kn @1 @Word16, kn @2 @Word16, kn @1000 @Word16, kn @32768 @Word16, kn @65535 @Word16 ]+      , testGroup "Word32"+          [ kn @0 @Word32, kn @1 @Word32, kn @2 @Word32, kn @7 @Word32, kn @65536 @Word32, kn @1000000007 @Word32+          , kn @(2 ^ 31) @Word32, kn @(2 ^ 32 - 1) @Word32+          ]+      , testGroup "Word64"+          [ kn @0 @Word64, kn @1 @Word64, kn @2 @Word64, kn @3 @Word64, kn @1000 @Word64+          , kn @(2 ^ 31) @Word64, kn @(2 ^ 32 - 1) @Word64, kn @(2 ^ 32) @Word64, kn @(2 ^ 32 + 1) @Word64+          , kn @(2 ^ 63) @Word64, kn @(2 ^ 63 + 1) @Word64, kn @14757395258967641293 @Word64, kn @(2 ^ 64 - 1) @Word64+          ]+      , testGroup "Word128"+          [ kn @0 @Word128, kn @1 @Word128, kn @2 @Word128, kn @3 @Word128, kn @1000 @Word128+          , kn @(2 ^ 63) @Word128, kn @(2 ^ 64 - 1) @Word128, kn @(2 ^ 64) @Word128, kn @(2 ^ 64 + 1) @Word128+          , kn @(10 ^ 30) @Word128, kn @(2 ^ 127) @Word128, kn @(2 ^ 128 - 1) @Word128+          ]+      ]+  , testGroup "long division"+      [ testProperty "divRem2By1" $+          forAll nonZeroGen $ \d@(NonZero dv) -> forAll interesting $ \hi' -> forAll interesting $ \lo ->+            let hi = hi' `mod` dv+                (q, r) = (toInteger hi * 2 ^ (64 :: Int) + toInteger lo) `quotRem` toInteger dv+            in divRem2By1 (newDivisor2By1 d) hi lo === (fromInteger q, fromInteger r)+      , testProperty "divisor2By1" $+          forAll nonZeroGen $ \d@(NonZero dv) -> divisor2By1 (newDivisor2By1 d) === dv+      , testProperty "longDivision" $+          forAll (listOf interesting) $ \limbs -> forAll nonZeroGen $ \d@(NonZero dv) ->+            let (qs, r, qs', r') = runST $ do+                  let num = primArrayFromList limbs+                  quotient <- newPrimArray (length limbs + 1)+                  setPrimArray quotient 0 (length limbs + 1) 0xdeadbeef+                  rem1 <- longDivision (newDivisor2By1 d) num quotient+                  q1 <- unsafeFreezePrimArray quotient+                  inPlace <- thawPrimArray num 0 (length limbs)+                  rem2 <- longDivisionInPlace (newDivisor2By1 d) inPlace+                  q2 <- unsafeFreezePrimArray inPlace+                  pure (primArrayToList q1, rem1, primArrayToList q2, rem2)+                (qI, rI) = fromLimbs limbs `quotRem` toInteger dv+            in (fromLimbs (take (length limbs) qs), toInteger r, drop (length limbs) qs, fromLimbs qs', toInteger r')+                 === (qI, rI, [0xdeadbeef], qI, rI)+      ]+  , testGroup "long multiplication"+      [ testProperty "multiply256By128UpperBits" $+          forAll interesting $ \aHi -> forAll interesting $ \aLo -> forAll interesting $ \b ->+            toInteger (multiply256By128UpperBits aHi aLo b)+              === ((toInteger aHi * 2 ^ (128 :: Int) + toInteger aLo) * toInteger b) `shiftR` 256+      , testProperty "longMultiply" $+          forAll (listOf interesting) $ \as -> forAll interesting $ \b -> forAll (vectorOf (length as) interesting) $ \ps ->+            let prod = runST $ do+                  p <- thawPrimArray (primArrayFromList (ps ++ [0, 0])) 0 (length as + 2)+                  longMultiply (primArrayFromList as) b p+                  primArrayToList <$> unsafeFreezePrimArray p+            in fromLimbs prod === fromLimbs ps + fromLimbs as * toInteger b+      , testCase "longMultiply overflow is an error" $ do+          let as = primArrayFromList [maxBound, maxBound :: Word64]+          r <- tryEvaluate $ runST $ do+            p <- newPrimArray 2+            setPrimArray p 0 2 maxBound+            longMultiply as maxBound p+            primArrayToList <$> unsafeFreezePrimArray p+          assertBool "expected an error" (not r)+      ]+  ]+  where+  tryEvaluate :: [Word64] -> IO Bool+  tryEvaluate xs = either (\(ErrorCall _) -> False) (const True) <$> try (evaluate (length xs)) -tests :: TestTree-tests =-  testGroup "Tests" [unitTests]+type Word' a = (Show a, Integral a, Bounded a, FiniteBits a, StrengthReduce a, Show (StrengthReduced a)) -instance Arbitrary Word128 where-  arbitrary = Word128 <$> arbitrary <*> arbitrary+fromLimbs :: [Word64] -> Integer+fromLimbs = foldr (\l acc -> acc * 2 ^ (64 :: Int) + toInteger l) 0 -deriving via (QC.NonZero a) instance (Num a, Eq a, Arbitrary a) => Arbitrary (NonZero a)+ref :: Integral a => a -> NonZero a -> (a, a)+ref n (NonZero d) = case toInteger n `quotRem` toInteger d of+  (q, r) -> (fromInteger q, fromInteger r) --- wrong-naiveMultiply256By128UpperBits :: Word128 -> Word128 -> Word128 -> Word128-naiveMultiply256By128UpperBits aHi aLo b =-  let-    a' = Word256-      { word256hi = word128Hi64 aHi-      , word256m1 = word128Lo64 aHi-      , word256m0 = word128Hi64 aLo-      , word256lo = word128Lo64 aLo-      }-    b' = Word256-      { word256hi = word128Hi64 b-      , word256m1 = word128Lo64 b-      , word256m0 = 0-      , word256lo = 0-      }-    Word256 h l _ _ = a' * b'-  in Word128 h l+interesting :: forall a. (Integral a, Bounded a, FiniteBits a) => Gen a+interesting = frequency+  [ (4, chooseBoundedIntegral (minBound, maxBound))+  , (2, fromInteger <$> choose (0, 1024))+  , (2, do k <- choose (0, bits - 1); o <- choose (-3, 3); pure (fromInteger (bit k + o)))+  , (1, (maxBound -) . fromInteger <$> choose (0, 1024))+  , (2, do k <- choose (1, bits); fromInteger <$> choose (0, bit k - 1))+  ]+  where+  bits = finiteBitSize (0 :: a) -unitTests :: TestTree-unitTests = testGroup "Unit tests"-  [ testGroup "Long multiplication"-      [ -- testProperty "multiply256By128UpperBits" $ equivalentOnArbitrary3 multiply256By128UpperBits naiveMultiply256By128UpperBits-      ]-  , testGroup "Long division"-      [ testGroup "StrengthReducedW64" -        [ testProperty "div64" $ \(a, b) -> ourDiv64 a b === theirDiv a b-        , testProperty "div32" $ \(a, b) -> ourDiv32 a b === theirDiv a b-        ]-      ]+nonZeroGen :: (Integral a, Bounded a, FiniteBits a) => Gen (NonZero a)+nonZeroGen = NonZero <$> interesting `suchThat` (/= 0)++edgeValues :: forall a. (Integral a, Bounded a, FiniteBits a) => [a]+edgeValues =+  [0, 1, 2, 3, 5, 7, 10, maxBound, maxBound - 1, maxBound `div` 2, maxBound `div` 3]+    ++ concat [[bit k - 1, bit k, bit k + 1] | k <- [1 .. finiteBitSize (0 :: a) - 1]]++allChecks :: Word' a => a -> NonZero a -> [(String, (a, a))]+allChecks n d =+  let sr = new d+  in [ ("divRem", divRem n sr)+     , ("divRemConst", divRemConst n sr)+     , ("divRemNonZero", divRemNonZero n d)+     , ("divRemNonZeroConst", divRemNonZeroConst n d)+     , ("div'/rem'", (div' n sr, rem' n sr))+     , ("divConst/remConst", (divConst n sr, remConst n sr))+     , ("divNonZero/remNonZero", (divNonZero n d, remNonZero n d))+     , ("divNonZeroConst/remNonZeroConst", (divNonZeroConst n d, remNonZeroConst n d))+     ]++checkPair :: Word' a => a -> NonZero a -> Assertion+checkPair n d = forM_ (allChecks n d) $ \(name, got) ->+  unless (got == ref n d) $+    assertFailure (name ++ " " ++ show n ++ " " ++ show d ++ ": got " ++ show got ++ ", expected " ++ show (ref n d))++widthTests :: forall a. Word' a => String -> TestTree+widthTests name = testGroup name+  [ testProperty "all variants agree with quotRem" $+      forAll interesting $ \n -> forAll nonZeroGen $ \d ->+        conjoin [counterexample f (got === ref n d) | (f, got) <- allChecks @a n d]+  , testProperty "divisor" $ forAll nonZeroGen $ \d@(NonZero dv) -> divisor (new @a d) === dv   ] -ourDiv64 :: NonZero Word64 -> NonZero Word64 -> Word64-ourDiv64 (NonZero dividend) divi =-  let sr = new64 divi-  in div64 dividend sr+exhaustiveWord8 :: Assertion+exhaustiveWord8 = forM_ [minBound .. maxBound :: Word8] $ \n ->+  forM_ [1 .. maxBound] $ \d -> checkPair n (NonZero d) -ourDiv32 :: NonZero Word32 -> NonZero Word32 -> Word32-ourDiv32 (NonZero dividend) divi =-  let sr = new StrengthReducedW32 divi-  in div' dividend sr+allDivisors :: forall a. Word' a => Assertion+allDivisors = forM_ [1 .. maxBound :: a] $ \d -> forM_ (edgeValues ++ [d - 1, d, d + 1, d * 2, d * 3 - 1]) $ \n ->+  checkPair n (NonZero d) -theirDiv :: Integral a => NonZero a -> NonZero a -> a-theirDiv (NonZero dividend) (NonZero divi) =-  dividend `div` divi+edgeCases :: forall a. Word' a => Assertion+edgeCases = forM_ (filter (/= 0) (edgeValues @a)) $ \d -> forM_ (edgeValues ++ [d - 1, d + 1, d * 2, d * 3 - 1]) $ \n ->+  checkPair n (NonZero d)++kd :: forall (d :: Nat) a. (KnownDivisor d a, KnownNat d, Integral a, Bounded a, FiniteBits a, Show a) => TestTree+kd = testProperty (show dv) $+  conjoin [check n | n <- edgeValues ++ [d - 1, d, d + 1, d * 2 - 1, d * 2]] .&&. forAll interesting check+  where+  dv = natVal (Proxy @d)+  d = fromIntegral dv :: a+  check n = counterexample (show n) $ (divRemK @d n, divK @d n, remK @d n) === (ref n (NonZero d), fst (ref n (NonZero d)), snd (ref n (NonZero d)))++kn :: forall (n :: Nat) a. (KnownNumerator n a, KnownNat n, Word' a) => TestTree+kn = testProperty (show nv) $+  conjoin [check (NonZero d) | d <- filter (/= 0) (edgeValues ++ [n - 1, n, n + 1, n `div` 2, n `div` 3])] .&&. forAll nonZeroGen check+  where+  nv = natVal (Proxy @n)+  n = fromIntegral nv :: a+  check d = counterexample (show d) $+    ( divRemN @n (new d), divN @n (new d), remN @n (new d)+    , divRemNonZeroN @n d, divNonZeroN @n d, remNonZeroN @n d+    ) === (ref n d, fst (ref n d), snd (ref n d), ref n d, fst (ref n d), snd (ref n d))