使用Haskell类型族或GADTs进行模数算术?

9
我经常需要在Haskell中执行模算术运算,其中模数通常很大且经常是质数(例如2000000011)。目前,我只使用像(modAdd m a b),(modMul m a b),(modDiv m a b)等函数。但这相当不方便,需要额外指定并携带一个附加参数,并在常规整数形式和单独的mod-形式中创建我的各种函数。
因此,创建一个新类似于以下内容的新类可能是一个好主意:
class Integral a => Mod a

m = 2000000011

instance Integral a => Num (Mod a) where
   ... defining (+), etc. as modulo m

然后可以使用常规函数执行常规算术,并定义有用的结构,例如

factorials :: [Mod Int]
factorials = 1:zipWith (*) factorials [1..]

但是这里有个问题:所有类型为 Mod Int 的值必须有相同的模数。然而,通常我需要在一个程序中使用多个模数(当然,始终只组合相同模数的值)。
我认为可以通过类似下面的方法解决这个问题,但我不理解得足够好以确保这一点:
class Integral a => Mod Nat a

其中Nat是一种以Peano方式编码模数的类型。这将是有利的:我可以拥有不同模数的值,类型检查器将防止我意外组合此值。

像这样的东西是否可行和有效?如果我尝试使用该模数,编译器或RTS是否会尝试实际构造巨大的(Succ(Succ(Succ ...重复2000000011次)),使解决方案实际上无用? RTS是否会在每个操作上尝试检查类型匹配?每个值的RTS表示是否会从本来可以只是未装箱的int被放大?

有更好的方法吗?

结论

感谢 cirdecdfeueruser5402tikhon-jelvis 提供的有用评论,我意识到(并不出乎意料)我不是第一个想到这个想法的人。特别地,Kiselyov和Shan最近发表了一篇论文,给出了一种实现方法,而tikhon-jelvis已经在Hackage上发布了一个名为(惊喜!)模数算术的解决方案,使用高级的ghc pragma提供了更好的语义。

未解决问题(对我来说)

我需要成为一个有用的助手来翻译文本。
幕后会发生什么?特别是,一个包含一百万个[Mod Int 2000000011]元素的列表是否会携带额外的一百万个2000000011副本?还是它编译成与单独携带模参数的一百万个Int列表相同的代码?后者更好。
性能补充说明
我在解决当前问题时运行了一些基准测试。第一次运行使用未装箱的10,000元素Int向量,并对其执行了10,000次操作:
   4,810,589,520 bytes allocated in the heap
         107,496 bytes copied during GC
       1,197,320 bytes maximum residency (1454 sample(s))
         734,960 bytes maximum slop
              10 MB total memory in use (0 MB lost due to fragmentation)

                                     Tot time (elapsed)  Avg pause  Max pause
  Gen  0      6905 colls,     0 par    0.109s   0.101s     0.0000s    0.0006s
  Gen  1      1454 colls,     0 par    0.812s   0.914s     0.0006s    0.0019s

  TASKS: 13 (1 bound, 12 peak workers (12 total), using -N11)

  SPARKS: 0 (0 converted, 0 overflowed, 0 dud, 0 GC'd, 0 fizzled)

  INIT    time    0.000s  (  0.001s elapsed)
  MUT     time    2.672s  (  2.597s elapsed)
  GC      time    0.922s  (  1.015s elapsed)
  EXIT    time    0.000s  (  0.001s elapsed)
  Total   time    3.594s  (  3.614s elapsed)

  Alloc rate    1,800,454,557 bytes per MUT second

  Productivity  74.3% of total user, 73.9% of total elapsed

第二次运行,我在一个未装箱的大小为10,000的向量(Mod Int 1000000007)上执行了相同的操作。这使得我的代码变得更简单,但时间大约需要3倍长(同时具有几乎相同的内存配置文件):

   4,810,911,824 bytes allocated in the heap
         107,128 bytes copied during GC
       1,199,408 bytes maximum residency (1453 sample(s))
         736,928 bytes maximum slop
              10 MB total memory in use (0 MB lost due to fragmentation)

                                     Tot time (elapsed)  Avg pause  Max pause
  Gen  0      6906 colls,     0 par    0.094s   0.107s     0.0000s    0.0007s
  Gen  1      1453 colls,     0 par    1.516s   1.750s     0.0012s    0.0035s

  TASKS: 13 (1 bound, 12 peak workers (12 total), using -N11)

  SPARKS: 0 (0 converted, 0 overflowed, 0 dud, 0 GC'd, 0 fizzled)

  INIT    time    0.000s  (  0.001s elapsed)
  MUT     time    8.562s  (  8.323s elapsed)
  GC      time    1.609s  (  1.857s elapsed)
  EXIT    time    0.000s  (  0.001s elapsed)
  Total   time   10.172s  ( 10.183s elapsed)

  Alloc rate    561,858,315 bytes per MUT second

  Productivity  84.2% of total user, 84.1% of total elapsed

我想知道为什么会发生这种情况,以及是否可以修复。尽管如此,我非常喜欢模数算术包,并将在性能不是绝对关键的情况下使用它。


7
你可能会对这篇论文Functional Pearl: Implicit Configuration以及相应的reflection包感兴趣。 - Cirdec
2
请注意,在“反射”中与“给定”有关的一切都被认为是不明智的;其余部分真的很酷。 - dfeuer
我现在正在阅读这篇论文。非常有趣。感谢Cirdec和dfeuer。这可能正是我一直在寻找的东西。我仍然困惑的主要问题是是否可能构建Kiselyov/Shan模数,而不考虑它们作为废弃参数的函数的奇怪性。 - CarlEdman
2个回答

4

较新版本的GHC内置了类型级数,这应该比使用Peano算术自己编写的更有效。您可以通过启用DataKinds来使用它们。作为奖励,您还将获得一些漂亮的语法:

factorials :: [Mod Int 20]

无论这是否有效取决于您如何实现 Mod 类型。在最一般的情况下,您可能希望在每次算术运算后进行 mod 。除非您在一个需要节省一些指令的热循环中,否则这应该没问题。(而在热循环内部,明确何时mod可能更好。)
我实际上在Hackage上的一个库中实现了这种类型:modular-arithmetic。它具有测试套件,但没有基准,因此我无法保证绝对性能,但它不会执行任何慢操作,对于我的目的来说足够快。(需要注意的是,涉及到小模数。)如果您尝试并遇到性能问题,我很乐意听取反馈以便我可以尝试解决问题。

1
非常好。看起来你已经做了我想做的事情。一个问题:为什么要定义一个新的“inv”函数,而不是实例化Fractional并使用recip?这也给你提供了fromRational和/,对于模数来说通常都很有意义。 - CarlEdman
1
@CarlEdman:啊,说得对。我只是没想到这个问题,但这似乎是有道理的,尽管可能有些行为不太直观?(fromRational感觉有点奇怪...)你可以打开一个GitHub issue 或发起一个 pull request 吗?谢谢! - Tikhon Jelvis

4

以下是使用 Data.Reflection 的一些有效代码:

{-# LANGUAGE Rank2Types #-}
{-# LANGUAGE FlexibleContexts #-}

import Data.Reflection
import Data.Proxy

data M a s = M a -- Note the phantom comes *after* the concrete

-- In `normalize` we're tying the knot to get the phantom types to align
-- note that reflect :: Reifies s a => forall proxy. proxy s -> a

normalize :: (Reifies s a, Integral a) => a -> M a s
normalize a = b where b = M (mod a (reflect b)) 

instance (Reifies s a, Integral a) => Num (M a s) where
  M a + M b = normalize (a + b)
  M a - M b = normalize (a - b)
  M a * M b = normalize (a * b)
  fromInteger n = normalize (fromInteger n)
  abs _     = error "abs not implemented"
  signum _  = error "sgn not implemented"

withModulus :: Integral a => a -> (forall s. Reifies s a => M a s) -> a
withModulus m ma = reify m (runM . asProxyOf ma)
  where asProxyOf :: f s -> Proxy s -> f s
        asProxyOf a _ = a

runM :: M a s -> a
runM (M a) = a

example :: (Reifies s a, Integral a) => M a s
example = normalize 3

example2 :: (Reifies s a, Integral a, Num (M a s)) => M a s
example2 = 3*3 + 5*5

mfactorial :: (Reifies s a, Integral a, Num (M a s)) => Int -> M a s
mfactorial n = product $ map fromIntegral [1..n]

test1 p n = withModulus p $ mfactorial n

madd :: (Reifies s Int, Num (M Int s)) => M Int s -> M Int s -> M Int s
madd a b = a + b

test2 :: Int -> Int -> Int -> Int
test2 p a b = withModulus p $ madd (fromIntegral a) (fromIntegral b)

1
为避免混淆,我强烈建议在使用启用forall关键字的扩展时使用ScopedTypeVariables。这里没有任何区别,但我不喜欢考虑它是开启还是关闭状态。 - dfeuer
1
实际上,在这里使用 ScopedTypeVariables 将非常有用,以避免需要打结。只需使用 normalize :: forall s a . (Reifies s a, Integral a) => a -> M a s 然后 normalize a = M (mod a (reflect (Proxy :: Proxy s))) - dfeuer
完全同意,当我第一次了解到类型变量不是有作用域时,我感到震惊。ScopedTypeVariables确实属于Haskell'。 - CarlEdman

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接