/*************** -*- Mode: MACSYMA; Package: MAXIMA -*-  ******************/

(kill(all), done);
done;

/* Examples from Beyer (1984), an early paper on algsys, which references
   an internal General Motors report by Morgan (1981).  Morgan (1983) has
   the same title and contains some of the problems.

   Beyer, William A., Solution of Simultaneous Polynomial Equations by
   Elimination in MACSYMA, Proceedings of the 1984 MACSYMA Users' Conference,
   p 110-120

   Morgan, Alexander P., A Method for Computing All Solutions to Systems of
   Polynomials Equations, ACM Transactions on Mathematical Software,
   9:1-17 (1983)
 */

/* Beyer example 1, p117 */
algsys([x^2+x*y-1,y^2+x-5],[x,y]);
[[x = -3.174087266251113,y = 2.859036144578313],
 [x = 2.152189020152884,y = -1.687545787545788],
 [x = 0.3937104979189148,y = 2.146226811046756],
 [x = -0.3718122830843519,y = -2.317717206132879]];

/* Beyer example 2, p117,  Morgan problem 1 */
algsys([4*x^3-3*x-y,x^2-y],[x,y]);
[[x = 1,y = 1],[x = -3/4,y = 9/16],[x = 0,y = 0]];

/* Beyer example 3, p117, Morgan problem 2 */
algsys([4*(x+y),4*(x+y)+(x-y)*((x-2)^2+y^2-1)],[x,y]);
[[x = -((sqrt(2)*%i-2)/2), y =   (sqrt(2)*%i-2)/2],
 [x =   (sqrt(2)*%i+2)/2,  y = -((sqrt(2)*%i+2)/2)],
 [x = 0,y = 0]];

/* Beyer example 4, p117. Don't understand notation in paper */

/* Beyer example 5, p118.  2 real and 6 complex solutions */
algsys([
 x^2+2*y^2-4,
 x^2+y^2+z-8,
 (x-1)^2+(2*y-sqrt(2))^2+(z-5)^2-4],
 [x,y,z]);
[[x = 0, y = sqrt(2), z = 6],
 [x = 2.998824877038494 - 0.879551236563151*%i,
  y = (- 1.656382492469821*%i) - 0.796198988104073,
  z = 2.637620128835544*%i + 1.890329867297412],
 [x = 0.879551236563151*%i + 2.998824877038494,
  y = 1.656382492469821*%i - 0.7961989881040729,
  z = 1.890329867297412 - 2.637620128835543*%i],
 [x = (- 1.511881394988993*%i) - 2.77266624385368,
  y = 1.573386231515891*%i - 1.332140330399878,
  z = 3.299053626354102 - 4.191942508596404*%i],
 [x = 1.511881394988993*%i - 2.77266624385368,
  y = (- 1.573386231515895*%i) - 1.332140330399874,
  z = 4.191942508596393*%i + 3.299053626354096],
 [x = (- 1.06051780109088*%i) - 1.226158633184815,
  y = 1.421232537317398 - 0.4574772330741239*%i,
  z = 5.81061650634849 - 1.30036305745376*%i],
 [x = 1.06051780109088*%i - 1.226158633184815,
  y = 0.4574772330741241*%i + 1.421232537317397,
  z = 1.300363057453758*%i + 5.810616506348491],
 [x = 2, y = 0, z = 4]];

/* Beyer example 6, p118 */
block([%rnum:0],
  algsys([x+10*y,z+w,(z-2*z)^2,(z-w)^2],[w,x,y,z]));
[[w = 0, x = %r1 , y = -(%r1/10), z = 0]];

/* Beyer example 7, p118.  Supposed to be two real solutions.
   This solution seems correct.  Typo in paper? */
algsys([
 x+y+z+w-1,
 x+y-z+w-3,
 x^2+y^2+z^2+w^2-4,
 (x-1)^2+y^2+z^2+w^2],
 [w,x,y,z]);
[[w = -((3^(3/2)*%i+1)/4), x = 5/2, y =   (3^(3/2)*%i-1)/4,  z = -1],
 [w =   (3^(3/2)*%i-1)/4,  x = 5/2, y = -((3^(3/2)*%i+1)/4), z = -1]];

/* Check Beyer example 7. Preprocess with poly_reduced_grobner(). OK */
algsys([z+1,2*x-5,4*y^2+2*y+7,2*y+2*w+1],[w,x,y,z]);
[[w = -((3^(3/2)*%i+1)/4), x = 5/2, y =   (3^(3/2)*%i-1)/4,  z = -1],
 [w =   (3^(3/2)*%i-1)/4,  x = 5/2, y = -((3^(3/2)*%i+1)/4), z = -1]];

/* Beyer example 8, p118. Two real solutions */
algsys([
 x1^2+x2^2+x3^2+x4-6,
 x1+x2+x5+x6-1,
 x1+3*x3+x6-3,
 x1+x2+x4-2,
 2*x3+x6,
 x5+x6],
 [x1,x2,x3,x4,x5,x6]);
[[x1 = 1,  x2 = 0,   x3 = 2,  x4 = 1, x5 = 4,   x6 = -4],
 [x1 = 5/3,x2 = -2/3,x3 = 4/3,x4 = 1, x5 = 8/3, x6 = -8/3]];

/* Beyer example 9, p118. Morgan problem 4
   Three real and two complex solutions */
algsys(
 [x1+x1+x2+x3+x4+x5-6=0,
  x2+x1+x2+x3+x4+x5-6=0,
  x3+x1+x2+x3+x4+x5-6=0,
  x4+x1+x2+x3+x4+x5-6=0,
  x1*x2*x3*x4*x5-1=0],
 [x1,x2,x3,x4,x5]);
