roots.cpp

   1 #include "stdafx.h"
   2 #include "defs.h"
   3
   4 extern void sort_stack(int);
   5 void roots(void);
   6 static void roots2(void);
   7 static void roots3(void);
   8 static void mini_solve(void);
   9
  10 #define POLY p1
  11 #define X p2
  12 #define A p3
  13 #define B p4
  14 #define C p5
  15 #define Y p6
  16
  17 void
  18 eval_roots(void)
  19 {
  20         // A == B -> A - B
  21
  22         p2 = cadr(p1);
  23
  24         if (car(p2) == symbol(SETQ) || car(p2) == symbol(TESTEQ)) {
  25                 push(cadr(p2));
  26                 eval();
  27                 push(caddr(p2));
  28                 eval();
  29                 subtract();
  30         } else {
  31                 push(p2);
  32                 eval();
  33                 p2 = pop();
  34                 if (car(p2) == symbol(SETQ) || car(p2) == symbol(TESTEQ)) {
  35                         push(cadr(p2));
  36                         eval();
  37                         push(caddr(p2));
  38                         eval();
  39                         subtract();
  40                 } else
  41                         push(p2);
  42         }
  43
  44         // 2nd arg, x
  45
  46         push(caddr(p1));
  47         eval();
  48         p2 = pop();
  49         if (p2 == symbol(NIL))
  50                 guess();
  51         else
  52                 push(p2);
  53
  54         p2 = pop();
  55         p1 = pop();
  56
  57         if (!ispoly(p1, p2))
  58                 stop("roots: 1st argument is not a polynomial");
  59
  60         push(p1);
  61         push(p2);
  62
  63         roots();
  64 }
  65
  66 void
  67 roots(void)
  68 {
  69         int h, i, n;
  70         h = tos - 2;
  71         roots2();
  72         n = tos - h;
  73         if (n == 0)
  74                 stop("roots: the polynomial is not factorable, try nroots");
  75         if (n == 1)
  76                 return;
  77         sort_stack(n);
  78         save();
  79         p1 = alloc_tensor(n);
  80         p1->u.tensor->ndim = 1;
  81         p1->u.tensor->dim[0] = n;
  82         for (i = 0; i < n; i++)
  83                 p1->u.tensor->elem[i] = stack[h + i];
  84         tos = h;
  85         push(p1);
  86         restore();
  87 }
  88
  89 static void
  90 roots2(void)
  91 {
  92         save();
  93
  94         p2 = pop();
  95         p1 = pop();
  96
  97         push(p1);
  98         push(p2);
  99         factorpoly();
 100
 101         p1 = pop();
 102
 103         if (car(p1) == symbol(MULTIPLY)) {
 104                 p1 = cdr(p1);
 105                 while (iscons(p1)) {
 106                         push(car(p1));
 107                         push(p2);
 108                         roots3();
 109                         p1 = cdr(p1);
 110                 }
 111         } else {
 112                 push(p1);
 113                 push(p2);
 114                 roots3();
 115         }
 116
 117         restore();
 118 }
 119
 120 static void
 121 roots3(void)
 122 {
 123         save();
 124         p2 = pop();
 125         p1 = pop();
 126         if (car(p1) == symbol(POWER) && ispoly(cadr(p1), p2) && isposint(caddr(p1))) {
 127                 push(cadr(p1));
 128                 push(p2);
 129                 mini_solve();
 130         } else if (ispoly(p1, p2)) {
 131                 push(p1);
 132                 push(p2);
 133                 mini_solve();
 134         }
 135         restore();
 136 }
 137
 138 //-----------------------------------------------------------------------------
 139 //
 140 //      Input:          stack[tos - 2]          polynomial
 141 //
 142 //                      stack[tos - 1]          dependent symbol
 143 //
 144 //      Output:         stack                   roots on stack
 145 //
 146 //                                              (input args are popped first)
 147 //
 148 //-----------------------------------------------------------------------------
 149
 150 static void
 151 mini_solve(void)
 152 {
 153         int n;
 154
 155         save();
 156
 157         X = pop();
 158         POLY = pop();
 159
 160         push(POLY);
 161         push(X);
 162
 163         n = coeff();
 164
 165         // AX + B, X = -B/A
 166
 167         if (n == 2) {
 168                 A = pop();
 169                 B = pop();
 170                 push(B);
 171                 push(A);
 172                 divide();
 173                 negate();
 174                 restore();
 175                 return;
 176         }
 177
 178         // AX^2 + BX + C, X = (-B +/- (B^2 - 4AC)^(1/2)) / (2A)
 179
 180         if (n == 3) {
 181                 A = pop();
 182                 B = pop();
 183                 C = pop();
 184                 push(B);
 185                 push(B);
 186                 multiply();
 187                 push_integer(4);
 188                 push(A);
 189                 multiply();
 190                 push(C);
 191                 multiply();
 192                 subtract();
 193                 push_rational(1, 2);
 194                 power();
 195                 Y = pop();
 196                 push(Y);                        // 1st root
 197                 push(B);
 198                 subtract();
 199                 push(A);
 200                 divide();
 201                 push_rational(1, 2);
 202                 multiply();
 203                 push(Y);                        // 2nd root
 204                 push(B);
 205                 add();
 206                 negate();
 207                 push(A);
 208                 divide();
 209                 push_rational(1, 2);
 210                 multiply();
 211                 restore();
 212                 return;
 213         }
 214
 215         tos -= n;
 216
 217         restore();
 218 }
 219
 220 #if SELFTEST
 221
 222 static char *s[] = {
 223
 224         "roots(x)",
 225         "0",
 226
 227         "roots(x^2)",
 228         "0",
 229
 230         "roots(x^3)",
 231         "0",
 232
 233         "roots(2 x)",
 234         "0",
 235
 236         "roots(2 x^2)",
 237         "0",
 238
 239         "roots(2 x^3)",
 240         "0",
 241
 242         "roots(6+11*x+6*x^2+x^3)",
 243         "(-3,-2,-1)",
 244
 245         "roots(a*x^2+b*x+c)",
 246         "(-b/(2*a)-(-4*a*c+b^2)^(1/2)/(2*a),-b/(2*a)+(-4*a*c+b^2)^(1/2)/(2*a))",
 247
 248         "roots(3+7*x+5*x^2+x^3)",
 249         "(-3,-1)",
 250
 251         "roots(x^3+x^2+x+1)",
 252         "(-1,-i,i)",
 253
 254         "roots(x^4+1)",
 255         "Stop: roots: the polynomial is not factorable, try nroots",
 256
 257         "roots(x^2==1)",
 258         "(-1,1)",
 259
 260         "roots(3 x + 12 == 24)",
 261         "4",
 262
 263         "y=roots(x^2+b*x+c/k)[1]",
 264         "",
 265
 266         "y^2+b*y+c/k",
 267         "0",
 268
 269         "y=roots(x^2+b*x+c/k)[2]",
 270         "",
 271
 272         "y^2+b*y+c/k",
 273         "0",
 274
 275         "y=roots(a*x^2+b*x+c/4)[1]",
 276         "",
 277
 278         "a*y^2+b*y+c/4",
 279         "0",
 280
 281         "y=roots(a*x^2+b*x+c/4)[2]",
 282         "",
 283
 284         "a*y^2+b*y+c/4",
 285         "0",
 286
 287         "y=quote(y)",
 288         "",
 289 };
 290
 291 void
 292 test_roots(void)
 293 {
 294         test(__FILE__, s, sizeof s / sizeof (char *));
 295 }
 296
 297 #endif