累乗剰余演算関数の整数オーバーフロー

Question

私は非常に集中的に使用しなければならないpowMod関数を書いています。

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 

：出発点は、カスタム pow関数である。出力を表現するのに十分な大きさと私は bigintで高価な計算を避けることができるので、

// 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を選びました

モジュロ（m）が乗数（a,b）よりも大きく、関数が負でない数値だけで動作すると仮定します。 ほとんどの場合、powModの機能は正しいです。ただし、a*bがint64の範囲を超えている可能性がありますが、(a*b)%mはではないmulMod機能に問題があります。たとえば、以下の ザ・オーバーフローの問題を示しています

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

はbigintタイプに頼ることなくint64オーバーフローを回避するための方法はありますか？

Answer 1

あなたは持っている問題はあなたのものすべてであります中間計算は、暗黙的にMOD 2 であり、それは一般・B MOD M =（・MOD 2 B）MOD M

真

ことはありません

これはあなたが計算しているものです。

ちょうど64ビットの数字を使って正しい計算を行う簡単な方法は考えられませんが、bigintsまでずっと進む必要はありません。 aとbは最大64ビットを持っている場合は、その完全な製品は、ほとんどの128ビットであり、あなたは、2つの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 

あなたが中間計算を行い、最終的にlongに変換するのはどうしてbigintを使うのではなく、bigintが正しい答えを与えるからです（mの範囲で無損失変換であることが保証されています）？私は、bigintsを使用することによるパフォーマンスの低下は、ほとんどのアプリケーションで受け入れられるべきだと思っています（自分の数学ルーチンを維持するという頭痛に比べて）。

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を返します。したがって、上記の等価式を使用すると良い結果は得られません。

良いニュースは、変更を1行に制限し、[a * b）％cがオーバーフローしないことを知っているので、64ビットに値を再入力する必要があります。計算。これにより、（実行パフォーマンスに関して）費用の部分を可能な限り小さく抑えることができます。

Answer 3

私はよくなく、（擬似コード）のようなもの

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がOKでない、まだしかし、あなたは、この結果を見ることができるので、

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

ヘルパー関数として使用されるべきであるF＃のを知りません結果が正方形になるとすぐに間違ってしまうので、uint64なので、nはuint32に収まる必要があります。

Answer 4

私はほとんど何も知りませんが、事実を適用することができます。

bが偶数の場合、Bの奇数とnようb = 2n + 1

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

と同様に、もし。あなたは明らかに、int64に適合する製品になるまで、aまたはnに必要な回数だけこれを繰り返すことができます。私は、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 

が、私は「unsigned long型ウォンを想定十分でしょうか？

編集：以下のアルゴリズムが少ない乗算を行うため、大BIGINTに速く上記1よりも3倍される - あなたは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