[[x1 = 0.916354556803995,
  x2 = 0.916354556803995,
  x3 = 0.916354556803995,
  x4 = 0.916354556803995,
  x5 = 1.418227215980025],
 [x1 = -0.5790430889759374,
  x2 = -0.5790430889759374,
  x3 = -0.5790430889759374,
  x4 = -0.5790430889759374,
  x5 = 8.895215311004785],
 [x1 = (-0.6100917318163317*%i)-0.06865574701986676,
  x2 = (-0.6100917318163317*%i)-0.06865574701986676,
  x3 = (-0.6100917318163317*%i)-0.06865574701986676,
  x4 = (-0.6100917318163317*%i)-0.06865574701986676,
  x5 = 3.050458659081659*%i+6.343278735099334],
 [x1 = 0.6100917318163317*%i-0.06865574701986676,
  x2 = 0.6100917318163318*%i-0.06865574701986678,
  x3 = 0.6100917318163317*%i-0.06865574701986676,
  x4 = 0.6100917318163317*%i-0.06865574701986676,
  x5 = 6.343278735099333-3.050458659081658*%i],
 [x1 = 1, x2 = 1, x3 = 1, x4 = 1, x5 = 1]];

/* Beyer example 10, p118. 4 real one-parameter solutions */
block([%rnum:0],
  algsys([x^2+y^2-1,x^2+y^2+z^2-5=0],[x,y,z]));
[[x = %r1, y = -sqrt(1-%r1^2), z =  2],
 [x = %r2, y =  sqrt(1-%r2^2), z =  2],
 [x = %r3, y = -sqrt(1-%r3^2), z = -2],
 [x = %r4, y =  sqrt(1-%r4^2), z = -2]];

/* Morgan problem 3.  Single solution (0,0,0,0) */
algsys([x+10*y,sqrt(5)*(z-w),(y-2*z)^2,sqrt(10)*(x-w)^2],[w,x,y,z]);
[[w = 0, x = 0, y = 0, z = 0]];

/* SF bug #3208 Can't solve this simple equation */
algsys([x^2+y^2=1,(x-2)^2+(y-3)^2=9],[x,y]);
[[x = -((3^(5/2)-10)/26),y = (2*3^(3/2)+15)/26],
 [x = (3^(5/2)+10)/26,y = -((2*3^(3/2)-15)/26)]];

/* Reduced from SF bug #3208 above */
algsys([676*y^2-20*3^(5/2)-333,676*y^2-4056*y+28*3^(7/2)+2007],[y]);
[[y = (2*3^(3/2)+15)/26]];

/* Reduced from SF bug #3208 above */
algsys([676*y^2+20*3^(5/2)-333,676*y^2-4056*y-28*3^(7/2)+2007],[y]);
[[y = -((2*3^(3/2)-15)/26)]];

/* from example(solve) */
algsys([4*x^2-y^2=12,x*y-x=2],[x,y]);
[[x = 2, y = 2],
 [x = 0.5202594388652008*%i - 0.1331240357358706,
  y = 0.07678378523787788 - 3.608003221870287*%i],
 [x = -0.5202594388652008*%i - 0.1331240357358706,
  y = 3.608003221870287*%i + 0.07678378523787788],
 [x = -1.733751846381093, y = -0.1535675710019696]];

/* failures from example(algsys) during testing */
block( [
     %rnum:0,
     f1:2*x*(1-l1)-2*(x-1)*l2,
     f2:l2-l1,
     f3:l1*(1-x**2-y),
     f4:l2*(y-(x-1)**2)
   ],
   algsys([f1,f2,f3,f4],[x,y,l1,l2]));
[[x = 0, y = %r1, l1 = 0, l2 = 0],
 [x = 1, y = 0,   l1 = 1, l2 = 1]];

block( [
   f1:x**2-y**2,
   f2:x**2-x+2*y**2-y-1
 ],
 algsys([f1,f2],[x,y]));
[[x = -(1/sqrt(3)), y =   1/sqrt(3)],
 [x =   1/sqrt(3),  y = -(1/sqrt(3))],
 [x = -(1/3),y = -(1/3)],
 [x = 1,y = 1]];

/* failure in example(charpoly) during testing */
 algsys([x2-2*x1=0,x2^2+x1^2=1],[x1,x2]);
 [[x1 = -(1/sqrt(5)), x2 = -(2/sqrt(5))],
  [x1 =   1/sqrt(5),  x2 =   2/sqrt(5)]];

block( [%rnum:0],
  algsys([g526+((sqrt(a^2+4)-a)/2+a)*g525,
          ((sqrt(a^2+4)-a)*g526)/2+g525],
	 [g525,g526]));
[[g525 = %r1,g526 = -((%r1*sqrt(a^2+4)+%r1*a)/2)]];

block( [%rnum:0],
  algsys([g625+(a-(sqrt(a^2+4)+a)/2)*g624,
           g624-((sqrt(a^2+4)+a)*g625)/2],
	   [g624,g625]));
[[g624 = %r1,g625 = (%r1*sqrt(a^2+4)-%r1*a)/2]];

/* Working example from bug #3150 */
algsys([x^2+y^2=1, x^2+(y-3)^2=9],[x,y]);
[[x = -sqrt(35)/6,y = 1/6],
 [x = sqrt(35)/6,y = 1/6]];

