| 1 | r""" |
| 2 | Interface to Maxima |
| 3 | |
| 4 | Maxima is a free GPL'd general purpose computer algebra system |
| 5 | whose development started in 1968 at MIT. It contains symbolic |
| 6 | manipulation algorithms, as well as implementations of special |
| 7 | functions, including elliptic functions and generalized |
| 8 | hypergeometric functions. Moreover, Maxima has implementations of |
| 9 | many functions relating to the invariant theory of the symmetric |
| 10 | group `S_n`. (However, the commands for group invariants, |
| 11 | and the corresponding Maxima documentation, are in French.) For many |
| 12 | links to Maxima documentation see |
| 13 | http://maxima.sourceforge.net/docs.shtml/. |
| 14 | |
| 15 | AUTHORS: |
| 16 | |
| 17 | - William Stein (2005-12): Initial version |
| 18 | |
| 19 | - David Joyner: Improved documentation |
| 20 | |
| 21 | - William Stein (2006-01-08): Fixed bug in parsing |
| 22 | |
| 23 | - William Stein (2006-02-22): comparisons (following suggestion of |
| 24 | David Joyner) |
| 25 | |
| 26 | - William Stein (2006-02-24): *greatly* improved robustness by adding |
| 27 | sequence numbers to IO bracketing in _eval_line |
| 28 | |
| 29 | If the string "error" (case insensitive) occurs in the output of |
| 30 | anything from Maxima, a RuntimeError exception is raised. |
| 31 | |
| 32 | EXAMPLES: We evaluate a very simple expression in Maxima. |
| 33 | |
| 34 | :: |
| 35 | |
| 36 | sage: maxima('3 * 5') |
| 37 | 15 |
| 38 | |
| 39 | We factor `x^5 - y^5` in Maxima in several different ways. |
| 40 | The first way yields a Maxima object. |
| 41 | |
| 42 | :: |
| 43 | |
| 44 | sage: F = maxima.factor('x^5 - y^5') |
| 45 | sage: F |
| 46 | -(y-x)*(y^4+x*y^3+x^2*y^2+x^3*y+x^4) |
| 47 | sage: type(F) |
| 48 | <class 'sage.interfaces.maxima_abstract.MaximaElement'> |
| 49 | |
| 50 | Note that Maxima objects can also be displayed using "ASCII art"; |
| 51 | to see a normal linear representation of any Maxima object x. Just |
| 52 | use the print command: use ``str(x)``. |
| 53 | |
| 54 | :: |
| 55 | |
| 56 | sage: print F |
| 57 | 4 3 2 2 3 4 |
| 58 | - (y - x) (y + x y + x y + x y + x ) |
| 59 | |
| 60 | You can always use ``repr(x)`` to obtain the linear |
| 61 | representation of an object. This can be useful for moving maxima |
| 62 | data to other systems. |
| 63 | |
| 64 | :: |
| 65 | |
| 66 | sage: repr(F) |
| 67 | '-(y-x)*(y^4+x*y^3+x^2*y^2+x^3*y+x^4)' |
| 68 | sage: F.str() |
| 69 | '-(y-x)*(y^4+x*y^3+x^2*y^2+x^3*y+x^4)' |
| 70 | |
| 71 | The ``maxima.eval`` command evaluates an expression in |
| 72 | maxima and returns the result as a *string* not a maxima object. |
| 73 | |
| 74 | :: |
| 75 | |
| 76 | sage: print maxima.eval('factor(x^5 - y^5)') |
| 77 | -(y-x)*(y^4+x*y^3+x^2*y^2+x^3*y+x^4) |
| 78 | |
| 79 | We can create the polynomial `f` as a Maxima polynomial, |
| 80 | then call the factor method on it. Notice that the notation |
| 81 | ``f.factor()`` is consistent with how the rest of Sage |
| 82 | works. |
| 83 | |
| 84 | :: |
| 85 | |
| 86 | sage: f = maxima('x^5 - y^5') |
| 87 | sage: f^2 |
| 88 | (x^5-y^5)^2 |
| 89 | sage: f.factor() |
| 90 | -(y-x)*(y^4+x*y^3+x^2*y^2+x^3*y+x^4) |
| 91 | |
| 92 | Control-C interruption works well with the maxima interface, |
| 93 | because of the excellent implementation of maxima. For example, try |
| 94 | the following sum but with a much bigger range, and hit control-C. |
| 95 | |
| 96 | :: |
| 97 | |
| 98 | sage: maxima('sum(1/x^2, x, 1, 10)') |
| 99 | 1968329/1270080 |
| 100 | |
| 101 | Tutorial |
| 102 | -------- |
| 103 | |
| 104 | We follow the tutorial at |
| 105 | http://maxima.sourceforge.net/docs/intromax/. |
| 106 | |
| 107 | :: |
| 108 | |
| 109 | sage: maxima('1/100 + 1/101') |
| 110 | 201/10100 |
| 111 | |
| 112 | :: |
| 113 | |
| 114 | sage: a = maxima('(1 + sqrt(2))^5'); a |
| 115 | (sqrt(2)+1)^5 |
| 116 | sage: a.expand() |
| 117 | 3*2^(7/2)+5*sqrt(2)+41 |
| 118 | |
| 119 | :: |
| 120 | |
| 121 | sage: a = maxima('(1 + sqrt(2))^5') |
| 122 | sage: float(a) |
| 123 | 82.012193308819747 |
| 124 | sage: a.numer() |
| 125 | 82.01219330881975 |
| 126 | |
| 127 | :: |
| 128 | |
| 129 | sage: maxima.eval('fpprec : 100') |
| 130 | '100' |
| 131 | sage: a.bfloat() |
| 132 | 8.20121933088197564152489730020812442785204843859314941221237124017312418754011041266612384955016056b1 |
| 133 | |
| 134 | :: |
| 135 | |
| 136 | sage: maxima('100!') |
| 137 | 93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000 |
| 138 | |
| 139 | :: |
| 140 | |
| 141 | sage: f = maxima('(x + 3*y + x^2*y)^3') |
| 142 | sage: f.expand() |
| 143 | x^6*y^3+9*x^4*y^3+27*x^2*y^3+27*y^3+3*x^5*y^2+18*x^3*y^2+27*x*y^2+3*x^4*y+9*x^2*y+x^3 |
| 144 | sage: f.subst('x=5/z') |
| 145 | (5/z+25*y/z^2+3*y)^3 |
| 146 | sage: g = f.subst('x=5/z') |
| 147 | sage: h = g.ratsimp(); h |
| 148 | (27*y^3*z^6+135*y^2*z^5+(675*y^3+225*y)*z^4+(2250*y^2+125)*z^3+(5625*y^3+1875*y)*z^2+9375*y^2*z+15625*y^3)/z^6 |
| 149 | sage: h.factor() |
| 150 | (3*y*z^2+5*z+25*y)^3/z^6 |
| 151 | |
| 152 | :: |
| 153 | |
| 154 | sage: eqn = maxima(['a+b*c=1', 'b-a*c=0', 'a+b=5']) |
| 155 | sage: s = eqn.solve('[a,b,c]'); s |
| 156 | [[a=(25*sqrt(79)*%i+25)/(6*sqrt(79)*%i-34),b=(5*sqrt(79)*%i+5)/(sqrt(79)*%i+11),c=(sqrt(79)*%i+1)/10],[a=(25*sqrt(79)*%i-25)/(6*sqrt(79)*%i+34),b=(5*sqrt(79)*%i-5)/(sqrt(79)*%i-11),c=-(sqrt(79)*%i-1)/10]] |
| 157 | |
| 158 | Here is an example of solving an algebraic equation:: |
| 159 | |
| 160 | sage: maxima('x^2+y^2=1').solve('y') |
| 161 | [y=-sqrt(1-x^2),y=sqrt(1-x^2)] |
| 162 | sage: maxima('x^2 + y^2 = (x^2 - y^2)/sqrt(x^2 + y^2)').solve('y') |
| 163 | [y=-sqrt((-y^2-x^2)*sqrt(y^2+x^2)+x^2),y=sqrt((-y^2-x^2)*sqrt(y^2+x^2)+x^2)] |
| 164 | |
| 165 | You can even nicely typeset the solution in latex:: |
| 166 | |
| 167 | sage: latex(s) |
| 168 | \left[ \left[ a={{25\,\sqrt{79}\,i+25}\over{6\,\sqrt{79}\,i-34}} , b={{5\,\sqrt{79}\,i+5}\over{\sqrt{79}\,i+11}} , c={{\sqrt{79}\,i+1 }\over{10}} \right] , \left[ a={{25\,\sqrt{79}\,i-25}\over{6\, \sqrt{79}\,i+34}} , b={{5\,\sqrt{79}\,i-5}\over{\sqrt{79}\,i-11}} , c=-{{\sqrt{79}\,i-1}\over{10}} \right] \right] |
| 169 | |
| 170 | To have the above appear onscreen via ``xdvi``, type |
| 171 | ``view(s)``. (TODO: For OS X should create pdf output |
| 172 | and use preview instead?) |
| 173 | |
| 174 | :: |
| 175 | |
| 176 | sage: e = maxima('sin(u + v) * cos(u)^3'); e |
| 177 | cos(u)^3*sin(v+u) |
| 178 | sage: f = e.trigexpand(); f |
| 179 | cos(u)^3*(cos(u)*sin(v)+sin(u)*cos(v)) |
| 180 | sage: f.trigreduce() |
| 181 | (sin(v+4*u)+sin(v-2*u))/8+(3*sin(v+2*u)+3*sin(v))/8 |
| 182 | sage: w = maxima('3 + k*%i') |
| 183 | sage: f = w^2 + maxima('%e')^w |
| 184 | sage: f.realpart() |
| 185 | %e^3*cos(k)-k^2+9 |
| 186 | |
| 187 | :: |
| 188 | |
| 189 | sage: f = maxima('x^3 * %e^(k*x) * sin(w*x)'); f |
| 190 | x^3*%e^(k*x)*sin(w*x) |
| 191 | sage: f.diff('x') |
| 192 | k*x^3*%e^(k*x)*sin(w*x)+3*x^2*%e^(k*x)*sin(w*x)+w*x^3*%e^(k*x)*cos(w*x) |
| 193 | sage: f.integrate('x') |
| 194 | (((k*w^6+3*k^3*w^4+3*k^5*w^2+k^7)*x^3+(3*w^6+3*k^2*w^4-3*k^4*w^2-3*k^6)*x^2+(-18*k*w^4-12*k^3*w^2+6*k^5)*x-6*w^4+36*k^2*w^2-6*k^4)*%e^(k*x)*sin(w*x)+((-w^7-3*k^2*w^5-3*k^4*w^3-k^6*w)*x^3+(6*k*w^5+12*k^3*w^3+6*k^5*w)*x^2+(6*w^5-12*k^2*w^3-18*k^4*w)*x-24*k*w^3+24*k^3*w)*%e^(k*x)*cos(w*x))/(w^8+4*k^2*w^6+6*k^4*w^4+4*k^6*w^2+k^8) |
| 195 | |
| 196 | :: |
| 197 | |
| 198 | sage: f = maxima('1/x^2') |
| 199 | sage: f.integrate('x', 1, 'inf') |
| 200 | 1 |
| 201 | sage: g = maxima('f/sinh(k*x)^4') |
| 202 | sage: g.taylor('x', 0, 3) |
| 203 | f/(k^4*x^4)-2*f/(3*k^2*x^2)+11*f/45-62*k^2*f*x^2/945 |
| 204 | |
| 205 | :: |
| 206 | |
| 207 | sage: maxima.taylor('asin(x)','x',0, 10) |
| 208 | x+x^3/6+3*x^5/40+5*x^7/112+35*x^9/1152 |
| 209 | |
| 210 | Examples involving matrices |
| 211 | --------------------------- |
| 212 | |
| 213 | We illustrate computing with the matrix whose `i,j` entry |
| 214 | is `i/j`, for `i,j=1,\ldots,4`. |
| 215 | |
| 216 | :: |
| 217 | |
| 218 | sage: f = maxima.eval('f[i,j] := i/j') |
| 219 | sage: A = maxima('genmatrix(f,4,4)'); A |
| 220 | matrix([1,1/2,1/3,1/4],[2,1,2/3,1/2],[3,3/2,1,3/4],[4,2,4/3,1]) |
| 221 | sage: A.determinant() |
| 222 | 0 |
| 223 | sage: A.echelon() |
| 224 | matrix([1,1/2,1/3,1/4],[0,0,0,0],[0,0,0,0],[0,0,0,0]) |
| 225 | sage: A.eigenvalues() |
| 226 | [[0,4],[3,1]] |
| 227 | sage: A.eigenvectors() |
| 228 | [[[0,4],[3,1]],[[[1,0,0,-4],[0,1,0,-2],[0,0,1,-4/3]],[[1,2,3,4]]]] |
| 229 | |
| 230 | We can also compute the echelon form in Sage:: |
| 231 | |
| 232 | sage: B = matrix(QQ, A) |
| 233 | sage: B.echelon_form() |
| 234 | [ 1 1/2 1/3 1/4] |
| 235 | [ 0 0 0 0] |
| 236 | [ 0 0 0 0] |
| 237 | [ 0 0 0 0] |
| 238 | sage: B.charpoly('x').factor() |
| 239 | (x - 4) * x^3 |
| 240 | |
| 241 | Laplace Transforms |
| 242 | ------------------ |
| 243 | |
| 244 | We illustrate Laplace transforms:: |
| 245 | |
| 246 | sage: _ = maxima.eval("f(t) := t*sin(t)") |
| 247 | sage: maxima("laplace(f(t),t,s)") |
| 248 | 2*s/(s^2+1)^2 |
| 249 | |
| 250 | :: |
| 251 | |
| 252 | sage: maxima("laplace(delta(t-3),t,s)") #Dirac delta function |
| 253 | %e^-(3*s) |
| 254 | |
| 255 | :: |
| 256 | |
| 257 | sage: _ = maxima.eval("f(t) := exp(t)*sin(t)") |
| 258 | sage: maxima("laplace(f(t),t,s)") |
| 259 | 1/(s^2-2*s+2) |
| 260 | |
| 261 | :: |
| 262 | |
| 263 | sage: _ = maxima.eval("f(t) := t^5*exp(t)*sin(t)") |
| 264 | sage: maxima("laplace(f(t),t,s)") |
| 265 | 360*(2*s-2)/(s^2-2*s+2)^4-480*(2*s-2)^3/(s^2-2*s+2)^5+120*(2*s-2)^5/(s^2-2*s+2)^6 |
| 266 | sage: print maxima("laplace(f(t),t,s)") |
| 267 | 3 5 |
| 268 | 360 (2 s - 2) 480 (2 s - 2) 120 (2 s - 2) |
| 269 | --------------- - --------------- + --------------- |
| 270 | 2 4 2 5 2 6 |
| 271 | (s - 2 s + 2) (s - 2 s + 2) (s - 2 s + 2) |
| 272 | |
| 273 | :: |
| 274 | |
| 275 | sage: maxima("laplace(diff(x(t),t),t,s)") |
| 276 | s*'laplace(x(t),t,s)-x(0) |
| 277 | |
| 278 | :: |
| 279 | |
| 280 | sage: maxima("laplace(diff(x(t),t,2),t,s)") |
| 281 | -?%at('diff(x(t),t,1),t=0)+s^2*'laplace(x(t),t,s)-x(0)*s |
| 282 | |
| 283 | It is difficult to read some of these without the 2d |
| 284 | representation:: |
| 285 | |
| 286 | sage: print maxima("laplace(diff(x(t),t,2),t,s)") |
| 287 | ! |
| 288 | d ! 2 |
| 289 | - -- (x(t))! + s laplace(x(t), t, s) - x(0) s |
| 290 | dt ! |
| 291 | !t = 0 |
| 292 | |
| 293 | Even better, use |
| 294 | ``view(maxima("laplace(diff(x(t),t,2),t,s)"))`` to see |
| 295 | a typeset version. |
| 296 | |
| 297 | Continued Fractions |
| 298 | ------------------- |
| 299 | |
| 300 | A continued fraction `a + 1/(b + 1/(c + \cdots))` is |
| 301 | represented in maxima by the list `[a, b, c, \ldots]`. |
| 302 | |
| 303 | :: |
| 304 | |
| 305 | sage: maxima("cf((1 + sqrt(5))/2)") |
| 306 | [1,1,1,1,2] |
| 307 | sage: maxima("cf ((1 + sqrt(341))/2)") |
| 308 | [9,1,2,1,2,1,17,1,2,1,2,1,17,1,2,1,2,1,17,2] |
| 309 | |
| 310 | Special examples |
| 311 | ---------------- |
| 312 | |
| 313 | In this section we illustrate calculations that would be awkward to |
| 314 | do (as far as I know) in non-symbolic computer algebra systems like |
| 315 | MAGMA or GAP. |
| 316 | |
| 317 | We compute the gcd of `2x^{n+4} - x^{n+2}` and |
| 318 | `4x^{n+1} + 3x^n` for arbitrary `n`. |
| 319 | |
| 320 | :: |
| 321 | |
| 322 | sage: f = maxima('2*x^(n+4) - x^(n+2)') |
| 323 | sage: g = maxima('4*x^(n+1) + 3*x^n') |
| 324 | sage: f.gcd(g) |
| 325 | x^n |
| 326 | |
| 327 | You can plot 3d graphs (via gnuplot):: |
| 328 | |
| 329 | sage: maxima('plot3d(x^2-y^2, [x,-2,2], [y,-2,2], [grid,12,12])') # not tested |
| 330 | [displays a 3 dimensional graph] |
| 331 | |
| 332 | You can formally evaluate sums (note the ``nusum`` |
| 333 | command):: |
| 334 | |
| 335 | sage: S = maxima('nusum(exp(1+2*i/n),i,1,n)') |
| 336 | sage: print S |
| 337 | 2/n + 3 2/n + 1 |
| 338 | %e %e |
| 339 | ----------------------- - ----------------------- |
| 340 | 1/n 1/n 1/n 1/n |
| 341 | (%e - 1) (%e + 1) (%e - 1) (%e + 1) |
| 342 | |
| 343 | We formally compute the limit as `n\to\infty` of |
| 344 | `2S/n` as follows:: |
| 345 | |
| 346 | sage: T = S*maxima('2/n') |
| 347 | sage: T.tlimit('n','inf') |
| 348 | %e^3-%e |
| 349 | |
| 350 | Miscellaneous |
| 351 | ------------- |
| 352 | |
| 353 | Obtaining digits of `\pi`:: |
| 354 | |
| 355 | sage: maxima.eval('fpprec : 100') |
| 356 | '100' |
| 357 | sage: maxima(pi).bfloat() |
| 358 | 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068b0 |
| 359 | |
| 360 | Defining functions in maxima:: |
| 361 | |
| 362 | sage: maxima.eval('fun[a] := a^2') |
| 363 | 'fun[a]:=a^2' |
| 364 | sage: maxima('fun[10]') |
| 365 | 100 |
| 366 | |
| 367 | Interactivity |
| 368 | ------------- |
| 369 | |
| 370 | Unfortunately maxima doesn't seem to have a non-interactive mode, |
| 371 | which is needed for the Sage interface. If any Sage call leads to |
| 372 | maxima interactively answering questions, then the questions can't be |
| 373 | answered and the maxima session may hang. See the discussion at |
| 374 | http://www.ma.utexas.edu/pipermail/maxima/2005/011061.html for some |
| 375 | ideas about how to fix this problem. An example that illustrates this |
| 376 | problem is ``maxima.eval('integrate (exp(a*x), x, 0, inf)')``. |
| 377 | |
| 378 | Latex Output |
| 379 | ------------ |
| 380 | |
| 381 | To TeX a maxima object do this:: |
| 382 | |
| 383 | sage: latex(maxima('sin(u) + sinh(v^2)')) |
| 384 | \sinh v^2+\sin u |
| 385 | |
| 386 | Here's another example:: |
| 387 | |
| 388 | sage: g = maxima('exp(3*%i*x)/(6*%i) + exp(%i*x)/(2*%i) + c') |
| 389 | sage: latex(g) |
| 390 | -{{i\,e^{3\,i\,x}}\over{6}}-{{i\,e^{i\,x}}\over{2}}+c |
| 391 | |
| 392 | Long Input |
| 393 | ---------- |
| 394 | |
| 395 | The MAXIMA interface reads in even very long input (using files) in |
| 396 | a robust manner, as long as you are creating a new object. |
| 397 | |
| 398 | .. note:: |
| 399 | |
| 400 | Using ``maxima.eval`` for long input is much less robust, and is |
| 401 | not recommended. |
| 402 | |
| 403 | :: |
| 404 | |
| 405 | sage: t = '"%s"'%10^10000 # ten thousand character string. |
| 406 | sage: a = maxima(t) |
| 407 | |
| 408 | TESTS: This working tests that a subtle bug has been fixed:: |
| 409 | |
| 410 | sage: f = maxima.function('x','gamma(x)') |
| 411 | sage: g = f(1/7) |
| 412 | sage: g |
| 413 | gamma(1/7) |
| 414 | sage: del f |
| 415 | sage: maxima(sin(x)) |
| 416 | sin(x) |
| 417 | |
| 418 | This tests to make sure we handle the case where Maxima asks if an |
| 419 | expression is positive or zero. |
| 420 | |
| 421 | :: |
| 422 | |
| 423 | sage: var('Ax,Bx,By') |
| 424 | (Ax, Bx, By) |
| 425 | sage: t = -Ax*sin(sqrt(Ax^2)/2)/(sqrt(Ax^2)*sqrt(By^2 + Bx^2)) |
| 426 | sage: t.limit(Ax=0,dir='above') |
| 427 | 0 |
| 428 | |
| 429 | A long complicated input expression:: |
| 430 | |
| 431 | sage: maxima._eval_line('((((((((((0) + ((1) / ((n0) ^ (0)))) + ((1) / ((n1) ^ (1)))) + ((1) / ((n2) ^ (2)))) + ((1) / ((n3) ^ (3)))) + ((1) / ((n4) ^ (4)))) + ((1) / ((n5) ^ (5)))) + ((1) / ((n6) ^ (6)))) + ((1) / ((n7) ^ (7)))) + ((1) / ((n8) ^ (8)))) + ((1) / ((n9) ^ (9)));') |
| 432 | '1/n9^9+1/n8^8+1/n7^7+1/n6^6+1/n5^5+1/n4^4+1/n3^3+1/n2^2+1/n1+1' |
| 433 | """ |
| 434 | |
| 435 | #***************************************************************************** |
| 436 | # Copyright (C) 2005 William Stein <wstein@gmail.com> |
| 437 | # |
| 438 | # Distributed under the terms of the GNU General Public License (GPL) |
| 439 | # |
| 440 | # This code is distributed in the hope that it will be useful, |
| 441 | # but WITHOUT ANY WARRANTY; without even the implied warranty of |
| 442 | # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
| 443 | # General Public License for more details. |
| 444 | # |
| 445 | # The full text of the GPL is available at: |
| 446 | # |
| 447 | # http://www.gnu.org/licenses/ |
| 448 | #***************************************************************************** |
| 449 | |
| 450 | from __future__ import with_statement |
| 451 | |
| 452 | import os, re, sys, subprocess |
| 453 | import pexpect |
| 454 | cygwin = os.uname()[0][:6]=="CYGWIN" |
| 455 | |
| 456 | from expect import Expect, ExpectElement, FunctionElement, ExpectFunction, gc_disabled, AsciiArtString |
| 457 | from pexpect import EOF |
| 458 | |
| 459 | from random import randrange |
| 460 | |
| 461 | ##import sage.rings.all |
| 462 | import sage.rings.complex_number |
| 463 | |
| 464 | from sage.misc.misc import verbose, DOT_SAGE, SAGE_ROOT |
| 465 | |
| 466 | from sage.misc.multireplace import multiple_replace |
| 467 | |
| 468 | COMMANDS_CACHE = '%s/maxima_commandlist_cache.sobj'%DOT_SAGE |
| 469 | |
| 470 | import sage.server.support |
| 471 | |
| 472 | import maxima_abstract |
| 473 | from maxima_abstract import MaximaFunctionElement, MaximaExpectFunction, MaximaElement, MaximaFunction, maxima_console |
| 474 | |
| 475 | # The Maxima "apropos" command, e.g., apropos(det) gives a list |
| 476 | # of all identifiers that begin in a certain way. This could |
| 477 | # maybe be useful somehow... (?) Also maxima has a lot for getting |
| 478 | # documentation from the system -- this could also be useful. |
| 479 | |
| 480 | #################################################################### |
| 481 | ## We begin here by initializing maxima in library more |
| 482 | from sage.libs.ecl import * |
| 483 | ecl_eval("(setf *load-verbose* NIL)") |
| 484 | ecl_eval("(require 'maxima)") |
| 485 | ecl_eval("(in-package :maxima)") |
| 486 | ecl_eval("(setq $nolabels t))") |
| 487 | ecl_eval("(defvar *MAXIMA-LANG-SUBDIR* NIL)") |
| 488 | ecl_eval("(set-pathnames)") |
| 489 | ecl_eval("(defun add-lineinfo (x) x)") |
| 490 | ecl_eval("""(defun retrieve (msg flag &aux (print? nil))(error (concatenate 'string "Maxima asks:" (meval (list '($string) msg)))))""") |
| 491 | ecl_eval('(defparameter *dev-null* (make-two-way-stream (make-concatenated-stream) (make-broadcast-stream)))') |
| 492 | ecl_eval('(defun principal nil (error "Divergent Integral"))') |
| 493 | |
| 494 | maxima_eval=ecl_eval(""" |
| 495 | (defun maxima-eval( form ) |
| 496 | (let ((result (catch 'macsyma-quit (cons 'maxima_eval (meval form))))) |
| 497 | ;(princ (list "result=" result)) |
| 498 | ;(terpri) |
| 499 | ;(princ (list "$error=" $error)) |
| 500 | ;(terpri) |
| 501 | (cond |
| 502 | ((and (consp result) (eq (car result) 'maxima_eval)) (cdr result)) |
| 503 | ((eq result 'maxima-error) (error (cadr $error))) |
| 504 | (t (error (concatenate 'string "Maxima condition. result:" (princ-to-string result) "$error:" (princ-to-string $error)))) |
| 505 | ) |
| 506 | ) |
| 507 | ) |
| 508 | """) |
| 509 | |
| 510 | #ecl_eval('(defun ask-evod (x even-odd)(error "Maxima asks a question"))') |
| 511 | #ecl_eval('(defun ask-integerp (x)(error "Maxima asks a question"))') |
| 512 | #ecl_eval('(defun ask-declare (x property)(error "Maxima asks a question"))') |
| 513 | #ecl_eval('(defun ask-prop (object property fun-or-number)(error "Maxima asks a question"))') |
| 514 | #ecl_eval('(defun asksign01 (a)(error "Maxima asks a question"))') |
| 515 | #ecl_eval('(defun asksign (x)(error "Maxima asks a question"))') |
| 516 | #ecl_eval('(defun asksign1 ($askexp)(error "Maxima asks a question"))') |
| 517 | #ecl_eval('(defun ask-greateq (x y)(error "Maxima asks a question"))') |
| 518 | #ecl_eval('(defun askinver (a)(error "Maxima asks a question"))') |
| 519 | #ecl_eval('(defun npask (exp)(error "Maxima asks a question"))') |
| 520 | |
| 521 | |
| 522 | import sage.symbolic.expression |
| 523 | from sage.symbolic.ring import is_SymbolicVariable |
| 524 | from sage.symbolic.ring import SR |
| 525 | |
| 526 | maxima_lib_instances = 0 |
| 527 | |
| 528 | maxprint=EclObject("$STRING") |
| 529 | meval=EclObject("MEVAL") |
| 530 | |
| 531 | def max_to_string(s): |
| 532 | return meval(EclObject([[maxprint],s])).python()[1:-1] |
| 533 | |
| 534 | class Maxima(maxima_abstract.Maxima): |
| 535 | """ |
| 536 | Interface to the Maxima interpreter. |
| 537 | """ |
| 538 | def __init__(self, script_subdirectory=None, logfile=None, server=None, |
| 539 | init_code = None): |
| 540 | """ |
| 541 | Create an instance of the Maxima interpreter. |
| 542 | |
| 543 | TESTS:: |
| 544 | |
| 545 | sage: maxima == loads(dumps(maxima)) |
| 546 | True |
| 547 | |
| 548 | We make sure labels are turned off (see trac 6816):: |
| 549 | |
| 550 | sage: 'nolabels:true' in maxima._Expect__init_code |
| 551 | True |
| 552 | """ |
| 553 | global maxima_lib_instances |
| 554 | if maxima_lib_instances > 0: |
| 555 | raise RuntimeError, "Maxima interface in library mode can only be instantiated once" |
| 556 | maxima_lib_instances += 1 |
| 557 | |
| 558 | #the single one gets defined in the Expect interface |
| 559 | self._eval_using_file_cutoff = 10**9 |
| 560 | #the double underscore one seems to by private to the Maxima interface |
| 561 | #this must have started as a typo by someone |
| 562 | self.__eval_using_file_cutoff = 10**9 |
| 563 | self._available_vars = [] |
| 564 | self._Expect__seq = 0 |
| 565 | self._session_number = 1 |
| 566 | self._Expect__name = "maxima" |
| 567 | |
| 568 | if True: |
| 569 | # display2d -- no ascii art output |
| 570 | # keepfloat -- don't automatically convert floats to rationals |
| 571 | init_code = ['display2d:false', 'domain: complex', 'keepfloat: true', 'load(to_poly_solver)', 'load(simplify_sum)'] |
| 572 | init_code.append('nolabels:true') |
| 573 | ecl_eval("(setf original-standard-output *standard-output*)") |
| 574 | ecl_eval("(setf *standard-output* *dev-null*)") |
| 575 | for l in init_code: |
| 576 | ecl_eval("#$%s$"%l) |
| 577 | ecl_eval("(setf *standard-output* original-standard-output)") |
| 578 | |
| 579 | def _function_class(self): |
| 580 | """ |
| 581 | EXAMPLES:: |
| 582 | |
| 583 | sage: maxima._function_class() |
| 584 | <class 'sage.interfaces.maxima_abstract.MaximaExpectFunction'> |
| 585 | """ |
| 586 | return MaximaExpectFunction |
| 587 | |
| 588 | def _start(self): |
| 589 | """ |
| 590 | Starts the Maxima interpreter. |
| 591 | |
| 592 | EXAMPLES:: |
| 593 | |
| 594 | sage: m = Maxima() |
| 595 | sage: m.is_running() |
| 596 | False |
| 597 | sage: m._start() |
| 598 | sage: m.is_running() |
| 599 | True |
| 600 | """ |
| 601 | ecl_eval(r"(defun tex-derivative (x l r) (tex (if $derivabbrev (tex-dabbrev x) (tex-d x '\partial)) l r lop rop ))") |
| 602 | |
| 603 | def __reduce__(self): |
| 604 | """ |
| 605 | EXAMPLES:: |
| 606 | |
| 607 | sage: maxima.__reduce__() |
| 608 | (<function reduce_load_Maxima at 0x...>, ()) |
| 609 | """ |
| 610 | return reduce_load_Maxima, tuple([]) |
| 611 | |
| 612 | def _sendline(self, str): |
| 613 | self._sendstr(str) |
| 614 | os.write(self._expect.child_fd, os.linesep) |
| 615 | |
| 616 | def _expect_expr(self, expr=None, timeout=None): |
| 617 | """ |
| 618 | EXAMPLES: |
| 619 | sage: a,b=var('a,b') |
| 620 | sage: integrate(1/(x^3 *(a+b*x)^(1/3)),x) |
| 621 | Traceback (most recent call last): |
| 622 | ... |
| 623 | TypeError: Computation failed since Maxima requested additional constraints (try the command 'assume(a>0)' before integral or limit evaluation, for example): |
| 624 | Is a positive or negative? |
| 625 | sage: assume(a>0) |
| 626 | sage: integrate(1/(x^3 *(a+b*x)^(1/3)),x) |
| 627 | 2/9*sqrt(3)*b^2*arctan(1/3*(2*(b*x + a)^(1/3) + a^(1/3))*sqrt(3)/a^(1/3))/a^(7/3) + 2/9*b^2*log((b*x + a)^(1/3) - a^(1/3))/a^(7/3) - 1/9*b^2*log((b*x + a)^(2/3) + (b*x + a)^(1/3)*a^(1/3) + a^(2/3))/a^(7/3) + 1/6*(4*(b*x + a)^(5/3)*b^2 - 7*(b*x + a)^(2/3)*a*b^2)/((b*x + a)^2*a^2 - 2*(b*x + a)*a^3 + a^4) |
| 628 | sage: var('x, n') |
| 629 | (x, n) |
| 630 | sage: integral(x^n,x) |
| 631 | Traceback (most recent call last): |
| 632 | ... |
| 633 | TypeError: Computation failed since Maxima requested additional constraints (try the command 'assume(n+1>0)' before integral or limit evaluation, for example): |
| 634 | Is n+1 zero or nonzero? |
| 635 | sage: assume(n+1>0) |
| 636 | sage: integral(x^n,x) |
| 637 | x^(n + 1)/(n + 1) |
| 638 | sage: forget() |
| 639 | """ |
| 640 | if expr is None: |
| 641 | expr = self._prompt_wait |
| 642 | if self._expect is None: |
| 643 | self._start() |
| 644 | try: |
| 645 | if timeout: |
| 646 | i = self._expect.expect(expr,timeout=timeout) |
| 647 | else: |
| 648 | i = self._expect.expect(expr) |
| 649 | if i > 0: |
| 650 | v = self._expect.before |
| 651 | |
| 652 | #We check to see if there is a "serious" error in Maxima. |
| 653 | #Note that this depends on the order of self._prompt_wait |
| 654 | if expr is self._prompt_wait and i > len(self._ask): |
| 655 | self.quit() |
| 656 | raise ValueError, "%s\nComputation failed due to a bug in Maxima -- NOTE: Maxima had to be restarted."%v |
| 657 | |
| 658 | j = v.find('Is ') |
| 659 | v = v[j:] |
| 660 | k = v.find(' ',4) |
| 661 | msg = "Computation failed since Maxima requested additional constraints (try the command 'assume(" + v[4:k] +">0)' before integral or limit evaluation, for example):\n" + v + self._ask[i-1] |
| 662 | self._sendline(";") |
| 663 | self._expect_expr() |
| 664 | raise ValueError, msg |
| 665 | except KeyboardInterrupt, msg: |
| 666 | #print self._expect.before |
| 667 | i = 0 |
| 668 | while True: |
| 669 | try: |
| 670 | print "Control-C pressed. Interrupting Maxima. Please wait a few seconds..." |
| 671 | self._sendstr('quit;\n'+chr(3)) |
| 672 | self._sendstr('quit;\n'+chr(3)) |
| 673 | self.interrupt() |
| 674 | self.interrupt() |
| 675 | except KeyboardInterrupt: |
| 676 | i += 1 |
| 677 | if i > 10: |
| 678 | break |
| 679 | pass |
| 680 | else: |
| 681 | break |
| 682 | raise KeyboardInterrupt, msg |
| 683 | |
| 684 | def _eval_line(self, line, allow_use_file=False, |
| 685 | wait_for_prompt=True, reformat=True, error_check=True): |
| 686 | return max_to_string(maxima_eval("#$%s$"%line)) |
| 687 | |
| 688 | |
| 689 | def _synchronize(self): |
| 690 | """ |
| 691 | Synchronize pexpect interface. |
| 692 | |
| 693 | This put a random integer (plus one!) into the output stream, then |
| 694 | waits for it, thus resynchronizing the stream. If the random |
| 695 | integer doesn't appear within 1 second, maxima is sent interrupt |
| 696 | signals. |
| 697 | |
| 698 | This way, even if you somehow left maxima in a busy state |
| 699 | computing, calling _synchronize gets everything fixed. |
| 700 | |
| 701 | EXAMPLES: This makes Maxima start a calculation:: |
| 702 | |
| 703 | sage: maxima._sendstr('1/1'*500) |
| 704 | |
| 705 | When you type this command, this synchronize command is implicitly |
| 706 | called, and the above running calculation is interrupted:: |
| 707 | |
| 708 | sage: maxima('2+2') |
| 709 | 4 |
| 710 | """ |
| 711 | marker = '__SAGE_SYNCHRO_MARKER_' |
| 712 | if self._expect is None: return |
| 713 | r = randrange(2147483647) |
| 714 | s = marker + str(r+1) |
| 715 | |
| 716 | # The 0; *is* necessary... it comes up in certain rare cases |
| 717 | # that are revealed by extensive testing. Don't delete it. -- william stein |
| 718 | cmd = '''0;sconcat("%s",(%s+1));\n'''%(marker,r) |
| 719 | self._sendstr(cmd) |
| 720 | try: |
| 721 | self._expect_expr(timeout=0.5) |
| 722 | if not s in self._before(): |
| 723 | self._expect_expr(s,timeout=0.5) |
| 724 | self._expect_expr(timeout=0.5) |
| 725 | except pexpect.TIMEOUT, msg: |
| 726 | self._interrupt() |
| 727 | except pexpect.EOF: |
| 728 | self._crash_msg() |
| 729 | self.quit() |
| 730 | |
| 731 | ########################################### |
| 732 | # Direct access to underlying lisp interpreter. |
| 733 | ########################################### |
| 734 | def lisp(self, cmd): |
| 735 | """ |
| 736 | Send a lisp command to maxima. |
| 737 | |
| 738 | .. note:: |
| 739 | |
| 740 | The output of this command is very raw - not pretty. |
| 741 | |
| 742 | EXAMPLES:: |
| 743 | |
| 744 | sage: maxima.lisp("(+ 2 17)") # random formatted output |
| 745 | :lisp (+ 2 17) |
| 746 | 19 |
| 747 | ( |
| 748 | """ |
| 749 | self._eval_line(':lisp %s\n""'%cmd, allow_use_file=False, wait_for_prompt=False, reformat=False, error_check=False) |
| 750 | self._expect_expr('(%i)') |
| 751 | return self._before() |
| 752 | |
| 753 | ########################################### |
| 754 | # Interactive help |
| 755 | ########################################### |
| 756 | def _command_runner(self, command, s, redirect=True): |
| 757 | """ |
| 758 | Run ``command`` in a new Maxima session and return its |
| 759 | output as an ``AsciiArtString``. |
| 760 | |
| 761 | If redirect is set to False, then the output of the command is not |
| 762 | returned as a string. Instead, it behaves like os.system. This is |
| 763 | used for interactive things like Maxima's demos. See maxima.demo? |
| 764 | |
| 765 | EXAMPLES:: |
| 766 | |
| 767 | sage: maxima._command_runner('describe', 'gcd') |
| 768 | -- Function: gcd (<p_1>, <p_2>, <x_1>, ...) |
| 769 | ... |
| 770 | """ |
| 771 | cmd = 'maxima --very-quiet -r "%s(%s);" '%(command, s) |
| 772 | if sage.server.support.EMBEDDED_MODE: |
| 773 | cmd += '< /dev/null' |
| 774 | |
| 775 | if redirect: |
| 776 | p = subprocess.Popen(cmd, shell=True, stdin=subprocess.PIPE, |
| 777 | stdout=subprocess.PIPE, stderr=subprocess.PIPE) |
| 778 | res = p.stdout.read() |
| 779 | # we are now getting three lines of commented verbosity |
| 780 | # every time Maxima starts, so we need to get rid of them |
| 781 | for _ in range(3): |
| 782 | res = res[res.find('\n')+1:] |
| 783 | return AsciiArtString(res) |
| 784 | else: |
| 785 | subprocess.Popen(cmd, shell=True) |
| 786 | |
| 787 | def _object_class(self): |
| 788 | """ |
| 789 | Return the Python class of Maxima elements. |
| 790 | |
| 791 | EXAMPLES:: |
| 792 | |
| 793 | sage: maxima._object_class() |
| 794 | <class 'sage.interfaces.maxima_abstract.MaximaElement'> |
| 795 | """ |
| 796 | return MaximaElement |
| 797 | |
| 798 | def _function_element_class(self): |
| 799 | """ |
| 800 | EXAMPLES:: |
| 801 | |
| 802 | sage: maxima._function_element_class() |
| 803 | <class 'sage.interfaces.maxima_abstract.MaximaFunctionElement'> |
| 804 | """ |
| 805 | return MaximaFunctionElement |
| 806 | |
| 807 | def function(self, args, defn, rep=None, latex=None): |
| 808 | """ |
| 809 | Return the Maxima function with given arguments and definition. |
| 810 | |
| 811 | INPUT: |
| 812 | |
| 813 | |
| 814 | - ``args`` - a string with variable names separated by |
| 815 | commas |
| 816 | |
| 817 | - ``defn`` - a string (or Maxima expression) that |
| 818 | defines a function of the arguments in Maxima. |
| 819 | |
| 820 | - ``rep`` - an optional string; if given, this is how |
| 821 | the function will print. |
| 822 | |
| 823 | |
| 824 | EXAMPLES:: |
| 825 | |
| 826 | sage: f = maxima.function('x', 'sin(x)') |
| 827 | sage: f(3.2) |
| 828 | -.05837414342758009 |
| 829 | sage: f = maxima.function('x,y', 'sin(x)+cos(y)') |
| 830 | sage: f(2,3.5) |
| 831 | sin(2)-.9364566872907963 |
| 832 | sage: f |
| 833 | sin(x)+cos(y) |
| 834 | |
| 835 | :: |
| 836 | |
| 837 | sage: g = f.integrate('z') |
| 838 | sage: g |
| 839 | (cos(y)+sin(x))*z |
| 840 | sage: g(1,2,3) |
| 841 | 3*(cos(2)+sin(1)) |
| 842 | |
| 843 | The function definition can be a maxima object:: |
| 844 | |
| 845 | sage: an_expr = maxima('sin(x)*gamma(x)') |
| 846 | sage: t = maxima.function('x', an_expr) |
| 847 | sage: t |
| 848 | gamma(x)*sin(x) |
| 849 | sage: t(2) |
| 850 | sin(2) |
| 851 | sage: float(t(2)) |
| 852 | 0.90929742682568171 |
| 853 | sage: loads(t.dumps()) |
| 854 | gamma(x)*sin(x) |
| 855 | """ |
| 856 | name = self._next_var_name() |
| 857 | if isinstance(defn, MaximaElement): |
| 858 | defn = defn.str() |
| 859 | elif not isinstance(defn, str): |
| 860 | defn = str(defn) |
| 861 | if isinstance(args, MaximaElement): |
| 862 | args = args.str() |
| 863 | elif not isinstance(args, str): |
| 864 | args = str(args) |
| 865 | cmd = '%s(%s) := %s'%(name, args, defn) |
| 866 | maxima._eval_line(cmd) |
| 867 | if rep is None: |
| 868 | rep = defn |
| 869 | f = MaximaFunction(self, name, rep, args, latex) |
| 870 | return f |
| 871 | |
| 872 | def set(self, var, value): |
| 873 | """ |
| 874 | Set the variable var to the given value. |
| 875 | |
| 876 | INPUT: |
| 877 | |
| 878 | |
| 879 | - ``var`` - string |
| 880 | |
| 881 | - ``value`` - string |
| 882 | |
| 883 | |
| 884 | EXAMPLES:: |
| 885 | |
| 886 | sage: maxima.set('xxxxx', '2') |
| 887 | sage: maxima.get('xxxxx') |
| 888 | '2' |
| 889 | """ |
| 890 | if not isinstance(value, str): |
| 891 | raise TypeError |
| 892 | cmd = '%s : %s$'%(var, value.rstrip(';')) |
| 893 | if len(cmd) > self.__eval_using_file_cutoff: |
| 894 | self._batch(cmd, batchload=True) |
| 895 | else: |
| 896 | self._eval_line(cmd) |
| 897 | #self._sendline(cmd) |
| 898 | #self._expect_expr() |
| 899 | #out = self._before() |
| 900 | #self._error_check(cmd, out) |
| 901 | |
| 902 | def clear(self, var): |
| 903 | """ |
| 904 | Clear the variable named var. |
| 905 | |
| 906 | EXAMPLES:: |
| 907 | |
| 908 | sage: maxima.set('xxxxx', '2') |
| 909 | sage: maxima.get('xxxxx') |
| 910 | '2' |
| 911 | sage: maxima.clear('xxxxx') |
| 912 | sage: maxima.get('xxxxx') |
| 913 | 'xxxxx' |
| 914 | """ |
| 915 | try: |
| 916 | self._expect.send('kill(%s)$'%var) |
| 917 | except (TypeError, AttributeError): |
| 918 | pass |
| 919 | |
| 920 | def get(self, var): |
| 921 | """ |
| 922 | Get the string value of the variable var. |
| 923 | |
| 924 | EXAMPLES:: |
| 925 | |
| 926 | sage: maxima.set('xxxxx', '2') |
| 927 | sage: maxima.get('xxxxx') |
| 928 | '2' |
| 929 | """ |
| 930 | s = self._eval_line('%s;'%var) |
| 931 | return s |
| 932 | |
| 933 | def version(self): |
| 934 | """ |
| 935 | Return the version of Maxima that Sage includes. |
| 936 | |
| 937 | EXAMPLES:: |
| 938 | |
| 939 | sage: maxima.version() |
| 940 | '5.19.1' |
| 941 | """ |
| 942 | return maxima_version() |
| 943 | |
| 944 | |
| 945 | ## def display2d(self, flag=True): |
| 946 | ## """ |
| 947 | ## Set the flag that determines whether Maxima objects are |
| 948 | ## printed using their 2-d ASCII art representation. When the |
| 949 | ## maxima interface starts the default is that objects are not |
| 950 | ## represented in 2-d. |
| 951 | |
| 952 | ## INPUT: |
| 953 | ## flag -- bool (default: True) |
| 954 | |
| 955 | ## EXAMPLES |
| 956 | ## sage: maxima('1/2') |
| 957 | ## 1/2 |
| 958 | ## sage: maxima.display2d(True) |
| 959 | ## sage: maxima('1/2') |
| 960 | ## 1 |
| 961 | ## - |
| 962 | ## 2 |
| 963 | ## sage: maxima.display2d(False) |
| 964 | ## """ |
| 965 | ## self._display2d = bool(flag) |
| 966 | |
| 967 | def is_MaximaElement(x): |
| 968 | """ |
| 969 | Returns True if x is of type MaximaElement. |
| 970 | |
| 971 | EXAMPLES:: |
| 972 | |
| 973 | sage: from sage.interfaces.maxima import is_MaximaElement |
| 974 | sage: m = maxima(1) |
| 975 | sage: is_MaximaElement(m) |
| 976 | True |
| 977 | sage: is_MaximaElement(1) |
| 978 | False |
| 979 | """ |
| 980 | return isinstance(x, MaximaElement) |
| 981 | |
| 982 | # An instance |
| 983 | maxima = Maxima() |
| 984 | |
| 985 | def reduce_load_Maxima(): |
| 986 | """ |
| 987 | EXAMPLES:: |
| 988 | |
| 989 | sage: from sage.interfaces.maxima import reduce_load_Maxima |
| 990 | sage: reduce_load_Maxima() |
| 991 | Maxima |
| 992 | """ |
| 993 | return maxima |
| 994 | |
| 995 | def maxima_version(): |
| 996 | """ |
| 997 | EXAMPLES:: |
| 998 | |
| 999 | sage: from sage.interfaces.maxima import maxima_version |
| 1000 | sage: maxima_version() |
| 1001 | '5.19.1' |
| 1002 | """ |
| 1003 | return os.popen('maxima --version').read().split()[-1] |
| 1004 | |
| 1005 | def __doctest_cleanup(): |
| 1006 | import sage.interfaces.quit |
| 1007 | sage.interfaces.quit.expect_quitall() |
| 1008 | |
| 1009 | |