# Uncomment the following line if you have Maple 9:
# macro( lattice=IntegerRelations:-LLL ):


macro(
	# main function:
	pyth_LLL=`algcurves/pyth_LLL`,

	# tools:
	isq=`algcurves/isq`,
	ifa=`algcurves/ifa`
):

# Find, if it exists, an integer point (x:y:z) on the
# conic a*x^2+b*y^2+c*z^2 = 0 where a,b,c are integers, not all 0.
#
# Mark van Hoeij, June 2001.
#
#
pyth_LLL:=proc(a::integer,b::integer,c::integer)
  local A,B,C,u,v,w,i,L,j,d,a1,a2,b1,b2,c1,c2;
if nargs=4 and args[4]=`sqrfree` then
	A,B,C := abs(a), abs(b), abs(c);

	u:=numtheory['msqrt'](-b/c,A);
	if u=FAIL then return FAIL fi;
	v:=numtheory['msqrt'](-c/a,B);
	if v=FAIL then return FAIL fi;
	w:=numtheory['msqrt'](-a/b,C);
	if w=FAIL then return FAIL fi;
	A,B,C := seq(round(evalf(10*sqrt(i),15)),i=[A,B,C]);

	# Consider the following system of equations:
	#
	# y = w*x mod C  }
	# z = u*y mod A  }   (**)
	# x = v*z mod B  }
	#
	# The solutions form a lattice L in Z^3 with det=A*B*C. The following
	# is a basis for L (the coordinates have been rescaled by multiplying
	# x by A, y by B and z by C).
	L :=[[A, B*w, C*(w*u - `if`(abs(b)=1,0,(w*u*v-1)/a/v mod b) * a )],
	     [0, B*c, C*(u*c -     (            u*c/a        mod b) * a )],
	     [0,   0, C*a*b]];

	L := lattice(L, integer);
	L := [seq([i[1]/A, i[2]/B, i[3]/C],i=L)];
	#
	# Now L is an LLL-reduced basis of the system (**) above.
	# Take some short vectors in the lattice L:
	#
	L := [ L[1], L[2],  L[2]+L[1], L[2]-L[1],
		     L[3],  L[3]+L[1], L[3]-L[1],
			    L[3]+L[2], L[3]-L[2],
		L[3]+L[2]+L[1], L[3]+L[2]-L[1],
		L[3]-L[2]+L[1], L[3]-L[2]-L[1]]:

	# Sort the short vectors:
	L:=sort([seq(`if`(member(0,i),NULL,i),i=L)], proc(u1,u2)
	  local l1,l2,i,j;
	  (l1,l2) := seq(length(i),i=[u1,u2]);
	  if l1=l2 then
	    (l1,l2) := seq(add(abs(j),j=i),i=[u1,u2])
	  fi;
	  if l1=l2 then
	    (l1,l2) := seq(length([i[1]/i[3],i[2]/i[3]]), i=[u1,u2])
	  fi;
	  evalb(l1 < l2)
	end);

 	for j in L do
	  if a*j[1]^2+b*j[2]^2+c*j[3]^2 = 0 then
		userinfo(2,'algcurves',`found point`,j,
		 `on conic`,a*_X^2+b*_Y^2+c*_Z^2);
		return j
	  fi
	od;
	userinfo(2,'algcurves',`pyth error args,L=`,[args],L);
	error FAIL
elif a=0 then [1,0,0]
elif b=0 then [0,1,0]
elif c=0 then [0,0,1]
elif (a>0 and b>0 and c>0) or (a<0 and b<0 and c<0) then FAIL
elif igcd(a,igcd(b,c))>1 then
  d:=igcd(a,igcd(b,c)); procname(a/d,b/d,c/d)
else
  if igcd(a,b)>1 then
    d:=igcd(a,b); v:=procname(a/d,b/d,c*d);
    if v<>FAIL then v:=[v[1],v[2],v[3]*d] fi
  elif igcd(a,c)>1 then
    d:=igcd(a,c); v:=procname(a/d,b*d,c/d);
    if v<>FAIL then v:=[v[1],v[2]*d,v[3]] fi
  elif igcd(b,c)>1 then
    d:=igcd(b,c); v:=procname(a*d,b/d,c/d);
    if v<>FAIL then v:=[v[1]*d,v[2],v[3]] fi
  else
    i := roots((j^2+a/b)*(j^2+b/c)*(j^2+c/a),j);
    if i<>[] then
	i:=i[1][1];
	j := denom(i), numer(i);
	if i^2=-a/b then return [j,0]
	elif i^2=-b/c then return [0,j]
	else return [numer(i), 0, denom(i)]
	fi
    fi;
    j := proc(u,v) local a;
	a:=abs(v);
	while type(a,even) do a:=a/2 od;
	numtheory['jacobi'](u,a)
    end:
    for i in [[-a*b,c], [-b*c,a], [-c*a,b]] do if j(op(i))<>1 then
	# Jacobi symbol shows non-existence of rational point.
	return FAIL
    fi od;
	a1,a2 := isq(a);
	b1,b2 := isq(b);
	c1,c2 := isq(c);
	v := procname(a1,b1,c1,`sqrfree`);
    if v<>FAIL then
	v := [v[1]*b2*c2, v[2]*c2*a2, v[3]*a2*b2];
	d:=igcd(igcd(v[1],v[2]),v[3]);
	v:=[v[1]/d,v[2]/d,v[3]/d]
    fi
  fi;
  v
fi
end:


# Square-free factorization.
isq := proc(a)
	local v,d,i,s,p1,p2,r,x,P;
	v := (proc(i) _Env_ifactor_easy := true; ifa(i) end)(a);
	s := v[1];
	d := indets(v,name);
	v := subs([seq(i=1,i=d)],v[2]);
	v := [seq(`if`(i[1]=1,NULL,i),i=v)];
	p1 := mul(``(i[1])^irem(i[2],2),i=v);
	p2 := mul(  (i[1])^iquo(i[2],2),i=v);
	if d<>{} then
		r := a/(s*expand(p1)*p2^2);
		d := roots(x^2-r,x);
		if d=[] then
			v := ifa(r)[2];
			if indets(v,name)<>{} then
				# Only possible if _Env_ifactor_easy = true
				error "too hard to ifactor, giving up"
			fi;
			p1 := p1 * mul(``(i[1])^irem(i[2],2),i=v);
			p2 := p2 * mul(  (i[1])^iquo(i[2],2),i=v)
		else
			p2 := p2*d[1][1]
		fi
	fi;
	P := expand(p1);
	ifactor(P) := p1;  # Because P will have to be factored again
			   # under the msolve command in procedure msqrt below.
	s*P, p2
end:

ifa := proc(n)
	local a;
	a := ifactors(n);
	if indets(a,'name')<>{} then
		# Clean up the mess caused by _Env_ifactor_easy := true
		forget(ifactor,n)
	fi;
	a
end: