您遇到的问题是所有中间计算都隐含mod 2
64,而通常事实并非如此。您正在计算的是:
a·b mod m = (a·b mod 264) mod m
我无法想到一种简单的方法来使用仅64位数字进行正确的计算,但您不必一直升级到大整数; 如果
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
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
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
可以给你正确的答案,为什么不直接使用bigint
进行中间计算,最后再转换为long(考虑到m的范围,这是一种无损转换)?我猜测相对于维护自己的数学例程而言,使用bigint
的性能惩罚应该是可以接受的。