/* Working example from bug #2736 */
algsys([x^2+y^2=1,(x-2)^2+(y-3)^2=10],[x,y]);
[[x = 1,y = 0],[x = -5/13,y = 12/13]];

/* bug #1106 only found 4 roots (of 6)
   Now crash with heap exhausted
 */
/*
block( [
  p1:-x*y^3+y^2+x^4-9*x/8,
  p2:y^4-x^3*y-9*y/8+x^2],
  algsys([p1,p2],[x,y]));
[[x = 1/2,y = 1],
 [x = 9/8,y = 9/8],
 [x = 1,y = 1/2],
 [x = 0,y = 0]]
 */

/* Above can be solved with
   load(grobner)$
   eqs:poly_reduced_grobner([p1,p2],[x,y])$
   algsys(eqs,[x,y]);
 */
algsys([-200704*y^8+437832*y^5-263633*y^2+53010*x,
        4096*y^10-10440*y^7+7073*y^4-729*y],[x,y]);
[[x = 0, y = 0],
 [x = 1, y = 1/2],
 [x = 9/8, y = 9/8],
 [x = 1/2, y = 1],
 [x = (sqrt(3)*%i-1)/4,     y = -((sqrt(3)*%i+1)/2)],
 [x = -((sqrt(3)*%i+1)/4),  y = (sqrt(3)*%i-1)/2],
 [x = (3^(5/2)*%i-9)/16,    y = -((3^(5/2)*%i+9)/16)],
 [x = -((3^(5/2)*%i+9)/16), y = (3^(5/2)*%i-9)/16],
 [x = (sqrt(3)*%i-1)/2,     y = -((sqrt(3)*%i+1)/4)],
 [x = -((sqrt(3)*%i+1)/2),  y = (sqrt(3)*%i-1)/4]];

/* bug #1266 no solution found */
algsys([x^2+x+1=0,x^2*y+1=0],[x,y]);
[[x =   (sqrt(3)*%i-1)/2,  y = -((sqrt(3)*%i-1)/2)],
 [x = -((sqrt(3)*%i+1)/2), y =   (sqrt(3)*%i+1)/2]];

/* bug #1302 no solution found.  Also bug #1469 */
algsys([(x-x0)^2+(y-y0)^2=r^2,(x-x1)^2+(y-y1)^2=r^2],[x,y]);
[]; /* wrong, but at least it returns */

/* bug #1369 algsys hangs */
/*
eq1 : a*x + c*y + d*y^2/2 = 0;
eq2 : b*x + e*x^2 + f*y - g*y^3 = h;
algsys([eq1,eq2],[x,y]);
*/

/* bug #1370 Can solve [F=0,F-G=0] but not [F,G]
   there are some work arounds in the report
 */
F:2*x^3*y^2+2*a*x^2*y^2+3*x^2*y+2*a*x*y+x+1;
2*x^3*y^2+2*a*x^2*y^2+3*x^2*y+2*a*x*y+x+1;
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;
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;
/* FIXME:  Testsuite issue.  Can't match form of solution
algsys([F=0,F-G=0],[x,y]);
[[x = -(((-a^3)+sqrt(a^2-4*a)*(a^2-2*a)+4*a^2)/(2*sqrt(a^2-4*a))),
  y = -(1/(a*sqrt(a^2-4*a)))],
 [x = -((a^3+sqrt(a^2-4*a)*(a^2-2*a)-4*a^2)/(2*sqrt(a^2-4*a))),
  y = 1/(a*sqrt(a^2-4*a))]];
*/

/* hangs
algsys([F=0,G=0],[x,y]);
[];
 */

/* Can solve with 
   eqs:poly_reduced_grobner([F,G],[x,y])
   algsys:(eqs,[x,y]);
 */
/* FIXME:  Testsuite issue.  Can't match form of solution
algsys([(a^4-4*a^3)*y^2-1,(a^4-4*a^3)*y+2*x+a^2-2*a],[x,y]);
[[x = -((a*sqrt(a^2-4*a)+a^2-2*a)/2), y =   sqrt(a^2-4*a)/(a^3-4*a^2)],
 [x =   (a*sqrt(a^2-4*a)-a^2+2*a)/2,  y = -(sqrt(a^2-4*a)/(a^3-4*a^2))]];
*/
 
/* Bug #1647 */
algsys([x^2+y^2-d^2, (x-h)^2+(y-k)^2-s^2], [x,y]);
[];  /* no solution found, but at least it returns */

/*
  Here is a workaround
  load(grobner);
  [x^2+y^2-d^2, (x-h)^2+(y-k)^2-s^2];
  poly_reduced_grobner(%, [x,y]);
  algsys(%, [x,y]);
 */
/* FIXME:  Testsuite issue.  Can't match form of solution 
algsys([2*k*y+2*h*x+s^2-k^2-h^2-d^2,
        (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
          +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]);
[[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
                           +2*d^2*h^2-d^4)
          +h*s^2-(h*k^2)-h^3-(d^2*h))
          /(2*k^2+2*h^2))),
  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)
          -k*s^2+k^3+(h^2+d^2)*k)
          /(2*k^2+2*h^2)],
 [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
                          +2*d^2*h^2-d^4)
          -h*s^2+h*k^2+h^3+d^2*h)
          /(2*k^2+2*h^2),
  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)
          +k*s^2-k^3+((-h^2)-d^2)*k)
          /(2*k^2+2*h^2)))]];
 */
 
/* SF bug #3210 ALGSYS does not return (regression)
 *   also #3212 ALGSYS throws error
 *
 * Note: %nrum is the counter for the %r variables introduced
 *       by algsys.  Set it to 0 for reproducible results
 */
