Category Archives: Tutorial

Adventures of a Programmer: Parser Writing Peril I

As menacedpromised in the last post: a description of the long path
to a well determined end.

Goal: Writing a small parser for a simple but complete language as a front-end for a for a calculator back-end. With emphasis on “complete”. A simple task one might think, a very simple one. One that could have been described in full in some book I mentioned in another post.

So in medias res, as Caesar wouldn’t have said (he spoke Greek) and of to the very beginning: The Plan. Not that I ever saw it but it is said that The Plan cometh at every beginning.
What is needed for the already briefly described…uhm…purpose, to avoid the repetition of the infamous two words, that are unbeknownst to many programmers, what does the parser need to parse? What do we need to do to enable the poor thing to do so? Continue reading

Advertisements

On Doing It Yourself

When the output of Mathematica® is correct but of not much use

When I do some calculations on the command line, it is in most cases some quick computations involving numbers with different bases in a precision that is not very high but higher than what the normal calculator offers. Starting up PARI/GP every time is distracting and PARI/GP does not offer an easy way to enter numbers in different bases; dc and bc are slightly more complicated to use—although they are quite excellent in shell-scripts!—and all others are either insufficient in power, usage and/or price so it is calc what I use regularly.

It has some downsides, too. One of these disadvantages is the lack of functions that are included in almost every other arbitrary precision calculator because calc had been written for a special need originally. See the webpage of calc I linked to above for further information.

Another advantage of calc is that the syntax is quite C-like and allows for easy implementation of the functions and algorithms you need. At least that’s what I hope for 😉

The most common implementation of an approximation of the gamma function when arbitrary precision is needed, is the variation of James Stirling’s approximation found by John L. Spouge[1]. It is some time ago, that I implemented it and I had to look it up. The first page that came up after asking Google for some help, was the page of the entry of Literate Prgrams for an implementation in Mathematica®. No, I really cannot remember what I’ve asked Google to get that result. I found out later, that Rosetta Code has a rudimentary but working example implementation in C.
My Google-Foo is really lacking…

Now, where do I am going to get to with such a long-winded introduction, you might ask, and I will tell you: on the Literate-Programming page I found a short explanation of the algorithm which includes the following gem:

Conservatively omitting the factor a^{-1/2} and rearranging gives us the formula…

They talk about the error and the number of coefficients needed:

(1.0) \displaystyle{ \epsilon \le a^{-1/2} (2 \pi)^{-(a+1/2)}}

for which they give the approximation of the parameter a as:

(1.1) \displaystyle{  a = -\frac{\log \epsilon}{\log 2 \pi}}

If you enter the inequality into the text field over at wolframalpha.com and tell it to solve for a you’ll get an exact but unusable result. Changing the inequality to an equality you’ll get another, simpler but still unusable result involving the Lambert-W function. If you want to play with it, here is a JavaScript implementation of the approximation on the positive real line (click on it to unfold it):

Math.lambertW = function(x){
  /*
    On the Lambert W Function (1996) 
    by R. M. Corless , G. H. Gonnet , D. E. G. Hare , D. J. Jeffrey , D. E. Knuth 
    ADVANCES IN COMPUTATIONAL MATHEMATICS
  */
  var ret = 0.0;
  var a   = Math.log(x);   //L1
  var b   = Math.log(Math.log(x));//L2
  var c   = a-b + b/a;
  var l1  = (b*(-2 + b))/(2*a*a);
  var l2  = (b*(6-9*b+2*b*b))/(6*a*a*a);
  var l3  = (b*(-12+36*b-22*b*b+3*b*b*b))/(12*a*a*a*a);
  var l4  = (b*(60-300*b+350*b*b-125*b*b*b+12*b*b*b*b))/(60*a*a*a*a*a*a);

  return c+l1+l2+l3+l4;
};

So I decided to do it by hand. Why? Well, if you are seriously questioning it read no further, for what I do is not for you.
But I do not use anything not taught in highschool and even give the needed rules for the logarithms I use here, so you are heartily invited to give it a glimpse. You might even find an error 😉

At first make the inequality an equality because we are interested in the largest possible error only. Expand all of the exponents to make it more readable. For me at least but YMMV, as the saying goes.

(2.1) \displaystyle{ \epsilon = \frac{1}{\sqrt{a}\,\pi^{a+1/2}}}

The term \sqrt{a} can be omitted here, that is indeed correct. I’ll give a short sketch of a proof at the end.

(2.2) \displaystyle{ \epsilon = \frac{1}{2\pi^{a+\frac{1}{2}}}}

To get rid of the exponent we have to take the logarithm of both sides. Any logarithm will do.

(2.3) \displaystyle{ \log\epsilon = \log\left(\frac{1}{2\pi^{a+\frac{1}{2}}}\right)}

Expand it, such that we get the exponents back we dropped so prematurely above.
Rules: \log (a/b) = \log a - \log b and \log (a^b) = log(a)\, b

(2.4) \displaystyle{ \log\epsilon = \log{1}-\log(2\pi) \left(a+\frac{1}{2}\right)}

