Пока наш главный математик оформляет ответ, вот вам несколько туповатый (из-за приспособления к Numba) код, который с помощью ускорения через Numba.njit считает искомое за 6 секунд. (Без Numba до 100 миллионов молотит за 1.5 минуты, что всё-равно сильно лучше, чем если считать в питоновских бесконечных int-ах.) Правда, оставляет некоторые сомнения в абсолютной точности результата, потому что подсчёт идёт в числах с плавающей точкой. В каком-то уникальном случае в конце числа может оказаться ...999999 и нужный результат, возможно, не найдётся. Хотя можно было бы добавить 1 в последний (а какой это?) разряд и тогда уже обрезать число. Ну и можно ещё более универсально написать и ещё ускорить, но вроде бы уже и так нормально.

import math
from numba import njit

@njit
def fact(n):
    f = 1
    i = 2
    flag1 = True
    flag2 = True
    flag3 = True
    flag4 = True
    flag5 = True
    while i <= n:
        f *= i
        while f >= 10:
            f /= 10
        if flag1 and math.floor(f * 10000) == 11533:
            print(f'{i}!', f)
            flag1 = False
        if flag2 and math.floor(f * 10000) == 51261:
            print(f'{i}!', f)
            flag2 = False
        if flag3 and math.floor(f * 100000) == 113131:
            print(f'{i}!', f)
            flag3 = False
        if flag4 and math.floor(f * 100000) == 131091:
            print(f'{i}!', f)
            flag4 = False
        if flag5 and math.floor(f * 100000000) == 151192033:
            print(f'{i}!', f)
            flag5 = False
        i += 1

fact(100_000_000)

Вывод:

16641! 1.153311564705438
46978! 5.126194982442562
266538! 1.3109139672369952
323172! 1.1313167610953019
37773304! 1.5119203331946358

Число всё время приводится к виду x.xxxxxx, чтобы потом не считать, на сколько разрядов нужно домножать перед сравнением. Поэтому и в выводе такой же вид.

P.S. В чистом Питоне я и имена выводил, но Numba очень привередлива, родные питоновские типы не любит, в частности словарь я сходу не смог победить. Поэтому вот так топорно.

P.P.S. На C++ наверное вообще за секунду всё посчитало бы.

Решение

Вычислим факториал с помощью интервальной арифметики. Интервал задаётся тройкой неотрицательных чисел

I(m, d, e) = [m·10^e, (m + d)·10^e].

Целому n соответствует интервал I(n, 0, 0) = [n·10⁰, (n + 0)·10⁰] = [n, n].

Умножение интервала на целое: I(m, d, e)·f = I(m·f, d·f, e).

Округление интервала на 10^r. Заметим что
⌊^m/_10^r⌋·10^r ≤ m
m + d ≤ ⌈^{(m + d)}/_10^r⌉·10^r

Откуда получается

I(m, d, e) ⊂ I(⌊^m/_10^r⌋, ⌈^{(m + d)}/_10^r⌉, e + r).

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

Округленный интервал содержит внутри исходный.

Программа читает желаемый префикс факториала и количество цифр в мантиссе. Затем она считает последовательные факториалы в виде округлённых интервалов. Например, расчет для пяти цифр в мантиссе:

Пяти цифр достаточно чтобы вычислить две старшие цифры 100!. Удивительная эффективность.

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

Для примеров из вопроса хватает двадцати знаков.

import time


class Interval:
    ''' [mantissa * 10 ** exponent, (mantissa + delta) * 10 ** exponent] '''

    def __init__(self, mantissa, delta=0, exponent=0):
        self.mantissa = mantissa
        self.delta = delta
        self.exponent = exponent

    def __mul__(self, f):
        return Interval(self.mantissa * f, self.delta * f, self.exponent)

    def rounded(self, m_digits):
        exp = max(0, len(str(self.mantissa + self.delta)) - m_digits)
        p = 10 ** exp
        m1 =  self.mantissa                       // p
        m2 = (self.mantissa + self.delta + p - 1) // p
        return Interval(m1, m2 - m1, self.exponent + exp)

    def __str__(self):
        if self.delta == 0:
            return f'{self.mantissa}e{self.exponent}'
        return f'({self.mantissa}+{self.delta})e{self.exponent}'


def main():
    p, m = input().split()
    m = int(m)
    c = 10 ** len(p)

    next_time = 0

    f = Interval(1)
    i = 1
    while True:
        f = (f * i).rounded(m)
        curr_time = time.time()
        if curr_time > next_time:
            print(i, f, end='\r')
            next_time = curr_time + 1
        if f.mantissa < f.delta * c:
            print(f'\n{i} {f}\nfailure')
            break
        if str(f.mantissa          ).startswith(p) and \
           str(f.mantissa + f.delta).startswith(p)     \
        :
            print(f'\n{i} {f}\n{f.mantissa}\n{f.mantissa + f.delta}\ndone')
            break
        i += 1


main()

$ echo 11533 20 | python factorial-prefix.py
16641 (11533115647054388190+7497)e63001
11533115647054388190
11533115647054395687
done

$ echo 51261 20 | python factorial-prefix.py
46978 (51261949824425729445+94135)e199057
51261949824425729445
51261949824425823580
done

$ echo 113131 20 | python factorial-prefix.py
323172 (11313167610953440853+142970)e1640127
11313167610953440853
11313167610953583823
done

$ echo 131091 20 | python factorial-prefix.py
266538 (13109139672370305808+136664)e1330399
13109139672370305808
13109139672370442472
done