block([%rnum:0], algsys([
 4*y*(y^2-28),
 8*y*(y^2-28)*z,
 -16*y*(y^2-23)*z,
 8*y*(y^2-10)*z,
 32*y*(z-y),
 32*y*(z+y),
 4*y*(y^2-28),
 8*y*(y^2-28)*z,
 -16*y*(y^2-23)*z,
 8*y*(y^2-10)*z,
 y*((4*y^2-112)*z^2-8*y^4+(185-4*x^2)*y^2),
 -2*y*((4*y^2-92)*z^2-2*y^4+(21-4*x^2)*y^2),
 y*((4*y^2-40)*z^2+(1-4*x^2)*y^2)],
[x,y,z]));
[[x=%r1,y=0,z=%r2]];

block([%rnum:0], algsys([
  -8*a*z,
  -3*(y-1)*y*z,-3*z*(z+y-1)*(z+y),
  2*(y-1)*y*z*(6*y*z-3*z-4*b*y^2-2*y^2+4*b*y+2*y),
  -8*a*z,
  -3*(y-1)*y*z,
  -4*z*((w^2-2*v*w+v^2+2*b-2*a)*z+8*a*y-4*a),
  -2*(y-1)*y*z*((6*y-3)*z+((-4*b)-2)*y^2+(4*b+2)*y),
  4*z
   *((2*w^2-4*v*w+2*v^2+2*b-1)*z^2
    +(((-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)
     *z-12*a*y^2+12*a*y-2*a),
  -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
                                 -8*x+3)
     *z^2
     +(((-32*b)-16)*y^3+(48*b+24)*y^2+((-16*b)-8)*y)*z+(8*b-8*a+4)*y^4
     +((-16*b)+16*a-8)*y^3+(8*b-8*a+4)*y^2),
  -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
                                +(16*w+16*v-16)*x-32*v*w+16*b+8)
                                *z^2
                               +((4*w^2-8*v*w+4*v^2+48*b-48*a+20)*y^2
                                +(((-8*w)-8*v+8)*x+16*v*w-48*b+48*a-24)
                                 *y+4*x^2-8*x+8*b-8*a+4)
                                *z+32*a*y^3-48*a*y^2+16*a*y),
  -2*z
    *(((4*w^2-8*v*w+4*v^2+2)*y+((-4*w)-4*v+4)*x+8*v*w-3)*z^3
     +(((-4*w^2)+8*v*w-4*v^2-24*b-8)*y^2
      +((8*w+8*v-8)*x-16*v*w+24*b+12)*y-4*x^2+8*x-4*b-2)
      *z^2+((16*b-16*a+8)*y^3+((-24*b)+24*a-12)*y^2+(8*b-8*a+4)*y)*z
     +4*a*y^4-8*a*y^3+4*a*y^2)],
[v,w,x,y,z]));
[[v=%r1,w=%r2,x=%r3,y=%r4,z=0]];

/* Test problems distilled from share/contrib/diffequations tests

   The odelin method for linear second order ODEs calls algys.
   Changes to algsys capability change the odelin results.
   A number of tests below were obtained by tracing algsys
   while solving ODES with odelin (or contrib_ode).
 */

/* Solved by odelin for Kamke ODE 2.265.  New solution 2016-09 */
 algsys(
 [y = -1,
  z = 1,
  (-w^2)+2*v*w-v^2+5,2*((w+v-1)*x+4*w^2-10*v*w+4*v^2-22),
  (-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,
  (-4*x^2)+((-8*w)-8*v+8)*x+8*x-4*w^2+24*v*w-4*v^2+48,
  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),
  -2*((-6*x^2)-(12*w+12*v-12)*x+((-4*w)-4*v+4)*x+12*x-10*w^2+52*v*w
     -10*v^2+100),
  (-13*x^2)-(24*w+24*v-24)*x+((-26*w)-26*v+26)*x+26*x-41*w^2+182*v*w
     -41*v^2+344],
 [v,w,x,y,z]);
[[v =   (sqrt(5)+5)/2,  w = -((sqrt(5)-5)/2), x = 3,  y = -1, z = 1],
 [v = -((sqrt(5)-5)/2), w =   (sqrt(5)+5)/2,  x = 3,  y = -1, z = 1],
 [v =   (sqrt(5)-3)/2,  w = -((sqrt(5)+3)/2), x = -1, y = -1, z = 1],
 [v = -((sqrt(5)+3)/2), w =   (sqrt(5)-3)/2,  x = -1, y = -1, z = 1],
 [v =   (sqrt(5)+1)/2,  w = -((sqrt(5)-1)/2), x = 3,  y = -1, z = 1],
 [v = -((sqrt(5)-1)/2), w =   (sqrt(5)+1)/2,  x = 3,  y = -1, z = 1],
 [v =   (sqrt(5)+1)/2,  w = -((sqrt(5)-1)/2), x = -1, y = -1, z = 1],
 [v = -((sqrt(5)-1)/2), w =   (sqrt(5)+1)/2,  x = -1, y = -1, z = 1]];

/* Solved by odelin for Murphy ODE 1.182.  New solution 2016-10 */
/* FIXME:  Testsuite issue.  Can't match form of solution 
algsys(
 [x = -2*sqrt(-a1*a3),
  4*sqrt(-a1*a3)*((-4*a1*a3*(2*y-z))-4*a0*a3*sqrt(-a1*a3)),
  8*(-a1*a3)^(3/2)*(z^2-2*z-a2^2+2*a2)],
  [x,y,z]);
[[x = -2*sqrt(-a1*a3),
  y = -((a0*sqrt(-a1*a3)+a1*a2-2*a1)/(2*a1)),
  z = 2-a2],
 [x = -2*sqrt(-a1*a3),
  y = -((a0*sqrt(-a1*a3)-a1*a2)/(2*a1)),
  z = a2]];
 */

