Re: real polynomials

[John Cremona <john.cremona@gmail.com>]

On 6 November 2015 at 16:15, Karim Belabas
<Karim.Belabas@math.u-bordeaux.fr> wrote:
> * John Cremona [2015-11-06 16:38]:
> [...]
>> I will give you a laugh by showing you my own version:
>>
>> GEN dbltorat(double x) // assumes x=a/2^e
>> {
>>   long e;
>>   GEN res = cgetg(3, t_FRAC);
>>   if (x==0)
>>     {
>>       gel(res,1) = cgeti(0);
>>       gel(res,2) = cgeti(1);
>>     }
>>   else
>>     {
>>       gel(res,1) = mantissa_real(dbltor(x), &e);
>>       gel(res,2) = cgeti(2<<e);
>>     }
>>   return res;
>> }
>>
>> (tested!  a little)
>
> There are a few problems indeed:

This is very instructive (possibly for others too), so thanks for
taking the time.


>
> 1) cgeti(0) / cgeti(1) / cgeti(2<<e) don't do what you think

No, of course I meant stoi().

>
> In fact, both 0 and 1 are illegal arguments corrupting the stack
> (immediately for 0, after other low-level commands such as setsigne() for 1) !
>
>
> stoi(0), stoi(1) will both work [ think signed_long-to-integer ].
> And so would utoi(0), utoi(1)   [ think unsigned_long-to-integer ]
>
> cgeti(n) is a low-level routine reserving n words for later assignment
> of an integer, it does not return a valid GEN in itself (and we must
> have n >= 3 ) !
>

OK, I did read section 4.4 but then misread it!


> 2) you probably meant  1 << e instead of 2 << e;

Indeed I did.  Or even 2**e (but not 1**e).

>
> 3) int2n(e) will be more robust than stoi(1 << e), which breaks for
>   e >= BITS_IN_LONG-1

I looked for such a function but did not find it.

>
> 4) t_FRAC members must be in normal form (lowest terms). Just dividing
> the numerator by the denominator will do that automatically:
>
>   num = mantissa_real(dbltor(x), &e);
>   den = int2n(e);
>   return gdiv(num, den);
>

For some reason I thought that calling gdiv on two t_INTs would give a
truncated answer.  In fact I wanted to make a t_INT out of the
mantissa and then right-shift it by e places, but the function which
does that does (I think?) truncate.

>
>
> Use high level constructors (such as gdiv / stoi) as much as possible, instead
> of low lever allocations (such as cgetg / cgeti); the latter is marginally more
> efficient but leads to (many) pitfalls...

That had been my intention!

>
> Don't bother with garbage collecting for first programs; if absolutely
> necessary (stack overflows), stick to the

I will not.  This is in fact a small program so gc can almost
certainly be ignored anyway.

>
>   pari_sp av = avma;
>   ...
>   return gerepilecopy(av, result);
>
> idiom. Good luck !
>

I will report on success here...

John

> Cheers,
>
>     K.B.
>
> P.S. There's a PARI/GP Atelier in Grenoble from January 11th-15th 2016.
> It's the ideal event to quickly learn the basics :-)
>
> --
> Karim Belabas, IMB (UMR 5251)  Tel: (+33) (0)5 40 00 26 17
> Universite de Bordeaux         Fax: (+33) (0)5 40 00 69 50
> 351, cours de la Liberation    http://www.math.u-bordeaux.fr/~kbelabas/
> F-33405 Talence (France)       http://pari.math.u-bordeaux.fr/  [PARI/GP]
> `
>