$ time echo 151192033 20 | python factorial-prefix.py
37773304 (15119203331949324965+22322287)e2698105593
15119203331949324965
15119203331971647252
done

real  1m44.529s
user  1m44.520s
sys   0m0.000s

Проверка

Убедимся что факториал вычислился верно. Тут две трудности: сам факториал вычисляется не быстро, ещё дольше он печатается. В Питоне печать длинных целых требует квадратичного времени от числа цифр в числе. Чтобы обойти это затруднение, перед печатью разделим факториал на степень десятки так, чтобы осталось около двадцати цифр.

Настин факториал требует около полутора часов для проверки, остальные не больше нескольких секунд:

import math

n = int(input())
f = math.factorial(n)
print(f // 10 ** (math.floor(math.log10(f)) - 20))

$ echo 16641 | python factorial-head.py
115331156470543919298

$ echo 46978 | python factorial-head.py
512619498244257763839

$ echo 323172 | python factorial-head.py
113131676109535123222

$ echo 266538 | python factorial-head.py
131091396723703741317

$ time echo 37773304 | python factorial-head.py
151192033319604853614

real  96m43.944s
user  96m35.068s
sys   0m2.688s

Основная идея: при вычислении старших знаков, меньшие не нужны.

Нам нужны первый n знаков, при вычислениях будем использовать только 2n старших разрядов (есть предположение, что это не всегда будет работать правильно, но для текущих данных верно):

n = 151192033

nsize = len(str(n))
nstr = str(n)
n2size = nsize * 2

fact = 1
i = 1

while str(fact)[:nsize] != nstr:
    i += 1
    fact *= i

    fact = int(str(fact)[:n2size])

print(i, fact)

Выполняется за ~40 секунд.

можно уменьшить количество вызовов функции str в цикле, перенеся условие выхода в сам цикл:

while True:
    i += 1
    fact *= i

    sfact = str(fact)
    if sfact[:nsize] == nstr:
        break
    fact = int(sfact[:n2size])

теперь ~28 секунд.

Реализовал тоже самое не с str, а с делением на степени десятки и вычисление для через math.log10:

pow10 = [10**(i + 1) for i in range(n2size)] + [1] * 100
pow10_size = 10**nsize
while True:
    i += 1
    fact *= i

    fact //= pow10[int(math.log10(fact)) - n2size]
    if fact // pow10_size == n:
        break

~26 секунд. Некоторые параметры просто взяты из ни откуда и если к примеру i станет слишком велик, то будет не очень хорошо

Замечательный ответ Danis заставил меня подумать над оптимизацией решения. Опять интервальная арифметика: после умножения границы интервала усекаются так, чтобы в них оставалось не более m старших цифр. Нижняя граница округляется вниз, верхняя – вверх. У усечённых границ выбираются старшие цифры, которые сравниваются с искомым префиксом. Процесс повторяется в цикле, пока не будет найден ответ, или пока не будет обнаружено, что m слишком мало, чтобы его найти.

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

Я покажу вычисления на примере: префикс 11533 длиной пять цифр и m = 10.

Перед итерацией 16640 интервал имеет вид

[4164975812, 4165002835].

После умножения:

[69305_19751_1680, 69305_64717_4400].

Чтобы усечь его до десяти цифр достаточно разделить границы на 10000. Левую округлить вниз, правую – вверх:

[69305_19751, 69305_64718].

Для проверки совпадения с префиксом обе границы усекаются до пяти цифр (это длина префикса 11533). Для этого делим на 100000:

69305, 69305.

Всё можно сделать арифметикой, если по левой границе интервала из некоторой таблицы извлечь два коэффициента: 10000 и 100000.

Таблица индексируется количеством бит в левой границе интервала. В числе 69_305_197_511_680 сорок шесть бит. Диапазон чисел из сорока шести бит:

[2⁴⁵, 2⁴⁶ - 1] = [35_184_372_088_832, 70_368_744_177_663].

Любое сорокашестибитовое имеет четырнадцать разрядов. Строка в таблице будет иметь вид:

46: 10000, 100000

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

[2⁴⁶, 2⁴⁷ - 1] = [70_368_744_177_664, 140_737_488_355_327].

В этом случае число сравнивается с 10¹⁴. Если число меньше, используется строка 47, иначе строка 48:

47: 10000, 100000, 10¹⁴
48: 100000, 100000

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

def n_digits(n):
    return 0 if n == 0 else 1 + n_digits(n // 10)


def main():
    n, m = map(int, input().split())

    ds = [n_digits(2 ** i // 2) for i in range((10 ** (m + 40)).bit_length())]
    ss = [10 ** max(0, d - m) for d in ds]
    ts = [10 ** d for d in ds]
    cs = [10 ** max(0, min(m, d) - n_digits(n)) for d in ds]

    low = 1
    high = 1
    i = 1
    while True:
        low *= i
        high *= i
        b = low.bit_length()
        j = b + (low >= ts[b])
        low //= ss[j]
        high = (high + ss[j] - 1) // ss[j]
        if n * (high - low) > low:
            print('failure')
            break
        if n == low // cs[j] == high // cs[j]:
            print(i)
            break
        i += 1


main()

Полминуты для Насти:

$ echo 11533 20 | python factorial-prefix.py
16641

$ echo 51261 20 | python factorial-prefix.py
46978

$ echo 113131 20 | python factorial-prefix.py
323172

$ echo 131091 20 | python factorial-prefix.py
266538

$ time -p echo 151192033 20 | python factorial-prefix.py
37773304
real 28.74
user 28.68
sys 0.00