/* Solved by odelin for Kamke ODE 1.105.  New solution 2016-10 */
/* FIXME:  Testsuite issue.  Can't match form of solution
algsys(
 [4*sqrt(-a*c)*((-4*a*c*(2*y-z))-4*a*sqrt(-a*c)*d),
  8*(-a*c)^(3/2)*(z^2-2*z-b^2+2*b)],
 [y,z]);
[[y = -((sqrt(-a*c)*d+(b-2)*c)/(2*c)), z = 2-b],
 [y = -((sqrt(-a*c)*d-b*c)/(2*c)),     z = b]];
 */
 
/* Solved by odelin for Kamke ODE 1.291.  New solution 2016-10 */
algsys([
 -2*(18*y^2-36*x*y+18*x^2-5),
 (72*y+72*x-72)*z+72*y^2-288*x*y+72*x^2+69,
 (-36*z^2)-(144*y+144*x-144)*z+72*z-36*y^2+360*x*y-36*x^2-203,
 72*z^2-((-72*y)-72*x+72)*z-144*z-144*x*y+159,(-36*z^2)+72*z-35],
 [x,y,z]);
[[x =   (sqrt(10)-3)/12,   y = -((sqrt(10)+3)/12),  z = 5/6],
 [x = -((sqrt(10)+3)/12),  y =   (sqrt(10)-3)/12,   z = 5/6],
 [x =   (sqrt(10)+13)/12,  y = -((sqrt(10)-13)/12), z = 5/6],
 [x = -((sqrt(10)-13)/12), y =   (sqrt(10)+13)/12,  z = 5/6],
 [x =   (sqrt(10)+15)/12,  y = -((sqrt(10)-15)/12), z = 7/6],
 [x = -((sqrt(10)-15)/12), y =   (sqrt(10)+15)/12,  z = 7/6],
 [x =   (sqrt(10)-1)/12,   y = -((sqrt(10)+1)/12),  z = 7/6],
 [x = -((sqrt(10)+1)/12),  y =   (sqrt(10)-1)/12,   z = 7/6]];

/* SF bug #2038
   solve(eqs,[x,y]) gives correct solution
   solve(subst(sqrt(3),q,eqs),[x,y]) yields []!
*/
algsys(eqs:[y^2+x^2-(400-400/(2/q+1))^2=0,
       (-y-400/(2/q+1)+400)^2+x^2-640000/(2/q+1)^2=0],[x,y]);
[[x = -(400*q*sqrt(-((q-2)/(q+2)))), y = -((400*q^2-800)/(q+2))],
 [x =   400*q*sqrt(-((q-2)/(q+2))),  y = -((400*q^2-800)/(q+2))]];

algsys(subst(q=sqrt(3),eqs),[x,y]);
[[x = 1200-800*sqrt(3),y = 400*sqrt(3)-800],
 [x = 800*sqrt(3)-1200,y = 400*sqrt(3)-800]];

algsys(subst(q=1-sqrt(3),eqs),[x,y]);
[[x = -((25*2^(9/2))/3^(1/4)), y = 800/sqrt(3)],
 [x =   (25*2^(9/2))/3^(1/4),  y = 800/sqrt(3)]];

/* SF bug #2736 (in comments)
   Solving with p=1/(1-cos(1)) gave incorrect solution */
eqs:[x*y=p,x-y=3]$
[x*y = p,x-y = 3];

algsys(subst(p=1/(1-cos(1)),eqs),[x,y])$
[[x =   -((sqrt(9*cos(1)^2-22*cos(1)+13)-3*cos(1)+3)/(2*cos(1)-2)),
  y =   -((sqrt(9*cos(1)^2-22*cos(1)+13)+3*cos(1)-3)/(2*cos(1)-2))],
 [x =     (sqrt(9*cos(1)^2-22*cos(1)+13)+3*cos(1)-3)/(2*cos(1)-2),
  y =     (sqrt(9*cos(1)^2-22*cos(1)+13)-3*cos(1)+3)/(2*cos(1)-2)]];

float(%)$
[[x =  3.603649840080336, y =  0.603649840080336],
 [x = -0.603649840080336, y = -3.603649840080336]];

float(algsys(subst(p=float(1/(1-cos(1))),eqs),[x,y]))$
[[x = -0.6036497963666444,y = -3.603649796366644],
 [x =  3.603649796366644, y =  0.6036497963666444]];

/* SF bug 2059 - variable order dependent
   No solution found for certain solution variable orders
 */
(eqs:[z-a=%i*b,z+a=c,a*z+3*(a-z)=13-12*%i,3*(z-a)+a*z=12*%i+13],
algsys(eqs,[a,b,c,z]));
[[a = (-2*%i)-3, b = 4, c = -6, z = 2*%i-3],
 [a = 3-2*%i,    b = 4, c =  6, z = 2*%i+3]];

/* Count the number of solutions for each permutation of variable */
map(length,
  map(lambda([u],algsys(eqs,u)), listify(permutations([a,b,c,z]))));
[2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2];

/* SF bug #453 algsys fails in simple case */
algsys([v^2-2*v+u^2+10*u+1,v^2-10*v+u^2+2*u+1],[v,u]);
[[v = -((sqrt(34)-6)/2), u =   (sqrt(2)*sqrt(17)-6)/2],
 [v =   (sqrt(34)+6)/2,  u = -((sqrt(2)*sqrt(17)+6)/2)]];

