[PURRS-devel] Bug found: it is in GiNaC's numer_denom

[Roberto Bagnara]

Hi there,

after following a dozen of wrong tracks (which required, among other
things, producing particular versions of GMP, CLN and GiNaC and a
patched version of valgrind to make sure we were not bitten by a
memory access bug) I believe I have found the cause of the long
standing bug we have in the system.

The problem is in GiNaC's numer_denom and a workaround is constituted
by the improvement to

     void Expr::numerator_denominator(Expr& x, Expr& y) const

I suggested last Friday.  It turns out that the problem in GiNaC's
numer_denom can be reproduced easily with ginsh.  Here is a sample
session where numer_denom is computed several times with the same
expression (using the up-arrow history recall mechanism):

$ ./ginsh
ginsh - GiNaC Interactive Shell (GiNaC V1.1.1)
   __,  _______  Copyright (C) 1999-2003 Johannes Gutenberg University Mainz,
  (__) *       | Germany.  This is free software with ABSOLUTELY NO WARRANTY.
   ._) i N a C | You are welcome to redistribute it under certain conditions.
<-------------' For details type `warranty;'.

Type ?? for a list of help topics.
 > 
numer_denom((-(4/79*I)*(-1/8+(1/8*I)*sqrt(79))*sqrt(79)-(4/79*I)*(1/8+(1/8*I)*sqrt(79))*sqrt(79))*d-(3/79*I)*(1/8-(1/8*I)*sqrt(79))^2*sqrt(79)+(3/79*I)*(1/8+(1/8*I)*sqrt(79))^2*sqrt(79)-(4/79*I)*(-3/32-(3/32*I)*sqrt(79))*sqrt(79)-(4/79*I)*(3/32-(3/32*I)*sqrt(79))*sqrt(79)+(-(4/79*I)*(-1/32-(1/32*I)*sqrt(79))*sqrt(79)-(4/79*I)*(1/32-(1/32*I)*sqrt(79))*sqrt(79)+(4/79*I)*(1/8-(1/8*I)*sqrt(79))^2*sqrt(79)-(4/79*I)*(1/8+(1/8*I)*sqrt(79))^2*sqrt(79))*c-3/4*(-(1/8-(1/8*I)*sqrt(79))^2-(1/4*I)*sqrt(79)+(1/8-(1/8*I)*sqrt(79))*(1/8+(1/8*I)*sqrt(79)))^(-1)*(1/8-(1/8*I)*sqrt(79))^3+3/4*(-(1/8-(1/8*I)*sqrt(79))^2-(1/4*I)*sqrt(79)+(1/8-(1/8*I)*sqrt(79))*(1/8+(1/8*I)*sqrt(79)))^(-1)*(1/8-(1/8*I)*sqrt(79))-3/4*((1/8+(1/8*I)*sqrt(79))^2-(1/4*I)*sqrt(79)+(-1/8+(1/8*I)*sqrt(79))*(1/8+(1/8*I)*sqrt(79)))^(-1)*(1/8+(1/8*I)*sqrt(79))+3/4*((1/8+(1/8*I)*sqrt(79))^2-(1/4*I)*sqrt(79)+(-1/8+(1/8*I)*sqrt(79))*(1/8+(1/8*I)*sqrt(79)))^(-1)*(1/8+(1/8*I)*sqrt(79))^3);
{-1264*d,-1264}
 > 
