I have some code which involves tiny numbers being put to the power of very large numbers. The function I'm looking at is

$\varphi = \omega(T) \left(1 - (1 - \epsilon)^{n_{e}(T)} \right)$

when $\epsilon $ is very small (~$10^{-16}$) and $n_{e}$ is large ($> 10^{10}$). Both $n_{e}$ and $\omega$ are known functions of $T$. At a known value ($T = 37$) the value of $\varphi$ is known, and we can in theory estimate $\epsilon$ by

$\epsilon = 1 - \left(1 - \frac{\varphi}{\omega} \right)^{\frac{1}{n_{e}}}$

However, when I put this into Mathematica I'm getting very strange behaviour; here's a MWE when $\varphi(37) = 0.2506$;

gammaex = 0.2506; 
omega[t_] := 2.43163218375*10^7*Exp[1700*(1/298.15 - 1/(273.15 + t))];

w[t_] := (3.414105049212413*10^12)/(omega[t]); 
v[t_] := Sqrt[661.6469313477045*(t + 273.15)];
ne[t_] := (v[t]*5.104757516005496*(10^7));

epsilon = SetAccuracy[ 1 - ( 1 - gammaex/w[37])^(1/ne[37]), 30];
test = w[37]*(1 - (1 - epsilon)^(ne[37]))

Now when I evaluate $\varphi(37)$ via the last line I should get the value 0.2506 but I do not; instead I get $\varphi = 0.289104$, over 15% off the true value. I thought maybe this was a precision problem, so I tried some other commands (Surd, Power...) and got the same wrong value. The output value for the epsilon is $\epsilon \approx 1.11 \times 10^{-16} $ from mathematica.

However I've evaluated this with WolframAlpha and got a markedly different value of $\epsilon \approx 9.61 \times 10^{-17}$, and as this link demonstrates the calculation seems to work with WolframAlpha, returning close to the expected value around 0.25*. Any ideas why the calculations are different, and how I can make Mathematica behave with such extreme values?

Incidentally, if I dump the WolframAlpha equation ( (112608)(1 - (1 - 9.6136910^-17)^(2.31246*10^10)) ) directly into Mathematica, it still comes up with the wrong value, so I assume it's a radical issue?