模数求幂函数中的整数溢出

Question

我正在写一个powMod函数，我必须使用相当密集的函数。其出发点是一个自定义pow功能：

// Compute power using multiplication and square. 
// pow (*) (^2) 1 x n = x^n 
let pow mul sq one x n = 
    let rec loop x' n' acc = 
     match n' with 
     | 0 -> acc 
     | _ -> let q = n'/2 
       let r = n'%2 
       let x2 = sq x' 
       if r = 0 then 
       loop x2 q acc 
       else 
       loop x2 q (mul x' acc) 
    loop x n one 

检查我的输入范围后，我选择了int64因为它足够大代表输出和我能避免昂贵的计算与bigint：

let mulMod m a b = (a*b)%m 
let squareMod m a = mulMod m a a 
let powMod m = pow (mulMod m) (squareMod m) 1L 

我假设模数（m）大于乘数（a,b），函数仅适用于非负数。 对于大多数情况，powMod函数是正确的;然而，问题在于mulMod函数，其中a*b可能超过int64范围，但(a*b)%m不是。 下面的例子演示了溢出问题：

let a = (pown 2L 40) - 1L 
let b = (pown 2L 32) - 1L 
let p = powMod a b 2 // p = -8589934591L -- wrong 

有什么办法避免int64溢出而不诉诸bigint类型？

Answer 1

，你有这样的问题：你的所有中间计算隐式模2 ，它不是通常正确的

A·b模m =（A·b模2 ）模m

这就是你正在计算的。

我想不出一个简单的方法来使用64位数进行正确的计算，但是您不必一路走到bigint;如果a和b最多有64位，则其全部产品具有至多128位，这样你就可以跟踪产品的两个64位整数（此处捆绑为一个自定义的结构）：

// bit width of a uint64, needed for mod calculation 
let width = 
    let rec loop w = function 
    | 0uL -> w 
    | n -> loop (w+1) (n >>> 1) 
    loop 0 

[<Struct; CustomComparison; CustomEquality>] 
type UInt128 = 
    val hi : uint64 
    val lo : uint64 
    new (hi,lo) = { lo = lo; hi = hi } 
    new (lo) = { lo = lo; hi = 0uL } 
    static member (+)(x:UInt128, y:UInt128) = 
     if x.lo > 0xffffffffuL - y.lo then 
      UInt128(x.hi + y.hi + 1uL, x.lo + y.lo) 
     else 
      UInt128(x.hi + y.hi, x.lo + y.lo) 
    static member (-)(x:UInt128, y:UInt128) = 
     if y.lo > x.lo then 
      UInt128(x.hi - y.hi - 1uL, x.lo - y.lo) 
     else 
      UInt128(x.hi - y.hi, x.lo - y.lo) 

    static member (*)(x:UInt128, y:UInt128) = 
     let a1 = ((x.lo &&& 0xffffffffuL) * (y.lo &&& 0xffffffffuL)) >>> 32 
     let a2 = (x.lo &&& 0xffffffffuL) * (y.lo >>> 32) 
     let a3 = (x.lo >>> 32) * (y.lo &&& 0xffffffffuL) 
     let sum = ((a1 + a2 + a3) >>> 32) + (x.lo >>> 32) * (y.lo >>> 32) 
     let sum = 
      if a2 > 0xffffffffffffffffuL - a1 || a1 + a2 > 0xffffffffffffffffuL - a3 then 
       0x100000000uL + sum 
      else 
       sum 
     UInt128(x.hi * y.lo + x.lo * y.hi + sum, x.lo * y.lo) 

    static member (>>>)(x:UInt128, n) = 
     UInt128(x.hi >>> n, x.lo >>> n) 

    static member (<<<)(x:UInt128, n) = 
     UInt128((x.hi <<< n) + (x.lo >>> (64 - n)), x.lo <<< n) 

    interface System.IComparable with 
     member x.CompareTo(y) = 
      match y with 
      | :? UInt128 as y -> 
       match x.hi.CompareTo(y.hi) with 
       | 0 -> x.lo.CompareTo(y.lo) 
       | n -> n 

    override x.Equals(y) = 
     match y with 
     | :? UInt128 as y -> x.hi = y.hi && x.lo = y.lo 
     | _ -> false 

    override x.GetHashCode() = x.hi.GetHashCode() + x.lo.GetHashCode() * 7 

    (* calculate mod via long-division *) 
    static member (%)(x:UInt128, d) = 
     let rec reduce (r:UInt128) d' = 
      if r.hi = 0uL then r.lo % d 
      else 
       let r' = if r < d' then r else r - d' 
       reduce r' (d' >>> 1) 
     let shift = width x.hi + (64 - width d) 
     reduce x (UInt128(0uL,d) <<< shift) 

