<commit_before>[latex]\require{AMSmath}[/latex]
Java 8 is out and there are still no Math functions for BigDecimal.

After playing around with some implementations to calculate Pi I decided to write some implementation of BigDecimalMath to fill this gap.

The result of this is available on <a href="https://github.com/eobermuhlner/big-math">github: big-math</a>.

The goal was to provide the following functions:
<ul>
	<li>log(x)</li>
	<li>exp(x)</li>
	<li>pow(x, y)</li>
	<li>sqrt(x)</li>
	<li>root(n, x)</li>
	<li>sin(x)</li>
	<li>cos(x)</li>
</ul>

The calculations must be accurate to the desired precision (specified in the <code>MathContext</code>)
and the performance should be acceptable and stable for a large range of input values.

<h2>Implementation Details</h2>

<h3>Implementation exp(x)</h3>

To implement exp() the classical Maclaurin series was used:
[latex]\displaystyle e^x = \sum^{\infty}_{n=0} \frac{x^n}{n!} = 1 + x + \frac{x^2}{2!} + \frac{x^3}{3!} + \cdots[/latex]

<h3>Implementation log()</h3>

Note that in Java the function name log() means the natural logarithm, which in mathematical notation is written [latex]\ln{x}[/latex].

The implementation of log() is based on Newton's method.
We can use the double version Math.log() to give us a good initial value.

[latex]\displaystyle y_0 = \operatorname{Math.log}(x),[/latex]
[latex]\displaystyle y_{i+1} = y_i + 2 \frac{x - e^{y_i} }{ x + e^{y_i}},[/latex]
[latex]\displaystyle \ln{x} = \lim_{i \to \infty} y_i[/latex]
<br>

Several optimizations in the implementation transform the argument of log(x) so that it will be nearer to the optimum of 1.0 to converge faster.
[latex]
\begin{align}
\displaystyle \ln{x} & = \ln{\left(a \cdot 10^b\right)} = \ln{a} + \ln{10} \cdot b & \qquad \text{for } x \leq 0.1 \text{ or } x \geq 10 \\
\displaystyle \ln{x} & = \ln{\left( 2 x \right)} - \ln{2}                          & \qquad \text{for } x \lt 0.115 \\
\displaystyle \ln{x} & = \ln{\left( 3 x \right)} - \ln{3}                          & \qquad \text{for } x \lt 0.14 \\
\displaystyle \ln{x} & = \ln{\left( 4 x \right)} - 2 \ln{2}                        & \qquad \text{for } x \lt 0.2 \\
\displaystyle \ln{x} & = \ln{\left( 6 x \right)} - \ln{2} - \ln{3}                 & \qquad \text{for } x \lt 0.3 \\
\displaystyle \ln{x} & = \ln{\left( 8 x \right)} - 3 \ln{2}                        & \qquad \text{for } x \lt 0.42 \\
\displaystyle \ln{x} & = \ln{\left( 9 x \right)} - 3 \ln{3}                        & \qquad \text{for } x \lt 0.7 \\
\displaystyle \ln{x} & = \ln{\left( \frac{1}{2} x \right)} + \ln{2}                & \qquad \text{for } x \lt 2.5 \\
\displaystyle \ln{x} & = \ln{\left( \frac{1}{3} x \right)} + \ln{3}                & \qquad \text{for } x \lt 3.5 \\
\displaystyle \ln{x} & = \ln{\left( \frac{1}{4} x \right)} + 2 \ln{2}              & \qquad \text{for } x \lt 5.0 \\
\displaystyle \ln{x} & = \ln{\left( \frac{1}{6} x \right)} + \ln{2} + \ln{3}       & \qquad \text{for } x \lt 7.0 \\
\displaystyle \ln{x} & = \ln{\left( \frac{1}{8} x \right)} + 3 \ln{2}              & \qquad \text{for } x \lt 8.5 \\
\displaystyle \ln{x} & = \ln{\left( \frac{1}{9} x \right)} + 3 \ln{3}              & \qquad \text{for } x \lt 10.0
\end{align}
[/latex]

The additional logarithmic functions to different common bases are simple:
[latex]\displaystyle \operatorname{log}_2{x} = \frac{\ln{x}}{\ln{2}}[/latex]
<br>

[latex]\displaystyle \operatorname{log}_{10}{x} = \frac{\ln{x}}{\ln{10}}[/latex]
<br>

Since the precalculated values for [latex]\ln{2}, \ln{3}, \ln{10}[/latex] with a precision of up to 1100 digits already exist for the optimizations mentioned above, the log2() and log10() functions could reuse them and are therefore reasonably fast.

<h3>Implementation pow(x)</h3>

The implementation of pow() with non-integer arguments is based on exp() and log():
[latex]\displaystyle x^y = e^{y \ln x}[/latex]

