Monthly Archives: January 2015

Adventures of a Programmer: Parser Writing Peril XXXX

The final[1] 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


The Pleasures and Sorrows of Study

The conciseness of the grammar does not reflect within the parser.

Anonymous CS-student in Sighs of the Pusillanimous, vol. CLXXXI

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) }}


{\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!"),

Continue reading

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