Блог пользователя Gerald

Автор Gerald, 14 лет назад, По-русски
Недавно я узнал, что есть способ перемножить два 64 битных числа по модулю третьего 64 битного числа без применения длинной арифметики. Вроде как такое перемножение можно сделать используя long double в 2 строчки кода. Прилично погуглив, я не обнаружил ничего путного на эту тему, хотя польза от такой операции весьма большая (например не прибегать к длинной арифметике в тесте миллера рабина на простоту большого 64 битного числа). Может кто-то знает как это можно сделать?  Заранее благодарен за помощь.
  • Проголосовать: нравится
  • +12
  • Проголосовать: не нравится

14 лет назад, # |
Rev. 2   Проголосовать: нравится +18 Проголосовать: не нравится
В long double действительно очень просто:

long long mul( long long a, long long b, long long m ) {
  long long q = (long long)((long double)a * (long double)b / (long double)m);
  long long r = a * b - q * m;
  return (r + 5 * m) % m;
}
  • 14 лет назад, # ^ |
    Rev. 2   Проголосовать: нравится 0 Проголосовать: не нравится
    Артем, а ты знаешь почему это работает? Кажется, что в long double должно не хватать точности чтобы перемножить 2 таких больших числа без потери информации.

    Да и вообще r = a*b - q*m, a*b - переполнится же?!
    • 14 лет назад, # ^ |
        Проголосовать: нравится +5 Проголосовать: не нравится
      Ну так получилось примерно частное. Если отнять, получится примерно остаток, разница с настоящим остатком будет (погрешность) * m. Код будет работать, если, например, погрешность по модулю не больше 5, и числа таковы, что 11 * m влезает в long long.
      • 14 лет назад, # ^ |
          Проголосовать: нравится 0 Проголосовать: не нравится
        А откуда 11 * m, или это опечатка, должно быть 5 * m?
        • 14 лет назад, # ^ |
            Проголосовать: нравится +3 Проголосовать: не нравится
          Мы предположили, что погрешность может быть +5. Тогда r + 5m + погрешность ~ 11m.
      • 14 лет назад, # ^ |
          Проголосовать: нравится 0 Проголосовать: не нравится
        Если я правильно понимаю, можно заместо (r + 5 * m) % m написать
        r = r % m;
        r = r < 0 ? r + m : r;
        так вроде надежнее будет, так ведь?
  • 14 лет назад, # ^ |
    Rev. 2   Проголосовать: нравится 0 Проголосовать: не нравится
    Тут, наверное, стоит учитывать, что по стандарту С да и С++ размер, а следовательно и точность long double определены просто как >=double. Конечно на современных платформах нормальные компиляторы используют тут 80-хбитную арифметику, но есть и другие. Например, компилятор от microsoft полагает long double = double, по крайней мере они так заявляют. А там кто его знает... И если они не врут, то там перемножение с long double (64-хбитном) будет работать только для относительно небольших чисел.
14 лет назад, # |
Rev. 2   Проголосовать: нравится +46 Проголосовать: не нравится
Общая идея состоит в двух равенствах:
(x+1)*y % mod = (x*y+y)%mod
2*x*y % mod = ((x*y)%mod *2 )%mod
Получаем следующую простую реализацию:
ll mulmod(ll a,ll b,ll m)
{
if (!a || !b) return 0;
if (a&1) return (b+mulmod(a-1,b,m))%mod;
return (mulmod(a/2,b)*2)%mod;

Это работает если 2м влазит в ЛЛ. Сложность  - двоичный логарифм А. И совсем несложно написать нерекурсивную реализацию, если немного подумать.
}

АПД: Прошу прощения, не заметил, что требовалось решение с лонг даблом. Но это тоже работает.
14 лет назад, # |
  Проголосовать: нравится +8 Проголосовать: не нравится
Если 2*x влезает в long long, можно сделать просто за 64 * 2 операции сложения/вычитания/сравнения. Правда, это наверное помедленней, чем хаки с long double, зато надёжно.
  • 14 лет назад, # ^ |
      Проголосовать: нравится 0 Проголосовать: не нравится
    Очень интерестно. То, что я написал выше, работает за такую же сложность, но использует взятие остатка. По всей видимости, это просто сравнение - вычитание М если получившееся число больше М?
    • 14 лет назад, # ^ |
        Проголосовать: нравится +13 Проголосовать: не нравится
      Да, тот же самый способ, но с вычитанием M (понятно, что так будет работать на порядок быстрее). Я просто не успел заметить ваш коммент - написал практически одновременно =)
14 лет назад, # |
  Проголосовать: нравится 0 Проголосовать: не нравится
Вот пример перемножения и возведения в степень по модулю - различие найдите сами (за пример спасибо cmd). Если необходимо, то в powmod можно использовать mulmod вместо умножения по модулю.

long long mulmod(long long a, long long b, long long MOD)
{
    if (b == 0) return 0;
    long long res = mulmod(a, b >> 1, MOD);
    res += res;
    res %= MOD;
    return (b & 1)? (a + res) % MOD : res;
}

long long powmod(long long a, long long b, long long MOD)
{
    if (b == 0) return 1;
    long long res = powmod(a, b >> 1, MOD);
    res *= res;
    res %= MOD;
    return (b & 1)? (a * res) % MOD : res;
}


»
9 лет назад, # |
Rev. 3   Проголосовать: нравится 0 Проголосовать: не нравится

Ignore this please :)

  • »
    »
    9 лет назад, # ^ |
    Rev. 2   Проголосовать: нравится +36 Проголосовать: не нравится

    Guys I don't know what's happening. When I saw that post, I didn't see any comments. But Gassa send me a link of this post. And I saw that there are already more comments. It's sad but true that I'm kinda fooled. Sorry for the inconvenience guys.

    • »
      »
      »
      9 лет назад, # ^ |
        Проголосовать: нравится 0 Проголосовать: не нравится

      That's because of wrong locale. The original post was written in the Russian version of the site, and all people who written comments were writing in Russian locale. But when you got here by the direct link, you've seen the English version of the thread that doesn't any comments in English, so you didn't see them.

      It's pretty often situtation when Russian version of the thread contains more comments than English version, and by default all those comments are hidden in English version of the site.

  • »
    »
    9 лет назад, # ^ |
      Проголосовать: нравится -20 Проголосовать: не нравится

    Из-за таких вот некропостеров я минуты 3 не мог понять, как красный, к тому же достаточно известный не знает достаточно простого алгоритма :D