Есть на LeetCode задачка. Medium уровня, на динамическое программирование, ничего особенного. Однако, если присмотреться внимательнее, она окажется интереснее, чем на первый взгляд. Кроме того, можно получить более быстрое решение, чем "официальное".

Условие (в адаптации для постсоветского читателя)

Недалеко от вокзала, есть маленькая столовая. Каждое утро в ней варят две одинаковых кастрюли супа: один с говядиной , другой, собственно, с котом . Порция супа составляет миллилитров.

• Для особо уважаемых клиентов наливают только суп с говядиной ( мл говядина);

• Для просто уважаемых на четверть разбавляют котом ( мл говядина / мл кот);

• Для обычных с улицы - пятьдесят на пятьдесят (мл говядина / мл кот);

• Для местных алкашей подают, в основном, кота, добавив чуть-чуть говядины для аромата (мл говядина / мл кот).

Клиенты каждой категории приходят в случайном порядке, каждая категория с вероятностью %.

Если какого-то вида супа для текущего клиента не хватило - ему подают сколько есть, не доливая ни того, ни другого. При этом заказ считается исполненным.

При заданном объёме кастрюли мл требуется определить, какова вероятность того, что говядина закончится раньше, и припозднившиеся гости будут жрать кушать только и исключительно несчастного кота. Кроме того, к ответу надо добавить вероятности того, что оба супа закончатся одновременно. Результат требуется получить с погрешностью не более .

Объем кастрюли может принимать значения от до .

Наблюдения и размышления

Прежде всего, заметим, что объем половника у них мл, поэтому, вместо того, чтобы считать в миллилитрах, будем считать в половниках. Тогда объём кастрюли будет половников. Если не кратно , то округляем в большую сторону — как указано в условии, остатки считаются за полную порцию.

В начальный момент у нас имеется половников говядины и половников кота.
После первого заказа может остаться:

1. , если гость особо уважаемый;

2. , если гость просто уважаемый;

3. , если гость обычный;

4. , если гость алкаш.

Или, после -го заказа:

1. , если гость особо уважаемый;

2. , если гость просто уважаемый;

3. , если гость обычный;

4. , если гость алкаш.

Тут мы, внезапно, замечаем, что общее количество супа с каждой порцией уменьшается ровно на половника, а разность увеличивается либо на 4, либо на 2, либо не меняется, либо уменьшается на 2. Учитывая, что с утра , она всегда чётная, как и сумма, которая в начальный момент .

Введём новые переменные: и .

В начальный момент .

После -го заказа эти величины меняются следующим образом:

Рассмотрим вероятность того, что после заказа разность равна . Если в этот момент , значит:

1. либо перед этим было , но пришёл особо уважаемый гость и разность изменилась на ;

2. либо было , но пришёл просто уважаемый гость и разность изменилась на ;

3. либо было , пришёл обычный гость и разность не изменилась;

4. либо было , но пришёл алкаш и разность изменилась на .

Все эти события имеют место с вероятностью . Тогда:

Таким образом, если мы создадим массив , зададим ,то сможем последовательно посчитать вероятности после каждого заказа.

Но в какой-то момент один (или оба) супа закончатся. Когда это произойдёт?

Возвращаясь от наших переменных обратно к количествам супов, получим:

Количество супа с говядиной становится равным при , а супа с котом при . То есть, как только будет достигнуто состояние, в котором , мы должны прекратить сервировку. Это значит, что соответствующие им вероятности , где , не внесут вклада в вероятности последующих событий и их надо обнулить. Если при этом первым закончился суп с говядиной, то есть , то мы должны сначала добавить эту вероятность к окончательному ответу.

А когда закончатся оба супа, тогда их полусумма станет равна . Таким образом, максимальное количество порций равно , при чем последняя может быть неполной.

Рассмотрим, например, .

Возможные состояния супов на каждом шаге отмечены кружочечком.

Границы серой области показывают уменьшение на каждом шаге. Если попало в серую область, значит после данного заказа какой-то из супов закончился. Если закончился суп с говядиной — кружок зелёный, если закончился суп с котом — красный, если оба — двухцветный.

В голубых кружочках указана искомая вероятность того, что суп с говядиной закончится раньше. Для каждого шага она равна сумме всех зелёных значений, лежащих ниже или на том же уровне (двухцветный входит с коэффициентом ).

Стрелками показано, какое состояние может получиться на следующем шаге (чтобы не загромождать, показаны только крайние случаи).

Отметим, что состояния, где один из супов закончился, на последующие не влияют, но если бы мы продолжали их строить, то они давали бы вклад только в серой области (отмечено прозрачной зеленой и красной заливкой) — если суп закончился, то обратно он не появится. Запомним этот факт, он нам понадобится позже.

