На чистом python, максимально просто, ваша функция выглядит так:
import math
_count_log10_e = math.log10(math.e) # 1./math.log(10)
assert 1. == _count_log10_e*math.log(10)
def count(n):
return 1 + int(math.lgamma(n + 1)*_count_log10_e)
assert count(5) == 3
assert count(50) == 65
assert count(500) == 1135
assert count(5000) == 16326
assert count(50000) == 213237
assert count(500000) == 2632342
assert count(5000000) == 31323382
assert count(50000000) == 363233781
for n in range(10):
assert count(n) == len(str(math.factorial(n)))
Подскажите какая нужна оптимизация!!!
Поскольку n! = Γ(n + 1), а число цифр это ⎣1 + log10 n!⎦, то math.lgamma() - она самая.
P.S.
Дело было ... и делать было ... Если "есть о чём поговорить", продолжаем разговор.
Ошибки точности
Внезапно, прям как зимний снегопад, оказалось, что для вычисления
длины факториала точности float бывает недостаточно. Причём, как
всегда, неожиданно рано. А в формуле длины факториала
1 + ⎣log₁₀ n!⎦ предполагается,
что log₁₀ n! не бывает целым, и точность промежуточных
вычислений должна позволять выяснить больше он ближайшего
целого или меньше.
- Начиная с
n >= 28_563_732, представление float (binary64
по IEC 60559) значения log₁₀ Γ(n + 1) может
оказываться целым (немного повезло, что первое округление в удачную
сторону, поэтому собственно ошибки count(n) начинаются немного
позднее, с n >= 46_464_264);
- Так же, лично для меня, оказалась неожиданностью, что у
CPython
своя собственная реализация гамма-функции и она несколько хуже,
чем стандартные BSD/macOS/Microsoft и даже glibc (мой
поверхностный поиск не дал информации о точности реализации math.lgamma,
но комментарии "it's probably not worth it", как бы, намекают).
Конечно, в условиях не сказано, какой должен быть диапазон, однако,
исходя из test.assert_equals(count(50000000),363233781) можно
использовать поправку
+ (n == 48_655_817) - (n in [46_464_264, 54_336_595, 54_528_830, 74_048_650]),
типа, с двойным запасом.
Впрочем, есть парочка, в некотором смысле, более честных, хотя и
заметно менее эффективных решений.
gmpy2.lgamma()
Поскольку на codewars.com
есть gmpy2 (жаль, что нет православного mpmath), то решение на нём
будет самым коротким среди тех, у которых диапазон много больше
50_000_000 (точность по умолчанию 107-бит, имитирует точность
компенсационной арифметики, а ля double-double, алгоритмы Кэхана,
Неймура и т.п.).
import gmpy2
def gcount(n, precision=107, check=True, **kwargs):
with gmpy2.context(precision=precision, **kwargs):
log10_e = gmpy2.log10(gmpy2.exp(1))
l10g = gmpy2.lgamma(n + 1)[0]*log10_e
i, f = gmpy2.modf(l10g)
if check and l10g:
l10ge = gmpy2.next_above(l10g) - gmpy2.next_below(l10g)
if f < l10ge or 1 - f < l10ge:
# Как вариант, рекурсивный вызов с повышенной точностью
raise ValueError(f"Не хватает точности {f, l10ge=}")
return 1 + i
decimal и ряд Стирлинга
Решение на "чистом" Python, честно покрывающее диапазон до
50_000_000 и далее, увы, немного длиннее предыдущего и
ещё немного медленнее (math.lgamma() ∼ 300 нс,
gmpy2.lgamma() ∼ 30 мкс, а используя decimal ∼ 200 мкс),
но, главное, оно тоже O(1). Точность decimal по
умолчанию - 28 знаков (≥ 93 бита), её хватит на многие и
многие миллиарды. 😉 В принципе и с 100⁵⁰⁰! может справится, хотя уже и на пределе.
import decimal
def dcount(n, check=True, **kwargs):
with decimal.localcontext(**kwargs):
log10_e = decimal.Decimal(1).exp().log10()
lg, lge = decimal_simple_lgamma(n + 1)
l10g = lg*log10_e
i = int(l10g)
if check and l10g:
l10ge = lge*log10_e
f = l10g - i
if f < l10ge or 1 - f < l10ge:
# Как вариант, рекурсивный вызов с повышенной точностью
raise ValueError(f"Не хватает точности {n, f, l10ge=}")
return 1 + i
def decimal_simple_lgamma(x) -> tuple:
"""Простая реализация логарифма гамма-функции для оценки
длины факториала. Использует только первый член ряда
Стирлинга (A046968, A046969), он же является оценкой
остаточного члена ряда сверху.
"""
n = decimal.Decimal(x) - 1
if 0 == n or 1 == n:
return (decimal.Decimal(0), decimal.Decimal(0))
pi = decimal.Decimal('3.14159265358979323846264338327950288419'
'7169399375105820974944592307816406286209')
lg = n*n.ln() - n + (2*pi*n).ln()/2
n12 = 12*n
in12 = 1/n12
err = in12 - 1/(n12 + 1)
return (lg + in12, abs(err) + 2*(lg.next_plus() - lg.next_minus()))
Double-double math.fsum() логарифмов
Забавно, что суммирование обычных логарифмов из math тоже вполне способно выдавать корректный
ответ для достаточно больших n. Ошибку суммы
можно оценить, грубо говоря, пусть у нас ln x рассчитывается
с точностью ± 0.5 ε, типа, дисперсия
σ² ≈ ε²/12:
СКО(Σ(i=1..n)( N(0,σ)×ln(i) )) ≈ σ × sqrt(∫(x=1..n) ln²(x) dx) =
= σ × sqrt(n×(ln²(n) - 2×ln(n) + 2) - 2)
Т.е. хотя и СКО (как и дисперсия) суммы растёт, относительная ошибка этой суммы даже немного снижается, типа sqrt(1/n).
Однако требования на необходимую относительная точность увеличивается, типа, 1/n² или 1/n³,
так что, да, где-то они встретятся, но
достаточно далеко.
В принципе, в коде ответа @FoxFox можно было бы просто заменить sum() на два sum() (или на два math.fsum() для <3.12), да ответы становятся гораздо, гораздо точнее, чем простой float.
Одна беда, оценка ошибки для ∑ log₁₀ i, что-то как-то не сходятся, то ли в библиотечной log10() засада (она реально, для этой цели, менее точная, чем log() или log2()), то ли руки кривые. Поэтому в макетике использована log(), для неё формулы по-проще, да и сравнение с гамма-функцией тривиально.
Вышеприведённая оценка ошибки, как и сумма, отлично работает до многих миллиардов, правда, вероятность параноидального ложноотрицательного завышения ошибки - 2..3 раза на миллиард. Возможно, в той формуле она просто немного.
import math
import numpy as np
def sum_ln_err(n):
return (math.ulp(1)/math.sqrt(12))*math.sqrt(
n*(math.log(n)**2 - 2*math.log(n) + 2) - 2)
class two_factorial_length_np:
_chunk_size = 20_000
_sum_not_double_double = (1e-37 != sum([1., 1e-37, -1.]))
def __init__(self, fsum=_sum_not_double_double):
self.fsum = fsum
self.n = 1
self.sum = np.zeros(2)
self.res = 1
def __call__(self, n, check=True):
if n <= 1:
return 1
if self.n == n: # TODO: check при повторном вызове игнорируется?
return self.res
assert self.n < n
while self.n < n:
c = min(n, self.n + self._chunk_size)
a = np.arange(self.n + 1, c + 1 + 3, dtype=float)
a[:-3] = np.log(a[:-3])
a[-3:-1] = self.sum
if self.fsum:
self.sum[0] = math.fsum(a[:-1])
a[-1] = -self.sum[0]
self.sum[1] = math.fsum(a)
else:
self.sum[0] = sum(a[:-1].tolist())
a[-1] = -self.sum[0]
self.sum[1] = sum(a.tolist())
assert self.sum[0] == self.sum[0] + self.sum[1]
self.n = c
l10g = two_mul(self.sum, _two_log10_e)
i = two_floor(l10g)
self.res = 1 + i
if check and l10g[0]:
f = two_add(l10g, -i) # TODO: случай i >= 2**53
err = sum_ln_err(self.n)*_two_log10_e[0]
if f[0] <= err or (1 - f[0]) <= err:
raise ValueError(f"Не хватает точности {f, err=}")
return self.res
_codewars_state = two_factorial_length_np()
def _codewars_count(n):
global _codewars_state
if n <= 1:
return 1
if n < _codewars_state.n:
_codewars_state = two_factorial_length_np()
return _codewars_state(n)
#----
# Функции double-double
# Половина нужна для умножения на константу, однако, другая
# половина требуется для контроля точности (вычитания) и
# небольшая функция two_floor()
import ctypes
import decimal
import sys
if sys.platform.startswith('win'):
libm = ctypes.cdll.msvcr120
elif sys.platform.startswith('darwin'):
libm = ctypes.CDLL("libc.dylib")
elif sys.platform.startswith('freebsd'):
libm = ctypes.CDLL("libm.so")
else:
libm = ctypes.CDLL("libm.so.6")
def quick_two_sum(a, b) -> tuple:
assert abs(a) >= abs(b)
hi = a + b
lo = b - (hi - a)
assert hi == hi + lo
return hi, lo
def two_sum(a, b) -> tuple:
if abs(a) >= abs(b):
return quick_two_sum(a, b)
return quick_two_sum(b, a)
def two_add(a: tuple, b) -> tuple:
if isinstance(b, tuple):
hi, lo = two_sum(a[0], b[0])
lh, ll = two_sum(a[1], b[1])
lo += lh
hi, lo = two_sum(hi, lo) # TODO:??? quick_two_sum(hi, lo)
lo += ll
else:
hi, lo = two_sum(a[0], b)
lo += a[1]
return two_sum(hi, lo) # TODO:??? quick_two_sum(hi, lo)
if 'fma' in math.__dict__:
two_fma = math.fma
else:
# Для <3.13 есть pyfma, libc (ctypes) или известный сплит-алгоритм
libm.fma.restype = ctypes.c_double
libm.fma.argtypes = (ctypes.c_double, ctypes.c_double, ctypes.c_double)
two_fma = libm.fma
def two_prod(a, b) -> tuple:
hi = a * b
lo = two_fma(a, b, -hi)
assert hi == hi + lo
return hi, lo
def two_mul(a: tuple, b: tuple) -> tuple:
assert a[0] == a[0] + a[1]
assert b[0] == b[0] + b[1]
hi, lo = two_prod(a[0], b[0])
lo += a[0] * b[1]
lo += a[1] * b[0]
return quick_two_sum(hi, lo)
def two_floor(a: tuple) -> int:
assert a[0] == a[0] + a[1]
i = math.floor(a[0])
if i == a[0]:
return i + math.floor(a[1])
return i
assert (two_prod(float.fromhex('0x1.0000000000001p+0'),
float.fromhex('0x1.0000000000001p+0'))
== (float.fromhex('0x1.0000000000002p+0'),
float.fromhex('0x1.0000000000000p-104')))
assert (two_mul(two_prod(float.fromhex('0x1.0000000000001p+0'),
float.fromhex('0x1.0000000000001p+0')),
two_prod(float.fromhex('0x1.fffffffffffffp-1'),
float.fromhex('0x1.fffffffffffffp-1')))
== (float.fromhex('0x1.0000000000001p+0'),
float.fromhex('-0x1.8000000000002p-105')))
assert 1 == two_floor((float.fromhex('0x1.0000000000001p+0'),
float.fromhex('-0x1.8000000000002p-105')))
assert 0 == two_floor((float.fromhex('0x1.0000000000000p+0'),
float.fromhex('-0x1.8000000000002p-105')))
with decimal.localcontext(prec = sys.float_info.dig*3) as ctx:
_two_log10_e = decimal.Decimal(1).exp().log10()
_two_log10_e_hi = float(_two_log10_e)
_two_log10_e_lo = float(_two_log10_e -
decimal.Decimal(_two_log10_e_hi))
assert _two_log10_e_hi == _two_log10_e_hi + _two_log10_e_lo
_two_log10_e = (_two_log10_e_hi, _two_log10_e_lo)
_two_ln_10 = decimal.Decimal(10).ln()
_two_ln_10_hi = float(_two_ln_10)
_two_ln_10_lo = float(_two_ln_10 -
decimal.Decimal(_two_ln_10_hi))
assert _two_ln_10_hi == _two_ln_10_hi + _two_ln_10_lo
_two_ln_10 = (_two_ln_10_hi, _two_ln_10_lo)
_test_one = two_mul(_two_log10_e, _two_ln_10)
assert 1 == _test_one[0]
assert 1 >= _test_one[1]/math.ulp(0.5)**2
Отдельное спасибо за ответы @StanislavVolodarskiy и @@CrazyElf.
Полный код с тестами, парой графиков и возможностью запуска в Colab: https://github.com/Serge3leo/temp-cola/blob/main/ruSO/1615333-Найти-длину-факториала/Варианты-и-график.ipynb