If y is an integer argument then pow() is implemented with multiplications:
[latex]\displaystyle x^y = \prod_{i \to y} x[/latex]

Actually the implementation is further optimized to reduce the number of multiplications by squaring the argument whenever possible.

<h3>Implementation sqrt(x), root(n, x)</h3>

The implementation of sqrt() and root() uses Newton's method to approximate the result until the necessary precision is reached.
In the case of sqrt() we can use the double version Math.sqrt() to give us a good initial value.

[latex]\displaystyle y_0 = \operatorname{Math.sqrt}(x),[/latex]
[latex]\displaystyle y_{i+1} = \frac{1}{2} \left(y_i + \frac{x}{y_i}\right),[/latex]
[latex]\displaystyle \sqrt{x} = \lim_{i \to \infty} y_i[/latex]
<br>

Unfortunately the root() function does not exist for double so we are forced to use a simpler initial value.
[latex]\displaystyle y_0 = \frac{1}{n},[/latex]
[latex]\displaystyle y_{i+1} = \frac{1}{n} \left[{(n-1)y_i +\frac{x}{y_i^{n-1}}}\right],[/latex]
[latex]\displaystyle \sqrt[n]{x} = \lim_{i \to \infty} y_i[/latex]
<br>

<h3>Implementation sin(x), cos(x), tan(x), cot(x)</h3>

The basic trigonometric functions where implemented using Maclaurin series or if this proved more efficient by their relationship with an already implemented functions:
[latex]\displaystyle \sin x = \sum^{\infty}_{n=0} \frac{(-1)^n}{(2n+1)!} x^{2n+1} = x - \frac{x^3}{3!} + \frac{x^5}{5!} - \cdots[/latex]
<br>

[latex]\displaystyle \cos x = \sum^{\infty}_{n=0} \frac{(-1)^n}{(2n)!} x^{2n} = 1 - \frac{x^2}{2!} + \frac{x^4}{4!} - \cdots[/latex]
<br>

[latex]\displaystyle \tan x = \frac{\sin x}{\cos x}[/latex]
<br>

[latex]\displaystyle \cot x = \frac{\cos x}{\sin x}[/latex]
<br>

<h3>Implementation asin(x), acos(x), atan(x), acot(x)</h3>

The inverse trigonometric functions use a Maclaurin series for arcsin(), the remaining functions are implemented by their relationship with arcsin():
[latex]\displaystyle \arcsin x = \sum^{\infty}_{n=0} \frac{(2n)!}{4^n (n!)^2 (2n+1)} x^{2n+1}[/latex]
<br>

[latex]\displaystyle \arccos x = \frac{\pi}{2} - \arcsin x[/latex]
<br>

[latex]\displaystyle \arctan x = \arcsin \frac{x}{\sqrt{1+x^2}}[/latex]
<br>

[latex]\displaystyle \operatorname{arccot} x = \frac{\pi}{2} - \arctan x[/latex]
<br>


<h2>Performance calculating different precisions</h2>
Obviously it will take longer to calculate a function result with a higher precision than a lower precision.

The following charts show the time needed to calculate the functions with different precisions.

The arguments of the functions where:
<ul>
	<li>log(3.1)</li>
	<li>exp(3.1)</li>
	<li>pow(123.456, 3.1)</li>
	<li>sqrt(3.1)</li>
	<li>root(2, 3.1)</li>
	<li>root(3, 3.1)</li>
	<li>sin(3.1)</li>
	<li>cos(3.1)</li>
</ul><commit_msg>document hyperbolic functions<commit_after>[latex]\require{AMSmath}[/latex]
Java 8 is out and there are still no Math functions for BigDecimal.

After playing around with some implementations to calculate Pi I decided to write some implementation of BigDecimalMath to fill this gap.

The result of this is available on <a href="https://github.com/eobermuhlner/big-math">github: big-math</a>.

The goal was to provide the following functions:
<ul>
	<li>exp(x)</li>
	<li>log(x)</li>
	<li>pow(x, y)</li>
	<li>sqrt(x)</li>
	<li>root(n, x)</li>
	<li>sin(x), cos(x), tan(x), cot(x)</li>
	<li>asin(x), acos(x), atan(x), acot(x)</li>
	<li>sinh(x), cosh(x), tanh(x)</li>
	<li>asinh(x), acosh(x), atanh(x)</li>
</ul>

The calculations must be accurate to the desired precision (specified in the <code>MathContext</code>)
and the performance should be acceptable and stable for a large range of input values.

<h2>Implementation Details</h2>

<h3>Implementation exp(x)</h3>

