tests/rtest_algsys.mac

   1 /*************** -*- Mode: MACSYMA; Package: MAXIMA -*-  ******************/
   2
   3 (kill(all), done);
   4 done;
   5
   6 /* Examples from Beyer (1984), an early paper on algsys, which references
   7    an internal General Motors report by Morgan (1981).  Morgan (1983) has
   8    the same title and contains some of the problems.
   9
  10    Beyer, William A., Solution of Simultaneous Polynomial Equations by
  11    Elimination in MACSYMA, Proceedings of the 1984 MACSYMA Users' Conference,
  12    p 110-120
  13
  14    Morgan, Alexander P., A Method for Computing All Solutions to Systems of
  15    Polynomials Equations, ACM Transactions on Mathematical Software,
  16    9:1-17 (1983)
  17  */
  18
  19 /* Beyer example 1, p117 */
  20 algsys([x^2+x*y-1,y^2+x-5],[x,y]);
  21 [[x = -3.174087266251113,y = 2.859036144578313],
  22  [x = 2.152189020152884,y = -1.687545787545788],
  23  [x = 0.3937104979189148,y = 2.146226811046756],
  24  [x = -0.3718122830843519,y = -2.317717206132879]];
  25
  26 /* Beyer example 2, p117,  Morgan problem 1 */
  27 algsys([4*x^3-3*x-y,x^2-y],[x,y]);
  28 [[x = 1,y = 1],[x = -3/4,y = 9/16],[x = 0,y = 0]];
  29
  30 /* Beyer example 3, p117, Morgan problem 2 */
  31 algsys([4*(x+y),4*(x+y)+(x-y)*((x-2)^2+y^2-1)],[x,y]);
  32 [[x = -((sqrt(2)*%i-2)/2), y =   (sqrt(2)*%i-2)/2],
  33  [x =   (sqrt(2)*%i+2)/2,  y = -((sqrt(2)*%i+2)/2)],
  34  [x = 0,y = 0]];
  35
  36 /* Beyer example 4, p117. Don't understand notation in paper */
  37
  38 /* Beyer example 5, p118.  2 real and 6 complex solutions */
  39 algsys([
  40  x^2+2*y^2-4,
  41  x^2+y^2+z-8,
  42  (x-1)^2+(2*y-sqrt(2))^2+(z-5)^2-4],
  43  [x,y,z]);
  44 [[x = 0, y = sqrt(2), z = 6],
  45  [x = 2.998824877038494 - 0.879551236563151*%i,
  46   y = (- 1.656382492469821*%i) - 0.796198988104073,
  47   z = 2.637620128835544*%i + 1.890329867297412],
  48  [x = 0.879551236563151*%i + 2.998824877038494,
  49   y = 1.656382492469821*%i - 0.7961989881040729,
  50   z = 1.890329867297412 - 2.637620128835543*%i],
  51  [x = (- 1.511881394988993*%i) - 2.77266624385368,
  52   y = 1.573386231515891*%i - 1.332140330399878,
  53   z = 3.299053626354102 - 4.191942508596404*%i],
  54  [x = 1.511881394988993*%i - 2.77266624385368,
  55   y = (- 1.573386231515895*%i) - 1.332140330399874,
  56   z = 4.191942508596393*%i + 3.299053626354096],
  57  [x = (- 1.06051780109088*%i) - 1.226158633184815,
  58   y = 1.421232537317398 - 0.4574772330741239*%i,
  59   z = 5.81061650634849 - 1.30036305745376*%i],
  60  [x = 1.06051780109088*%i - 1.226158633184815,
  61   y = 0.4574772330741241*%i + 1.421232537317397,
  62   z = 1.300363057453758*%i + 5.810616506348491],
  63  [x = 2, y = 0, z = 4]];
  64
  65 /* Beyer example 6, p118 */
  66 block([%rnum:0],
  67   algsys([x+10*y,z+w,(z-2*z)^2,(z-w)^2],[w,x,y,z]));
  68 [[w = 0, x = %r1 , y = -(%r1/10), z = 0]];
  69
  70 /* Beyer example 7, p118.  Supposed to be two real solutions.
  71    This solution seems correct.  Typo in paper? */
  72 algsys([
  73  x+y+z+w-1,
  74  x+y-z+w-3,
  75  x^2+y^2+z^2+w^2-4,
  76  (x-1)^2+y^2+z^2+w^2],
  77  [w,x,y,z]);
  78 [[w = -((3^(3/2)*%i+1)/4), x = 5/2, y =   (3^(3/2)*%i-1)/4,  z = -1],
  79  [w =   (3^(3/2)*%i-1)/4,  x = 5/2, y = -((3^(3/2)*%i+1)/4), z = -1]];
  80
  81 /* Check Beyer example 7. Preprocess with poly_reduced_grobner(). OK */
  82 algsys([z+1,2*x-5,4*y^2+2*y+7,2*y+2*w+1],[w,x,y,z]);
  83 [[w = -((3^(3/2)*%i+1)/4), x = 5/2, y =   (3^(3/2)*%i-1)/4,  z = -1],
  84  [w =   (3^(3/2)*%i-1)/4,  x = 5/2, y = -((3^(3/2)*%i+1)/4), z = -1]];
  85
  86 /* Beyer example 8, p118. Two real solutions */
  87 algsys([
  88  x1^2+x2^2+x3^2+x4-6,
  89  x1+x2+x5+x6-1,
  90  x1+3*x3+x6-3,
  91  x1+x2+x4-2,
  92  2*x3+x6,
  93  x5+x6],
  94  [x1,x2,x3,x4,x5,x6]);
  95 [[x1 = 1,  x2 = 0,   x3 = 2,  x4 = 1, x5 = 4,   x6 = -4],
  96  [x1 = 5/3,x2 = -2/3,x3 = 4/3,x4 = 1, x5 = 8/3, x6 = -8/3]];
  97
  98 /* Beyer example 9, p118. Morgan problem 4
  99    Three real and two complex solutions */
 100 algsys(
 101  [x1+x1+x2+x3+x4+x5-6=0,
 102   x2+x1+x2+x3+x4+x5-6=0,
 103   x3+x1+x2+x3+x4+x5-6=0,
 104   x4+x1+x2+x3+x4+x5-6=0,
 105   x1*x2*x3*x4*x5-1=0],
 106  [x1,x2,x3,x4,x5]);
 107 [[x1 = 0.916354556803995,
 108   x2 = 0.916354556803995,
 109   x3 = 0.916354556803995,
 110   x4 = 0.916354556803995,
 111   x5 = 1.418227215980025],
 112  [x1 = -0.5790430889759374,
 113   x2 = -0.5790430889759374,
 114   x3 = -0.5790430889759374,
 115   x4 = -0.5790430889759374,
 116   x5 = 8.895215311004785],
 117  [x1 = (-0.6100917318163317*%i)-0.06865574701986676,
 118   x2 = (-0.6100917318163317*%i)-0.06865574701986676,
 119   x3 = (-0.6100917318163317*%i)-0.06865574701986676,
 120   x4 = (-0.6100917318163317*%i)-0.06865574701986676,
 121   x5 = 3.050458659081659*%i+6.343278735099334],
 122  [x1 = 0.6100917318163317*%i-0.06865574701986676,
 123   x2 = 0.6100917318163318*%i-0.06865574701986678,
 124   x3 = 0.6100917318163317*%i-0.06865574701986676,
 125   x4 = 0.6100917318163317*%i-0.06865574701986676,
 126   x5 = 6.343278735099333-3.050458659081658*%i],
 127  [x1 = 1, x2 = 1, x3 = 1, x4 = 1, x5 = 1]];
 128
 129 /* Beyer example 10, p118. 4 real one-parameter solutions */
 130 block([%rnum:0],
 131   algsys([x^2+y^2-1,x^2+y^2+z^2-5=0],[x,y,z]));
 132 [[x = %r1, y = -sqrt(1-%r1^2), z =  2],
 133  [x = %r2, y =  sqrt(1-%r2^2), z =  2],
 134  [x = %r3, y = -sqrt(1-%r3^2), z = -2],
 135  [x = %r4, y =  sqrt(1-%r4^2), z = -2]];
 136
 137 /* Morgan problem 3.  Single solution (0,0,0,0) */
 138 algsys([x+10*y,sqrt(5)*(z-w),(y-2*z)^2,sqrt(10)*(x-w)^2],[w,x,y,z]);
 139 [[w = 0, x = 0, y = 0, z = 0]];
 140
 141 /* SF bug #3208 Can't solve this simple equation */
 142 algsys([x^2+y^2=1,(x-2)^2+(y-3)^2=9],[x,y]);
 143 [[x = -((3^(5/2)-10)/26),y = (2*3^(3/2)+15)/26],
 144  [x = (3^(5/2)+10)/26,y = -((2*3^(3/2)-15)/26)]];
 145
 146 /* Reduced from SF bug #3208 above */
 147 algsys([676*y^2-20*3^(5/2)-333,676*y^2-4056*y+28*3^(7/2)+2007],[y]);
 148 [[y = (2*3^(3/2)+15)/26]];
 149
 150 /* Reduced from SF bug #3208 above */
 151 algsys([676*y^2+20*3^(5/2)-333,676*y^2-4056*y-28*3^(7/2)+2007],[y]);
 152 [[y = -((2*3^(3/2)-15)/26)]];
 153
 154 /* from example(solve) */
 155 algsys([4*x^2-y^2=12,x*y-x=2],[x,y]);
 156 [[x = 2, y = 2],
 157  [x = 0.5202594388652008*%i - 0.1331240357358706,
 158   y = 0.07678378523787788 - 3.608003221870287*%i],
 159  [x = -0.5202594388652008*%i - 0.1331240357358706,
 160   y = 3.608003221870287*%i + 0.07678378523787788],
 161  [x = -1.733751846381093, y = -0.1535675710019696]];
 162
 163 /* failures from example(algsys) during testing */
 164 block( [
 165      %rnum:0,
 166      f1:2*x*(1-l1)-2*(x-1)*l2,
 167      f2:l2-l1,
 168      f3:l1*(1-x**2-y),
 169      f4:l2*(y-(x-1)**2)
 170    ],
 171    algsys([f1,f2,f3,f4],[x,y,l1,l2]));
 172 [[x = 0, y = %r1, l1 = 0, l2 = 0],
 173  [x = 1, y = 0,   l1 = 1, l2 = 1]];
 174
 175 block( [
 176    f1:x**2-y**2,
 177    f2:x**2-x+2*y**2-y-1
 178  ],
 179  algsys([f1,f2],[x,y]));
 180 [[x = -(1/sqrt(3)), y =   1/sqrt(3)],
 181  [x =   1/sqrt(3),  y = -(1/sqrt(3))],
 182  [x = -(1/3),y = -(1/3)],
 183  [x = 1,y = 1]];
 184
 185 /* failure in example(charpoly) during testing */
 186  algsys([x2-2*x1=0,x2^2+x1^2=1],[x1,x2]);
 187  [[x1 = -(1/sqrt(5)), x2 = -(2/sqrt(5))],
 188   [x1 =   1/sqrt(5),  x2 =   2/sqrt(5)]];
 189
 190 block( [%rnum:0],
 191   algsys([g526+((sqrt(a^2+4)-a)/2+a)*g525,
 192           ((sqrt(a^2+4)-a)*g526)/2+g525],
 193          [g525,g526]));
 194 [[g525 = %r1,g526 = -((%r1*sqrt(a^2+4)+%r1*a)/2)]];
 195
 196 block( [%rnum:0],
 197   algsys([g625+(a-(sqrt(a^2+4)+a)/2)*g624,
 198            g624-((sqrt(a^2+4)+a)*g625)/2],
 199            [g624,g625]));
 200 [[g624 = %r1,g625 = (%r1*sqrt(a^2+4)-%r1*a)/2]];
 201
 202 /* Working example from bug #3150 */
 203 algsys([x^2+y^2=1, x^2+(y-3)^2=9],[x,y]);
 204 [[x = -sqrt(35)/6,y = 1/6],
 205  [x = sqrt(35)/6,y = 1/6]];
 206
 207 /* Working example from bug #2736 */
 208 algsys([x^2+y^2=1,(x-2)^2+(y-3)^2=10],[x,y]);
 209 [[x = 1,y = 0],[x = -5/13,y = 12/13]];
 210
 211 /* bug #1106 only found 4 roots (of 6)
 212    Now crash with heap exhausted
 213  */
 214 /*
 215 block( [
 216   p1:-x*y^3+y^2+x^4-9*x/8,
 217   p2:y^4-x^3*y-9*y/8+x^2],
 218   algsys([p1,p2],[x,y]));
 219 [[x = 1/2,y = 1],
 220  [x = 9/8,y = 9/8],
 221  [x = 1,y = 1/2],
 222  [x = 0,y = 0]]
 223  */
 224
 225 /* Above can be solved with
 226    load(grobner)$
 227    eqs:poly_reduced_grobner([p1,p2],[x,y])$
 228    algsys(eqs,[x,y]);
 229  */
 230 algsys([-200704*y^8+437832*y^5-263633*y^2+53010*x,
 231         4096*y^10-10440*y^7+7073*y^4-729*y],[x,y]);
 232 [[x = 0, y = 0],
 233  [x = 1, y = 1/2],
 234  [x = 9/8, y = 9/8],
 235  [x = 1/2, y = 1],
 236  [x = (sqrt(3)*%i-1)/4,     y = -((sqrt(3)*%i+1)/2)],
 237  [x = -((sqrt(3)*%i+1)/4),  y = (sqrt(3)*%i-1)/2],
 238  [x = (3^(5/2)*%i-9)/16,    y = -((3^(5/2)*%i+9)/16)],
 239  [x = -((3^(5/2)*%i+9)/16), y = (3^(5/2)*%i-9)/16],
 240  [x = (sqrt(3)*%i-1)/2,     y = -((sqrt(3)*%i+1)/4)],
 241  [x = -((sqrt(3)*%i+1)/2),  y = (sqrt(3)*%i-1)/4]];
 242
 243 /* bug #1266 no solution found */
 244 algsys([x^2+x+1=0,x^2*y+1=0],[x,y]);
 245 [[x =   (sqrt(3)*%i-1)/2,  y = -((sqrt(3)*%i-1)/2)],
 246  [x = -((sqrt(3)*%i+1)/2), y =   (sqrt(3)*%i+1)/2]];
 247
 248 /* bug #1302 no solution found.  Also bug #1469 */
 249 algsys([(x-x0)^2+(y-y0)^2=r^2,(x-x1)^2+(y-y1)^2=r^2],[x,y]);
 250 []; /* wrong, but at least it returns */
 251
 252 /* bug #1369 algsys hangs */
 253 /*
 254 eq1 : a*x + c*y + d*y^2/2 = 0;
 255 eq2 : b*x + e*x^2 + f*y - g*y^3 = h;
 256 algsys([eq1,eq2],[x,y]);
 257 */
 258
 259 /* bug #1370 Can solve [F=0,F-G=0] but not [F,G]
 260    there are some work arounds in the report
 261  */
 262 F:2*x^3*y^2+2*a*x^2*y^2+3*x^2*y+2*a*x*y+x+1;
 263 2*x^3*y^2+2*a*x^2*y^2+3*x^2*y+2*a*x*y+x+1;
 264 G:x^3*y^2+2*a*x^2*y^2+a^2*x*y^2+x^2*y+2*a*x*y+a^2*y+1;
 265 x^3*y^2+2*a*x^2*y^2+a^2*x*y^2+x^2*y+2*a*x*y+a^2*y+1;
 266 /* FIXME:  Testsuite issue.  Can't match form of solution
 267 algsys([F=0,F-G=0],[x,y]);
 268 [[x = -(((-a^3)+sqrt(a^2-4*a)*(a^2-2*a)+4*a^2)/(2*sqrt(a^2-4*a))),
 269   y = -(1/(a*sqrt(a^2-4*a)))],
 270  [x = -((a^3+sqrt(a^2-4*a)*(a^2-2*a)-4*a^2)/(2*sqrt(a^2-4*a))),
 271   y = 1/(a*sqrt(a^2-4*a))]];
 272 */
 273
 274 /* hangs
 275 algsys([F=0,G=0],[x,y]);
 276 [];
 277  */
 278
 279 /* Can solve with
 280    eqs:poly_reduced_grobner([F,G],[x,y])
 281    algsys:(eqs,[x,y]);
 282  */
 283 /* FIXME:  Testsuite issue.  Can't match form of solution
 284 algsys([(a^4-4*a^3)*y^2-1,(a^4-4*a^3)*y+2*x+a^2-2*a],[x,y]);
 285 [[x = -((a*sqrt(a^2-4*a)+a^2-2*a)/2), y =   sqrt(a^2-4*a)/(a^3-4*a^2)],
 286  [x =   (a*sqrt(a^2-4*a)-a^2+2*a)/2,  y = -(sqrt(a^2-4*a)/(a^3-4*a^2))]];
 287 */
 288
 289 /* Bug #1647 */
 290 algsys([x^2+y^2-d^2, (x-h)^2+(y-k)^2-s^2], [x,y]);
 291 [];  /* no solution found, but at least it returns */
 292
 293 /*
 294   Here is a workaround
 295   load(grobner);
 296   [x^2+y^2-d^2, (x-h)^2+(y-k)^2-s^2];
 297   poly_reduced_grobner(%, [x,y]);
 298   algsys(%, [x,y]);
 299  */
 300 /* FIXME:  Testsuite issue.  Can't match form of solution
 301 algsys([2*k*y+2*h*x+s^2-k^2-h^2-d^2,
 302         (4*k^2+4*h^2)*y^2+(k*(4*s^2-4*d^2)-4*k^3-4*h^2*k)*y+s^4-2*d^2*s^2
 303           +h^2*((-2*s^2)+2*k^2-2*d^2)+k^2*(2*d^2-2*s^2)+k^4+h^4+d^4],[x,y]);
 304 [[x = -(((k*sqrt((-s^4)+(2*k^2+2*h^2+2*d^2)*s^2-k^4+(2*d^2-2*h^2)*k^2-h^4
 305                            +2*d^2*h^2-d^4)
 306           +h*s^2-(h*k^2)-h^3-(d^2*h))
 307           /(2*k^2+2*h^2))),
 308   y = (h*sqrt((-s^2)+2*d*s+k^2+h^2-d^2)*sqrt(s^2+2*d*s-k^2-h^2+d^2)
 309           -k*s^2+k^3+(h^2+d^2)*k)
 310           /(2*k^2+2*h^2)],
 311  [x = (k*sqrt((-s^4)+(2*k^2+2*h^2+2*d^2)*s^2-k^4+(2*d^2-2*h^2)*k^2-h^4
 312                           +2*d^2*h^2-d^4)
 313           -h*s^2+h*k^2+h^3+d^2*h)
 314           /(2*k^2+2*h^2),
 315   y = -(((h*sqrt((-s^2)+2*d*s+k^2+h^2-d^2)*sqrt(s^2+2*d*s-k^2-h^2+d^2)
 316           +k*s^2-k^3+((-h^2)-d^2)*k)
 317           /(2*k^2+2*h^2)))]];
 318  */
 319
 320 /* SF bug #3210 ALGSYS does not return (regression)
 321  *   also #3212 ALGSYS throws error
 322  *
 323  * Note: %nrum is the counter for the %r variables introduced
 324  *       by algsys.  Set it to 0 for reproducible results
 325  */
 326 block([%rnum:0], algsys([
 327  4*y*(y^2-28),
 328  8*y*(y^2-28)*z,
 329  -16*y*(y^2-23)*z,
 330  8*y*(y^2-10)*z,
 331  32*y*(z-y),
 332  32*y*(z+y),
 333  4*y*(y^2-28),
 334  8*y*(y^2-28)*z,
 335  -16*y*(y^2-23)*z,
 336  8*y*(y^2-10)*z,
 337  y*((4*y^2-112)*z^2-8*y^4+(185-4*x^2)*y^2),
 338  -2*y*((4*y^2-92)*z^2-2*y^4+(21-4*x^2)*y^2),
 339  y*((4*y^2-40)*z^2+(1-4*x^2)*y^2)],
 340 [x,y,z]));
 341 [[x=%r1,y=0,z=%r2]];
 342
 343 block([%rnum:0], algsys([
 344   -8*a*z,
 345   -3*(y-1)*y*z,-3*z*(z+y-1)*(z+y),
 346   2*(y-1)*y*z*(6*y*z-3*z-4*b*y^2-2*y^2+4*b*y+2*y),
 347   -8*a*z,
 348   -3*(y-1)*y*z,
 349   -4*z*((w^2-2*v*w+v^2+2*b-2*a)*z+8*a*y-4*a),
 350   -2*(y-1)*y*z*((6*y-3)*z+((-4*b)-2)*y^2+(4*b+2)*y),
 351   4*z
 352    *((2*w^2-4*v*w+2*v^2+2*b-1)*z^2
 353     +(((-2*w^2)+4*v*w-2*v^2-8*b+8*a-2)*y+(2*w+2*v-2)*x-4*v*w+4*b-4*a+2)
 354      *z-12*a*y^2+12*a*y-2*a),
 355   -z*(((4*w^2-8*v*w+4*v^2+14)*y^2+(((-8*w)-8*v+8)*x+16*v*w-18)*y+4*x^2
 356                                  -8*x+3)
 357      *z^2
 358      +(((-32*b)-16)*y^3+(48*b+24)*y^2+((-16*b)-8)*y)*z+(8*b-8*a+4)*y^4
 359      +((-16*b)+16*a-8)*y^3+(8*b-8*a+4)*y^2),
 360   -z*((4*w^2-8*v*w+4*v^2-1)*z^3+(((-16*w^2)+32*v*w-16*v^2-32*b)*y
 361                                 +(16*w+16*v-16)*x-32*v*w+16*b+8)
 362                                 *z^2
 363                                +((4*w^2-8*v*w+4*v^2+48*b-48*a+20)*y^2
 364                                 +(((-8*w)-8*v+8)*x+16*v*w-48*b+48*a-24)
 365                                  *y+4*x^2-8*x+8*b-8*a+4)
 366                                 *z+32*a*y^3-48*a*y^2+16*a*y),
 367   -2*z
 368     *(((4*w^2-8*v*w+4*v^2+2)*y+((-4*w)-4*v+4)*x+8*v*w-3)*z^3
 369      +(((-4*w^2)+8*v*w-4*v^2-24*b-8)*y^2
 370       +((8*w+8*v-8)*x-16*v*w+24*b+12)*y-4*x^2+8*x-4*b-2)
 371       *z^2+((16*b-16*a+8)*y^3+((-24*b)+24*a-12)*y^2+(8*b-8*a+4)*y)*z
 372      +4*a*y^4-8*a*y^3+4*a*y^2)],
 373 [v,w,x,y,z]));
 374 [[v=%r1,w=%r2,x=%r3,y=%r4,z=0]];
 375
 376 /* Test problems distilled from share/contrib/diffequations tests
 377
 378    The odelin method for linear second order ODEs calls algys.
 379    Changes to algsys capability change the odelin results.
 380    A number of tests below were obtained by tracing algsys
 381    while solving ODES with odelin (or contrib_ode).
 382  */
 383
 384 /* Solved by odelin for Kamke ODE 2.265.  New solution 2016-09 */
 385  algsys(
 386  [y = -1,
 387   z = 1,
 388   (-w^2)+2*v*w-v^2+5,2*((w+v-1)*x+4*w^2-10*v*w+4*v^2-22),
 389   (-x^2)-(12*w+12*v-12)*x+((-2*w)-2*v+2)*x+2*x-26*w^2+80*v*w-26*v^2+161,
 390   (-4*x^2)+((-8*w)-8*v+8)*x+8*x-4*w^2+24*v*w-4*v^2+48,
 391   2*(3*x^2+(13*w+13*v-13)*x+(6*w+6*v-6)*x-6*x+22*w^2-82*v*w+22*v^2-157),
 392   -2*((-6*x^2)-(12*w+12*v-12)*x+((-4*w)-4*v+4)*x+12*x-10*w^2+52*v*w
 393      -10*v^2+100),
 394   (-13*x^2)-(24*w+24*v-24)*x+((-26*w)-26*v+26)*x+26*x-41*w^2+182*v*w
 395      -41*v^2+344],
 396  [v,w,x,y,z]);
 397 [[v =   (sqrt(5)+5)/2,  w = -((sqrt(5)-5)/2), x = 3,  y = -1, z = 1],
 398  [v = -((sqrt(5)-5)/2), w =   (sqrt(5)+5)/2,  x = 3,  y = -1, z = 1],
 399  [v =   (sqrt(5)-3)/2,  w = -((sqrt(5)+3)/2), x = -1, y = -1, z = 1],
 400  [v = -((sqrt(5)+3)/2), w =   (sqrt(5)-3)/2,  x = -1, y = -1, z = 1],
 401  [v =   (sqrt(5)+1)/2,  w = -((sqrt(5)-1)/2), x = 3,  y = -1, z = 1],
 402  [v = -((sqrt(5)-1)/2), w =   (sqrt(5)+1)/2,  x = 3,  y = -1, z = 1],
 403  [v =   (sqrt(5)+1)/2,  w = -((sqrt(5)-1)/2), x = -1, y = -1, z = 1],
 404  [v = -((sqrt(5)-1)/2), w =   (sqrt(5)+1)/2,  x = -1, y = -1, z = 1]];
 405
 406 /* Solved by odelin for Murphy ODE 1.182.  New solution 2016-10 */
 407 /* FIXME:  Testsuite issue.  Can't match form of solution
 408 algsys(
 409  [x = -2*sqrt(-a1*a3),
 410   4*sqrt(-a1*a3)*((-4*a1*a3*(2*y-z))-4*a0*a3*sqrt(-a1*a3)),
 411   8*(-a1*a3)^(3/2)*(z^2-2*z-a2^2+2*a2)],
 412   [x,y,z]);
 413 [[x = -2*sqrt(-a1*a3),
 414   y = -((a0*sqrt(-a1*a3)+a1*a2-2*a1)/(2*a1)),
 415   z = 2-a2],
 416  [x = -2*sqrt(-a1*a3),
 417   y = -((a0*sqrt(-a1*a3)-a1*a2)/(2*a1)),
 418   z = a2]];
 419  */
 420
 421 /* Solved by odelin for Kamke ODE 1.105.  New solution 2016-10 */
 422 /* FIXME:  Testsuite issue.  Can't match form of solution
 423 algsys(
 424  [4*sqrt(-a*c)*((-4*a*c*(2*y-z))-4*a*sqrt(-a*c)*d),
 425   8*(-a*c)^(3/2)*(z^2-2*z-b^2+2*b)],
 426  [y,z]);
 427 [[y = -((sqrt(-a*c)*d+(b-2)*c)/(2*c)), z = 2-b],
 428  [y = -((sqrt(-a*c)*d-b*c)/(2*c)),     z = b]];
 429  */
 430
 431 /* Solved by odelin for Kamke ODE 1.291.  New solution 2016-10 */
 432 algsys([
 433  -2*(18*y^2-36*x*y+18*x^2-5),
 434  (72*y+72*x-72)*z+72*y^2-288*x*y+72*x^2+69,
 435  (-36*z^2)-(144*y+144*x-144)*z+72*z-36*y^2+360*x*y-36*x^2-203,
 436  72*z^2-((-72*y)-72*x+72)*z-144*z-144*x*y+159,(-36*z^2)+72*z-35],
 437  [x,y,z]);
 438 [[x =   (sqrt(10)-3)/12,   y = -((sqrt(10)+3)/12),  z = 5/6],
 439  [x = -((sqrt(10)+3)/12),  y =   (sqrt(10)-3)/12,   z = 5/6],
 440  [x =   (sqrt(10)+13)/12,  y = -((sqrt(10)-13)/12), z = 5/6],
 441  [x = -((sqrt(10)-13)/12), y =   (sqrt(10)+13)/12,  z = 5/6],
 442  [x =   (sqrt(10)+15)/12,  y = -((sqrt(10)-15)/12), z = 7/6],
 443  [x = -((sqrt(10)-15)/12), y =   (sqrt(10)+15)/12,  z = 7/6],
 444  [x =   (sqrt(10)-1)/12,   y = -((sqrt(10)+1)/12),  z = 7/6],
 445  [x = -((sqrt(10)+1)/12),  y =   (sqrt(10)-1)/12,   z = 7/6]];
 446
 447 /* SF bug #2038
 448    solve(eqs,[x,y]) gives correct solution
 449    solve(subst(sqrt(3),q,eqs),[x,y]) yields []!
 450 */
 451 algsys(eqs:[y^2+x^2-(400-400/(2/q+1))^2=0,
 452        (-y-400/(2/q+1)+400)^2+x^2-640000/(2/q+1)^2=0],[x,y]);
 453 [[x = -(400*q*sqrt(-((q-2)/(q+2)))), y = -((400*q^2-800)/(q+2))],
 454  [x =   400*q*sqrt(-((q-2)/(q+2))),  y = -((400*q^2-800)/(q+2))]];
 455
 456 algsys(subst(q=sqrt(3),eqs),[x,y]);
 457 [[x = 1200-800*sqrt(3),y = 400*sqrt(3)-800],
 458  [x = 800*sqrt(3)-1200,y = 400*sqrt(3)-800]];
 459
 460 algsys(subst(q=1-sqrt(3),eqs),[x,y]);
 461 [[x = -((25*2^(9/2))/3^(1/4)), y = 800/sqrt(3)],
 462  [x =   (25*2^(9/2))/3^(1/4),  y = 800/sqrt(3)]];
 463
 464 /* SF bug #2736 (in comments)
 465    Solving with p=1/(1-cos(1)) gave incorrect solution */
 466 eqs:[x*y=p,x-y=3]$
 467 [x*y = p,x-y = 3];
 468
 469 algsys(subst(p=1/(1-cos(1)),eqs),[x,y])$
 470 [[x =   -((sqrt(9*cos(1)^2-22*cos(1)+13)-3*cos(1)+3)/(2*cos(1)-2)),
 471   y =   -((sqrt(9*cos(1)^2-22*cos(1)+13)+3*cos(1)-3)/(2*cos(1)-2))],
 472  [x =     (sqrt(9*cos(1)^2-22*cos(1)+13)+3*cos(1)-3)/(2*cos(1)-2),
 473   y =     (sqrt(9*cos(1)^2-22*cos(1)+13)-3*cos(1)+3)/(2*cos(1)-2)]];
 474
 475 float(%)$
 476 [[x =  3.603649840080336, y =  0.603649840080336],
 477  [x = -0.603649840080336, y = -3.603649840080336]];
 478
 479 float(algsys(subst(p=float(1/(1-cos(1))),eqs),[x,y]))$
 480 [[x = -0.6036497963666444,y = -3.603649796366644],
 481  [x =  3.603649796366644, y =  0.6036497963666444]];
 482
 483 /* SF bug 2059 - variable order dependent
 484    No solution found for certain solution variable orders
 485  */
 486 (eqs:[z-a=%i*b,z+a=c,a*z+3*(a-z)=13-12*%i,3*(z-a)+a*z=12*%i+13],
 487 algsys(eqs,[a,b,c,z]));
 488 [[a = (-2*%i)-3, b = 4, c = -6, z = 2*%i-3],
 489  [a = 3-2*%i,    b = 4, c =  6, z = 2*%i+3]];
 490
 491 /* Count the number of solutions for each permutation of variable */
 492 map(length,
 493   map(lambda([u],algsys(eqs,u)), listify(permutations([a,b,c,z]))));
 494 [2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2];
 495
 496 /* SF bug #453 algsys fails in simple case */
 497 algsys([v^2-2*v+u^2+10*u+1,v^2-10*v+u^2+2*u+1],[v,u]);
 498 [[v = -((sqrt(34)-6)/2), u =   (sqrt(2)*sqrt(17)-6)/2],
 499  [v =   (sqrt(34)+6)/2,  u = -((sqrt(2)*sqrt(17)+6)/2)]];
 500
 501 algsys([2*v^2-12*v+1,v^2-2*v+u^2+10*u+1],[v,u]);
 502
 503 [[v = -((sqrt(2)*sqrt(17)-6)/2), u = -((sqrt(34)+14)/2)],
 504  [v =   (sqrt(2)*sqrt(17)+6)/2,  u =   (sqrt(34)-14)/2],
 505  [v =   (sqrt(2)*sqrt(17)+6)/2,  u = -((sqrt(34)+6)/2)],
 506  [v = -((sqrt(2)*sqrt(17)-6)/2), u =   (sqrt(34)-6)/2]];
 507
 508 /* Solved by odelin for Kamke ODE 2.184. */
 509 algsys(
 510 [w^2-a^2*w,2*w*x-2*a^2*x,
 511  w*(2*y+2*a*c+b^2)-3*a^2*y+x^2-2*a*b*x+A*a^2,
 512  w*(2*z+4*b*c)-4*a^2*z+2*x*y-4*a*b*y+2*A*a*b,
 513  x*(2*z+2*b*c)-6*a*b*z+y^2+((-2*a*c)-b^2)*y+3*c^2*w+2*A*a*c+A*b^2,
 514  2*y*z+((-4*a*c)-2*b^2)*z+2*c^2*x+2*A*b*c,z^2-2*b*c*z+c^2*y+A*c^2],
 515  [w,x,y,z]);
 516 [[w = a^2,
 517   x = -((a*sqrt((-4*a*c)+b^2-4*A)-3*a*b)/2),
 518   y = -((b*sqrt((-4*a*c)+b^2-4*A)-2*a*c-b^2)/2),
 519   z = -((c*sqrt((-4*a*c)+b^2-4*A)-b*c)/2)],
 520  [w = a^2,
 521   x = (a*sqrt((-4*a*c)+b^2-4*A)+3*a*b)/2,
 522   y = (b*sqrt((-4*a*c)+b^2-4*A)+2*a*c+b^2)/2,
 523   z = (c*sqrt((-4*a*c)+b^2-4*A)+b*c)/2]];
 524
 525 /* Solved by odelin for Kamke ODE 2.119. */
 526  algsys([x^2-x-b^2-b,2*x*y-2*y,y^2],[x,y]);
 527 [[x = -b, y = 0], [x = b+1, y = 0]];
 528
 529 /* Solved by odelin for Kamke ODE 2.331. */
 530 algsys([x^2-16*x+48,2*x*y-32*y-192,y^2+64*y+64*x+192],[x,y]);
 531 [[x = 12,y = -24],[x = 4,y = -8]];
 532
 533 /* Solved by odelin for Kamke ODE 2.382. */
 534 algsys([
 535  w^2-w,
 536  2*w*x-2*x,w*(2*y+b^2+4*a*b+a^2)-3*y+x^2+(2*b+2*a)*x-c,
 537  w*(2*z-4*a*b^2-4*a^2*b)-4*z+2*x*y+(4*b+4*a)*y+(2*b+2*a)*c,
 538  x*(2*z-2*a*b^2-2*a^2*b)+(6*b+6*a)*z+y^2+((-b^2)-4*a*b-a^2)*y
 539    +3*a^2*b^2*w+((-b^2)-4*a*b-a^2)*c,
 540  2*y*z+((-2*b^2)-8*a*b-2*a^2)*z+2*a^2*b^2*x+(2*a*b^2+2*a^2*b)*c,
 541  z^2+(2*a*b^2+2*a^2*b)*z+a^2*b^2*y-a^2*b^2*c],
 542 [w,x,y,z]);
 543 [[w = 1,
 544  x = -((sqrt(4*c+b^2-(2*a*b)+a^2)+3*b+3*a)/2),
 545  y = ((b+a)*sqrt(4*c+b^2-(2*a)*b+a^2)+b^2+4*a*b+a^2)/2,
 546  z = -((a*b*sqrt(4*c+b^2-(2*a*b)+a^2)+a*b^2+a^2*b)/2)],
 547 [w = 1,
 548  x = (sqrt(4*c+b^2-(2*a*b)+a^2)-3*b-3*a)/2,
 549  y = -(((b+a)*sqrt(4*c+b^2-2*a*b+a^2)-b^2-4*a*b-a^2)/2),
 550  z = (a*b*sqrt(4*c+b^2-2*a*b+a^2)-a*b^2-a^2*b)/2]];
 551
 552 /* Solved by odelin for Kamke ODE 2.383. */
 553 algsys([
 554   w^2-4*w,
 555   2*w*x-8*x,
 556   w*(2*y+4*b^2+16*a*b+4*a^2)-12*y+x^2+(8*b+8*a)*x+b*(8*a*b^2-16*a^2*b+8*a^3)
 557                             +(4-4*a^2)*b^2+b^2*((-4*b^2)+8*a*b-4*a^2)
 558                             +(8*a^3-8*a)*b-4*a^4+4*a^2,
 559   w*(2*z-16*a*b^2-16*a^2*b)-16*z+2*x*y+(16*b+16*a)*y
 560                            +b*((-16*a*b^3)+16*a^2*b^2+16*a^3*b-16*a^4)
 561                            +b^2*(8*b^3-8*a*b^2-8*a^2*b+8*a^3)+(8*a^2-8)*b^3
 562                            +(8*a-8*a^3)*b^2+(8*a^2-8*a^4)*b+8*a^5-8*a^3,
 563   x*(2*z-8*a*b^2-8*a^2*b)+(24*b+24*a)*z+y^2+((-4*b^2)-16*a*b-4*a^2)*y
 564                          +12*a^2*b^2*w
 565                          +b*(8*a*b^4+16*a^2*b^3-48*a^3*b^2+16*a^4*b+8*a^5)
 566                          +(4-4*a^2)*b^4
 567                          +b^2*((-4*b^4)-8*a*b^3+24*a^2*b^2-8*a^3*b-4*a^4)
 568                          +(8*a-8*a^3)*b^3+(24*a^4-24*a^2)*b^2+(8*a^3-8*a^5)*b
 569                          -4*a^6+4*a^4,
 570   2*y*z+((-8*b^2)-32*a*b-8*a^2)*z+8*a^2*b^2*x
 571        +b*((-16*a^2*b^4)+16*a^3*b^3+16*a^4*b^2-16*a^5*b)
 572        +b^2*(8*a*b^4-8*a^2*b^3-8*a^3*b^2+8*a^4*b)+(8*a^3-8*a)*b^4
 573        +(8*a^2-8*a^4)*b^3+(8*a^3-8*a^5)*b^2+(8*a^6-8*a^4)*b,
 574   z^2+(8*a*b^2+8*a^2*b)*z+4*a^2*b^2*y+b*(8*a^3*b^4-16*a^4*b^3+8*a^5*b^2)
 575      +b^2*((-4*a^2*b^4)+8*a^3*b^3-4*a^4*b^2)+(4*a^2-4*a^4)*b^4
 576      +(8*a^5-8*a^3)*b^3+(4*a^4-4*a^6)*b^2],
 577  [w,x,y,z]);
 578 [[w = 4,
 579   x = 2*b^2+((-4*a)-6)*b+2*a^2-6*a,
 580   y = (-2*b^3)+(2*a+2)*b^2+(2*a^2+8*a)*b-2*a^3+2*a^2,
 581   z = 2*a*b^3+((-4*a^2)-2*a)*b^2+(2*a^3-2*a^2)*b],
 582  [w = 4,
 583   x = (-2*b^2)+(4*a-6)*b-2*a^2-6*a,
 584   y = 2*b^3+(2-2*a)*b^2+(8*a-2*a^2)*b+2*a^3+2*a^2,
 585   z = (-2*a*b^3)+(4*a^2-2*a)*b^2+((-2*a^3)-2*a^2)*b]];
 586
 587 /* Eugene L. Allgower, Kurt Georg and Rick Miranda, The Method of
 588    Resultants for Computing Real Solutions of Polynomial Systems,
 589    SIAM Journal on Numerical Analysis, Vol 29 No 3, Jun 1992, pp 831-844
 590  */
 591
 592 /* Example 1, p840, qe 12.  Robot arm coordinates.
 593    Real solutions have y = +/- 1.6237
 594  */
 595 (p1:x^2*y^2+(2/q2)*x^2*y+(2/q2)*(1-b)*x*y^2+x^2+y^2-(2/q2)*(1+b)*x-(2/q2)*y+1,
 596  p2:(1-b-q1)*x^2*y^2+(-1-b-q1)*x^2-4*x*y+(-1+b-q1)*y^2+(1+b-q1),
 597  algsys(subst([b=1/2,q1=4/5,q2=2/5],[p1,p2]),[x,y]));
 598 [[x = %i,y = -%i],
 599  [x = -%i,y = %i],
 600  [x = %i,y = -3*%i],
 601  [x = -%i,y = 3*%i],
 602  [x =   (sqrt(319)+8)/17,  y = -(sqrt(29)/sqrt(11))],
 603  [x = -((sqrt(319)-8)/17), y =   sqrt(29)/sqrt(11)]];
 604 float([%[5],%[6]]);
 605 [[x = 1.521210064675985, y = -1.623688281771977],
 606  [x = -0.58003359408775, y =  1.623688281771977]];
 607
 608 /* Example 2, p 841.  Intersection of a sphere with two paraboloids,
 609    Two real solutions (not found by maxima 5.38.0)
 610    y = z = (sqrt(5)-1)/2, x = +/- sqrt(z-y^2)
 611  */
 612 algsys([x^2+y^2+z^2-1,z-x^2-y^2,y-x^2-z^2],[x,y,z]);
 613 [[x =  sqrt(2)*%i,         y =   (sqrt(5)-1)/2,  z = -((sqrt(5)+1)/2)],
 614  [x =  sqrt(2)*%i,         y = -((sqrt(5)+1)/2), z =   (sqrt(5)-1)/2],
 615  [x = -sqrt(2)*%i,         y =   (sqrt(5)-1)/2,  z = -((sqrt(5)+1)/2)],
 616  [x = -sqrt(2)*%i,         y = -((sqrt(5)+1)/2), z =   (sqrt(5)-1)/2],
 617  [x =  sqrt((-sqrt(5))-2), y = -((sqrt(5)+1)/2), z = -((sqrt(5)+1)/2)],
 618  [x = -sqrt((-sqrt(5))-2), y = -((sqrt(5)+1)/2), z = -((sqrt(5)+1)/2)],
 619  [x =  sqrt(sqrt(5)-2),    y =   (sqrt(5)-1)/2,  z =   (sqrt(5)-1)/2],
 620  [x = -sqrt(sqrt(5)-2),    y =   (sqrt(5)-1)/2,  z =   (sqrt(5)-1)/2]];
 621
 622 /* Two examples from the maxima mailing list 2016-10-15 */
 623
 624 algsys(eqs:[x+2*y-z=6,2*x+y*z-z^2=-1,3*x-y+2*z^2=3],[x,y,z]);
 625 [[x = (27*sqrt(13)*%i-41)/14,
 626   y = -((17*sqrt(13)*%i-87)/14),
 627   z = -((sqrt(13)*%i-7)/2)],
 628  [x = -((27*sqrt(13)*%i+41)/14),
 629   y = (17*sqrt(13)*%i+87)/14,
 630   z = (sqrt(13)*%i+7)/2],
 631  [x = 1, y = 2, z = -1]];
 632
 633 block([%rnum:0],
 634   algsys(eqs:[b*c+a^2-1=0,b*d+a*b=0,c*d+a*c=0,d^2+b*c-1=0],[a,b,c,d]));
 635 [[a = -sqrt(1-%r1*%r2), b = %r1, c = %r2, d =  sqrt(1-%r1*%r2)],
 636  [a =  sqrt(1-%r3*%r4), b = %r3, c = %r4, d = -sqrt(1-%r3*%r4)],
 637  [a =  1, b = %r5, c = 0, d = -1],
 638  [a = -1, b = %r6, c = 0, d =  1],
 639  [a = -1, b = 0,   c = 0, d = -1],
 640  [a =  1, b = 0,   c = 0, d =  1]];
 641
 642 /* Example from Stavros Macrakis, maxima mailing list, 2017-07-31
 643    Not solved by 5.38.0, solved by 5.39.0 */
 644 block([%rnum:0],
 645  algsys(eqs:[b*c+a^2,b*d+a*b,c*d+a*c,d^2+b*c],[a,b,c,d]));
 646 [[a = -sqrt(-%r1*%r2), b = %r1, c = %r2, d =  sqrt(-%r1*%r2)],
 647  [a =  sqrt(-%r3*%r4), b = %r3, c = %r4, d = -sqrt(-%r3*%r4)],
 648  [a =  0,              b = %r5, c = 0,   d = 0]];
 649
 650
 651 /* SF bug 3321 - apparent regression from 5.38.1
 652
 653    In fact, correct solution depends on the signs of constants g and h.
 654    Version 5.38.1 returns wrong solution when g>0 and h<0.
 655    Would be better if algsys asked for signs, but ...
 656  */
 657
 658 eqs:[y-x=0, g*x*y-h=0, z+(x+1)/y-x-1=0];
 659 [y-x=0, g*x*y-h=0, z+(x+1)/y-x-1=0];
 660
 661 (assume(h>0),0);
 662 0;
 663 s:algsys(eqs,[x,y,z]);
 664 [[x = sqrt(h)/sqrt(g),
 665   y = sqrt(h)/sqrt(g),
 666   z = (h-g)/(sqrt(g)*sqrt(h))],
 667  [x = -(sqrt(h)/sqrt(g)),
 668   y = -(sqrt(h)/sqrt(g)),
 669   z = -((h-g)/(sqrt(g)*sqrt(h)))]];
 670 ratsimp(subst(s[1],eqs));
 671 [0=0, 0=0, 0=0];
 672 ratsimp(subst(s[2],eqs));
 673 [0=0, 0=0, 0=0];
 674 (forget(h>0),0);
 675 0;
 676
 677 (assume(h<0),0);
 678 0;
 679 s:algsys(eqs,[x,y,z]);
 680 [[x = -((%i*sqrt(-h))/sqrt(g)),
 681   y = -((%i*sqrt(-h))/sqrt(g)),
 682   z = -((sqrt(-h)*(%i*h-%i*g))/(sqrt(g)*h))],
 683  [x = (%i*sqrt(-h))/sqrt(g),
 684   y = (%i*sqrt(-h))/sqrt(g),
 685   z = (sqrt(-h)*(%i*h-%i*g))/(sqrt(g)*h)]];
 686 ratsimp(subst(s[1],eqs));
 687 [0=0, 0=0, 0=0];
 688 ratsimp(subst(s[2],eqs));
 689 [0=0, 0=0, 0=0];
 690 (forget(h<0),0);
 691 0;
 692
 693 /* SF bug 3375.  Failure to calculate eigenvalues
 694    Solution is ugly, so just check it contains %r1 */
 695 block([%rnum:0], s:algsys(eqs:[
 696   y/10-(((-(sqrt(3)*%i)/2)-1/2)*((sqrt(37)*%i)/(8*3^(3/2))+1/8)^(1/3)
 697    +((sqrt(3)*%i)/2-1/2)/(3*((sqrt(37)*%i)/(8*3^(3/2))+1/8)^(1/3)))*z,
 698   x/20-(((-(sqrt(3)*%i)/2)-1/2)*((sqrt(37)*%i)/(8*3^(3/2))+1/8)^(1/3)
 699    +((sqrt(3)*%i)/2-1/2)/(3*((sqrt(37)*%i)/(8*3^(3/2))+1/8)^(1/3)))*y,
 700   50*z+20*y+((-((-(sqrt(3)*%i)/2)-1/2)*((sqrt(37)*%i)/(8*3^(3/2))+1/8)^(1/3))
 701    -((sqrt(3)*%i)/2-1/2)/(3*((sqrt(37)*%i)/(8*3^(3/2))+1/8)^(1/3)))*x],
 702   [x,y,z]),0);
 703 0;
 704 /* Solution of homogeneous eqns must have free parameter */
 705 freeof(%r1,s);
 706 false;
 707
 708 /* SF bug 380.  Missed solutions for algsys([a*b-c*d],[a,b,c,d]) */
 709 block([%rnum:0],algsys([a*b-c*d],[a,b,c,d]));
 710 [[a = (%r2*%r3)/%r1,b = %r1,c = %r2,d = %r3],
 711  [a = %r4,b = 0,c = 0,d = %r5],
 712  [a = %r6,b = 0,c = %r7,d = 0]];
 713
 714 block([%rnum:0],algsys([a*b*c-d*e*f],[a,b,c,d,e,f]));
 715 [[a = (%r3*%r4*%r5)/(%r1*%r2),b = %r1,c = %r2,d = %r3,e = %r4,f = %r5],
 716  [a = %r6,b = 0,c = %r7,d = 0,e = %r8,f = %r9],
 717  [a = %r10,b = %r11,c = 0,d = 0,e = %r12,f = %r13],
 718  [a = %r14,b = 0,c = %r15,d = %r16,e = 0,f = %r17],
 719  [a = %r18,b = %r19,c = 0,d = %r20,e = 0,f = %r21],
 720  [a = %r22,b = 0,c = %r23,d = %r24,e = %r25,f = 0],
 721  [a = %r26,b = %r27,c = 0,d = %r28,e = %r29,f = 0]];
 722
 723 block([%rnum:0],algsys([a^2*b^2-c*d],[a,b,c,d]));
 724 [[a = sqrt(%r2*%r3)/%r1,b = %r1,c = %r2,d = %r3],
 725  [a = -sqrt(%r5*%r6)/%r4,b = %r4,c = %r5,d = %r6],
 726  [a = %r7,b = 0,c = 0,d = %r8],
 727  [a = %r9,b = 0,c = %r10,d = 0]];
 728
 729 block([%rnum:0],algsys([a*b^2*c-d*e],[a,b,c,d,e]));
 730 [[a = (%r3*%r4)/(%r1^2*%r2),b = %r1,c = %r2,d = %r3,e = %r4],
 731  [a = %r5,b = 0,c = %r6,d = 0,e = %r7],
 732  [a = %r8,b = %r9,c = 0,d = 0,e = %r10],
 733  [a = %r11,b = 0,c = %r12,d = %r13,e = 0],
 734  [a = %r14,b = %r15,c = 0,d = %r16,e = 0]];
 735
 736 /* Reduced from rtest_to_poly_solve.mac
 737    Third solution found after bug #380 fixed */
 738 block([%rnum:0],
 739   eqs:[a5*b1=0,a1*b2=r3,b2+a1=r2,a2*b2+a3*b1=r4,r4+a3*b2+a5=d5],
 740   algsys(eqs,[a1,a2,a3,a5,b1,b2,r2,r3,r4]))$
 741 [[a1=%r1,a2=%r2,a3=%r3,a5=d5+((-%r3)-%r2)*%r4,b1=0,b2=%r4,
 742   r2=%r4+%r1,r3=%r1*%r4,r4=%r2*%r4],
 743  [a1=%r5,a2=(d5-%r6*%r8-%r6*%r7)/%r8,a3=%r6,a5=0,b1=%r7,
 744   b2=%r8,r2=%r8+%r5,r3=%r5*%r8,r4=d5-%r6*%r8],
 745  [a1=%r9,a2=%r10,a3=d5/%r11,a5=0,b1=%r11,b2=0,r2=%r9,
 746   r3=0,r4=d5]];
 747
 748
 749 (kill(eqt,eqs,F,G,p,p1,p2,s),done)$
 750 done;