3

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.

Emil Eriksson
  • 2,110
  • 1
  • 21
  • 31
  • @KellyBundy Sorry, perhaps pow(x, 1/k) is more accurate. Here's a plot in Desmos where you can play around with different values of k: https://www.desmos.com/calculator/vg5eini1o1 – Emil Eriksson Jun 20 '23 at 19:47

2 Answers2

5

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).

ftorre
  • 501
  • 1
  • 9
  • This is perfect. This is also valid for `k = 1` so no need for special casing. I initially tried to do this but I guess I was too lazy to try harder to simplify it. – Emil Eriksson Jun 20 '23 at 20:49
1

As already answered by ftorre. Simplifying the equations into one would reduce the chance of numerical errors. This problem is, however, caused by the limited precision in the floating-point representation of the variables. Floating-point representation says about the way numbers are stored in actual memory. The default precision for float type variables in Python is double precision (meaning 8+8=16 values after the decimal). This can vary with the specific Python implementation, and when multiple operations are cascaded, the final result often comes out unexpected as above.

A more robust solution to your problem would be to use the decimal module, which lets you choose the precision. Here's an example of your code using the module:

from decimal import Decimal, getcontext

def sym_gamma(x, k):
    if k == 1.0:
        return x

    precision = 30
    getcontext().prec = precision

    a = Decimal(1.0) / Decimal(k) - Decimal(1.0)
    b = Decimal(1.0) / a
    c = Decimal(k) + Decimal(1.0) / Decimal(k) - Decimal(2.0)

    x = Decimal(x)

    return float(Decimal(1.0) / (a - c * x) - b)

This code also preserves your values of a, b, and c. You can set the precision to your desired value (I have preset it to 30)

Dev Bhuyan
  • 541
  • 5
  • 19
  • 1
    If you don't define `precision` inside the function and do `for precision in range(14, 36): print(precision, sym_gamma(0.5, 1.00000001))`, you can see how increasing it affects the result. – Kelly Bundy Jun 20 '23 at 21:09
  • 1
    And `x = Decimal(str(x))` and `k = Decimal(str(k))` help further (at least if those numbers are entered in decimal, like `1.00000001`, since it undoes the damage of the evaluation to `float`). – Kelly Bundy Jun 20 '23 at 21:10