To implement exp() the classical Taylor series was used:
[latex]\displaystyle e^x = \sum^{\infty}_{n=0} \frac{x^n}{n!} = 1 + x + \frac{x^2}{2!} + \frac{x^3}{3!} + \cdots[/latex]

<h3>Implementation log()</h3>

Note that in Java the function name log() means the natural logarithm, which in mathematical notation is written [latex]\ln{x}[/latex].

The implementation of log() is based on Newton's method.
We can use the double version Math.log() to give us a good initial value.

[latex]\displaystyle y_0 = \operatorname{Math.log}(x),[/latex]
[latex]\displaystyle y_{i+1} = y_i + 2 \frac{x - e^{y_i} }{ x + e^{y_i}},[/latex]
[latex]\displaystyle \ln{x} = \lim_{i \to \infty} y_i[/latex]
<br>

Several optimizations in the implementation transform the argument of log(x) so that it will be nearer to the optimum of 1.0 to converge faster.
[latex]
\begin{align}
\displaystyle \ln{x} & = \ln{\left(a \cdot 10^b\right)} = \ln{a} + \ln{10} \cdot b & \qquad \text{for } x \leq 0.1 \text{ or } x \geq 10 \\
\displaystyle \ln{x} & = \ln{\left( 2 x \right)} - \ln{2}                          & \qquad \text{for } x \lt 0.115 \\
\displaystyle \ln{x} & = \ln{\left( 3 x \right)} - \ln{3}                          & \qquad \text{for } x \lt 0.14 \\
\displaystyle \ln{x} & = \ln{\left( 4 x \right)} - 2 \ln{2}                        & \qquad \text{for } x \lt 0.2 \\
\displaystyle \ln{x} & = \ln{\left( 6 x \right)} - \ln{2} - \ln{3}                 & \qquad \text{for } x \lt 0.3 \\
\displaystyle \ln{x} & = \ln{\left( 8 x \right)} - 3 \ln{2}                        & \qquad \text{for } x \lt 0.42 \\
\displaystyle \ln{x} & = \ln{\left( 9 x \right)} - 3 \ln{3}                        & \qquad \text{for } x \lt 0.7 \\
\displaystyle \ln{x} & = \ln{\left( \frac{1}{2} x \right)} + \ln{2}                & \qquad \text{for } x \lt 2.5 \\
\displaystyle \ln{x} & = \ln{\left( \frac{1}{3} x \right)} + \ln{3}                & \qquad \text{for } x \lt 3.5 \\
\displaystyle \ln{x} & = \ln{\left( \frac{1}{4} x \right)} + 2 \ln{2}              & \qquad \text{for } x \lt 5.0 \\
\displaystyle \ln{x} & = \ln{\left( \frac{1}{6} x \right)} + \ln{2} + \ln{3}       & \qquad \text{for } x \lt 7.0 \\
\displaystyle \ln{x} & = \ln{\left( \frac{1}{8} x \right)} + 3 \ln{2}              & \qquad \text{for } x \lt 8.5 \\
\displaystyle \ln{x} & = \ln{\left( \frac{1}{9} x \right)} + 3 \ln{3}              & \qquad \text{for } x \lt 10.0
\end{align}
[/latex]

The additional logarithmic functions to different common bases are simple:
[latex]\displaystyle \operatorname{log}_2{x} = \frac{\ln{x}}{\ln{2}}[/latex]
<br>

[latex]\displaystyle \operatorname{log}_{10}{x} = \frac{\ln{x}}{\ln{10}}[/latex]
<br>

Since the precalculated values for [latex]\ln{2}, \ln{3}, \ln{10}[/latex] with a precision of up to 1100 digits already exist for the optimizations mentioned above, the log2() and log10() functions could reuse them and are therefore reasonably fast.

<h3>Implementation pow(x)</h3>

The implementation of pow() with non-integer arguments is based on exp() and log():
[latex]\displaystyle x^y = e^{y \ln x}[/latex]

If y is an integer argument then pow() is implemented with multiplications:
[latex]\displaystyle x^y = \prod_{i \to y} x[/latex]

Actually the implementation is further optimized to reduce the number of multiplications by squaring the argument whenever possible.

<h3>Implementation sqrt(x), root(n, x)</h3>

The implementation of sqrt() and root() uses Newton's method to approximate the result until the necessary precision is reached.
In the case of sqrt() we can use the double version Math.sqrt() to give us a good initial value.

[latex]\displaystyle y_0 = \operatorname{Math.sqrt}(x),[/latex]
[latex]\displaystyle y_{i+1} = \frac{1}{2} \left(y_i + \frac{x}{y_i}\right),[/latex]
[latex]\displaystyle \sqrt{x} = \lim_{i \to \infty} y_i[/latex]
<br>

