Factorial calculation giving wrong answer with float or decimal number

Question

this code is working fine to calculate factorial,but giving wrong answer with float or decimal number

    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 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 {
        // https://rosettacode.org/wiki/Gamma_function
        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()
    }

For the code below, the test is failing after running nine the code

@Test
    fun sqrt() {

        val result =  NumberFormatter.factorial(BigDecimal(3.03))
        assertEquals(6.23120891259, result.toFloat())
    }

got this expected:<6.23120891259> but was:<5.8062997>

Could anyone help with this. Thanks

Answer 1

The code has precision loss due to the conversion between BigDecimal and double. you can update the gammaLanczos to

fun gammaLanczos(x: BigDecimal): BigDecimal {
// Lanczos approximation parameters
   val p = arrayOf(
       676.5203681218851,
      -1259.1392167224028,
      771.3234287776531,
      -176.6150291621406,
      12.507343278686905,
      -0.13857109526572012,
      9.984369578019572e-6,
      1.5056327351493116e-7
  )
  val g = 7.0
  var z = x.toDouble() - 1.0

  var a = 0.99999999999980993
  for (i in p.indices) {
      a += p[i] / (z + i + 1)
  }

  val t = z + g + 0.5
  val sqrtTwoPi = sqrt(2.0 * PI)
  val firstPart = sqrtTwoPi * t.pow(z + 0.5) * exp(-t)
  val result = firstPart * a

  return BigDecimal(result, MathContext.DECIMAL64)
}

fun BigDecimal.pow(exponent: Double): BigDecimal {
   return this.toDouble().pow(exponent).toBigDecimal(MathContext.DECIMAL64)
}

Answer 2

An alternative approach would be to use the Apache commons math library. The version three docs for Gamma are here.

Sample code to calculate 3.03!

Gamma.gamma(1 + 3.03)
    .also { check(it == 6.231208912587357 ) }

Version four has been in the works for a while. If you are interested, the current Java source is here.