Adventures of a Programmer: Parser Writing Peril VI

Yep, you guessed right, still no grammar 😉
But I implemented a square root for the big rationals from the last post. It is not very elegant, especially too ignorant about the exact value of the precision which it puts too high, but it works.

#ifndef MPI_LOG2
#   define MPI_LOG2 3.321928094887362347870319429489390175864831393024580612054756
#endif
#include <math.h>
static int mpq_epsilon;
void mpq_set_epsilon(int eps){
   if(eps <= 0){
     mpq_epsilon = 0;
   }
   else{
     /* eps assumed to be precision in number of decimal digits */
     mpq_epsilon = (int)((double)eps*MPI_LOG2) + (int)log2((double)eps);
   }
}
int mpq_get_epsilon(void){
   return mpq_epsilon;
}
/* sqrt() meant as a test for epsilon */
int mpq_sqrt(mp_rat *a, mp_rat *c){
  if(mpq_isneg(a)){
    return MP_VAL;
  }
  mpq_copy(a,c);
  printf("epsilon = %d\n",mpq_get_epsilon());
  /*
    This is too brutal.
    There are better ways to do it.
  */
  int eps  =  mpq_get_epsilon();
  mp_mul_2d(&c->numerator,  eps,&c->numerator);
  mp_mul_2d(&c->denominator,eps,&c->denominator);

  mp_sqrt(&c->numerator,&c->numerator);
  mp_sqrt(&c->denominator,&c->denominator);
  mpq_reduce(c);
  return MP_OKAY;
}

The angstzuschlag of log2(eps) (the logarithm base two of eps which holds the number of decimal digits of precision) is purely arbitrary.
The problem with this approach is that any perfect square will result in unnecessary inexact results when eps is not zero.
Changing the bit-shifting to multiplying with exponents of ten will get rid of that problem, lets…seeeeee…
Segmentation Fault!
Oops!
What the…
Ah, the function mpq_reduce() uses mpq_set_int() which itself uses mpq_reduce() which resulted in a nice crash. The first problem is easy to solve:

#define mp_isone(a) (((a)->used == 1 && (a)->dp[0] == 1) ? MP_YES : MP_NO)
#define MPQ_DIVISION_BY_ZERO -997
int mpq_reduce(mp_rat *a){
  mp_int gcd;
  int e;
  /* check if num is zero -> return 0 if it is*/
  if(mp_iszero(&a->numerator)){
    mp_set_int(&a->numerator,0);
    mp_set_int(&a->denominator,1);
    if( (e = mpq_set_sign(a,MP_ZPOS) ) != MP_OKAY){
      return e;
    }
    return MP_OKAY;
  }
  /* check if den is 1 -> return if it is*/
  if(mp_isone(&a->denominator)){
    mpq_normalize_sign(a);
    return MP_OKAY;
  }
  /* check if den is zero -> return undefined if it is*/
  if(mp_iszero(&a->denominator)){
    mp_set_int(&a->numerator,0);
    mp_set_int(&a->denominator,0);

    return MPQ_DIVISION_BY_ZERO;
  }
  /* check if den = num -> return 1 */
  if(mpq_isone(a)){
    mpq_normalize_sign(a);
    mp_set_int(&a->numerator,1);
    mp_set_int(&a->denominator,1);

    return MP_OKAY; 
  }
  mp_init(&gcd);
  /* compute gcd(num,den) */
  if( (e = mp_gcd(&a->numerator,&a->denominator,&gcd) ) != MP_OKAY){
    return e;
  }
  /* set num to num/gcd */
  if( (e = mp_div(&a->numerator,&gcd,&a->numerator,NULL) ) != MP_OKAY){
    return e;
  }
  /* set den to den/gcd */
  if( (e = mp_div(&a->denominator,&gcd,&a->denominator,NULL) ) != MP_OKAY){
    return e;
  }
  /* not to forget: */
  mpq_normalize_sign(a);
  return MP_OKAY;
}

The second problem is not so easy to solve because it shouldn’t have segfaulted before it hit the stack’s recursion limit. I leave that for later.
Where have I been…ah, yes, the problem with perfect squares. Using decimals as the multiplier doesn’t help at all, just give it a try.

  //mp_mul_2d(&c->numerator,  eps,&c->numerator);
  //mp_mul_2d(&c->denominator,eps,&c->denominator);
  mp_int tenpow;
  mp_init(&tenpow);
  mp_set_int(&tenpow,10);
  mp_expt_d(&tenpow,eps,&tenpow);
  mp_mul(&c->numerator, &tenpow ,&c->numerator);
  mp_mul(&c->denominator,&tenpow,&c->denominator);
  ...
  mp_clear(&tenpow);

It makes no difference if you expand the fraction with a power of ten or a power of two if you expand both sides of the fraction.
The actual precision gained with this method is also too low. The number of digits n of a square root of a number x with N is about N/2 depending on the number system, exactly half with base two. So we need twice the precision to get a square root of normal precision. We are using shifts, so simply double eps.
Note to self for further optimization: subtract number of digits of min(nbits(num),nbits(den)) from eps.

  /* check if eps<<1 is still smaller than the datatype used for eps!*/
  mp_mul_2d(&c->numerator,  eps<<1,&c->numerator);
  mp_mul_2d(&c->denominator,eps<<1,&c->denominator);

And a funny thing happens: it gives the exact result now.
Why?
Well, I do have a printf() giving me the actual value. I set epsilon to 50 which gives a result of 171 here. Which is odd, because doubling it is even curioser. Ok, ok, I’ll stop here 😉
The reason is that the number of digits of the square root of a number in base two is exactly half of the number of digits of the radicand. If the number of digits of the radicand is odd it cannot work it out exactly.
A hypothesis easily tested 😉
Two fractions given: 529/1849 ((23/43)^2), 3481/11664 ((59/108)^2), 1764/1849 ((42/43)^2), and 44521/46656 ((211/216)^2). Uh, no, that’s illegible. Let’s make a table.
With eps = 50

a/b (a/b)^2 odd eps even eps
23/43 529/1849 1258…143/235…312 23/43
59/108 3481/11664 322…847/590…109 59/108
42/43 1764/1849 229…653/235…312 42/43
211/216 44521/46656 115…927/118…219 211/216

Well, that’s what most scientist call a “failed chance to publish”. Because not many journals accept experiments with a negative outcome. That causes many already failed experiments to get duplicated because the experimentalist had no chance to know that it doesn’t work in the first playe. That heap of a male cow’s excrements wastes a lot of money and hours.
So, if you have such a failed experiment you can publish it at e.g.: arxive.org. But check if you will not get punished for it. Yes, that might happen!
But I digress.
So it is the even multiplier that counts?
Why, you ask?
Yes, that’s an interesting question.

int mpq_sqrt(mp_rat *a, mp_rat *c){
  int eps;
  if(mpq_isneg(a)){
    return MP_VAL;
  }

  /* place shortcut for zero and one  here */
  /* place shortcut for integers here */

  mpq_copy(a,c);
  eps  =  mpq_get_epsilon();
  if((eps & 0x1)){
    eps++;
  }
  eps <<= 1;

  /* reduce size of eps r.t. to min(digits(num),digits(den)) here */

  mp_mul_2d(&c->numerator,  eps,&c->numerator);
  mp_mul_2d(&c->denominator,eps,&c->denominator);
  mp_sqrt(&c->numerator,&c->numerator);
  mp_sqrt(&c->denominator,&c->denominator);
  mpq_reduce(c);
  return MP_OKAY;
}
Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s