let mulMod m a b = 
    UInt128(a) * UInt128(b) % m 

(* squareMod, powMod basically as before: *) 
let squareMod m a = mulMod m a a 
let powMod m = pow (mulMod m) (squareMod m) 1uL 

let a = (pown 2uL 40) - 1uL 
let b = (pown 2uL 32) - 1uL 
let p = powMod a b 2 

话虽如此，因为bigint s会给你正确的答案，为什么不使用bigint s做中间计算，并在最后转换为long（这保证是无损转换给定m的范围）？我怀疑对于大多数应用程序来说，使用bigint的性能损失应该可以接受（与维护自己的数学例程相比）。

Answer 2

根据Wikipedia以下公式是等效的。你的代码使用第一个，改变第二个应该解决溢出问题。

c = (a x b) mod(m) 
c = (a x (b mod(m))) mod(m) 

希望这会有所帮助。

根据您的意见如下 - 如果< = m和b < = m，m> sqrt（maxint64），那么我不确定没有更大的存储空间就不可能有解决方案。对于m的大数值，b mod m将返回b，所以使用上面的等价公式不会有好处。

好消息是，您应该能够限制对单行的更改并将值重新输入回64位[因为我们知道（a * b）％c应该不会溢出]，然后继续计算。这将开销（在执行性能方面）限制为尽可能小的一段代码。

Answer 3

我不知道F＃是很好，但类似的信息（伪）

pmod a b n // a^b % n 
pmod a 0 n = 1 
pmod a 1 n = a%n 
pmod a b n = match b%2 
    | 0 -> ((pmod (a) (b/2) n)^2) % n 
    | 1 -> ((pmod (a) (b-1) n) * a) % n 

PMOD是不正常不过也应作为辅助功能，使

PowMod a b n = pmod (a%n) b n 

你可以看到这个结果将错误结果的平方将磨碎uint64所以n必须在uint32

Answer 4

我对f＃几乎一无所知，但我认为你可以应用的事实：

若b奇数和n使得b = 2n + 1

a * b mod(m) = 2 * a * n + a mod(m) 
      = 2 * (a*n mod(m)) + a mod(m) 

并且类似地，如果b为偶数。您可以根据需要在a或n上多次重复此操作，直到您获得可以放入int64的产品。如果m> maxint64/2，我认为仍然有可能发生溢出。

Answer 5

不回答你的问题，但是写的函数这种方式将使其更通用，方便使用，它也似乎更高效：

let inline pow x n = 
    let zero = LanguagePrimitives.GenericZero 
    let rec loop x acc = function 
     | n when n = zero -> acc 
     | n -> 
      let q = n >>> 1 
      let acc = if n = (q <<< 1) then acc else x * acc 
      loop (x * x) acc q 
    loop x LanguagePrimitives.GenericOne n;; 

for x = 0 to 1000000 do 
    pow 3UL 31UL |> ignore 

而且，我想无符号长亿韩元”够了吗？

编辑：下面的算法是因为执行少的乘法比一个快3倍以上的大BIGINT - 可以帮助你赞成BIGINT的选择：

let inline pow2 x n = 
    let zero = LanguagePrimitives.GenericZero 
    let one = LanguagePrimitives.GenericOne 
    let rec loop x data = function 
     | c when c <<< 1 <= n -> 
      let c = c <<< 1 
      let x = x * x 
      loop x (Map.add -c x data) c 
     | c -> reduce x data (n - c) 
    and reduce acc data = function 
     | c when c = zero -> acc 
     | c -> 
      let next, value = data |> Seq.pick (fun (KeyValue (n, v)) -> if -n <= c then Some (-n, v) else None) 
      reduce (acc * value) data (c - next) 
    loop x (Map [-1, x]) one;; 

for x = 10000 downto 9000 do 
    pow2 7I x |> ignore