improve qfsigndicho + more tests

libre/Qfmodife.gp

@@ -94,7 +94,6 @@ DEBUGLEVEL_myqfsolve=0;
         \\choix du vecteur
         Y=vector(n,i,random(bor)-bor>>1)~;
         H=Y~*Qt*Y;
-        print1(count," ");
         if(DEBUGLEVEL_myqfsolve>0,
             print("Essai de completion:",count);
             print("Y=",Y);
@@ -132,6 +131,7 @@ DEBUGLEVEL_myqfsolve=0;
             );\\fin if "ispseudoprime"
         );\\fin du if prim
     );\\fin Boucle principale (while)
+    print("count = ", count);
 
     \\On a un "bon" vecteur Y, donc on complète la matrice
     for(i=1,n,

programme/Pinceau.gp

@@ -1,26 +1,28 @@
 \\ Ce programme renvoie deux entiers u et v tels que uA+vB soit équilibré.
 dichoqfsign(A,B)=
-{ my(n = #A, pol, sol, s, m, s1, a, b, q, t=1, c, d, u, v);
-  \\ test de m=1
-  if( qfsign(A+B)[1] >= floor(n/2), return(1));
-  \\ test de -1
-  if( qfsign(-A+B)[1] >= floor(n/2), return(-1));
-  pol = matdet(A*x+B);
+{ my(n = #A, pol, sol, s, m, s1, t, a, b, c, d, q);
+
+  pol = matdet(A*'x + B);
   s = polrootsreal(pol);
   m = #s;
   s1 = max( abs(floor(s[1])) , abs(ceil(s[m]))) + 1;
-  [a,b] = [-s1,s1];
+  [a,b] = [-s1, s1 + 1];
+  c = qfsign(a*A+B)[1];
+  d = n - c;
   \\dichotomie
-  while( t==0 || q != ceil(n/2),   
-    c = qfsign(a*A+B)[1];
+  while (1,
+    t = (a+b) / 2;
     d = qfsign(t*A+B)[1];
     if( (ceil(n/2) - c) * (ceil(n/2) - d) < 0, b = t, a = t);
-    t = (a+b)/2;
-    q = qfsign(t*A+B)[1]
+    if (t != 0 && abs(d - n/2) <= 1, break);
+  );
+  q = denominator(t);
+  for (i = 1, oo,
+    my(t2 = bestappr(t, 2^i));
+    if (denominator(t2) > q, break);
+    if (abs(qfsign(t2*A + B)[1]) - n/2 <= 1, t = t2; break);
   );
-  u = numerator(t);
-  v= denominator(t);
-  return([u,v]);
+  return ([numerator(t), denominator(t)]);
 }
 
 WittPinceau(A,B,r)=

programme/Qfn.gp

@@ -45,7 +45,7 @@ Qfn(A,B,r)=
     );
     cpt++;
   );    
-  print("compteur ", cpt );
+  if (mm < 5, print("compteur = ", cpt));
   if (type(X) == "t_MAT", X = X[,1]);
   return(X~ * P[1..mm,]);
 }

test.in

@@ -10,6 +10,7 @@ default(parisize, "100M");
 
 \r libre/Qfmodife.gp
 \r programme/Qfn.gp
+\r programme/Pinceau.gp
 \r programme/groupe_orthogonal.gp
 \r programme/Changements_de_bases.gp
 \r qf13.gp
@@ -25,3 +26,12 @@ a = Alea(13,10^5);
 b = Alea(13,10^5);
 
 Qfn(a,b, 5)
+
+{
+for(i=1,50,
+  print(i,": ");
+  a = Alea(13,10^3);
+  b = Alea(13,10^3);
+  Qfn(a,b, 5)
+)
+}

test.out

@@ -29,27 +29,14 @@ Chargement de gramschmidt.gp
 - delta_diag(Qdiag,i)
 
 Chargement de minim_p.gp
-1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 2
-9 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54
- 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 
-80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 
-104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 
-123 124 125 126 compteur 0
-time = 35 ms.
+count = 126
+time = 39 ms.
 [3119900795401839, 1280939358528426, 782263412998367, -1149573585411029, -28
 87941316713948, 792878712740508, -3155975643213335, 7866565031850939/2, -521
 2658107497085/2, 252813867674597/2, -1065691787770711/2, -3360313638892447/2
 , -3752350536797747]
-1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 2
-9 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54
- 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 
-80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 
-104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 
-123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 
-142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 
-161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 
-180 181 182 183 184 185 186 187 188 189 compteur 0
-time = 1,472 ms.
+count = 189
+time = 1,492 ms.
 [301153599971557349161330412933556411633558600263976918969038064486487359836
 9200697513646237373148808195467377602229674073496665514632108014573090091542
 3239887278869067923359604887982108120471659836826112503477400624512212225557
@@ -100,3 +87,104 @@ time = 1,472 ms.
 3106242449187279939067399649845620284466317962213192696544391119500460373933
 3016117693109837142506251514633418571262138985416547528427623789791289300327
 840568156543265038974329369692479978557306786450289779056915188292711633444]
+1: 
+count = 335
+2: 
+count = 310
+3: 
+count = 297
+4: 
+count = 469
+5: 
+count = 250
+6: 
+count = 465
+7: 
+count = 332
+8: 
+count = 78
+9: 
+count = 337
+10: 
+count = 279
+11: 
+count = 456
+12: 
+count = 324
+13: 
+count = 170
+14: 
+count = 28
+15: 
+count = 1120
+16: 
+count = 2795
+17: 
+count = 268
+18: 
+count = 506
+19: 
+count = 18
+20: 
+count = 704
+21: 
+count = 438
+22: 
+count = 106
+23: 
+count = 124
+24: 
+count = 916
+25: 
+count = 75
+26: 
+count = 286
+27: 
+count = 32
+28: 
+count = 83
+29: 
+count = 70
+30: 
+count = 823
+31: 
+count = 522
+32: 
+count = 702
+33: 
+count = 8
+34: 
+count = 354
+35: 
+count = 11
+36: 
+count = 347
+37: 
+count = 7
+38: 
+count = 929
+39: 
+count = 858
+40: 
+count = 183
+41: 
+count = 284
+42: 
+count = 168
+43: 
+count = 1128
+44: 
+count = 399
+45: 
+count = 17
+46: 
+count = 202
+47: 
+count = 485
+48: 
+count = 41
+49: 
+count = 326
+50: 
+count = 186
+time = 27,597 ms.