# Multi-Precision Floats: Radix Conversion with Large Exponents

[UPDATE] bugfix in code listing
The method for the radix conversion of floating point numbers described more theoretically here and more practically here has one large drawback: albeit always correct it is abysmally slow for large magnitudes.

Is there a faster way to do it? Yes, there is! Continue reading

Advertisements

# Adventures of a Programmer: Parser Writing Peril XXXX

The final grammar together with the lexer. The whole file with a lot of comments is on Github. It is currently in one file for convenience but needs to be split into a lexer (Flex) and a parser (Bison) later.

The file has a lot of comments already, not much to say here, but nevertheless… Continue reading

# Computation of the Arcsin/Arccos

The MacLaurin (Taylor series at x=0) series for the arcsin ${\displaystyle{ \arcsin(x) = \sum_{k=0}^{\infty} \frac{(2k-1)!{\kern-1.25pt}!}{(2k)!{\kern-1.25pt}!} \frac{x^{2k+1}}{2k+1} }}$

is a very slowly converging one, especially close to the branchpoints at -1 and 1.
The MacLaurin series for arccos is the same as for arcsin because of ${\displaystyle{ \arccos x = \frac{\pi}{2} - \arcsin x }}$

so it has the very same problem.

The series for the arctan on the other side is simpler to compute.
Another good reason: I have it already implemented 😉

The arcsin is related to the atan by ${\displaystyle{ \arcsin x = 2\arctan\biggl(\frac{x}{1 + \sqrt{1 - x^2} } \biggr) }}$

which seems to be of no help at all because the fraction does not get very small and the series for atan is also quite slow near -1 and 1. We can use another relation for this range. ${\displaystyle{ \arcsin x = \mathop{\mathrm{sgn}}(x)\arctan\Biggl(\sqrt{\frac{x^2}{1 - x^2 }} \Biggr) }}$

Seems not of much use but we can expand the root ${\displaystyle{ \arcsin x = \mathop{\mathrm{sgn}}(x)\arctan\biggl(\frac{x}{\sqrt{1 - x^2 }} \biggr) }}$

The closer x comes to the branchpoint the larger the fraction, but with ${\displaystyle{ \arctan x = \frac{\pi}{2} - \arctan \frac{1}{x} }}$

the actual value gets close to zero when x goes close to one and $\arctan 0 = 0$ so $\arcsin 1 = \frac{\pi}{2}$ which is correct and as intended. The point where $\frac{x}{\sqrt{1 - x^2 }} = 1$ is at $x = \frac{1}{\sqrt{2}}$ and at this point $\arctan 1 = \arcsin \frac{1}{\sqrt{2}} = \frac{\pi}{4}$.

My cutoff points will be $\arcsin = 2\arctan \bigl(\frac{x}{1 + \sqrt{1 - x^2} } \bigr)$ for $|x| \le 0.8$ and $\arctan\bigl(\frac{x}{\sqrt{1 - x^2 }} \bigr)$ for $|x| > 0.8$.

No code yet, just twiddling with math, sorry 😉

For $x \in\mathbb{R} \land |x| > 1$ the real part of the result is always $\frac{\pi}{2}$ and the imaginary part $\mathop{\mathrm{sgn}}(z) \log\bigl(x + \sqrt{x^2-1}\bigr)$. This comes from one of the definitions ${\displaystyle{ \arcsin x = -i\log\biggl(ix + \sqrt{1 - x^2 }\biggr) }}$

I will use this definition for complex arguments.

# Closed Form for n-th digit Generalization in Benford’s Law

The sum as described in the Wikipedia article has a closed form. It took me some time but after some twiddling with the sum itself I decided to exponentiate the whole thing to get rid of the logarithm and had a simple product I was able to handle.

With $d$ the digit and $b$ the base: ${\displaystyle{ \prod_{n=b^{k-2}}^{b^{k-1}-1}\left(1+\frac{1}{d+b n}\right) }}$

This can be expressed as a fraction of two rising factorials (Pochhammer symbol). ${\displaystyle{ \frac{\Bigl(\frac{b^k+d b+b}{b^2}\Bigr)_{b^{k-2}\left( b - 1\right)}}{\Bigl(\frac{b^k+d b}{b^2}\Bigr)_{b^{k-2}\left( b - 1\right) }} }}$

The rising factorial $(x)_n$ can be expressed as a fraction of two Gamma functions ${\displaystyle{ (x)_n = \frac{\Gamma(x + n)}{\Gamma(x)} }}$

With which we get ${\displaystyle{\frac{\Gamma\bigl(b^{k-2}+\frac{d}{b}\bigr) \Gamma\bigl(b^{k-1}+\frac{d}{b}+\frac{1}{b}\bigr)}{\Gamma\bigl(b^{k-2}+\frac{d}{b}+\frac{1}{b}\bigr) \Gamma\bigl(b^{k-1}+\frac{d}{b}\bigr)} }}$

