4837.sagews

#An algorithm which transforms any Diophantine equation
#into an equivalent system of equations of the forms
#xi = 1, xi + xj = xk, xi * xj = xk
#Apoloniusz Tyszka, Krzysztof Molenda, Maciej Sporysz
#
#Example for D(x1, x2) = x1*x2 - 1
#
#author of code: Krzysztof Molenda ([email protected])
#code for SAGE 5.0.1 (http://www.sagemath.org/)


W.<x1, x2> = PolynomialRing(ZZ)
M = 1;
d_1 = 1; d_2 = 1;

mo = [ x1^(i1)*x2^(i2) for i1 in [0..d_1] for i2 in [0..d_2] ];
def pol( mon, a ):
       return a[0]*mo[0]+a[1]*mon[1]+a[2]*mon[2]+a[3]*mon[3];
t = [x1, x2];
t.extend([ pol(mo, [j1, j2, j3, j4]) for j1 in [-M..M] for j2 in [-M..M] for j3 in [-M..M] for j4 in [-M..M] 
           if pol(mo, [j1, j2, j3, j4]) != x1 and pol(mo, [j1, j2, j3, j4]) != x2 ]);
print 'card(T) = {n}'.format(n=len(t));
print "T=" + str( t[:3] + [...] + t[len(t)-2:] );
#print full t to the file
o=open('t.txt', 'w'); o.write( "t=\n" + str(t) + "\n\ncard(t) = " + str(len(t)) ); o.close();

#Remark:
#The index method does a linear search, and stops at the first matching item. If no matching item is found, it raises a ValueError exception.
#Lists in Phyton are numerated from 0, thus index()+1
print 'T[{idx}]={eq}'.format( idx = t.index(1)+1, eq='1' );
print 'T[{idx}]={eq}'.format( idx = t.index(x1*x2-1)+1, eq='x1*x2-1' );
print 'T[{idx}]={eq}'.format( idx = t.index(x1*x2)+1, eq='x1*x2' );

h1 = [ (i+1, j+1, k+1) for i in range(len(t)) for j in range(len(t)) for k in range(len(t)) if t[i]+t[j]==t[k] ];
print 'card(H1) = {n}'.format(n=len(h1));
print "H1=" + str( h1[:2] + [...] + h1[len(h1)-2:] );
#print full h1 to the file
o=open('h1.txt', 'w'); o.write( "h1=\n" + str(h1) + "\n\ncard(h1) = " + str(len(h1)) ); o.close();

h2 = [ (i+1, j+1, k+1) for i in range(len(t)) for j in range(len(t)) for k in range(len(t)) if t[i]*t[j]==t[k] ];
print 'card(H2) = {n}'.format(n=len(h2));
print "H2=" + str( h2[:2] + [...] + h2[len(h2)-2:] );
#print full h2 to the file
o=open('h2.txt', 'w'); o.write( "h2=\n" + str(h2) + "\n\ncard(h2) = " + str(len(h2)) ); o.close();

card(T) = 81 T=[x1, x2, -x1*x2 - x1 - x2 - 1, Ellipsis, x1 + x2 + 1, x1*x2 + x1 + x2 + 1] T[68]=1 T[17]=x1*x2-1 T[44]=x1*x2 card(H1) = 2401 H1=[(1, 2, 53), (1, 3, 6), Ellipsis, (81, 42, 80), (81, 43, 81)] card(H2) = 549 H2=[(1, 2, 44), (1, 7, 39), Ellipsis, (81, 43, 43), (81, 68, 81)]