如何提升sym_gamma函数的数值稳定性？解决k趋近1时的异常

原sym_gamma函数意图实现类似pow(x, 1/k)的功能，同时满足关于直线y=1-x对称、区间[0,1]两端无0或1斜率的特性，但当k接近1时会出现数值不稳定（比如sym_gamma(0.5, 1.00000001)返回0.0，预期应接近0.5），且原函数依赖k=1的特殊分支，无法平滑过渡。

通过代数化简原函数的表达式，可以得到一个完全等价、且在k=1处自然连续的形式，彻底避免数值不稳定：

def sym_gamma(x, k):
    return k * x / (1 - (1 - k) * x)

化简推导 #

原函数的表达式可逐步化简：

1. 定义a = 1/k - 1 = (1 - k)/k，b = 1/a = k/(1 - k)，c = k + 1/k - 2 = (k-1)²/k
2. 将a、b、c代入原返回式：

   1/(a - c*x) - b = 1/[(1-k)/k - (k-1)²x/k] - k/(1-k)
   

3. 提取公因子(1-k)/k到分母外，通分后分子分母的(1-k)项可抵消，最终化简为k*x / (1 - (1 - k)*x)。

特性验证 #

• k=1时的连续性：当k=1时，表达式自然变为1*x/(1 - 0*x) = x，无需特殊分支，完全符合预期。
• k接近1时的稳定性：以k=1.00000001、x=0.5为例，计算得(1.00000001*0.5)/(1 - (-0.00000001)*0.5) ≈ 0.500000005/1.000000005 ≈ 0.4999999975，与预期的0.5几乎一致。
• 对称性：若y = sym_gamma(x, k)，则sym_gamma(1-y, 1/k) = 1 - x，满足关于直线y=1-x对称的要求。
• 两端斜率：x=0处导数为k，x=1处导数为1/k，均不为0或1（k≠0,1时），符合设计需求。

内容的提问来源于stack exchange，提问作者Emil Sahlén