algsys([2*v^2-12*v+1,v^2-2*v+u^2+10*u+1],[v,u]);

[[v = -((sqrt(2)*sqrt(17)-6)/2), u = -((sqrt(34)+14)/2)],
 [v =   (sqrt(2)*sqrt(17)+6)/2,  u =   (sqrt(34)-14)/2],
 [v =   (sqrt(2)*sqrt(17)+6)/2,  u = -((sqrt(34)+6)/2)],
 [v = -((sqrt(2)*sqrt(17)-6)/2), u =   (sqrt(34)-6)/2]];

/* Solved by odelin for Kamke ODE 2.184. */
algsys(
[w^2-a^2*w,2*w*x-2*a^2*x,
 w*(2*y+2*a*c+b^2)-3*a^2*y+x^2-2*a*b*x+A*a^2,
 w*(2*z+4*b*c)-4*a^2*z+2*x*y-4*a*b*y+2*A*a*b,
 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,
 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],
 [w,x,y,z]);
[[w = a^2,
  x = -((a*sqrt((-4*a*c)+b^2-4*A)-3*a*b)/2),
  y = -((b*sqrt((-4*a*c)+b^2-4*A)-2*a*c-b^2)/2),
  z = -((c*sqrt((-4*a*c)+b^2-4*A)-b*c)/2)],
 [w = a^2,
  x = (a*sqrt((-4*a*c)+b^2-4*A)+3*a*b)/2,
  y = (b*sqrt((-4*a*c)+b^2-4*A)+2*a*c+b^2)/2,
  z = (c*sqrt((-4*a*c)+b^2-4*A)+b*c)/2]];

/* Solved by odelin for Kamke ODE 2.119. */
 algsys([x^2-x-b^2-b,2*x*y-2*y,y^2],[x,y]);
[[x = -b, y = 0], [x = b+1, y = 0]];

/* Solved by odelin for Kamke ODE 2.331. */
algsys([x^2-16*x+48,2*x*y-32*y-192,y^2+64*y+64*x+192],[x,y]);
[[x = 12,y = -24],[x = 4,y = -8]];

/* Solved by odelin for Kamke ODE 2.382. */
algsys([
 w^2-w,
 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,
 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,
 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
   +3*a^2*b^2*w+((-b^2)-4*a*b-a^2)*c,
 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,
 z^2+(2*a*b^2+2*a^2*b)*z+a^2*b^2*y-a^2*b^2*c],
[w,x,y,z]);
[[w = 1,
 x = -((sqrt(4*c+b^2-(2*a*b)+a^2)+3*b+3*a)/2),
 y = ((b+a)*sqrt(4*c+b^2-(2*a)*b+a^2)+b^2+4*a*b+a^2)/2,
 z = -((a*b*sqrt(4*c+b^2-(2*a*b)+a^2)+a*b^2+a^2*b)/2)],
[w = 1,
 x = (sqrt(4*c+b^2-(2*a*b)+a^2)-3*b-3*a)/2,
 y = -(((b+a)*sqrt(4*c+b^2-2*a*b+a^2)-b^2-4*a*b-a^2)/2),
 z = (a*b*sqrt(4*c+b^2-2*a*b+a^2)-a*b^2-a^2*b)/2]];

/* Solved by odelin for Kamke ODE 2.383. */
algsys([
  w^2-4*w,
  2*w*x-8*x,
  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)
                            +(4-4*a^2)*b^2+b^2*((-4*b^2)+8*a*b-4*a^2)
                            +(8*a^3-8*a)*b-4*a^4+4*a^2,
  w*(2*z-16*a*b^2-16*a^2*b)-16*z+2*x*y+(16*b+16*a)*y
                           +b*((-16*a*b^3)+16*a^2*b^2+16*a^3*b-16*a^4)
                           +b^2*(8*b^3-8*a*b^2-8*a^2*b+8*a^3)+(8*a^2-8)*b^3
                           +(8*a-8*a^3)*b^2+(8*a^2-8*a^4)*b+8*a^5-8*a^3,
  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
                         +12*a^2*b^2*w
                         +b*(8*a*b^4+16*a^2*b^3-48*a^3*b^2+16*a^4*b+8*a^5)
                         +(4-4*a^2)*b^4
                         +b^2*((-4*b^4)-8*a*b^3+24*a^2*b^2-8*a^3*b-4*a^4)
                         +(8*a-8*a^3)*b^3+(24*a^4-24*a^2)*b^2+(8*a^3-8*a^5)*b
                         -4*a^6+4*a^4,
  2*y*z+((-8*b^2)-32*a*b-8*a^2)*z+8*a^2*b^2*x
       +b*((-16*a^2*b^4)+16*a^3*b^3+16*a^4*b^2-16*a^5*b)
       +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
       +(8*a^2-8*a^4)*b^3+(8*a^3-8*a^5)*b^2+(8*a^6-8*a^4)*b,
  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)
     +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
     +(8*a^5-8*a^3)*b^3+(4*a^4-4*a^6)*b^2],
 [w,x,y,z]);