numer_denom((-(4/79*I)*(-1/8+(1/8*I)*sqrt(79))*sqrt(79)-(4/79*I)*(1/8+(1/8*I)*sqrt(79))*sqrt(79))*d-(3/79*I)*(1/8-(1/8*I)*sqrt(79))^2*sqrt(79)+(3/79*I)*(1/8+(1/8*I)*sqrt(79))^2*sqrt(79)-(4/79*I)*(-3/32-(3/32*I)*sqrt(79))*sqrt(79)-(4/79*I)*(3/32-(3/32*I)*sqrt(79))*sqrt(79)+(-(4/79*I)*(-1/32-(1/32*I)*sqrt(79))*sqrt(79)-(4/79*I)*(1/32-(1/32*I)*sqrt(79))*sqrt(79)+(4/79*I)*(1/8-(1/8*I)*sqrt(79))^2*sqrt(79)-(4/79*I)*(1/8+(1/8*I)*sqrt(79))^2*sqrt(79))*c-3/4*(-(1/8-(1/8*I)*sqrt(79))^2-(1/4*I)*sqrt(79)+(1/8-(1/8*I)*sqrt(79))*(1/8+(1/8*I)*sqrt(79)))^(-1)*(1/8-(1/8*I)*sqrt(79))^3+3/4*(-(1/8-(1/8*I)*sqrt(79))^2-(1/4*I)*sqrt(79)+(1/8-(1/8*I)*sqrt(79))*(1/8+(1/8*I)*sqrt(79)))^(-1)*(1/8-(1/8*I)*sqrt(79))-3/4*((1/8+(1/8*I)*sqrt(79))^2-(1/4*I)*sqrt(79)+(-1/8+(1/8*I)*sqrt(79))*(1/8+(1/8*I)*sqrt(79)))^(-1)*(1/8+(1/8*I)*sqrt(79))+3/4*((1/8+(1/8*I)*sqrt(79))^2-(1/4*I)*sqrt(79)+(-1/8+(1/8*I)*sqrt(79))*(1/8+(1/8*I)*sqrt(79)))^(-1)*(1/8+(1/8*I)*sqrt(79))^3);
{-1264*d,-1264}
 > 
numer_denom((-(4/79*I)*(-1/8+(1/8*I)*sqrt(79))*sqrt(79)-(4/79*I)*(1/8+(1/8*I)*sqrt(79))*sqrt(79))*d-(3/79*I)*(1/8-(1/8*I)*sqrt(79))^2*sqrt(79)+(3/79*I)*(1/8+(1/8*I)*sqrt(79))^2*sqrt(79)-(4/79*I)*(-3/32-(3/32*I)*sqrt(79))*sqrt(79)-(4/79*I)*(3/32-(3/32*I)*sqrt(79))*sqrt(79)+(-(4/79*I)*(-1/32-(1/32*I)*sqrt(79))*sqrt(79)-(4/79*I)*(1/32-(1/32*I)*sqrt(79))*sqrt(79)+(4/79*I)*(1/8-(1/8*I)*sqrt(79))^2*sqrt(79)-(4/79*I)*(1/8+(1/8*I)*sqrt(79))^2*sqrt(79))*c-3/4*(-(1/8-(1/8*I)*sqrt(79))^2-(1/4*I)*sqrt(79)+(1/8-(1/8*I)*sqrt(79))*(1/8+(1/8*I)*sqrt(79)))^(-1)*(1/8-(1/8*I)*sqrt(79))^3+3/4*(-(1/8-(1/8*I)*sqrt(79))^2-(1/4*I)*sqrt(79)+(1/8-(1/8*I)*sqrt(79))*(1/8+(1/8*I)*sqrt(79)))^(-1)*(1/8-(1/8*I)*sqrt(79))-3/4*((1/8+(1/8*I)*sqrt(79))^2-(1/4*I)*sqrt(79)+(-1/8+(1/8*I)*sqrt(79))*(1/8+(1/8*I)*sqrt(79)))^(-1)*(1/8+(1/8*I)*sqrt(79))+3/4*((1/8+(1/8*I)*sqrt(79))^2-(1/4*I)*sqrt(79)+(-1/8+(1/8*I)*sqrt(79))*(1/8+(1/8*I)*sqrt(79)))^(-1)*(1/8+(1/8*I)*sqrt(79))^3);
{1264*d,1264}
 > 