Рассмотрим также случай для чётного .

Мы видим, что для нечетного и следующего за ним четного количество раздач не меняется.

Алгоритм №1. Динамическое программирование

Таким образом стандартное решение данной задачи в неоптимизированном случае выглядит следующим образом.

1. Находим начальное

2. Создаем массив вероятностей , состоящий из строк, столбцов. Обозначим индекс элемента соответствующего . Пусть индексация начинается с 0, тогда . Строка соответствует началу дня, когда порций съедено, строка — когда оба супа кончились;

3. Задаём , все остальные элементы в этой строке равны ;

4. Инициализируем переменную, в которой будем хранить результат ;

5. Инициализируем переменную, в которой будем хранить текущее значение полусуммы ;

6. Запускаем цикл по строкам от до включительно;

7. Вычисляем на текущем шаге: . Если получилась меньше , то присваиваем (любой остаток считается полноценной порцией);

8. Запускаем цикл по всем столбцам от до ;

9. Для каждого элемента строки считаем вероятность:

   Если индекс переполняется, соответствующее слагаемое не учитываем;

10. Если оба супа закончились одновременно , прибавляем половину к ;

11. Иначе если только суп с котом закончился , зануляем его ;

12. Иначе если закончился суп с говядиной , прибавляем полученную вероятность к ответу и зануляем ;

13. Переходим к следующему элементу и повторяем 9-13;

14. Переходим к следующей строке и повторяем 9-13;

15. Ответ находится в .

Вот программка на C++. Разумеется, можно сделать красиво, можно заметить, что массив избыточен, что нам надо хранить только предыдущую и следующую строчки и т.п.