[[w = 4,
  x = 2*b^2+((-4*a)-6)*b+2*a^2-6*a,
  y = (-2*b^3)+(2*a+2)*b^2+(2*a^2+8*a)*b-2*a^3+2*a^2,
  z = 2*a*b^3+((-4*a^2)-2*a)*b^2+(2*a^3-2*a^2)*b],
 [w = 4,
  x = (-2*b^2)+(4*a-6)*b-2*a^2-6*a,
  y = 2*b^3+(2-2*a)*b^2+(8*a-2*a^2)*b+2*a^3+2*a^2,
  z = (-2*a*b^3)+(4*a^2-2*a)*b^2+((-2*a^3)-2*a^2)*b]];

/* Eugene L. Allgower, Kurt Georg and Rick Miranda, The Method of
   Resultants for Computing Real Solutions of Polynomial Systems, 
   SIAM Journal on Numerical Analysis, Vol 29 No 3, Jun 1992, pp 831-844
 */

/* Example 1, p840, qe 12.  Robot arm coordinates.
   Real solutions have y = +/- 1.6237 
 */
(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,
 p2:(1-b-q1)*x^2*y^2+(-1-b-q1)*x^2-4*x*y+(-1+b-q1)*y^2+(1+b-q1),
 algsys(subst([b=1/2,q1=4/5,q2=2/5],[p1,p2]),[x,y]));
[[x = %i,y = -%i],
 [x = -%i,y = %i],
 [x = %i,y = -3*%i],
 [x = -%i,y = 3*%i],
 [x =   (sqrt(319)+8)/17,  y = -(sqrt(29)/sqrt(11))],
 [x = -((sqrt(319)-8)/17), y =   sqrt(29)/sqrt(11)]];
float([%[5],%[6]]);
[[x = 1.521210064675985, y = -1.623688281771977],
 [x = -0.58003359408775, y =  1.623688281771977]];

/* Example 2, p 841.  Intersection of a sphere with two paraboloids,
   Two real solutions (not found by maxima 5.38.0)
   y = z = (sqrt(5)-1)/2, x = +/- sqrt(z-y^2)
 */
algsys([x^2+y^2+z^2-1,z-x^2-y^2,y-x^2-z^2],[x,y,z]);
[[x =  sqrt(2)*%i,         y =   (sqrt(5)-1)/2,  z = -((sqrt(5)+1)/2)],
 [x =  sqrt(2)*%i,         y = -((sqrt(5)+1)/2), z =   (sqrt(5)-1)/2],
 [x = -sqrt(2)*%i,         y =   (sqrt(5)-1)/2,  z = -((sqrt(5)+1)/2)],
 [x = -sqrt(2)*%i,         y = -((sqrt(5)+1)/2), z =   (sqrt(5)-1)/2],
 [x =  sqrt((-sqrt(5))-2), y = -((sqrt(5)+1)/2), z = -((sqrt(5)+1)/2)],
 [x = -sqrt((-sqrt(5))-2), y = -((sqrt(5)+1)/2), z = -((sqrt(5)+1)/2)],
 [x =  sqrt(sqrt(5)-2),    y =   (sqrt(5)-1)/2,  z =   (sqrt(5)-1)/2],
 [x = -sqrt(sqrt(5)-2),    y =   (sqrt(5)-1)/2,  z =   (sqrt(5)-1)/2]];

/* Two examples from the maxima mailing list 2016-10-15 */

algsys(eqs:[x+2*y-z=6,2*x+y*z-z^2=-1,3*x-y+2*z^2=3],[x,y,z]);
[[x = (27*sqrt(13)*%i-41)/14,
  y = -((17*sqrt(13)*%i-87)/14),
  z = -((sqrt(13)*%i-7)/2)],
 [x = -((27*sqrt(13)*%i+41)/14),
  y = (17*sqrt(13)*%i+87)/14,
  z = (sqrt(13)*%i+7)/2],
 [x = 1, y = 2, z = -1]];

block([%rnum:0],
  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]));
[[a = -sqrt(1-%r1*%r2), b = %r1, c = %r2, d =  sqrt(1-%r1*%r2)],
 [a =  sqrt(1-%r3*%r4), b = %r3, c = %r4, d = -sqrt(1-%r3*%r4)],
 [a =  1, b = %r5, c = 0, d = -1],
 [a = -1, b = %r6, c = 0, d =  1],
 [a = -1, b = 0,   c = 0, d = -1],
 [a =  1, b = 0,   c = 0, d =  1]];

/* Example from Stavros Macrakis, maxima mailing list, 2017-07-31
   Not solved by 5.38.0, solved by 5.39.0 */
block([%rnum:0],
 algsys(eqs:[b*c+a^2,b*d+a*b,c*d+a*c,d^2+b*c],[a,b,c,d]));
[[a = -sqrt(-%r1*%r2), b = %r1, c = %r2, d =  sqrt(-%r1*%r2)],
 [a =  sqrt(-%r3*%r4), b = %r3, c = %r4, d = -sqrt(-%r3*%r4)],
 [a =  0,              b = %r5, c = 0,   d = 0]];


/* SF bug 3321 - apparent regression from 5.38.1

   In fact, correct solution depends on the signs of constants g and h.
   Version 5.38.1 returns wrong solution when g>0 and h<0.
   Would be better if algsys asked for signs, but ...
 */

eqs:[y-x=0, g*x*y-h=0, z+(x+1)/y-x-1=0];
[y-x=0, g*x*y-h=0, z+(x+1)/y-x-1=0];

(assume(h>0),0);
0;
s:algsys(eqs,[x,y,z]);
[[x = sqrt(h)/sqrt(g),
  y = sqrt(h)/sqrt(g),
  z = (h-g)/(sqrt(g)*sqrt(h))],
 [x = -(sqrt(h)/sqrt(g)),
  y = -(sqrt(h)/sqrt(g)),
  z = -((h-g)/(sqrt(g)*sqrt(h)))]];
