Это не ответ, а расширенный комментарий.
Вы уверены, что вам нужна "точная" формула?
Я тут прикинул, разложив функцию возведения в степень в ряд Тейлора.
Ваша функция представляется таким рядом:

Здесь в скобках стоит обобщённый биномиальный коэффициент: 
Это разложение сходится при |s-1| < 1
Я сделал прототип на Python
P=2.0
Q=1.0
T=1.0
def C(s,t):
return P - Q/(s-1) + (Q + Q/(s-1))/pow(s,t/T)
# Обобщённый биномиальный коэффициент
def binom(a,n):
res = 1
for k in range(1,n+1):
res *= (a-k+1)/k
return res
# Вычисление разложением в ряд. По умолчанию ряд вычисляется до 4 степени
def C2(s,t, n=4):
d = s - 1
f = 0
for i in range(1,n+1):
f += binom(-t/T, i)*(d**(i-1))
return P + Q + Q*(1+d)*f
# Сравнение исходной формулы и разложения в ряд до 20 степени (оверкил)
# Перебираем epsilon от 0.1 до 1e-14
for p in range(1, 15):
s = 1+pow(10,-p)
c1 = C(s,0.5)
c2 = C2(s,0.5, n=20)
print(f"e = 1e-{p}: s = {s}, c1 = {c1}, c2 = {c2}, abs(c2-c1) = {abs(c2-c1)}")
Вот что получилось
e = 1e-1: c1 = 2.4880884817015154, c2 = 2.4880884817015154, abs(c2-c1) = 0.0
e = 1e-2: c1 = 2.4987562112089137, c2 = 2.4987562112089026, abs(c2-c1) = 1.1102230246251565e-14
e = 1e-3: c1 = 2.4998750624608874, c2 = 2.4998750624609647, abs(c2-c1) = 7.72715225139109e-14
e = 1e-4: c1 = 2.4999875006251386, c2 = 2.499987500624961, abs(c2-c1) = 1.7763568394002505e-13
e = 1e-5: c1 = 2.4999987500195857, c2 = 2.4999987500062497, abs(c2-c1) = 1.333599897179738e-11
e = 1e-6: c1 = 2.499999874853529, c2 = 2.4999998750000625, abs(c2-c1) = 1.4653345203896606e-10
e = 1e-7: c1 = 2.499999986961484, c2 = 2.4999999875000007, abs(c2-c1) = 5.385167867189011e-10
e = 1e-8: c1 = 2.5, c2 = 2.49999999875, abs(c2-c1) = 1.2500001034254637e-09
e = 1e-9: c1 = 2.5, c2 = 2.499999999875, abs(c2-c1) = 1.2500001034254637e-10
e = 1e-10: c1 = 2.5, c2 = 2.4999999999874998, abs(c2-c1) = 1.2500223078859563e-11
e = 1e-11: c1 = 2.5, c2 = 2.49999999999875, abs(c2-c1) = 1.2501111257279263e-12
e = 1e-12: c1 = 2.5, c2 = 2.4999999999998748, abs(c2-c1) = 1.2523315717771766e-13
e = 1e-13: c1 = 2.5, c2 = 2.4999999999999876, abs(c2-c1) = 1.2434497875801753e-14
e = 1e-14: c1 = 2.515625, c2 = 2.4999999999999987, abs(c2-c1) = 0.015625000000001332
Видно, что почти для всех эпсилонов расхождение между исходной формулой и разложением в ряд не хуже чем 1e-9
Заметная разница только для 1e-14
, но это и так предел точности представления вещественных чисел в double
.
Если для вас точность девять знаков после запятой достаточна, то не парьтесь с приблизительной формулой.
UPD: сравнение с приближённой формулой
Сравнил приближённую формулу C(s->1, t) = P + Q*(1 - t/T)
с результатом разложения в ряд. В выводе ниже c2
результат разложения в ряд, c3
- значение приближенной формулы.
e = 1e-1: c3 = 2.5, c2 = 2.4880884817015154, abs(c2-c3) = 0.01191151829848458
e = 1e-2: c3 = 2.5, c2 = 2.4987562112089026, abs(c2-c3) = 0.0012437887910974332
e = 1e-3: c3 = 2.5, c2 = 2.4998750624609647, abs(c2-c3) = 0.00012493753903530802
e = 1e-4: c3 = 2.5, c2 = 2.499987500624961, abs(c2-c3) = 1.2499375038999005e-05
e = 1e-5: c3 = 2.5, c2 = 2.4999987500062497, abs(c2-c3) = 1.2499937502852276e-06
e = 1e-6: c3 = 2.5, c2 = 2.4999998750000625, abs(c2-c3) = 1.2499993751191596e-07
e = 1e-7: c3 = 2.5, c2 = 2.4999999875000007, abs(c2-c3) = 1.2499999257897798e-08
e = 1e-8: c3 = 2.5, c2 = 2.49999999875, abs(c2-c3) = 1.2500001034254637e-09
e = 1e-9: c3 = 2.5, c2 = 2.499999999875, abs(c2-c3) = 1.2500001034254637e-10
e = 1e-10: c3 = 2.5, c2 = 2.4999999999874998, abs(c2-c3) = 1.2500223078859563e-11
e = 1e-11: c3 = 2.5, c2 = 2.49999999999875, abs(c2-c3) = 1.2501111257279263e-12
e = 1e-12: c3 = 2.5, c2 = 2.4999999999998748, abs(c2-c3) = 1.2523315717771766e-13
e = 1e-13: c3 = 2.5, c2 = 2.4999999999999876, abs(c2-c3) = 1.2434497875801753e-14
e = 1e-14: c3 = 2.5, c2 = 2.4999999999999987, abs(c2-c3) = 1.3322676295501878e-15
Получается, что ошибки округления начинают сказываться только для epsilon порядка 1e-14
, а при выбранной вами точности 0.00001
ошибка приближенной формулы на пять порядков хуже, чем у точной. Поэтому я могу лишь повторить "посыл" - не беспокойтесь о потери точности из-за малости отклонения от 1.
s==1
, а вся полуосьs>0
. Для значенийabs(s-1) <= EPSILON
вычисляю по формуле (2), но в окрестностиEPSILON
всё равно в знаменателе получаются близкие к нулю числа и вопрос точности остаётся актуален.s2 -> (s - 1)
,p - q/s2 + (q + q/s2)/((1+s2)^(t/T))
,lim s2 -> 0
=>p - q/s2 + (q + q/s2)/(1+s2*(t/T))