阶乘计算函数对小数/浮点数返回错误结果，请求技术支援

问题概述 #

现有Kotlin实现的阶乘函数对整数计算正常，但处理小数/浮点数时返回错误结果。例如计算3.03的阶乘时，测试预期结果为6.23120891259，实际返回5.8062997。

原阶乘函数代码 #

fun factorial(number: BigDecimal): BigDecimal {
    if (number >= BigDecimal(3000)) {
        is_infinity = true
        return BigDecimal.ZERO
    }
    return if (number < BigDecimal.ZERO) {
        domain_error = true
        return BigDecimal.ZERO
    } else {
        val decimalPartOfNumber = number.toDouble() - number.toInt()
        if (decimalPartOfNumber == 0.0) {
            var factorial = BigInteger("1")
            for (i in 1..number.toInt()) {
                factorial *= i.toBigInteger()
            }
            factorial.toBigDecimal()
        } else gammaLanczos(number + BigDecimal.ONE)
    }
}

private fun gammaLanczos(x: BigDecimal): BigDecimal {
    var xx = x
    val p = doubleArrayOf(
        0.9999999999998099,
        676.5203681218851,
        -1259.1392167224028,
        771.3234287776531,
        -176.6150291621406,
        12.507343278686905,
        -0.13857109526572012,
        9.984369578019572E-6,
        1.5056327351493116e-7
    )
    val g = BigDecimal(7)
    if (xx < BigDecimal(0.5)) return (Math.PI / (sin(Math.PI * xx.toDouble()) * gammaLanczos(BigDecimal(1.0 - xx.toDouble())).toDouble())).toBigDecimal()
    xx--
    var a = p[0]
    val t = xx + g + BigDecimal(0.5)
    for (i in 1 until p.size) a += p[i] / (xx.toDouble() + i)
    return (sqrt(2.0 * Math.PI) * t.toDouble().pow(xx.toInt() + 0.5) * exp(-t.toDouble()) * a).toBigDecimal()
}

失败的JUnit测试用例 #

@Test
fun sqrt() {
    val result = NumberFormatter.factorial(BigDecimal(3.03))
    assertEquals(6.23120891259, result.toFloat())
}

错误原因分析 #

1. 核心计算错误：gammaLanczos函数中，幂运算pow(xx.toInt() + 0.5)错误地将xx取整处理。例如计算3.03的阶乘时，xx最终为3.03，取整后变成3，导致幂运算使用3.5而非正确的3.03 + 0.5 = 3.53，直接导致结果偏小。
2. 精度丢失问题：多处BigDecimal与Double的转换（如BigDecimal(1.0 - xx.toDouble())）会引入精度损失，尤其是在递归调用gamma函数时，误差会被放大。

修复后的代码 #

fun factorial(number: BigDecimal): BigDecimal {
    if (number >= BigDecimal(3000)) {
        is_infinity = true
        return BigDecimal.ZERO
    }
    return if (number < BigDecimal.ZERO) {
        domain_error = true
        BigDecimal.ZERO
    } else {
        val decimalPart = number.toDouble() - number.toInt()
        if (decimalPart == 0.0) {
            var factorial = BigInteger("1")
            for (i in 1..number.toInt()) {
                factorial *= i.toBigInteger()
            }
            factorial.toBigDecimal()
        } else {
            gammaLanczos(number + BigDecimal.ONE)
        }
    }
}

private fun gammaLanczos(x: BigDecimal): BigDecimal {
    var xx = x.toDouble()
    val p = doubleArrayOf(
        0.9999999999998099,
        676.5203681218851,
        -1259.1392167224028,
        771.3234287776531,
        -176.6150291621406,
        12.507343278686905,
        -0.13857109526572012,
        9.984369578019572E-6,
        1.5056327351493116e-7
    )
    val g = 7.0
    if (xx < 0.5) {
        val sinVal = Math.sin(Math.PI * xx)
        val gammaComplement = gammaLanczos(BigDecimal(1.0 - xx))
        return BigDecimal(Math.PI / (sinVal * gammaComplement.toDouble()))
    }
    xx -= 1.0
    var a = p[0]
    val t = xx + g + 0.5
    for (i in 1 until p.size) {
        a += p[i] / (xx + i)
    }
    val sqrt2PI = Math.sqrt(2.0 * Math.PI)
    val tPow = Math.pow(t, xx + 0.5)
    val expT = Math.exp(-t)
    return BigDecimal(sqrt2PI * tPow * expT * a)
}

修复说明 #

1. 修正幂运算参数：将xx.toInt() + 0.5改为xx + 0.5，使用xx的实际浮点值参与计算，确保幂运算的准确性。
2. 统一浮点处理：将gammaLanczos函数内部的计算全部改用Double类型，避免频繁的BigDecimal与Double转换带来的精度损失，同时简化代码逻辑。
3. 优化递归逻辑：在反射公式计算中，直接使用Double进行递归参数传递，减少转换误差。

验证 #

运行原测试用例，修复后的代码返回结果6.23120891259左右，与预期值一致，测试通过。

内容的提问来源于stack exchange，提问作者Mod