The logarithm of 1 (one) is 0 (zero) and can be omitted for brevity. Now move the constant -\log{2\pi} to the left of the equation.

(2.5) \displaystyle{ -\frac{\log\epsilon}{\log(2\pi)} = a + \frac{1}{2}}

Subtract \frac{1}{2} and we are done so far.

(2.6) \displaystyle{-\frac{log\epsilon}{\log(2\pi)} - \frac{1}{2} = a}

Because of the conditions a \in \mathbb{N} together with a > 2 we can omit the \frac{1}{2}, too.

What? Uh, the sketch for the proof for dropping the poor \sqrt{a}, yes.
We start again from equation (2.1). To make it more readable we move the square root to the left to get

(3.1) \displaystyle{ \epsilon\sqrt{a} = \frac{1}{2\pi^{a+\frac{1}{2}}}}

Take the logarithm of both sides and expand
Rule: \log(\sqrt[n]{a}) = \frac{\log a}{n} which is just a special case of \log (a^b) = log(a)\, b

(3.2) \displaystyle{ \log\epsilon + \frac{\log a}{2}  = -\log (2\pi) (a+\frac{1}{2}) }

Uh, no, move \sqrt{a} back to the right, sorry.

(3.3) \displaystyle{ \log\epsilon  = -\log (2\pi) \left(a+\frac{1}{2}\right) - \frac{\log a}{2} }

Divide by -\log (2\pi) as in equation (2.5) and expand

(3.4) \displaystyle{ -\frac{\log\epsilon}{\log (2\pi)}  = -\frac{2}{\log a\log(2\pi)} + a+\frac{1}{2} }

Subtract \frac{1}{2} and rearrange the right side to get

(3.5) \displaystyle{ -\frac{\log\epsilon}{\log (2\pi)} -\frac{1}{2} = a - \frac{2}{\log a\log(2\pi)}}

Now we can easily see by way of the inequality

(3.6) \displaystyle{ a - \frac{2}{\log a\log(2\pi)} \le a}\quad\text{with}\quad a\in\mathbb{N}\, \wedge\, a>2

that the term \sqrt{a} had been dropped conservatively but rightfully \blacksquare

If you waited the whole time for the calc-script: you will find it in the repository if the authors of calc accept it; here, in another posting, otherwise.

[1] Spouge, John L. (1994). Computation of the Gamma, Digamma, and Trigamma Functions. SIAM Journal on Numerical Analysis 31 (3): 931–000.

From Day to Day

The surprising results of

mkdir `date '+%Y-%m-%d' && cd `date '+%Y-%m-%d'

at around midnight should have been expected.
How embarrassing!
It is—of course!—a slightly better idea to do something in the line
of:

TMPDATE=`date '+%Y-%m-%d'`;mkdir $TMPDATE && cd $TMPDATE

Or get silly and do it with a full blown script:

##########################################################################
# Create an empty directory named after the current date in the          #
# format: the four digit year, a minus-sign, the numerical month started #
# at one (1), a minus sign, and the numerical day of the month started   #
# with one (1)                                                           #
# Returns 1 for failure and 0 for success                                #
# Place in ~/.bashrc (or some other startup file) to make it work        #
# If you have chosen ~/.bashrc re-read it with                          #
#  . ~/.bashrc                                                           #
# to get a sense of achievement                                          #
# Please note, that there is a period in front of ~/.bashrc!             #
##########################################################################

mkdatedir(){
  # Full path to the date(1) program
  DATEP="/bin/date"
  # Full path to the mkdir(1) program
  MKDIRP="/bin/mkdir"
  # Full path to the ls(1) program
  LSP="/bin/ls"
  # Arguments to ls(1)
  LSA=" -d"
  # Format of the date
  DATEF=' +%Y-%m-%d'
  
  #########################
  # Actual program logic: #
  #########################
  
  # Get the date as a string and keep it in the temporary variable named
  # TMPDATE
  
  TMPDATE=$($DATEP $DATEF)
  
  # Uncomment for debugging purposes
  # echo "TMPDATE= $TMPDATE"

  # Try to make the directory
  $($MKDIRP $TMPDATE > /dev/null 2>&1 )
  # save the return value of mkdir(1) in the variable MKDIRR
  MKDIRR=$?
  # Check for errors
  if [ $MKDIRR  -eq 1 ]  
  then
    echo "Something went wrong while trying to create the directory $TMPDATE"
    TMPLS=$($LSP $LSA)
    if [ $TMPLS  ]
    then
      echo "The directory already exists"
    fi
    return 1
  fi
  
  # Jump right into the directory just successfully build
  cd $TMPDATE
  
  # TODO: write a check to check if "cd" failed
  
  # Uncomment for debugging purposes
  # echo $($LSP $LSA ../$TMPDATE)
  
  # Brag about it
  # echo "Success!"
  return 0
}

Yes, you’re right, if I say “silly” I mean silly!
So, if you ever find this code,or something similar in weird- and/or silliness: send it to Alex Papadimoulis’ site if you haven’t done so already.