double soupServings(int n) {
        int s0 = (n-1)/25+1;
        
        vector<vector<double>> P((s0+1)/2+1, vector<double>(2*s0+1,0.));

        int d0 =s0;
        P[0][d0] = 1.;

        double ans = 0.;
        int s = s0;

        for (int i = 1; i< P.size(); ++i) {
            s -= 2;
            if (s<0) s = 0;
            for (int j = 0; j< P.front().size(); ++j) {
                
                if (j-2 >= 0) P[i][j] += P[i-1][j-2]*0.25;
                if (j-1 >= 0) P[i][j] += P[i-1][j-1]*0.25;
                P[i][j] += P[i-1][j]*0.25;
                if (j+1 < P.front().size()) P[i][j] += P[i-1][j+1]*0.25;

                if (j==d0 && s <= 0) {
                    ans += P[i][j]*0.5;
                    P[i][j] = 0;
                }
                else if (j-d0 >= s) {
                    ans += P[i][j];
                    P[i][j] = 0;
                }
                else if (j-d0 <= -s) P[i][j] = 0;
            }
        }

        return ans;

И все бы хорошо, но временная сложность получается , также как и затраты памяти. А допустимое значение может быть аж . То есть нам надо будет несколько эксабайт памяти и очень много времени.

Однако, мы можем заметить, что с ростом ответ приближается к , и при отличается от менее, чем на . А требуемая точность и есть . То есть мы можем просто добавить в начале строчку if (n>5000) return 1.;. При мы имеем , то есть получается массив , а если оптимизировать и того меньше. Именно такой тип решения и предлагается в качестве "официального" на LeetCode.

Наблюдения и размышления (продолжение)

Вышеизложенное всем известно и никому не интересно. А что, если придумать решение по времени и по затратам памяти? Давайте придумаем.

Прежде всего вспомним, что область влияния состояний, которые соответствуют закончившемуся супу, ограничена серыми областями на рисунке. Что будет, если мы не будем обнулять эти вероятности? Раз они остаются в серой зоне, мы можем их оставить "на потом", и сложить уже только на последнем шаге.

Мы получили весьма интригующую картину. Вероятности оказались симметричны, более того, получившееся дерево одинаково для всех , разница только в том, где его обрубать.

К сожалению, ответ не совпадает. И мы видим, почему. Состояние с на предпоследнем шаге влияет на правую серую область и на центральное значение, а также состояние с на предпоследнем шаге влияет на центральное.

Тогда давайте уменьшим число шагов на 1, а последний сделаем "ручками". Для нечетных надо дополнительно обработать только один белый кружок .

Для четных ручками будем обрабатывать аж 3 кружка:

Ура! Получилось.

Но что дальше? Нам все равно придется строить пирамиду. Или нет?

Избавившись от граничный условий, мы получили следующую задачу. Есть целое число , в начальный момент . На каждом шаге к нему добавляется , , , с одинаковой вероятностью, обозначим это событие как . Найти вероятность того, что на шаге получится .

Если бы у нас было, к примеру, только два возможных исхода, мы получили бы биномиальное распределение https://ru.wikipedia.org/wiki/Биномиальное_распределение.

Но у нас не два, а четыре варианта для каждого шага. И тут мы вспоминаем, что четыре - это дважды два. Что если вместо каждого броска с четырьмя исходами, делать два броска с двумя исходами? Вместо того, чтобы выбирать одно из четырех возможных слагаемых , выберем сначала одно из двух , а потом еще одно из другой пары .

Выпишем все возможные исходы:

1. выбрали и : , вероятность ;

2. выбрали и : , вероятность ;

3. выбрали и : , вероятность ;

4. выбрали и : , вероятность .

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

Какова вероятность, что после пар бросков ?

Пусть в паре раз выпала двойка, вероятность этого равна .

В свою очередь, в паре раз выпала , вероятность .

Чтобы получить необходимо, чтобы . Просуммируем вероятности:

Выразим через :

Видим, что не может быть меньше, чем , и больше, чем .

Формулу получили, осталось применить её к нашей задаче.

Посчитаем вероятность того, что на шаге супа с котом больше или столько же, сколько супа с говядиной . Почему не наоборот? Потому что это число меньше, а значит, при одинаковой относительной погрешности вычисления будем иметь меньшую абсолютную погрешность. Кстати отметим, что минимальное возможное значение на каждом шаге равно . Оно соответствует случаю, когда каждый раз выпадала .

Вернемся к исходной форме выражения:

не может быть больше, чем , и в то же время . Если внешнее суммирование делать по , то его пределы будут от до . Тогда:

Однако, если мы применим данную формулу на предпоследнем шаге, результат будет включать вероятности, которые мы хотели обрабатывать "ручками", для случая , а для четных также и . К тому же, нам будут нужны сами эти вероятности отдельно, а также вероятность для четных .

Выпишем формулы для них:

Посмотрев на рисунки 4 и 5 и подумав, получим окончательное выражение для искомой вероятности того, что суп с говядиной закончится раньше. Вспоминаем, что максимальное количество порций , значит на предпоследнем шаге :

1. для нечетных s:

2. для четных s:

Теперь бы все это посчитать.

Для начала, нам нужны биномиальные коэффициенты. Поскольку у нас фиксировано, можно рассчитать их заранее. Проницательно отметив, что

будем считать только половину из них. Но если велико, то и биномиальные коэффициенты большие, может случиться переполнение. К счастью, нам нужны не сами коэффициенты, а вероятности, то есть:

Тут возникает обратная проблема: — очень маленькая величина.

А так как факториалы считаются последовательно, следующий выражается через предыдущий, то умножение на даст .Чтобы решить эту проблему, будем делить на степени двойки постепенно. Сначала на 1, потом на 4, потом на 16 и т.д. Почему именно так? Потому что мы считаем только половину из биномиальных коэффициентов, соответственно, к концу дойдет до нужной степени. А самое главное, наибольшее значение как раз в середине:

Параллельно будем считать сумму:

Чуть раньше, вы наверное, удивлялись, как же так, ведь содержит две суммы, почему же я заявляю время ? А потому, что суммы по считается здесь один раз и за линейное время. Очень удобно: посчитал коэффициент и сразу прибавил.

Кстати, вы помните, что ?

Алгоритм №2.1. Расчет вероятностей биномиального распределения

1. Создаём массивы для хранения биномиальных коэффициентов и сумм размером каждый;

2. Инициализируем переменную , в которой будем хранить текущий биномиальный коэффициент деленный на какую-то (сами не знаем, какую, но можем сочинить) степень двойки;

3. Записываем первый коэффициент (помним, что целочисленные степени двойки считаются очень просто);

4. То же самое записываем в первую ячейку массива сумм ;

5. Запускаем цикл по от до включительно;

6. Модифицируем значение:

   Вы можете легко получить это выражение, если разделите на ;

7. Записываем биномиальный коэффициент с учетом уже использованных степеней двойки: ;

8. Вычисляем текущую сумму: ;

9. Повторяем 6-9

Удобно сделать отдельный класс, который при инициализации будет считать коэффициенты (точнее, вероятности), а потом по запросу выдавать, в зависимости от того, из какой половины — посчитанной или симметричной — требуется коэффициент:

class Binom {
        vector<double> binomials;
        vector<double> sums;
        int i;

    public:
        Binom(int _i) : i(_i), binomials(_i/2 + 1), sums(_i/2 + 1, 0) {
            double val = 1.;

            binomials.front() = val*exp2(-i);
            sums.front() = binomials.front();
            for (int k = 1; k < binomials.size(); k++) {
                val *= (i-k+1.) / k;
                binomials[k] = val*exp2(-(i-2*k+2));
                sums[k] = sums[k-1] + binomials[k];
                val /= 4;
            }
        }

        double get_k(size_t k) {
            if (k < 0) return 0.;
            if (k < binomials.size()) return binomials[k];
            if (k <= i) return binomials[i - k];
            return 0;
        }

        double get_sum(int k) {
            if (k < 0) return 0.;
            if (k < sums.size()) return sums[k];
            if (k < i) return 1. - sums[i - k - 1];
            return 1.;
        }
    };

Отметим, что здесь мы возвращаем коэффициент 0, если выходит за допустимые пределы, что сделает подсчет сумм удобнее. Сумма, соответственно, 0 для , 1 для , но это добавили просто для согласованности, реально таких значений не будет.

Наконец, начинаем считать. Только предварительно переобозначим индекс суммирования в выражении для . Вместо:

напишем:

И суммировать все выражения будем до , все равно для неправильных значений коэффициенты будут 0.

Да, вы ведь еще помните, что ?

Алгоритм №2.2. O(n)/O(n)

1. Находим начальное ;

2. Находим индекс предпоследнего шага ;

3. Считаем вероятности биномиального распределения и их суммы;

4. Задаем начальные значения сумм ;

5. Запускаем цикл по от до включительно;

6. Вычисляем сумму по в :

7. Добавляем слагаемые к искомым вероятностям:

8. Увеличиваем на и повторяем 6-8;

9. если нечетное, выдать ответ:

10. если четное, выдать ответ:

double soupServings(int n) {
    if (n == 0) return 0.5;
    if (n > 5000) return 1.;
                               
    int s = (n - 1) / 25 + 1;
    int i = (s - 1) / 2;

    Binom bin(i);

    double ps = 0., p0 = 0., p_1 = 0., p1 = 0.;

    for (int q = 0; q*2 <= i; q++) {
        ps += (1. - bin.get_sum(q*2-1)) * bin.get_k(q);
        p0 += bin.get_k(2*q) * bin.get_k(q);
        p_1 += bin.get_k(2*q + 1) * bin.get_k(q);
        p1 += bin.get_k(2*q + 1) * bin.get_k(q + 1);
    }

    double ans = 1. - ps + p0*0.625;
    if (!(s&1)) ans +=  p_1*0.375 - p1*0.125;
    return ans;
}

Решение по времени и по затратам памяти готово.

Кстати, if (n > 5000) return 1.; в начале оставляем все равно. Мы ведь не дурные, чтобы каждый раз заново вычислять единицу, правда?.

Алгоритм №3. O(n)/O(1)

И тут кто-то вспомнил, я ведь обещал по использованию памяти. Неужели наврал?

Нет, и у нас уже есть всё, что нужно. Мы просто будем вычислять биномиальные коэффициенты по мере надобности, параллельно вычислению сумм. Выглядит это примерно так:

double soupServings(int n) {
    if (n == 0) return 0.5;
    if (n > 5000) return 1.;
    
    int s = (n - 1) / 25 + 1;
    int i = (s+1) / 2 - 1;

    double ps = 1., ksum = 1.;
    double p0 = 1., p_1 = i, p1 = 0.;

    double bin_q = 1.0, bin_2q = i;
    double iP1 = i + 1;

    for (int q = 1; q*2 <= i; q++) {
        // Renormalization
        bin_q /= 4;
        bin_2q /= 4;

        ksum /= 4; // Includes bin_2q
        ps /= 4;   // Includes bin_q
        p0 /= 16;  // Includes bin_q * bin_2q
        p_1 /= 16; // Includes bin_q * bin_2q
        p1 /= 16;  // Includes bin_q * bin_2q

        // Calculation
        bin_q *= iP1/q - 1.;
        p1 += bin_q * bin_2q;
        ksum += bin_2q;
        ps += bin_q * (1. - ksum * exp2(2*q - i));

        bin_2q *= iP1/(2*q) - 1.;
        ksum += bin_2q;
        p0 += bin_q * bin_2q;

        bin_2q *= iP1/(2*q+1) - 1;
        p_1 += bin_q * bin_2q;
    }

    
    if (i & 1) {
        int q = (i+1) / 2;
        bin_q *= iP1/q - 1.;
        p1 += bin_q * bin_2q;

        ps /= 2;
        p0 /= 4;
        p_1 /= 4;
        p1 /= 4;
    }
    
	double ans = 1. - ps + p0*0.625;
    if (!(s&1)) ans +=  p_1*0.375 - p1*0.125;
    return ans;

}

Не спрашивайте, как, но оно работает ?.

Заключение

Что в итоге? Решили очень полезную для владельцев тошниловок задачу. Порисовали картинки. Поиграли с биномиальным распределением. Научились считать отношения больших чисел без переполнений. Немного повыпендривались.