Trac 15767: explain application of Sturm's theorem

src/ext/pari/simon/ell.gp

@@ -879,6 +879,37 @@ if( DEBUGLEVEL_ell >= 4, print("    end of nflocallysoluble"));
 
 \\ places reelles 
 
+/*
+  Peter Bruin, August 2014: code adapted to fix precision problems.
+
+  We want to count the number of real roots of a polynomial of the
+  form x^4 + s*x^2 + a*x + b.  The action of complex conjugation on
+  the complex roots determines the sign of the discriminant Delta:
+
+    cycle type () or (1 2)(3 4), i.e. 4 or 0 real roots: Delta > 0
+    cycle type (1 2), i.e. exactly 2 real roots: Delta < 0
+
+  We may assume Delta > 0 and have to decide whether there are 0 or 4
+  real roots.  We apply Sturm's theorem to the following sequence
+  (with t = 2*s^3 - 8*b*s + 9*a^2):
+
+    p0 = x^4 + s*x^2 + a*x + b
+    p1 = 4*x^3 + 2*s*x + a
+    p2 = x*p1 - 4*p0
+       = -2*s*x^2 - 3*a*x - 4*b
+    p3 = (3*a - 2*s*x)*p2 - s^2*p1
+       = -t*x - a*(s^2 + 12*b)
+    p4 = Delta
+
+  Note that the case of 4 real roots can only occur if both s and t
+  are non-zero, otherwise the Sturm sequence is shorter than this one.
+  By Sturm's theorem, the number of roots equals
+
+    (sign changes in [1, -4, -2*s, t, Delta]) -
+    (sign changes in [1, 4, -2*s, -t, Delta]).
+
+  Hence the only way to have 4 real roots is if s < 0 and t < 0.
+*/
   if( nf.r1,
     c = pollead(pol);
     pol /= c;
@@ -889,7 +920,6 @@ if( DEBUGLEVEL_ell >= 4, print("    end of nflocallysoluble"));
     for( i = 1, nf.r1,
       if( nfrealsign(nf,c,i) > 0, next);
       if( nfrealsign(nf,Delta,i) < 0, next);
-      \\ Sturm => 4 roots if s < 0 and t < 0, otherwise 0 roots
       if( nfrealsign(nf,s,i) < 0 && nfrealsign(nf,t,i) < 0, next);
 if( DEBUGLEVEL_ell >= 2, print("  not ELS at infinity"));
 if( DEBUGLEVEL_ell >= 4, print("    end of nflocallysoluble"));