The following function is supposed to work similarly to pow(x, 1/k) but to be symmetric around the line y = 1 - x as well as not having a 0 or 1 slope at either end of [0, 1]:
def sym_gamma(x, k):
if k == 1.0:
return x
a = 1.0 / k - 1.0
b = 1.0 / a
c = k + 1.0 / k - 2.0;
return 1.0 / (a - c * x) - b
As can be seen, it is not defined when k = 1 so when that is the case, I simply return x. However, this special case handling is not enough since the function also behaves poorly when x is not equal to but very close to 1.0. For example sym_gamma(0.5, 1.00000001) yields 0.0 while it's supposed to return something very close to 0.5.
How can achieve the same thing without the poor stability? I know that I can introduce a tolerance with respect to k equaling 1.0 but it feels like a hack and I would also want to make sure that the function is perfectly smooth with regards to k.
Simplifying your expression seems to help with the precision. Numerical errors tends to accumulate in each operation. Thus, reducing the number of operation will reduce the chance of numerical errors.
We can notice that:
a = (1 - k) / k
b = k / (1 - k)
c = (1 - k) ** 2 / k
a - c * x = (1 - k) * (1 + x*k - x) / k
1.0 / (a - c * x) - b = x*k / (1 - x * (1 - k))
Then you can simply rewrite your method:
def sym_gamma(x, k):
return x*k / (1 - x * (1 - k))
Instead of performing several division, only one division is computed. This method returns 0.5000000025 for sym_gamma(0.5, 1.00000001).