John Cremona on Fri, 06 Nov 2015 17:51:27 +0100 |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
Re: real polynomials |
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] > ` >