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