numer_denom((-(4/79*I)*(-1/8+(1/8*I)*sqrt(79))*sqrt(79)-(4/79*I)*(1/8+(1/8*I)*sqrt(79))*sqrt(79))*d-(3/79*I)*(1/8-(1/8*I)*sqrt(79))^2*sqrt(79)+(3/79*I)*(1/8+(1/8*I)*sqrt(79))^2*sqrt(79)-(4/79*I)*(-3/32-(3/32*I)*sqrt(79))*sqrt(79)-(4/79*I)*(3/32-(3/32*I)*sqrt(79))*sqrt(79)+(-(4/79*I)*(-1/32-(1/32*I)*sqrt(79))*sqrt(79)-(4/79*I)*(1/32-(1/32*I)*sqrt(79))*sqrt(79)+(4/79*I)*(1/8-(1/8*I)*sqrt(79))^2*sqrt(79)-(4/79*I)*(1/8+(1/8*I)*sqrt(79))^2*sqrt(79))*c-3/4*(-(1/8-(1/8*I)*sqrt(79))^2-(1/4*I)*sqrt(79)+(1/8-(1/8*I)*sqrt(79))*(1/8+(1/8*I)*sqrt(79)))^(-1)*(1/8-(1/8*I)*sqrt(79))^3+3/4*(-(1/8-(1/8*I)*sqrt(79))^2-(1/4*I)*sqrt(79)+(1/8-(1/8*I)*sqrt(79))*(1/8+(1/8*I)*sqrt(79)))^(-1)*(1/8-(1/8*I)*sqrt(79))-3/4*((1/8+(1/8*I)*sqrt(79))^2-(1/4*I)*sqrt(79)+(-1/8+(1/8*I)*sqrt(79))*(1/8+(1/8*I)*sqrt(79)))^(-1)*(1/8+(1/8*I)*sqrt(79))+3/4*((1/8+(1/8*I)*sqrt(79))^2-(1/4*I)*sqrt(79)+(-1/8+(1/8*I)*sqrt(79))*(1/8+(1/8*I)*sqrt(79)))^(-1)*(1/8+(1/8*I)*sqrt(79))^3);
{-1264*d,-1264}
 > 
numer_denom((-(4/79*I)*(-1/8+(1/8*I)*sqrt(79))*sqrt(79)-(4/79*I)*(1/8+(1/8*I)*sqrt(79))*sqrt(79))*d-(3/79*I)*(1/8-(1/8*I)*sqrt(79))^2*sqrt(79)+(3/79*I)*(1/8+(1/8*I)*sqrt(79))^2*sqrt(79)-(4/79*I)*(-3/32-(3/32*I)*sqrt(79))*sqrt(79)-(4/79*I)*(3/32-(3/32*I)*sqrt(79))*sqrt(79)+(-(4/79*I)*(-1/32-(1/32*I)*sqrt(79))*sqrt(79)-(4/79*I)*(1/32-(1/32*I)*sqrt(79))*sqrt(79)+(4/79*I)*(1/8-(1/8*I)*sqrt(79))^2*sqrt(79)-(4/79*I)*(1/8+(1/8*I)*sqrt(79))^2*sqrt(79))*c-3/4*(-(1/8-(1/8*I)*sqrt(79))^2-(1/4*I)*sqrt(79)+(1/8-(1/8*I)*sqrt(79))*(1/8+(1/8*I)*sqrt(79)))^(-1)*(1/8-(1/8*I)*sqrt(79))^3+3/4*(-(1/8-(1/8*I)*sqrt(79))^2-(1/4*I)*sqrt(79)+(1/8-(1/8*I)*sqrt(79))*(1/8+(1/8*I)*sqrt(79)))^(-1)*(1/8-(1/8*I)*sqrt(79))-3/4*((1/8+(1/8*I)*sqrt(79))^2-(1/4*I)*sqrt(79)+(-1/8+(1/8*I)*sqrt(79))*(1/8+(1/8*I)*sqrt(79)))^(-1)*(1/8+(1/8*I)*sqrt(79))+3/4*((1/8+(1/8*I)*sqrt(79))^2-(1/4*I)*sqrt(79)+(-1/8+(1/8*I)*sqrt(79))*(1/8+(1/8*I)*sqrt(79)))^(-1)*(1/8+(1/8*I)*sqrt(79))^3);
{-1264*d,-1264}
 > 
