[GiNaC-list] gcd: wrong sign of cofactor (heur_gcd is broken?)

[Sheplyakov Alexei]

Hello!

This simple program fails:

#include <iostream>
#include <stdexcept>
#include <cassert>
#include <ginac/ginac.h>
using namespace std;
using namespace GiNaC;

int main(int argc, char** argv)
{
	symbol x("x");
	symbol y("y");
		ex a = pow(x,2) - pow(y,2);
		ex b = pow(x,4) - pow(y, 4);
		ex ca, cb;
		cout << "GCD(" << a << ", " << b << ") = ";
		ex ab_gcd =  gcd(a, b, &ca, &cb, false);
		cout << ab_gcd << endl;
		assert((a-ca*ab_gcd).expand().is_zero());
		assert((b-cb*ab_gcd).expand().is_zero());
		return 0;
}

Documentation (normal.cpp:1256) says 
\begin{quote}
/*
 * Compute GCD (Greatest Common Divisor) of multivariate 
 * polynomials a(X) and b(X) in Z[X]. Optionally also compute 
 * the cofactors of a and b, defined by a = ca * gcd(a, b)
 * and b = cb * gcd(a, b).
 */
\end{quote}

so this behaviour is probably a bug.

I've got a patch to fix it (see attachment #1), but I'm not sure if
it is correct.

Relevant part of gdb session is attached too.

P.S.

I use GiNaC 1.3.1 from CVS, g++-3.4 from Debian testing.


--
Best regards,
  Alexei

-------------- next part --------------
Index: normal.cpp
===================================================================
RCS file: /home/cvs/GiNaC/ginac/normal.cpp,v
retrieving revision 1.94.2.3
diff -r1.94.2.3 normal.cpp
1250,1253d1249
< 				ex lc = g.lcoeff(x);
< 				if (is_exactly_a<numeric>(lc) && ex_to<numeric>(lc).is_negative())
< 					return -g;
< 				else
-------------- next part --------------
    // Apply evaluation homomorphism and calculate GCD     |#1  0xb7efcb23 in GiNaC::gcd (a=@0xbfffec40,               
    ex cp, cq;                                             |    b=@0xbfffec50, ca=0xbfffebd0, cb=0xbfffebf0,           
    ex gamma = heur_gcd(p.subs(x == xi, subs_options::no_pa|    check_args=false)                                      
ttern), q.subs(x == xi, subs_options::no_pattern), &cp, &cq|    at ../../ginac-1.3.1-orig/ginac/normal.cpp:1462        
, var+1).expand();                                         |#2  0x08049cd2 in main (argc=1, argv=0xbfffed24)           
    if (!is_exactly_a<fail>(gamma)) {                      |    at gcd_bugreport.cpp:16                                
                                                           |(gdb) call a.dbgprint()                                    
      // Reconstruct polynomial from GCD of mapped polynomi|-y^2+x^2                                                   
als                                                        |(gdb) call b.dbgprint()                                    
      ex g = interpolate(gamma, xi, x, maxdeg);            |-y^4+x^4                                                   
                                                           |(gdb) call gc.value.debug_print()                          
      // Remove integer content                            |(cl_I) 1                                                   
      g /= g.integer_content();                            |(gdb) call g.dbgprint()                                    
                                                           |-y^2+x^2                                                   
      // If the calculated polynomial divides both p and q,|(gdb) call lc.dbgprint()                                   
 this is the GCD                                           |-1                                                         
      ex dummy;                                            |(gdb) call ca->dbgprint()                                  
      if (divide_in_z(p, g, ca ? *ca : dummy, var) && divid|1                                                          
e_in_z(q, g, cb ? *cb : dummy, var)) {                     |(gdb) call cb->dbgprint()                                  
        g *= gc;                                           |y^2+x^2                                                    
        // I think one needs just return g; here [-A.S.]   |(gdb)                                                      
        ex lc = g.lcoeff(x);                               |                                                           
        if (is_exactly_a<numeric>(lc) && ex_to<numeric>(lc)|                                                      
.is_negative())                                            |                                                           
          return -g;                                       |[No File] [RO]------------------------- 0,0-1 -------   All
        else                                               |int main(int argc, char** argv)                            
          return g;                                        |{                                                          
      }                                                    |  symbol x("x");                                           
    }                                                      |  symbol y("y");                                           
                                                           |    ex a = pow(x,2) - pow(y,2);                            
    // Next evaluation point                               |    ex b = pow(x,4) - pow(y, 4);                           
    xi = iquo(xi * isqrt(isqrt(xi)) * numeric(73794), numer|    ex ca, cb;                                             
ic(27011));                                                |    cout << "GCD(" << a << ", " << b << ") = ";            
  }                                                        |    ex ab_gcd =  gcd(a, b, &ca, &cb, false);               
  return (new fail())->setflag(status_flags::dynallocated);|    cout << ab_gcd << endl;                                
}                                                          |    cout << "a cofactor: " << ca << endl;                  
                                                           |    cout << "b cofactor: " << cb << endl; 
                                                           |