By setting $d = b-1$ we can already see that the result will get smaller the larger $b$ gets with the obvious limit ${\displaystyle{ \lim_{b\to\infty} \frac{\Gamma\bigl(b^{k-2}+\frac{d}{b}\bigr) \Gamma\bigl(b^{k-1}+\frac{d}{b}+\frac{1}{b}\bigr)}{\Gamma\bigl(b^{k-2}+\frac{d}{b}+\frac{1}{b}\bigr) \Gamma\bigl(b^{k-1}+\frac{d}{b}\bigr)} = 0 \quad\text{with }d = b-1 \land d\in\mathbb{N} }}$

The arguments of the Gamma functions get quite large which makes the results hard to handle but we already need the logarithm of the product, that means we can use the logarithm of the Gamma function and save one step (we are on the real-line, so I used the equality $\log\Gamma(x) = \log\left(\Gamma\left(x\right)\right)$ to avoid an even denser parentheses fence than it already is) ${\displaystyle{ \Biggl(\log\Gamma\biggl(b^{k-2}+\frac{d}{b}\biggr) + \log\Gamma\biggl(b^{k-1}+\frac{d}{b}+\frac{1}{b}\biggr)\Biggr) - \Biggl(\log\Gamma\biggl(b^{k-2}+\frac{d}{b}+\frac{1}{b}\biggr) + \log\Gamma\biggl(b^{k-1}+\frac{d}{b}\biggr)\Biggr) }}$

We can simplify a bit with $\frac{d}{b} + \frac{1}{b} = 1$ because $d = b-1$. Furthermore $b^{k-2} + \frac{d}{b} = \frac{b^{k-1} + d}{b}$ which can save us a multiplication here and there. So far we got ${\displaystyle{ \Biggl(\log\Gamma\biggl(\frac{b^{k-1}+d}{b}\biggr) + \log\Gamma\biggl(b^{k-1}+1\biggr)\Biggr) - \Biggl(\log\Gamma\biggl(b^{k-2}+1\biggr) + \log\Gamma\biggl(b^{k-1}+\frac{d}{b}\biggr)\Biggr) }}$

or ${\displaystyle{ \Biggl(\log\Gamma\biggl(\frac{b^{k-1}+d}{b}\biggr) + \log\Gamma\biggl(b^{k-1}+1\biggr)\Biggr) - \Biggl(\log\Gamma\biggl(\frac{b^{k-1}+d + 1}{b}\biggr) + \log\Gamma\biggl(b^{k-1}+\frac{d}{b}\biggr)\Biggr) }}$

In the special case I used it for $b$ was a power of two such that none of the exponentiations had to be calculated, the numbers could be build in O(1).

Calling the whole thing $f(d,b)$ we can compute the probability that the nth digit is $b-1$ with $\frac{f(d,b)}{ \log(b)}$. For general $0 \le d < b$ use the unsimplified version above for $f(d,b)$.

Beware: the numbers get really large, you will need a slightly higher amount of precision for the loggamma function if you calculate with a large base. I had 226 as the base and had trouble to compute the values for places higher than the 20th in PARI/GP (it said *** lngamma: division by zero). It worked in calc with a precision of 300 decimal digits but it is sloooow. Who wrote that lame gamma function? Who? Me? Myself? Really? Ouch! 😉

# JavaScript: Make a Real typeof

The problems with JavaScript’s typeof operator are manifold. It is so restricted, that its only use is in testing for undefined. The main culprit cause these restrictions is the loose typing together with the automatic type coercion: "2" + 2 results in 22, "2" == 2 is true, t = 2;typeof t returns "number" and t = new Number(2);typeof t returns "object" to name just a few.
The worst case is probably:

// somewhere on top of the code
var a = 2;
// some thousand lines and/or several scripts later
var b = new Number(2);
if(a === b){
// CPR: Cardiopulmonary resuscitation
console.log("Go on with CPR!");
} else {
console.log("He's dead, Jim!"),
}


# Conversion of Binary Floating Point Numbers—Second Try

The conversion of a string to a binary encoded floating point number is not easy, as shown in my last post. The information given there is a bit too much on the theoretical side. I’ll try to change that here with some real code: String to Bigfloat and Bigfloat to String. Continue reading

# Factoring an Integer in JavaScript with a Wheel

Update: bugfix in code

Factoring an integer with a wheel is a known algorithm for which many examples exist for many programming languages. Google might be of some help to find them.
Most of them are either way too simple or way too complicated; I would even say that all of them fall into one of these two categories but I do not know all of them.
One example that could fit both categories at once is the program factor from the GNU-coreutils which uses a wheel to do its work. I did a JavaScript version of it if that language suits you more. The program itself is quite simple yet hard to understand.
For example: what are these numbers at the start and where do they come from? Continue reading