14

たとえば、BigInteger ライブラリを使用せずにビルドする方法はあり(853467 * 21660421200929) % 100000000000007ますか (各数値は 64 ビット整数に収まりますが、乗算結果は収まらないことに注意してください)。

この解決策は非効率的です：

int64_t mulmod(int64_t a, int64_t b, int64_t m) {
    if (b < a)
        std::swap(a, b);
    int64_t res = 0;
    for (int64_t i = 0; i < a; i++) {
        res += b;
        res %= m;
    }
    return res;
}
4

7 に答える 7

23

ロシアの農民の乗算を使用する必要があります。倍増を繰り返してすべての値を計算し、 のth ビットが設定さ(b*2^i)%mれている場合はそれらを追加します。ia

uint64_t mulmod(uint64_t a, uint64_t b, uint64_t m) {
    int64_t res = 0;
    while (a != 0) {
        if (a & 1) res = (res + b) % m;
        a >>= 1;
        b = (b << 1) % m;
    }
    return res;
}

O(log(a))時間ではなく時間がかかるため、アルゴリズムが改善されますO(a)。

警告: 符号なしで、m63 ビット以下の場合にのみ機能します。

于 2012-08-29T04:42:02.077 に答える
20

Keith Randall の答えmは良いですが、彼が言ったように、 63 ビット以下の場合にのみ機能するという警告があります。

以下に、2 つの利点がある変更を示します。

  1. m64ビットでも動作します。
  2. 一部のプロセッサではコストがかかる可能性があるモジュロ演算を使用する必要はありません。

res -= m(および行は、期待される結果を得るために 64 ビットの符号なし整数オーバーフローに依存していることに注意してくださいtemp_b -= m。符号なし整数オーバーフローは C および C++ で明確に定義されているため、これで問題ありません。このため、符号なし整数型を使用することが重要です。 )

uint64_t mulmod(uint64_t a, uint64_t b, uint64_t m) {
    uint64_t res = 0;
    uint64_t temp_b;

    /* Only needed if b may be >= m */
    if (b >= m) {
        if (m > UINT64_MAX / 2u)
            b -= m;
        else
            b %= m;
    }

    while (a != 0) {
        if (a & 1) {
            /* Add b to res, modulo m, without overflow */
            if (b >= m - res) /* Equiv to if (res + b >= m), without overflow */
                res -= m;
            res += b;
        }
        a >>= 1;

        /* Double b, modulo m */
        temp_b = b;
        if (b >= m - b)       /* Equiv to if (2 * b >= m), without overflow */
            temp_b -= m;
        b += temp_b;
    }
    return res;
}
于 2013-09-08T04:02:27.767 に答える
5

どちらの方法でもうまくいきます。最初のものはあなたのものと同じですが、私はあなたの番号をULLを明示するために変更しました。2つ目は、アセンブラ表記を使用します。これは、より高速に動作するはずです。すでに述べたモンゴメリリダクションのように、暗号化で使用されるアルゴリズム（RSAおよびRSAベースの暗号化がほとんどだと思います）もありますが、それらの実装には時間がかかると思います。

#include <algorithm>
#include <iostream>

__uint64_t mulmod1(__uint64_t a, __uint64_t b, __uint64_t m) {
  if (b < a)
    std::swap(a, b);
  __uint64_t res = 0;
  for (__uint64_t i = 0; i < a; i++) {
    res += b;
    res %= m;
  }
  return res;
}

__uint64_t mulmod2(__uint64_t a, __uint64_t b, __uint64_t m) {
  __uint64_t r;
  __asm__
  ( "mulq %2\n\t"
      "divq %3"
      : "=&d" (r), "+%a" (a)
      : "rm" (b), "rm" (m)
      : "cc"
  );
  return r;
}

int main() {
  using namespace std;
  __uint64_t a = 853467ULL;
  __uint64_t b = 21660421200929ULL;
  __uint64_t c = 100000000000007ULL;

  cout << mulmod1(a, b, c) << endl;
  cout << mulmod2(a, b, c) << endl;
  return 0;
}
于 2012-08-28T22:45:13.807 に答える
4

反復倍増アルゴリズムの改良点は、オーバーフローなしで一度に計算できるビット数をチェックすることです。両方の引数に対して早期終了チェックを行うことができます - N が素数でないという (ありそうもない?) イベントを高速化します。

たとえば、100000000000007 == 0x00005af3107a4007 の場合、各反復ごとに 16 (または 17) ビットを計算できます。この例では、実際の反復回数は 3 回です。

// just a conceptual routine
int get_leading_zeroes(uint64_t n)
{
   int a=0;
   while ((n & 0x8000000000000000) == 0) { a++; n<<=1; }
   return a;
}

uint64_t mulmod(uint64_t a, uint64_t b, uint64_t n)
{
     uint64_t result = 0;
     int N = get_leading_zeroes(n);
     uint64_t mask = (1<<N) - 1;
     a %= n;
     b %= n;  // Make sure all values are originally in the proper range?
     // n is not necessarily a prime -- so both a & b can end up being zero
     while (a>0 && b>0)
     {
         result = (result + (b & mask) * a) % n;  // no overflow
         b>>=N;
         a = (a << N) % n;
     }
     return result;
}
于 2014-02-20T08:08:21.140 に答える
3

乗算を加算に分割する何かを試すことができます。

// compute (a * b) % m:

unsigned int multmod(unsigned int a, unsigned int b, unsigned int m)
{
    unsigned int result = 0;

    a %= m;
    b %= m;

    while (b)
    {
        if (b % 2 != 0)
        {
            result = (result + a) % m;
        }

        a = (a * 2) % m;
        b /= 2;
    }

    return result;
}
于 2012-08-28T22:30:55.923 に答える
1

アルゴリズムの改善を提案できます。

実際には、毎回加算し、各反復後にモジュロを実行するa * bことにより、反復的に計算します。b毎回追加することをお勧めしますがb * x、オーバーフローしないxように決定されます。b * x

int64_t mulmod(int64_t a, int64_t b, int64_t m)
{
    a %= m;
    b %= m;

    int64_t x = 1;
    int64_t bx = b;

    while (x < a)
    {
        int64_t bb = bx * 2;
        if (bb <= bx)
            break; // overflow

        x *= 2;
        bx = bb;
    }

    int64_t ans = 0;

    for (; x < a; a -= x)
        ans = (ans + bx) % m;

    return (ans + a*b) % m;
}
于 2012-08-28T22:46:18.270 に答える