Support for trigonometric functions

[breandan]

Describe the solution you'd like: Support for trigonometric functions, e.g. sin cos tan. Can be implemented using the MacLaurin series.

Describe alternatives you've considered: https://github.com/eobermuhlner/kotlin-big-math

[hakanai]

I had a look at what this might look like in Kotlin and got a mostly-working implementation for sin(x), but the one place I'm not 100% sure of is rounding.

Because it's Kotlin, using lazy sequences seemed like an appropriate way to deal with infinite series.

The general gist:

abstract class SeriesCalculator {
    fun calculate(x: BigDecimal, decimalMode: DecimalMode): BigDecimal {
        val oneMoreDigitDecimalMode = decimalMode.copy(decimalPrecision = decimalMode.decimalPrecision + 1)
        val epsilon = BigDecimal.ONE.moveDecimalPoint(-oneMoreDigitDecimalMode.decimalPrecision)

        val powerSequence = createPowerSequence(x, decimalMode)
        // TODO: This sequence can be cached, but I'm not entirely sure how just yet
        val factorSequence = createFactorSequence(decimalMode)

        return powerSequence.zip(factorSequence)
            .map { (a, b) -> a.multiply(b, oneMoreDigitDecimalMode) }
            .takeWhile { step -> step.abs() > epsilon }
            .fold(BigDecimal.ZERO) { acc, step -> acc.add(step, oneMoreDigitDecimalMode) }
            .roundSignificand(decimalMode)
    }

    /**
     * Implemented by subclasses to create an appropriate sequence of powers for the series.
     */
    abstract fun createPowerSequence(x: BigDecimal, decimalMode: DecimalMode): Sequence<BigDecimal>

    /**
     * Implemented by subclasses to create an appropriate sequence of factors for the series.
     */
    abstract fun createFactorSequence(decimalMode: DecimalMode): Sequence<BigDecimal>

    /**
     * Utility function for subclasses to use which produces a sequence of all powers of `x`,
     * starting at 0.
     */
    protected fun createAllPowersSequence(x: BigDecimal): Sequence<BigDecimal> {
        return generateSequence(BigDecimal.ONE) { n -> n * x }
    }

    /**
     * Utility function for subclasses to use which produces a sequence of all factorials,
     * starting at 0! = 1, then 1! = 1, 2! = 2, 3! = 6, etc.
     */
    protected fun createAllFactorialsSequence(): Sequence<BigDecimal> {
        return sequence {
            var n = BigDecimal.ONE
            var i = 0
            while (true) {
                yield(n)
                i++
                n *= i
            }
        }
    }
}

And then for sin(x):

object SinCalculator : SeriesCalculator() {
    override fun createPowerSequence(x: BigDecimal, decimalMode: DecimalMode): Sequence<BigDecimal> {
        return createAllPowersSequence(x)
            .filterIndexed { i, _ -> i % 2 != 0 }
    }

    override fun createFactorSequence(decimalMode: DecimalMode): Sequence<BigDecimal> {
        // XXX: No constant for -1 and there's probably a better way to do this sign flipping anyway
        var sign = BigDecimal.ZERO - BigDecimal.ONE
        return createAllFactorialsSequence()
            .filterIndexed { i, _ -> i % 2 != 0 }
            .map { n ->
                sign = -sign
                sign * BigDecimal.ONE.divide(n, decimalMode)
            }
    }
}

Test:

class BigDecimalTrigTest {

    @Test
    fun sinTest() {
        val decimalMode = DecimalMode(decimalPrecision = 20, roundingMode = RoundingMode.ROUND_HALF_TO_EVEN)

        assertEquals(
            BigDecimal.parseStringWithMode("0", decimalMode),
            BigDecimal.parseStringWithMode("0", decimalMode).sin()
        )

        assertEquals(
            BigDecimal.parseStringWithMode("0.19866933079506121545", decimalMode),
            BigDecimal.parseStringWithMode("0.2", decimalMode).sin()
        )

        assertEquals(
            BigDecimal.parseStringWithMode("0.38941834230865049166", decimalMode),
            BigDecimal.parseStringWithMode("0.4", decimalMode).sin()
        )

        assertEquals(
            BigDecimal.parseStringWithMode("0.56464247339503535720", decimalMode),
            BigDecimal.parseStringWithMode("0.6", decimalMode).sin()
        )
    }
}

The result is basically correct:

java.lang.AssertionError: expected:<1.9866933079506121545E-1> but was:<1.9866933079506121546E-1>

The expected values were determined using bc, so I'm sure the difference is going to be something subtle about what rounding mode bc is using vs. what I chose for the test.

I also wasn't entirely sure how to properly deal with rounding mode in the code itself so I ended up making it do the maths with one digit more precision than necessary, and then rounding at the end. Whether one digit extra is sufficient, I have no idea. Without that, the actual result turned out to be 40 digits of precision. :)

So open issues with my attempt:

1. Caching the factor sequence would be ideal because it never changes
2. Sort out precision / rounding mode handling
3. Figure out what to do if the initial number doesn't specify the rounding mode (add a DecimalMode param to sin itself?)

[ionspin]

Hi @hakanai , thanks for effort, hopefully I'll have some time over weekend to look into this more, but no guarantees. If you make more progress and get a pull request up that would be great!

[hakanai]

I was trying to figure out whether Kotlin had a thread-safe way to cache the values, but it really seems like the core Kotlin class library pretends that multithreading doesn't exist, and to get any kind of synchronisation you have to depend on coroutines or similar. I don't think they considered the possibility that library code might want to be thread-safe while not using multithreading itself!

And since this library is currently dependency-free, it may be fine to just forget about caching the sequence for now. The need for DecimalMode makes it complicated anyway. That other library gets to avoid the problem because they have a BigRational class which doesn't need precision.

The code to sum the sequence can still be shared between multiple sequence classes. I wrote cos and exp and the story is similar for those.

At this point the main issue is figuring out why the expected results differ from the actual results even if I increase the precision. It could be that the extra precision we need to use when computing the sum is much higher than just one more decimal place.

[hakanai]

I created #235 with exp, sin, cos, tan, sinh, cosh, tanh, which is about half of what you'd want.

From some light reading, it looks like the vast majority of what's left depends on having ln and sqrt.