numer_denom((-(4/79*I)*(-1/8+(1/8*I)*sqrt(79))*sqrt(79)-(4/79*I)*(1/8+(1/8*I)*sqrt(79))*sqrt(79))*d-(3/79*I)*(1/8-(1/8*I)*sqrt(79))^2*sqrt(79)+(3/79*I)*(1/8+(1/8*I)*sqrt(79))^2*sqrt(79)-(4/79*I)*(-3/32-(3/32*I)*sqrt(79))*sqrt(79)-(4/79*I)*(3/32-(3/32*I)*sqrt(79))*sqrt(79)+(-(4/79*I)*(-1/32-(1/32*I)*sqrt(79))*sqrt(79)-(4/79*I)*(1/32-(1/32*I)*sqrt(79))*sqrt(79)+(4/79*I)*(1/8-(1/8*I)*sqrt(79))^2*sqrt(79)-(4/79*I)*(1/8+(1/8*I)*sqrt(79))^2*sqrt(79))*c-3/4*(-(1/8-(1/8*I)*sqrt(79))^2-(1/4*I)*sqrt(79)+(1/8-(1/8*I)*sqrt(79))*(1/8+(1/8*I)*sqrt(79)))^(-1)*(1/8-(1/8*I)*sqrt(79))^3+3/4*(-(1/8-(1/8*I)*sqrt(79))^2-(1/4*I)*sqrt(79)+(1/8-(1/8*I)*sqrt(79))*(1/8+(1/8*I)*sqrt(79)))^(-1)*(1/8-(1/8*I)*sqrt(79))-3/4*((1/8+(1/8*I)*sqrt(79))^2-(1/4*I)*sqrt(79)+(-1/8+(1/8*I)*sqrt(79))*(1/8+(1/8*I)*sqrt(79)))^(-1)*(1/8+(1/8*I)*sqrt(79))+3/4*((1/8+(1/8*I)*sqrt(79))^2-(1/4*I)*sqrt(79)+(-1/8+(1/8*I)*sqrt(79))*(1/8+(1/8*I)*sqrt(79)))^(-1)*(1/8+(1/8*I)*sqrt(79))^3);
{1264*d,1264}
 > 
numer_denom((-(4/79*I)*(-1/8+(1/8*I)*sqrt(79))*sqrt(79)-(4/79*I)*(1/8+(1/8*I)*sqrt(79))*sqrt(79))*d-(3/79*I)*(1/8-(1/8*I)*sqrt(79))^2*sqrt(79)+(3/79*I)*(1/8+(1/8*I)*sqrt(79))^2*sqrt(79)-(4/79*I)*(-3/32-(3/32*I)*sqrt(79))*sqrt(79)-(4/79*I)*(3/32-(3/32*I)*sqrt(79))*sqrt(79)+(-(4/79*I)*(-1/32-(1/32*I)*sqrt(79))*sqrt(79)-(4/79*I)*(1/32-(1/32*I)*sqrt(79))*sqrt(79)+(4/79*I)*(1/8-(1/8*I)*sqrt(79))^2*sqrt(79)-(4/79*I)*(1/8+(1/8*I)*sqrt(79))^2*sqrt(79))*c-3/4*(-(1/8-(1/8*I)*sqrt(79))^2-(1/4*I)*sqrt(79)+(1/8-(1/8*I)*sqrt(79))*(1/8+(1/8*I)*sqrt(79)))^(-1)*(1/8-(1/8*I)*sqrt(79))^3+3/4*(-(1/8-(1/8*I)*sqrt(79))^2-(1/4*I)*sqrt(79)+(1/8-(1/8*I)*sqrt(79))*(1/8+(1/8*I)*sqrt(79)))^(-1)*(1/8-(1/8*I)*sqrt(79))-3/4*((1/8+(1/8*I)*sqrt(79))^2-(1/4*I)*sqrt(79)+(-1/8+(1/8*I)*sqrt(79))*(1/8+(1/8*I)*sqrt(79)))^(-1)*(1/8+(1/8*I)*sqrt(79))+3/4*((1/8+(1/8*I)*sqrt(79))^2-(1/4*I)*sqrt(79)+(-1/8+(1/8*I)*sqrt(79))*(1/8+(1/8*I)*sqrt(79)))^(-1)*(1/8+(1/8*I)*sqrt(79))^3);
{-1264*d,-1264}
 > 