ratsimp(subst(s[1],eqs));
[0=0, 0=0, 0=0];
ratsimp(subst(s[2],eqs));
[0=0, 0=0, 0=0];
(forget(h>0),0);
0;

(assume(h<0),0);
0;
s:algsys(eqs,[x,y,z]);
[[x = -((%i*sqrt(-h))/sqrt(g)),
  y = -((%i*sqrt(-h))/sqrt(g)),
  z = -((sqrt(-h)*(%i*h-%i*g))/(sqrt(g)*h))],
 [x = (%i*sqrt(-h))/sqrt(g),
  y = (%i*sqrt(-h))/sqrt(g),
  z = (sqrt(-h)*(%i*h-%i*g))/(sqrt(g)*h)]];
ratsimp(subst(s[1],eqs));
[0=0, 0=0, 0=0];
ratsimp(subst(s[2],eqs));
[0=0, 0=0, 0=0];
(forget(h<0),0);
0;

/* SF bug 3375.  Failure to calculate eigenvalues
   Solution is ugly, so just check it contains %r1 */
block([%rnum:0], s:algsys(eqs:[
  y/10-(((-(sqrt(3)*%i)/2)-1/2)*((sqrt(37)*%i)/(8*3^(3/2))+1/8)^(1/3)
   +((sqrt(3)*%i)/2-1/2)/(3*((sqrt(37)*%i)/(8*3^(3/2))+1/8)^(1/3)))*z,
  x/20-(((-(sqrt(3)*%i)/2)-1/2)*((sqrt(37)*%i)/(8*3^(3/2))+1/8)^(1/3)
   +((sqrt(3)*%i)/2-1/2)/(3*((sqrt(37)*%i)/(8*3^(3/2))+1/8)^(1/3)))*y,
  50*z+20*y+((-((-(sqrt(3)*%i)/2)-1/2)*((sqrt(37)*%i)/(8*3^(3/2))+1/8)^(1/3))
   -((sqrt(3)*%i)/2-1/2)/(3*((sqrt(37)*%i)/(8*3^(3/2))+1/8)^(1/3)))*x],
  [x,y,z]),0);
0;
/* Solution of homogeneous eqns must have free parameter */ 
freeof(%r1,s);
false;

/* SF bug 380.  Missed solutions for algsys([a*b-c*d],[a,b,c,d]) */
block([%rnum:0],algsys([a*b-c*d],[a,b,c,d]));
[[a = (%r2*%r3)/%r1,b = %r1,c = %r2,d = %r3],
 [a = %r4,b = 0,c = 0,d = %r5],
 [a = %r6,b = 0,c = %r7,d = 0]];

block([%rnum:0],algsys([a*b*c-d*e*f],[a,b,c,d,e,f]));
[[a = (%r3*%r4*%r5)/(%r1*%r2),b = %r1,c = %r2,d = %r3,e = %r4,f = %r5],
 [a = %r6,b = 0,c = %r7,d = 0,e = %r8,f = %r9],
 [a = %r10,b = %r11,c = 0,d = 0,e = %r12,f = %r13],
 [a = %r14,b = 0,c = %r15,d = %r16,e = 0,f = %r17],
 [a = %r18,b = %r19,c = 0,d = %r20,e = 0,f = %r21],
 [a = %r22,b = 0,c = %r23,d = %r24,e = %r25,f = 0],
 [a = %r26,b = %r27,c = 0,d = %r28,e = %r29,f = 0]];

block([%rnum:0],algsys([a^2*b^2-c*d],[a,b,c,d]));
[[a = sqrt(%r2*%r3)/%r1,b = %r1,c = %r2,d = %r3],
 [a = -sqrt(%r5*%r6)/%r4,b = %r4,c = %r5,d = %r6],
 [a = %r7,b = 0,c = 0,d = %r8],
 [a = %r9,b = 0,c = %r10,d = 0]];

block([%rnum:0],algsys([a*b^2*c-d*e],[a,b,c,d,e]));
[[a = (%r3*%r4)/(%r1^2*%r2),b = %r1,c = %r2,d = %r3,e = %r4],
 [a = %r5,b = 0,c = %r6,d = 0,e = %r7],
 [a = %r8,b = %r9,c = 0,d = 0,e = %r10],
 [a = %r11,b = 0,c = %r12,d = %r13,e = 0],
 [a = %r14,b = %r15,c = 0,d = %r16,e = 0]];

/* Reduced from rtest_to_poly_solve.mac
   Third solution found after bug #380 fixed */
block([%rnum:0],
  eqs:[a5*b1=0,a1*b2=r3,b2+a1=r2,a2*b2+a3*b1=r4,r4+a3*b2+a5=d5],
  algsys(eqs,[a1,a2,a3,a5,b1,b2,r2,r3,r4]))$
[[a1=%r1,a2=%r2,a3=%r3,a5=d5+((-%r3)-%r2)*%r4,b1=0,b2=%r4,
  r2=%r4+%r1,r3=%r1*%r4,r4=%r2*%r4],
 [a1=%r5,a2=(d5-%r6*%r8-%r6*%r7)/%r8,a3=%r6,a5=0,b1=%r7,
  b2=%r8,r2=%r8+%r5,r3=%r5*%r8,r4=d5-%r6*%r8],
 [a1=%r9,a2=%r10,a3=d5/%r11,a5=0,b1=%r11,b2=0,r2=%r9,
  r3=0,r4=d5]];


(kill(eqt,eqs,F,G,p,p1,p2,s),done)$
done;