Кончилась длина long long в с++

Question

#include<bits/stdc++.h>
using namespace std;
long long fac_1(int k1){
if(k1==0){
    return 1;
}
else{
    return k1*fac_1(k1-1);
}
}
int main(){
long long n,k,a,b;
cin >> n>>k;
a=n+k-1;
b=n+k-2;
for(int y=1;n-1>y;y++){
    a=a*b;
    b--;
}
cout <<(a/fac_1(n-1));
}

Есть такой код, суть его в том чтобы определить сколько способами можно положить k предметов в n коробок - C(n+k-1,k). Проблема в том, что n+k<=60 и при значениях близким к 30 long long кончается, можно ли это как-то пофиксить?

Answer 1

Вы пытаетесь посчитать выражение C(m,n) = n! / ((n-m)! * m!) Очевидно, что факториал растет быстро и напрямую немного плохо.

Но давайте посмотрим на эту часть n!/m! (у нас n>m). Если расписать, то это (1*2*...*(n-1)*n) / (1*2*...*(m-1)*m) - а это просто (m+1) * (m+2) * ... * n. Если m и n близкие числа, то это легко.

Теперь получается вот так ((m+1) * (m+2) * ... * n) / (1*2*....*(n-m)). Обычно это уже можно посчитать и для 99% олимпиадных задач достаточно. И это легко закодировать.

long long int C(int n, int m)
{
    long long int a = 1;
    for (int i = m+1; i<=n; i++) a = a*i;
    long long int b = 1;
    for (int i = 1; i <= (n-m); i++) b = b*i;
    return a/b;
}

C(35,25)=183579396
C(40,30)=847660528

но вот C(50,30) оно уже не хочет...

И тут упираемся, то можно сделать очередной трюк - просто множить "верхние числа" и потихоньку делить "на подходящее нижнее". Под подходящим имеется ввиду такое, на которое делится без остатка.

Но если пойти ещё дальше, то можно разложить на множители верхнее и нижнее произведение и потом поудалять немного. Хотя перед тем как удалять, можно просто сверить два списка, там вполне могут быть уже одинаковые числа.

я сел и накодил такое решение. Оно немного странное, но зато рабочее:)

#include <set>
#include <algorithm>
#include <numeric>

// раскладывает на множители. Да, она супер неэффективная
// но у нас числа до 100 и этого более чем достаточно.
// для 10 она вернет { 2 5 }
// для 12 => { 2 2 3 }
// но можно взять и просто руками захардкодить
std::multiset<int> divid(int n)
{
    std::multiset<int> r;
    for (int i = 2; i < n; i++) {
        while (n%i == 0) {
            n = n / i;
            r.insert(i);
        }
    }
    if (n != 1) r.insert(n);
    return r;
}

// эта функция получает два множества и просто вычитает одно с другого.
// тут куча странных функций, но лучше читать документацию:)
// да, вместо мультисета можно было бы использовать map<int,int>
// ключ - делитель, значение - кол-во, но это на любителя
unsigned long long int C1(int n, int m)
{
    std::multiset<int> a;
    for (int i = m+1; i<=n; i++) { auto x = divid(i); a.insert(x.begin(), x.end());}
    std::multiset<int> b;
    for (int i = 1; i <= (n-m); i++) { auto x = divid(i); b.insert(x.begin(), x.end());}
    std::multiset<int> r;
    std::set_difference(a.begin(),a.end(), b.begin(), b.end(), std::inserter(r,r.begin()));
    return std::accumulate(r.begin(), r.end(), 1ull, std::multiplies<unsigned long long int>());  
}

C1(50,30) = 47129212243960
C1(60, 20) = 4191844505805495

Answer 2

Воспользуйтесь соотношениями:

C(k, k) = 1
C(n, k) = n * C(n - 1, k) / (n - k)

Получается что какое количество дробей ни возьми произведение всегда целое:

k + 1   k + 2   k + 3         k + i
----- * ----- * ----- * ... * ----- * ...
  1       2       3             i

Что приводит к такой процедуре:

unsigned long long cnk(unsigned long long n, unsigned long long k) {
    unsigned long long p = 1;
    for (unsigned long long j = k + 1; j <= n; ++j) {
        p = p * j / (j - k);
    }
    return p;
}

Переполнения не будет пока n <= 62. Это решение можно улучшить до оптимального: если результат помещается в разрядную сетку мы его вычислим.

Единственная проблема кода выше - это переполнение произведения p * j. Избавимся от него. Пусть известно что a * b / c - целое, тогда

a * b / c = (a / НОД(a, c)) * (b / (c / НОД(a, c)))

В правой части все частные целые и нет переполнения (мы предполагаем что результат левой части в тип помещается).

Получается такое решение:

unsigned long long gcd(unsigned long long a, unsigned long long b) {
    while (b > 0) {
        unsigned long long c = a % b;
        a = b;
        b = c;
    }
    return a;
}

// a * b / c
unsigned long long mul_div(unsigned long long a, unsigned long long b, unsigned long long c) {
    unsigned long long ac = gcd(a, c);
    return (a / ac) * (b / (c / ac));
}

unsigned long long cnk(unsigned long long n, unsigned long long k) {
    unsigned long long p = 1;
    for (unsigned long long j = k + 1; j <= n; ++j) {
        p = mul_div(p, j, j - k);
    }
    return p;
}

Оно может вычислить Cnk для любых комбинаций n и k для которых результат помещается в разрядную сетку. Например любые n <= 67:

cnk(67, 33) = 14226520737620288370

Последняя вещь которую нужно сделать - уменьшить число итераций если n очень велико, а k мало. В таком случае будем заменять k на n - k:

unsigned long long cnk(unsigned long long n, unsigned long long k) {
    if (k < n - k) {
        k = n - k;
    }
    unsigned long long p = 1;
    for (unsigned long long j = k + 1; j <= n; ++j) {
        p = mul_div(p, j, j - k);
    }
    return p;
}

Максимальные n для небольших k:

cnk(18446744073709551614,  0) =                    1
cnk(18446744073709551614,  1) = 18446744073709551614
cnk(          6074001000,  2) = 18446744070963499500
cnk(             4801280,  3) = 18446738006366306560
cnk(              145056,  4) = 18446483332847246040
cnk(               18580,  5) = 18442234518422931216
cnk(                4868,  6) = 18426048264474788224
cnk(                1913,  7) = 18399302838933135756
cnk(                 967,  8) = 18419736117819661560
cnk(                 577,  9) = 18348006354228436600
cnk(                 386, 10) = 17991165343481265936

Есть пара сочетаний n и k для которых программа не работает, хотя результат помещается в разрядную сетку. Это cnk(2^64 - 1, 1) = 2^64 - 1 и cnk(2^64 - 1, 0) = 1. Для них счётчик цика переполняется и цикл не завершается. Нет в мире совершенства. :)

БЛОГ НА HUSL

Кончилась длина long long в с++

Ответы (2 шт):