numer_denom((-(4/79*I)*(-1/8+(1/8*I)*sqrt(79))*sqrt(79)-(4/79*I)*(1/8+(1/8*I)*sqrt(79))*sqrt(79))*d-(3/79*I)*(1/8-(1/8*I)*sqrt(79))^2*sqrt(79)+(3/79*I)*(1/8+(1/8*I)*sqrt(79))^2*sqrt(79)-(4/79*I)*(-3/32-(3/32*I)*sqrt(79))*sqrt(79)-(4/79*I)*(3/32-(3/32*I)*sqrt(79))*sqrt(79)+(-(4/79*I)*(-1/32-(1/32*I)*sqrt(79))*sqrt(79)-(4/79*I)*(1/32-(1/32*I)*sqrt(79))*sqrt(79)+(4/79*I)*(1/8-(1/8*I)*sqrt(79))^2*sqrt(79)-(4/79*I)*(1/8+(1/8*I)*sqrt(79))^2*sqrt(79))*c-3/4*(-(1/8-(1/8*I)*sqrt(79))^2-(1/4*I)*sqrt(79)+(1/8-(1/8*I)*sqrt(79))*(1/8+(1/8*I)*sqrt(79)))^(-1)*(1/8-(1/8*I)*sqrt(79))^3+3/4*(-(1/8-(1/8*I)*sqrt(79))^2-(1/4*I)*sqrt(79)+(1/8-(1/8*I)*sqrt(79))*(1/8+(1/8*I)*sqrt(79)))^(-1)*(1/8-(1/8*I)*sqrt(79))-3/4*((1/8+(1/8*I)*sqrt(79))^2-(1/4*I)*sqrt(79)+(-1/8+(1/8*I)*sqrt(79))*(1/8+(1/8*I)*sqrt(79)))^(-1)*(1/8+(1/8*I)*sqrt(79))+3/4*((1/8+(1/8*I)*sqrt(79))^2-(1/4*I)*sqrt(79)+(-1/8+(1/8*I)*sqrt(79))*(1/8+(1/8*I)*sqrt(79)))^(-1)*(1/8+(1/8*I)*sqrt(79))^3);
{1264*d,1264}
 > 
numer_denom((-(4/79*I)*(-1/8+(1/8*I)*sqrt(79))*sqrt(79)-(4/79*I)*(1/8+(1/8*I)*sqrt(79))*sqrt(79))*d-(3/79*I)*(1/8-(1/8*I)*sqrt(79))^2*sqrt(79)+(3/79*I)*(1/8+(1/8*I)*sqrt(79))^2*sqrt(79)-(4/79*I)*(-3/32-(3/32*I)*sqrt(79))*sqrt(79)-(4/79*I)*(3/32-(3/32*I)*sqrt(79))*sqrt(79)+(-(4/79*I)*(-1/32-(1/32*I)*sqrt(79))*sqrt(79)-(4/79*I)*(1/32-(1/32*I)*sqrt(79))*sqrt(79)+(4/79*I)*(1/8-(1/8*I)*sqrt(79))^2*sqrt(79)-(4/79*I)*(1/8+(1/8*I)*sqrt(79))^2*sqrt(79))*c-3/4*(-(1/8-(1/8*I)*sqrt(79))^2-(1/4*I)*sqrt(79)+(1/8-(1/8*I)*sqrt(79))*(1/8+(1/8*I)*sqrt(79)))^(-1)*(1/8-(1/8*I)*sqrt(79))^3+3/4*(-(1/8-(1/8*I)*sqrt(79))^2-(1/4*I)*sqrt(79)+(1/8-(1/8*I)*sqrt(79))*(1/8+(1/8*I)*sqrt(79)))^(-1)*(1/8-(1/8*I)*sqrt(79))-3/4*((1/8+(1/8*I)*sqrt(79))^2-(1/4*I)*sqrt(79)+(-1/8+(1/8*I)*sqrt(79))*(1/8+(1/8*I)*sqrt(79)))^(-1)*(1/8+(1/8*I)*sqrt(79))+3/4*((1/8+(1/8*I)*sqrt(79))^2-(1/4*I)*sqrt(79)+(-1/8+(1/8*I)*sqrt(79))*(1/8+(1/8*I)*sqrt(79)))^(-1)*(1/8+(1/8*I)*sqrt(79))^3);
{1264*d,1264}
 > 