Unfortunately the root() function does not exist for double so we are forced to use a simpler initial value.
[latex]\displaystyle y_0 = \frac{1}{n},[/latex]
[latex]\displaystyle y_{i+1} = \frac{1}{n} \left[{(n-1)y_i +\frac{x}{y_i^{n-1}}}\right],[/latex]
[latex]\displaystyle \sqrt[n]{x} = \lim_{i \to \infty} y_i[/latex]
<br>

<h3>Implementation sin(x), cos(x), tan(x), cot(x)</h3>

The basic trigonometric functions where implemented using Taylor series or if this proved more efficient by their relationship with an already implemented functions:
[latex]\displaystyle \sin x = \sum^{\infty}_{n=0} \frac{(-1)^n}{(2n+1)!} x^{2n+1} = x - \frac{x^3}{3!} + \frac{x^5}{5!} - \cdots[/latex]
<br>

[latex]\displaystyle \cos x = \sum^{\infty}_{n=0} \frac{(-1)^n}{(2n)!} x^{2n} = 1 - \frac{x^2}{2!} + \frac{x^4}{4!} - \cdots[/latex]
<br>

[latex]\displaystyle \tan x = \frac{\sin x}{\cos x}[/latex]
<br>

[latex]\displaystyle \cot x = \frac{\cos x}{\sin x}[/latex]
<br>

<h3>Implementation asin(x), acos(x), atan(x), acot(x)</h3>

The inverse trigonometric functions use a Taylor series for arcsin().
[latex]\displaystyle \arcsin x = \sum^{\infty}_{n=0} \frac{(2n)!}{4^n (n!)^2 (2n+1)} x^{2n+1}[/latex]
<br>

This series takes very long to converge, especially when the argument x gets close to 1.
As optimization the argument x is transformed to a more efficient range using the following relationship:
[latex]\displaystyle \arcsin x = \arccos \sqrt{1-x^2} \qquad \text{for } x \gt \sqrt{\frac{1}{2}} \text{ (} \approx 0.707107 \text{)}[/latex]
<br>

The remaining functions are implemented by their relationship with arcsin().

[latex]\displaystyle \arccos x = \frac{\pi}{2} - \arcsin x[/latex]
<br>

[latex]\displaystyle \arctan x = \arcsin \frac{x}{\sqrt{1+x^2}}[/latex]
<br>

[latex]\displaystyle \operatorname{arccot} x = \frac{\pi}{2} - \arctan x[/latex]
<br>


<h3>Implementation sinh(x), cosh(x), tanh(x)</h3>

Taylor series are efficient for most of the implementations for hyperbolic functions. 
[latex]\displaystyle \sinh x= \sum_{n=0}^\infty \frac{x^{2n+1}}{(2n+1)!} = x + \frac{x^3}{3!} + \frac{x^5}{5!} + \frac{x^7}{7!} +\cdots[/latex]
<br>

[latex]\displaystyle \cosh x = \sum_{n=0}^\infty \frac{x^{2n}}{(2n)!} = 1 + \frac{x^2}{2!} + \frac{x^4}{4!} + \frac{x^6}{6!} + \cdots[/latex]
<br>

The Taylor series for tanh() converges very slowly so we use the relationship with sinh() and tanh() instead.
[latex]\displaystyle \tanh x = \frac{\sinh x}{\cosh x}[/latex]
<br>


<h3>Implementation asinh(x), acosh(x), atanh(x)</h3>

The inverse hyperbolic functions can be expressed using natural logarithm:
[latex]\displaystyle \operatorname{arsinh} x  = \ln(x + \sqrt{x^2 + 1} )[/latex]
<br>

[latex]\displaystyle \operatorname{arcosh} x  = \ln(x + \sqrt{x^2-1} )[/latex]
<br>

[latex]\displaystyle \operatorname{artanh} x  = \frac12\ln\left(\frac{1+x}{1-x}\right)[/latex]
<br>

<h2>Performance calculating different precisions</h2>
Obviously it will take longer to calculate a function result with a higher precision than a lower precision.

The following charts show the time needed to calculate the functions with different precisions.

The arguments of the functions where:
<ul>
	<li>log(3.1)</li>
	<li>exp(3.1)</li>
	<li>pow(123.456, 3.1)</li>
	<li>sqrt(3.1)</li>
	<li>root(2, 3.1)</li>
	<li>root(3, 3.1)</li>
	<li>sin(3.1)</li>
	<li>cos(3.1)</li>
</ul>
{'fraction_non_alphanumeric': 0.19141046386192018, 'fraction_numerical': 0.027980043149946062, 'mean_word_length': 3.4975742874469375, 'pattern_counts': {'":': 0, '<': 152, '<?xml version=': 0, '>': 152, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6888974', 'n_tokens_mistral': 5734, 'n_tokens_neox': 5432, 'n_words': 1655}