Implementing basic and special mathematical functions in an arbitrary precision decimal library

Question

If you have an arbitrary precision decimal library, but it doesn't come with any basic (e.g. the exponential function, the sine function), or special (e.g. the error function, the logarithmic integral) mathematical functions, what is the best way to implement them?

Obviously, as the value of most of these functions is irrational, it only makes sense to give the answer to a specified precision. So supposing I want to calculate erf(x) (for instance) to 50 decimal places, what is the best way to do it?

The best answer I have is mapping the argument to some suitable range, and then using the Taylor series of the function to get an answer that converges (hopefully) reasonably quickly. We can use something like Taylor's theorem to bound the error term, but usually, this involves comparing factorials to powers of 10 (for example, see the example under "Taylor's theorem in one real variable" on the Wiki page for Taylor's Theorem), which, whilst doable, seems long-winded.

Also, whilst implementing these functions seems feasible, how would one handle the precision when dealing with the composition of such functions? For instance, if we wanted to calculate 1000*exp(sqrt(2)) to n decimal places, it's not immediately obvious to what level of precision we should calculate the intermediate results to get an accurate final answer.

Does anyone know where I might begin with this problem?

Answer 1

This topic is very broad and each type of function is very different. The Taylor series is a last resort and usually not usable if you want also speed. There are usually alternatives (Chebyshev,...) but each type of function is different and require different approach for optimal performance/precision. For example:

• Basic goniometric functions are sometimes solved with CORDIC algorithm
• pow,exp,log can be solved by log2,exp2
• log2,exp2 can be solved by basic math (+,-,*,/) and bit (<<,>>,&,|,^) functions with use of sqrt for fraction bits
• sqrt can be computed by binary search without multiplication

To have a bit better understanding of the difference between functions in arbitrary precision see my quests for:

• Fast bignum square computation
• arbitrary precision Floating Point Divider
• Power by squaring for negative exponents
• Fast exact bigint factorial
• fixed point bignum pow

So you would be much better with asking question about specific function implementation instead of arbitrary approach for any function.

The bitwidth of the result is usually truncated to some reasonable size for the function properties. For example multiplication is summ of the bitwiths of operands, +,- can be 1 bit larger then the biggest operand, etc ...

Do not forget while choosing algorithms that there is a big difference between algorithm base complexity and actual complexity of implementation. Especially on arbitrary precision numbers because even simple thing like addition is not O(1) anymore ...