numer_denom((-(4/79*I)*(-1/8+(1/8*I)*sqrt(79))*sqrt(79)-(4/79*I)*(1/8+(1/8*I)*sqrt(79))*sqrt(79))*d-(3/79*I)*(1/8-(1/8*I)*sqrt(79))^2*sqrt(79)+(3/79*I)*(1/8+(1/8*I)*sqrt(79))^2*sqrt(79)-(4/79*I)*(-3/32-(3/32*I)*sqrt(79))*sqrt(79)-(4/79*I)*(3/32-(3/32*I)*sqrt(79))*sqrt(79)+(-(4/79*I)*(-1/32-(1/32*I)*sqrt(79))*sqrt(79)-(4/79*I)*(1/32-(1/32*I)*sqrt(79))*sqrt(79)+(4/79*I)*(1/8-(1/8*I)*sqrt(79))^2*sqrt(79)-(4/79*I)*(1/8+(1/8*I)*sqrt(79))^2*sqrt(79))*c-3/4*(-(1/8-(1/8*I)*sqrt(79))^2-(1/4*I)*sqrt(79)+(1/8-(1/8*I)*sqrt(79))*(1/8+(1/8*I)*sqrt(79)))^(-1)*(1/8-(1/8*I)*sqrt(79))^3+3/4*(-(1/8-(1/8*I)*sqrt(79))^2-(1/4*I)*sqrt(79)+(1/8-(1/8*I)*sqrt(79))*(1/8+(1/8*I)*sqrt(79)))^(-1)*(1/8-(1/8*I)*sqrt(79))-3/4*((1/8+(1/8*I)*sqrt(79))^2-(1/4*I)*sqrt(79)+(-1/8+(1/8*I)*sqrt(79))*(1/8+(1/8*I)*sqrt(79)))^(-1)*(1/8+(1/8*I)*sqrt(79))+3/4*((1/8+(1/8*I)*sqrt(79))^2-(1/4*I)*sqrt(79)+(-1/8+(1/8*I)*sqrt(79))*(1/8+(1/8*I)*sqrt(79)))^(-1)*(1/8+(1/8*I)*sqrt(79))^3);
{-1264*d,-1264}
 >

Even disregarding the fact that the numerator and denominator returned
are trivially not coprime, numer_denom does not behave like a mathematical
function.  With the same input, it produces the sequence of answers

{-1264*d,-1264}
{-1264*d,-1264}
{1264*d,1264}
{-1264*d,-1264}
{-1264*d,-1264}
{1264*d,1264}
{-1264*d,-1264}
{1264*d,1264}
{1264*d,1264}
{-1264*d,-1264}

So, in the old implementation of Expr::numerator_denominator()
(which was silly performance-wise, but perfectly legal nonetheless),
in

   x =  Base::numer_denom().op(0);
   y =  Base::numer_denom().op(1);

we could end up with a pair (x, y) not corresponding to the original
expression.  I guess this may well be considered a bug in GiNaC
(even though, given our previous experiences, I wouldn't bet upon it).
Ciao

      Roberto

-- 
Prof. Roberto Bagnara
Computer Science Group
Department of Mathematics, University of Parma, Italy
http://www.cs.unipr.it/~bagnara/
mailto:bagnara at cs.unipr.it