Coverage for local/lib/python2.7/site-packages/sage/calculus/calculus.py: 64%

(sqrt(5) + sqrt(3) + sqrt(2))^8 - 40*(sqrt(5) + sqrt(3) + sqrt(2))^6 + 352*(sqrt(5) + sqrt(3) + sqrt(2))^4 - 960*(sqrt(5) + sqrt(3) + sqrt(2))^2 + 576

sage: f(a).expand()

sage: a = sin(pi/7)

sage: f = a.minpoly(algorithm='numerical'); f

x^6 - 7/4*x^4 + 7/8*x^2 - 7/64

sage: f(a).horner(a).numerical_approx(100)

0.00000000000000000000000000000

The degree must be high enough (default tops out at 24).

sage: a = sqrt(3) + sqrt(2)

sage: a.minpoly(algorithm='numerical', bits=100, degree=3)

Traceback (most recent call last):

...

ValueError: Could not find minimal polynomial (100 bits, degree 3).

sage: a.minpoly(algorithm='numerical', bits=100, degree=10)

x^4 - 10*x^2 + 1

sage: cos(pi/33).minpoly(algorithm='algebraic')

x^10 + 1/2*x^9 - 5/2*x^8 - 5/4*x^7 + 17/8*x^6 + 17/16*x^5 - 43/64*x^4 - 43/128*x^3 + 3/64*x^2 + 3/128*x + 1/1024

sage: cos(pi/33).minpoly(algorithm='numerical')

x^10 + 1/2*x^9 - 5/2*x^8 - 5/4*x^7 + 17/8*x^6 + 17/16*x^5 - 43/64*x^4 - 43/128*x^3 + 3/64*x^2 + 3/128*x + 1/1024

Sometimes it fails, as it must given that some numbers aren't algebraic::

sage: sin(1).minpoly(algorithm='numerical')

Traceback (most recent call last):

...

ValueError: Could not find minimal polynomial (1000 bits, degree 24).

.. note::

Of course, failure to produce a minimal polynomial does not

necessarily indicate that this number is transcendental.

"""

if algorithm is None or algorithm.startswith('numeric'):

bits_list = [bits] if bits else [100,200,500,1000]

degree_list = [degree] if degree else [2,4,8,12,24]

for bits in bits_list:

a = ex.numerical_approx(bits)

check_bits = int(1.25 * bits + 80)

aa = ex.numerical_approx(check_bits)

for degree in degree_list:

f = QQ[var](algdep(a, degree)) # TODO: use the known_bits parameter?

# If indeed we have found a minimal polynomial,

# it should be accurate to a much higher precision.

error = abs(f(aa))

dx = ~RR(Integer(1) << (check_bits - degree - 2))

expected_error = abs(f.derivative()(CC(aa))) * dx

if error < expected_error:

# Degree might have been an over-estimate,

# factor because we want (irreducible) minpoly.

ff = f.factor()

for g, e in ff:

lead = g.leading_coefficient()

if lead != 1:

g = g / lead

expected_error = abs(g.derivative()(CC(aa))) * dx

error = abs(g(aa))

if error < expected_error:

# See if we can prove equality exactly

if g(ex).simplify_trig().canonicalize_radical() == 0:

return g

# Otherwise fall back to numerical guess

elif epsilon and error < epsilon:

return g

elif algorithm is not None:

raise NotImplementedError("Could not prove minimal polynomial %s (epsilon %s)" % (g, RR(error).str(no_sci=False)))

if algorithm is not None:

raise ValueError("Could not find minimal polynomial (%s bits, degree %s)." % (bits, degree))

if algorithm is None or algorithm == 'algebraic':

from sage.rings.all import QQbar

return QQ[var](QQbar(ex).minpoly())

raise ValueError("Unknown algorithm: %s" % algorithm)

###################################################################

# limits

###################################################################

def limit(ex, dir=None, taylor=False, algorithm='maxima', **argv):

r"""

Return the limit as the variable `v` approaches `a`

from the given direction.

expr.limit(x = a)

expr.limit(x = a, dir='+')

INPUT:

- ``dir`` - (default: None); dir may have the value

'plus' (or '+' or 'right' or 'above') for a limit from above,

'minus' (or '-' or 'left' or 'below') for a limit from below, or may be omitted

(implying a two-sided limit is to be computed).

- ``taylor`` - (default: False); if True, use Taylor

series, which allows more limits to be computed (but may also

crash in some obscure cases due to bugs in Maxima).

- ``**argv`` - 1 named parameter

.. note::

The output may also use 'und' (undefined), 'ind'

(indefinite but bounded), and 'infinity' (complex

infinity).

EXAMPLES::

sage: x = var('x')

sage: f = (1+1/x)^x

sage: f.limit(x = oo)

sage: f.limit(x = 5)

7776/3125

sage: f.limit(x = 1.2)

2.06961575467...

sage: f.limit(x = I, taylor=True)

(-I + 1)^I

sage: f(x=1.2)

2.0696157546720...

sage: f(x=I)

(-I + 1)^I

sage: CDF(f(x=I))

2.0628722350809046 + 0.7450070621797239*I

sage: CDF(f.limit(x = I))

2.0628722350809046 + 0.7450070621797239*I

Notice that Maxima may ask for more information::

sage: var('a')

sage: limit(x^a,x=0)

Traceback (most recent call last):

...

ValueError: Computation failed since Maxima requested additional

constraints; using the 'assume' command before evaluation

*may* help (example of legal syntax is 'assume(a>0)', see

`assume?` for more details)

Is a positive, negative or zero?

With this example, Maxima is looking for a LOT of information::

sage: assume(a>0)

sage: limit(x^a,x=0)

Traceback (most recent call last):

...

ValueError: Computation failed since Maxima requested additional

constraints; using the 'assume' command before evaluation *may* help

(example of legal syntax is 'assume(a>0)', see `assume?` for

more details)

Is a an integer?

sage: assume(a,'integer')

sage: limit(x^a,x=0)

Traceback (most recent call last):

...

ValueError: Computation failed since Maxima requested additional

constraints; using the 'assume' command before evaluation *may* help

(example of legal syntax is 'assume(a>0)', see `assume?` for

more details)

Is a an even number?

sage: assume(a,'even')

sage: limit(x^a,x=0)

sage: forget()

More examples::

sage: limit(x*log(x), x = 0, dir='+')

sage: lim((x+1)^(1/x), x = 0)

sage: lim(e^x/x, x = oo)

+Infinity

sage: lim(e^x/x, x = -oo)

sage: lim(-e^x/x, x = oo)

-Infinity

sage: lim((cos(x))/(x^2), x = 0)

+Infinity

sage: lim(sqrt(x^2+1) - x, x = oo)

sage: lim(x^2/(sec(x)-1), x=0)

sage: lim(cos(x)/(cos(x)-1), x=0)

-Infinity

sage: lim(x*sin(1/x), x=0)

sage: limit(e^(-1/x), x=0, dir='right')

sage: limit(e^(-1/x), x=0, dir='left')

+Infinity

sage: f = log(log(x))/log(x)

sage: forget(); assume(x<-2); lim(f, x=0, taylor=True)

sage: forget()

Here ind means "indefinite but bounded"::

sage: lim(sin(1/x), x = 0)

ind

TESTS::

sage: lim(x^2, x=2, dir='nugget')

Traceback (most recent call last):

...

ValueError: dir must be one of None, 'plus', '+', 'above', 'right',

'minus', '-', 'below', 'left'

We check that :trac:`3718` is fixed, so that

Maxima gives correct limits for the floor function::

sage: limit(floor(x), x=0, dir='-')

-1

sage: limit(floor(x), x=0, dir='+')

sage: limit(floor(x), x=0)

und

Maxima gives the right answer here, too, showing

that :trac:`4142` is fixed::

sage: f = sqrt(1-x^2)

sage: g = diff(f, x); g

-x/sqrt(-x^2 + 1)

sage: limit(g, x=1, dir='-')

-Infinity

sage: limit(1/x, x=0)

Infinity

sage: limit(1/x, x=0, dir='+')

+Infinity

sage: limit(1/x, x=0, dir='-')

-Infinity

Check that :trac:`8942` is fixed::

sage: f(x) = (cos(pi/4-x) - tan(x)) / (1 - sin(pi/4+x))

sage: limit(f(x), x = pi/4, dir='minus')

+Infinity

sage: limit(f(x), x = pi/4, dir='plus')

-Infinity

sage: limit(f(x), x = pi/4)

Infinity

Check that :trac:`12708` is fixed::

sage: limit(tanh(x),x=0)

Check that :trac:`15386` is fixed::

sage: n = var('n')

sage: assume(n>0)

sage: sequence = -(3*n^2 + 1)*(-1)^n/sqrt(n^5 + 8*n^3 + 8)

sage: limit(sequence, n=infinity)

Check if :trac:`23048` is fixed::

sage: (1/(x-3)).limit(x=3, dir='below')

-Infinity

"""

if not isinstance(ex, Expression):

ex = SR(ex)

if len(argv) != 1:

raise ValueError("call the limit function like this, e.g. limit(expr, x=2).")

else:

k, = argv.keys()

v = var(k)

a = argv[k]

if taylor and algorithm == 'maxima':

algorithm = 'maxima_taylor'

if dir not in [None, 'plus', '+', 'right', 'minus', '-', 'left',

'above', 'below']:

raise ValueError("dir must be one of None, 'plus', '+', 'above', 'right', 'minus', '-', 'below', 'left'")

if algorithm == 'maxima':

if dir is None:

l = maxima.sr_limit(ex, v, a)

elif dir in ['plus', '+', 'right', 'above']:

l = maxima.sr_limit(ex, v, a, 'plus')

elif dir in ['minus', '-', 'left', 'below']:

l = maxima.sr_limit(ex, v, a, 'minus')

elif algorithm == 'maxima_taylor':

if dir is None:

l = maxima.sr_tlimit(ex, v, a)

elif dir in ['plus', '+', 'right', 'above']:

l = maxima.sr_tlimit(ex, v, a, 'plus')

elif dir in ['minus', '-', 'left', 'below']:

l = maxima.sr_tlimit(ex, v, a, 'minus')

elif algorithm == 'sympy':

if dir is None:

import sympy

l = sympy.limit(ex._sympy_(), v._sympy_(), a._sympy_())

else:

raise NotImplementedError("sympy does not support one-sided limits")

#return l.sage()

return ex.parent()(l)

# lim is alias for limit

lim = limit

###################################################################

# Laplace transform

###################################################################

def laplace(ex, t, s, algorithm='maxima'):

r"""

Return the Laplace transform with respect to the variable `t` and

transform parameter `s`, if possible.

If this function cannot find a solution, a formal function is returned.

The function that is returned may be viewed as a function of `s`.

DEFINITION:

The Laplace transform of a function `f(t)`, defined for all real numbers

`t \geq 0`, is the function `F(s)` defined by

.. MATH::

F(s) = \int_{0}^{\infty} e^{-st} f(t) dt.

INPUT:

- ``ex`` - a symbolic expression

- ``t`` - independent variable

- ``s`` - transform parameter

- ``algorithm`` - (default: ``'maxima'``) one of

- ``'maxima'`` - use Maxima (the default)

- ``'sympy'`` - use SymPy

- ``'giac'`` - use Giac

NOTES:

- The ``'sympy'`` algorithm returns the tuple (`F`, `a`, ``cond``) where `F` is the Laplace

transform of `f(t)`, `Re(s)>a` is the half-plane of convergence, and cond

are auxiliary convergence conditions.

.. SEEALSO::

:func:`inverse_laplace`

EXAMPLES:

We compute a few Laplace transforms::

sage: var('x, s, z, t, t0')

(x, s, z, t, t0)

sage: sin(x).laplace(x, s)

1/(s^2 + 1)

sage: (z + exp(x)).laplace(x, s)

z/s + 1/(s - 1)

sage: log(t/t0).laplace(t, s)

-(euler_gamma + log(s) + log(t0))/s

We do a formal calculation::

sage: f = function('f')(x)

sage: g = f.diff(x); g

diff(f(x), x)

sage: g.laplace(x, s)

s*laplace(f(x), x, s) - f(0)

A BATTLE BETWEEN the X-women and the Y-men (by David

Joyner): Solve

.. MATH::

x' = -16y, x(0)=270, y' = -x + 1, y(0) = 90.

This models a fight between two sides, the "X-women" and the

"Y-men", where the X-women have 270 initially and the Y-men have

90, but the Y-men are better at fighting, because of the higher

factor of "-16" vs "-1", and also get an occasional reinforcement,

because of the "+1" term.

sage: var('t')

sage: t = var('t')

sage: x = function('x')(t)

sage: y = function('y')(t)

sage: de1 = x.diff(t) + 16*y

sage: de2 = y.diff(t) + x - 1

sage: de1.laplace(t, s)

s*laplace(x(t), t, s) + 16*laplace(y(t), t, s) - x(0)

sage: de2.laplace(t, s)

s*laplace(y(t), t, s) - 1/s + laplace(x(t), t, s) - y(0)

Next we form the augmented matrix of the above system::

sage: A = matrix([[s, 16, 270],[1, s, 90+1/s]])

sage: E = A.echelon_form()

sage: xt = E[0,2].inverse_laplace(s,t)

sage: yt = E[1,2].inverse_laplace(s,t)

sage: xt

-91/2*e^(4*t) + 629/2*e^(-4*t) + 1

sage: yt

91/8*e^(4*t) + 629/8*e^(-4*t)

sage: p1 = plot(xt,0,1/2,rgbcolor=(1,0,0))

sage: p2 = plot(yt,0,1/2,rgbcolor=(0,1,0))

sage: (p1+p2).save(os.path.join(SAGE_TMP, "de_plot.png"))

Another example::

sage: var('a,s,t')

(a, s, t)

sage: f = exp (2*t + a) * sin(t) * t; f

t*e^(a + 2*t)*sin(t)

sage: L = laplace(f, t, s); L

2*(s - 2)*e^a/(s^2 - 4*s + 5)^2

sage: inverse_laplace(L, s, t)

t*e^(a + 2*t)*sin(t)

Unable to compute solution with Maxima::

sage: laplace(heaviside(t-1), t, s)

laplace(heaviside(t - 1), t, s)

Heaviside step function can be handled with different interfaces.

Try with giac::

sage: laplace(heaviside(t-1), t, s, algorithm='giac')

e^(-s)/s

Try with SymPy::

sage: laplace(heaviside(t-1), t, s, algorithm='sympy')

(e^(-s)/s, 0, True)

TESTS:

Testing Giac::

sage: var('t, s')

(t, s)

sage: laplace(5*cos(3*t-2)*heaviside(t-2), t, s, algorithm='giac')

5*(s*cos(4)*e^(-2*s) - 3*e^(-2*s)*sin(4))/(s^2 + 9)

Check unevaluated expression from Giac (it is locale-dependent, see

:trac:`22833`)::

sage: var('n')

sage: laplace(t^n, t, s, algorithm='giac')

laplace(t^n, t, s)

Testing SymPy::

sage: laplace(t^n, t, s, algorithm='sympy')

(s^(-n)*gamma(n + 1)/s, 0, -re(n) < 1)

Testing Maxima::

sage: laplace(t^n, t, s, algorithm='maxima')

s^(-n - 1)*gamma(n + 1)

Testing expression that is not parsed from SymPy to Sage::

sage: laplace(cos(t^2), t, s, algorithm='sympy')

Traceback (most recent call last):

...

AttributeError: Unable to convert SymPy result (=sqrt(pi)*(sqrt(2)*sin(s**2/4)*fresnelc(sqrt(2)*s/(2*sqrt(pi))) -

sqrt(2)*cos(s**2/4)*fresnels(sqrt(2)*s/(2*sqrt(pi))) + cos(s**2/4 + pi/4))/2) into Sage

"""

if not isinstance(ex, (Expression, Function)):

ex = SR(ex)

if algorithm == 'maxima':

return ex.parent()(ex._maxima_().laplace(var(t), var(s)))

elif algorithm == 'sympy':

ex_sy, t, s = [expr._sympy_() for expr in (ex, t, s)]

from sympy import laplace_transform

from sage.interfaces.sympy import sympy_init

sympy_init()

result = laplace_transform(ex_sy, t, s)

if isinstance(result, tuple):

try:

(result, a, cond) = result

return result._sage_(), a, cond

except AttributeError:

raise AttributeError("Unable to convert SymPy result (={}) into"

" Sage".format(result))

elif 'LaplaceTransform' in format(result):

return dummy_laplace(ex, t, s)

else:

return result

elif algorithm == 'giac':

from sage.interfaces.giac import giac

try:

result = giac.laplace(ex, t, s)

except TypeError:

raise ValueError("Giac cannot make sense of: %s" % ex_gi)

if 'integrate' in format(result) or 'integration' in format(result):

return dummy_laplace(ex, t, s)

else:

return result.sage()

else:

raise ValueError("Unknown algorithm: %s" % algorithm)

def inverse_laplace(ex, s, t, algorithm='maxima'):

r"""

Return the inverse Laplace transform with respect to the variable `t` and

transform parameter `s`, if possible.

If this function cannot find a solution, a formal function is returned.

The function that is returned may be viewed as a function of `t`.

DEFINITION:

The inverse Laplace transform of a function `F(s)` is the function

`f(t)`, defined by

.. MATH::

F(s) = \frac{1}{2\pi i} \int_{\gamma-i\infty}^{\gamma + i\infty} e^{st} F(s) dt,

where `\gamma` is chosen so that the contour path of

integration is in the region of convergence of `F(s)`.

INPUT:

- ``ex`` - a symbolic expression

- ``s`` - transform parameter

- ``t`` - independent variable

- ``algorithm`` - (default: ``'maxima'``) one of

- ``'maxima'`` - use Maxima (the default)

- ``'sympy'`` - use SymPy

- ``'giac'`` - use Giac

.. SEEALSO::

:func:`laplace`

EXAMPLES::

sage: var('w, m')

(w, m)

sage: f = (1/(w^2+10)).inverse_laplace(w, m); f

1/10*sqrt(10)*sin(sqrt(10)*m)

sage: laplace(f, m, w)

1/(w^2 + 10)

sage: f(t) = t*cos(t)

sage: s = var('s')

sage: L = laplace(f, t, s); L

t |--> 2*s^2/(s^2 + 1)^2 - 1/(s^2 + 1)

sage: inverse_laplace(L, s, t)

t |--> t*cos(t)

sage: inverse_laplace(1/(s^3+1), s, t)

1/3*(sqrt(3)*sin(1/2*sqrt(3)*t) - cos(1/2*sqrt(3)*t))*e^(1/2*t) + 1/3*e^(-t)

No explicit inverse Laplace transform, so one is returned formally a

function ``ilt``::

sage: inverse_laplace(cos(s), s, t)

ilt(cos(s), s, t)

Transform an expression involving a time-shift, via SymPy::

sage: inverse_laplace(1/s^2*exp(-s), s, t, algorithm='sympy')

-(log(e^(-t)) + 1)*heaviside(t - 1)

The same instance with Giac::

sage: inverse_laplace(1/s^2*exp(-s), s, t, algorithm='giac')

(t - 1)*heaviside(t - 1)

Transform a rational expression::

sage: inverse_laplace((2*s^2*exp(-2*s) - exp(-s))/(s^3+1), s, t, algorithm='giac')

-1/3*(sqrt(3)*e^(1/2*t - 1/2)*sin(1/2*sqrt(3)*(t - 1)) - cos(1/2*sqrt(3)*(t - 1))*e^(1/2*t - 1/2) +

e^(-t + 1))*heaviside(t - 1) + 2/3*(2*cos(1/2*sqrt(3)*(t - 2))*e^(1/2*t - 1) + e^(-t + 2))*heaviside(t - 2)

Dirac delta function can also be handled::

sage: inverse_laplace(1, s, t, algorithm='giac')

dirac_delta(t)

TESTS:

Testing unevaluated expression from Maxima::

sage: var('t, s')

(t, s)

sage: inverse_laplace(exp(-s)/s, s, t)

ilt(e^(-s)/s, s, t)

Testing Giac::

sage: inverse_laplace(exp(-s)/s, s, t, algorithm='giac')

heaviside(t - 1)

Testing SymPy::

sage: inverse_laplace(exp(-s)/s, s, t, algorithm='sympy')

heaviside(t - 1)

Testing unevaluated expression from Giac::

sage: n = var('n')

sage: inverse_laplace(1/s^n, s, t, algorithm='giac')

ilt(1/(s^n), t, s)

Try with Maxima::

sage: inverse_laplace(1/s^n, s, t, algorithm='maxima')

ilt(1/(s^n), s, t)

Try with SymPy::

sage: inverse_laplace(1/s^n, s, t, algorithm='sympy')

t^(n - 1)*heaviside(t)/gamma(n)

Testing unevaluated expression from SymPy::

sage: inverse_laplace(cos(s), s, t, algorithm='sympy')

ilt(cos(s), t, s)

Testing the same with Giac::

sage: inverse_laplace(cos(s), s, t, algorithm='giac')

ilt(cos(s), t, s)

"""

if not isinstance(ex, Expression):

ex = SR(ex)

if algorithm == 'maxima':

return ex.parent()(ex._maxima_().ilt(var(s), var(t)))

elif algorithm == 'sympy':

ex_sy, s, t = [expr._sympy_() for expr in (ex, s, t)]

from sympy import inverse_laplace_transform

from sage.interfaces.sympy import sympy_init

sympy_init()

result = inverse_laplace_transform(ex_sy, s, t)

try:

return result._sage_()

except AttributeError:

if 'InverseLaplaceTransform' in format(result):

return dummy_inverse_laplace(ex, t, s)

else:

raise AttributeError("Unable to convert SymPy result (={}) into"

" Sage".format(result))

elif algorithm == 'giac':

from sage.interfaces.giac import giac

try:

result = giac.invlaplace(ex, s, t)

except TypeError:

raise ValueError("Giac cannot make sense of: %s" % ex)

if 'ilaplace' in format(result):

return dummy_inverse_laplace(ex, t, s)

else:

return result.sage()

else:

raise ValueError("Unknown algorithm: %s" % algorithm)

###################################################################

# symbolic evaluation "at" a point

###################################################################

def at(ex, *args, **kwds):

"""

Parses ``at`` formulations from other systems, such as Maxima.

Replaces evaluation 'at' a point with substitution method of

a symbolic expression.

EXAMPLES:

We do not import ``at`` at the top level, but we can use it

as a synonym for substitution if we import it::

sage: g = x^3-3

sage: from sage.calculus.calculus import at

sage: at(g, x=1)

-2

sage: g.subs(x=1)

-2

We find a formal Taylor expansion::

sage: h,x = var('h,x')

sage: u = function('u')

sage: u(x + h)

u(h + x)

sage: diff(u(x+h), x)

D[0](u)(h + x)

sage: taylor(u(x+h),h,0,4)

1/24*h^4*diff(u(x), x, x, x, x) + 1/6*h^3*diff(u(x), x, x, x) + 1/2*h^2*diff(u(x), x, x) + h*diff(u(x), x) + u(x)

We compute a Laplace transform::

sage: var('s,t')

(s, t)

sage: f=function('f')(t)

sage: f.diff(t,2)

diff(f(t), t, t)

sage: f.diff(t,2).laplace(t,s)

s^2*laplace(f(t), t, s) - s*f(0) - D[0](f)(0)

We can also accept a non-keyword list of expression substitutions,

like Maxima does (:trac:`12796`)::

sage: from sage.calculus.calculus import at

sage: f = function('f')

sage: at(f(x), [x == 1])

f(1)

TESTS:

Our one non-keyword argument must be a list::

sage: from sage.calculus.calculus import at

sage: f = function('f')

sage: at(f(x), x == 1)

Traceback (most recent call last):

...

TypeError: at can take at most one argument, which must be a list

We should convert our first argument to a symbolic expression::

sage: from sage.calculus.calculus import at

sage: at(int(1), x=1)

"""

if not isinstance(ex, (Expression, Function)):

ex = SR(ex)

kwds={ (k[10:] if k[:10] == "_SAGE_VAR_" else k):v for k,v in six.iteritems(kwds)}

if len(args) == 1 and isinstance(args[0],list):

for c in args[0]:

kwds[str(c.lhs())]=c.rhs()

else:

if len(args) !=0:

raise TypeError("at can take at most one argument, which must be a list")

return ex.subs(**kwds)

def dummy_diff(*args):

"""

This function is called when 'diff' appears in a Maxima string.

EXAMPLES::

sage: from sage.calculus.calculus import dummy_diff

sage: x,y = var('x,y')

sage: dummy_diff(sin(x*y), x, SR(2), y, SR(1))

-x*y^2*cos(x*y) - 2*y*sin(x*y)

Here the function is used implicitly::

sage: a = var('a')

sage: f = function('cr')(a)

sage: g = f.diff(a); g

diff(cr(a), a)

"""

f = args[0]

args = list(args[1:])

for i in range(1, len(args), 2):

args[i] = Integer(args[i])

return f.diff(*args)

def dummy_integrate(*args):

"""

This function is called to create formal wrappers of integrals that

Maxima can't compute:

EXAMPLES::

sage: from sage.calculus.calculus import dummy_integrate

sage: f = function('f')

sage: dummy_integrate(f(x), x)

integrate(f(x), x)

sage: a,b = var('a,b')

sage: dummy_integrate(f(x), x, a, b)

integrate(f(x), x, a, b)

"""

if len(args) == 4:

return definite_integral(*args, hold=True)

else:

return indefinite_integral(*args, hold=True)

def dummy_laplace(*args):

"""

This function is called to create formal wrappers of laplace transforms

that Maxima can't compute:

EXAMPLES::

sage: from sage.calculus.calculus import dummy_laplace

sage: s,t = var('s,t')

sage: f = function('f')

sage: dummy_laplace(f(t),t,s)

laplace(f(t), t, s)

"""

return _laplace(args[0], var(repr(args[1])), var(repr(args[2])))

def dummy_inverse_laplace(*args):

"""

This function is called to create formal wrappers of inverse laplace

transforms that Maxima can't compute:

EXAMPLES::

sage: from sage.calculus.calculus import dummy_inverse_laplace

sage: s,t = var('s,t')

sage: F = function('F')

sage: dummy_inverse_laplace(F(s),s,t)

ilt(F(s), s, t)

"""

return _inverse_laplace(args[0], var(repr(args[1])), var(repr(args[2])))

#######################################################

# Helper functions for printing latex expression

#######################################################

def _laplace_latex_(self, *args):

r"""

Return LaTeX expression for Laplace transform of a symbolic function.

EXAMPLES::

sage: from sage.calculus.calculus import _laplace_latex_

sage: var('s,t')

(s, t)

sage: f = function('f')(t)

sage: _laplace_latex_(0,f,t,s)

'\\mathcal{L}\\left(f\\left(t\\right), t, s\\right)'

sage: latex(laplace(f, t, s))

\mathcal{L}\left(f\left(t\right), t, s\right)

"""

return "\\mathcal{L}\\left(%s\\right)"%(', '.join([latex(x) for x in args]))

def _inverse_laplace_latex_(self, *args):

r"""

Return LaTeX expression for inverse Laplace transform

of a symbolic function.

EXAMPLES::

sage: from sage.calculus.calculus import _inverse_laplace_latex_

sage: var('s,t')

(s, t)

sage: F = function('F')(s)

sage: _inverse_laplace_latex_(0,F,s,t)

'\\mathcal{L}^{-1}\\left(F\\left(s\\right), s, t\\right)'

sage: latex(inverse_laplace(F,s,t))

\mathcal{L}^{-1}\left(F\left(s\right), s, t\right)

"""

return "\\mathcal{L}^{-1}\\left(%s\\right)"%(', '.join([latex(x) for x in args]))

# Return un-evaluated expression as instances of SFunction class

_laplace = function_factory('laplace', print_latex_func=_laplace_latex_)

_inverse_laplace = function_factory('ilt',

print_latex_func=_inverse_laplace_latex_)

######################################i################

#######################################################

# Conversion dict for special maxima objects

# c,k1,k2 are from ode2()

symtable = {'%pi':'pi', '%e': 'e', '%i':'I', '%gamma':'euler_gamma',\

'%c' : '_C', '%k1' : '_K1', '%k2' : '_K2',

'e':'_e', 'i':'_i', 'I':'_I'}

import re

import six

maxima_tick = re.compile("'[a-z|A-Z|0-9|_]*")

maxima_qp = re.compile("\?\%[a-z|A-Z|0-9|_]*") # e.g., ?%jacobi_cd

maxima_var = re.compile("[a-z|A-Z|0-9|_\%]*") # e.g., %jacobi_cd

sci_not = re.compile("(-?(?:0|[1-9]\d*))(\.\d+)?([eE][-+]\d+)")

polylog_ex = re.compile('li\[([^\[\]]*)\]\(')

maxima_polygamma = re.compile("psi\[([^\[\]]*)\]\(") # matches psi[n]( where n is a number

maxima_hyper = re.compile("\%f\[\d+,\d+\]") # matches %f[m,n]

def symbolic_expression_from_maxima_string(x, equals_sub=False, maxima=maxima):

"""

Given a string representation of a Maxima expression, parse it and

return the corresponding Sage symbolic expression.

INPUT:

- ``x`` - a string

- ``equals_sub`` - (default: False) if True, replace

'=' by '==' in self

- ``maxima`` - (default: the calculus package's

Maxima) the Maxima interpreter to use.

EXAMPLES::

sage: from sage.calculus.calculus import symbolic_expression_from_maxima_string as sefms

sage: sefms('x^%e + %e^%pi + %i + sin(0)')

x^e + e^pi + I

sage: f = function('f')(x)

sage: sefms('?%at(f(x),x=2)#1')

f(2) != 1

sage: a = sage.calculus.calculus.maxima("x#0"); a

x#0

sage: a.sage()

x != 0

TESTS:

:trac:`8459` fixed::

sage: maxima('3*li[2](u)+8*li[33](exp(u))').sage()

3*dilog(u) + 8*polylog(33, e^u)

Check if :trac:`8345` is fixed::

sage: assume(x,'complex')

sage: t = x.conjugate()

sage: latex(t)

\overline{x}

sage: latex(t._maxima_()._sage_())

\overline{x}

Check that we can understand maxima's not-equals (:trac:`8969`)::

sage: from sage.calculus.calculus import symbolic_expression_from_maxima_string as sefms

sage: sefms("x!=3") == (factorial(x) == 3)

True

sage: sefms("x # 3") == SR(x != 3)

True

sage: solve([x != 5], x)

#0: solve_rat_ineq(ineq=_SAGE_VAR_x # 5)

[[x - 5 != 0]]

sage: solve([2*x==3, x != 5], x)

[[x == (3/2), (-7/2) != 0]]

Make sure that we don't accidentally pick up variables in the maxima namespace (:trac:`8734`)::

sage: sage.calculus.calculus.maxima('my_new_var : 2')

sage: var('my_new_var').full_simplify()

my_new_var

ODE solution constants are treated differently (:trac:`16007`)::

sage: from sage.calculus.calculus import symbolic_expression_from_maxima_string as sefms

sage: sefms('%k1*x + %k2*y + %c')

_K1*x + _K2*y + _C

Check that some hypothetical variables don't end up as special constants (:trac:`6882`)::

sage: from sage.calculus.calculus import symbolic_expression_from_maxima_string as sefms

sage: sefms('%i')^2

-1

sage: ln(sefms('%e'))

sage: sefms('i')^2

_i^2

sage: sefms('I')^2

_I^2

sage: sefms('ln(e)')

ln(_e)

sage: sefms('%inf')

+Infinity

"""

global _syms

syms = symbol_table.get('maxima', {}).copy()

if not len(x):

raise RuntimeError("invalid symbolic expression -- ''")

maxima.set('_tmp_', x)

# This is inefficient since it so rarely is needed:

#r = maxima._eval_line('listofvars(_tmp_);')[1:-1]

s = maxima._eval_line('_tmp_;')

formal_functions = maxima_tick.findall(s)

if len(formal_functions):

for X in formal_functions:

try:

syms[X[1:]] = _syms[X[1:]]

except KeyError:

syms[X[1:]] = function_factory(X[1:])

# You might think there is a potential very subtle bug if 'foo

# is in a string literal -- but string literals should *never*

# ever be part of a symbolic expression.

s = s.replace("'","")

delayed_functions = maxima_qp.findall(s)

if len(delayed_functions):

for X in delayed_functions:

if X == '?%at': # we will replace Maxima's "at" with symbolic evaluation, not an SFunction

pass

else:

syms[X[2:]] = function_factory(X[2:])

s = s.replace("?%", "")

s = maxima_hyper.sub('hypergeometric', s)

# Look up every variable in the symtable keys and fill a replacement list.

cursor = 0

l = []

for m in maxima_var.finditer(s):

if m.group(0) in symtable:

l.append(s[cursor:m.start()])

l.append(symtable.get(m.group(0)))

cursor = m.end()

if cursor > 0:

l.append(s[cursor:])

s = "".join(l)

s = s.replace("%","")

s = s.replace("#","!=") # a lot of this code should be refactored somewhere...

#we apply the square-bracket replacing patterns repeatedly

#to ensure that nested brackets get handled (from inside to out)

while True:

olds = s

s = polylog_ex.sub('polylog(\\1,', s)

s = maxima_polygamma.sub('psi(\g<1>,', s) # this replaces psi[n](foo) with psi(n,foo), ensuring that derivatives of the digamma function are parsed properly below

if s == olds: break

if equals_sub:

s = s.replace('=','==')

# unfortunately, this will turn != into !==, which we correct

s = s.replace("!==", "!=")

#replace %union from to_poly_solve with a list

if s[0:5]=='union':

s = s[5:]

s = s[s.find("(")+1:s.rfind(")")]

s = "[" + s + "]" # turn it into a string that looks like a list

#replace %solve from to_poly_solve with the expressions

if s[0:5]=='solve':

s = s[5:]

s = s[s.find("(")+1:s.find("]")+1]

#replace all instances of Maxima's scientific notation

#with regular notation

search = sci_not.search(s)

while not search is None:

(start, end) = search.span()

r = create_RealNumber(s[start:end]).str(no_sci=2, truncate=True)

s = s.replace(s[start:end], r)

search = sci_not.search(s)

# have to do this here, otherwise maxima_tick catches it

syms['diff'] = dummy_diff

syms['integrate'] = dummy_integrate

syms['laplace'] = dummy_laplace

syms['ilt'] = dummy_inverse_laplace

syms['at'] = at

global is_simplified

try:

# use a global flag so all expressions obtained via

# evaluation of maxima code are assumed pre-simplified

is_simplified = True

_syms = symbol_table['functions'].copy()

try:

global _augmented_syms

_augmented_syms = syms

return SRM_parser.parse_sequence(s)

finally:

_augmented_syms = {}

except SyntaxError:

raise TypeError("unable to make sense of Maxima expression '%s' in Sage"%s)

finally:

is_simplified = False

# Comma format options for Maxima

def mapped_opts(v):

"""

Used internally when creating a string of options to pass to

Maxima.

INPUT:

- ``v`` - an object

OUTPUT: a string.

The main use of this is to turn Python bools into lower case

strings.

EXAMPLES::

sage: sage.calculus.calculus.mapped_opts(True)

'true'

sage: sage.calculus.calculus.mapped_opts(False)

'false'

sage: sage.calculus.calculus.mapped_opts('bar')

'bar'

"""

if isinstance(v, bool):

return str(v).lower()

return str(v)

def maxima_options(**kwds):

"""

Used internally to create a string of options to pass to Maxima.

EXAMPLES::

sage: sage.calculus.calculus.maxima_options(an_option=True, another=False, foo='bar')

'an_option=true,foo=bar,another=false'

"""

return ','.join(['%s=%s'%(key,mapped_opts(val)) for key, val in six.iteritems(kwds)])

# Parser for symbolic ring elements

# We keep two dictionaries syms_cur and syms_default to keep the current symbol

# table and the state of the table at startup respectively. These are used by

# the restore() function (see sage.misc.reset).

# The dictionary _syms is used as a lookup table for the system function

# registry by _find_func() below. It gets updated by

# symbolic_expression_from_string() before calling the parser.

_syms = syms_cur = symbol_table.get('functions', {})

syms_default = dict(syms_cur)

# This dictionary is used to pass a lookup table other than the system registry

# to the parser. A global variable is necessary since the parser calls the

# _find_var() and _find_func() functions below without extra arguments.

_augmented_syms = {}

def _find_var(name):

"""

Function to pass to Parser for constructing

variables from strings. For internal use.

EXAMPLES::

sage: y = var('y')

sage: sage.calculus.calculus._find_var('y')

sage: sage.calculus.calculus._find_var('I')

"""

try:

res = _augmented_syms[name]

except KeyError:

pass

else:

# _augmented_syms might contain entries pointing to functions if

# previous computations polluted the maxima workspace

if not isinstance(res, Function):

return res

try:

return SR.symbols[name]

except KeyError:

pass

# try to find the name in the global namespace

# needed for identifiers like 'e', etc.

import sage.all

try:

return SR(sage.all.__dict__[name])

except (KeyError, TypeError):

return var(name)

def _find_func(name, create_when_missing = True):

"""

Function to pass to Parser for constructing

functions from strings. For internal use.

EXAMPLES::

sage: sage.calculus.calculus._find_func('limit')

limit

sage: sage.calculus.calculus._find_func('zeta_zeros')

zeta_zeros

sage: f(x)=sin(x)

sage: sage.calculus.calculus._find_func('f')

sage: sage.calculus.calculus._find_func('g', create_when_missing=False)

sage: s = sage.calculus.calculus._find_func('sin')

sage: s(0)

"""

try:

func = _augmented_syms.get(name)

if func is None:

func = _syms[name]

if not isinstance(func, Expression):

return func

except KeyError:

pass

import sage.all

try:

func = SR(sage.all.__dict__[name])

if not isinstance(func, Expression):

return func

except (KeyError, TypeError):

if create_when_missing:

return function_factory(name)

else:

return None

SR_parser = Parser(make_int = lambda x: SR(Integer(x)),

make_float = lambda x: SR(RealDoubleElement(x)),

make_var = _find_var,

make_function = _find_func)

def symbolic_expression_from_string(s, syms=None, accept_sequence=False):

"""

Given a string, (attempt to) parse it and return the

corresponding Sage symbolic expression. Normally used

to return Maxima output to the user.

INPUT:

- ``s`` - a string

- ``syms`` - (default: None) dictionary of

strings to be regarded as symbols or functions

- ``accept_sequence`` - (default: False) controls whether

to allow a (possibly nested) set of lists and tuples

as input

EXAMPLES::

sage: y = var('y')

sage: sage.calculus.calculus.symbolic_expression_from_string('[sin(0)*x^2,3*spam+e^pi]',syms={'spam':y},accept_sequence=True)

[0, 3*y + e^pi]

"""

global _syms

_syms = symbol_table['functions'].copy()

parse_func = SR_parser.parse_sequence if accept_sequence else SR_parser.parse_expression

if syms is None:

return parse_func(s)

else:

try:

global _augmented_syms

_augmented_syms = syms

return parse_func(s)

finally:

_augmented_syms = {}

def _find_Mvar(name):

"""

Function to pass to Parser for constructing

variables from strings. For internal use.

EXAMPLES::

sage: y = var('y')

sage: sage.calculus.calculus._find_var('y')

sage: sage.calculus.calculus._find_var('I')

"""

if name[:10] == "_SAGE_VAR_":

return var(name[10:])

res = _augmented_syms.get(name)

if res is not None and not isinstance(res, Function):

return res

# try to find the name in the global namespace

# needed for identifiers like 'e', etc.

import sage.all

try:

return SR(sage.all.__dict__[name])

except (KeyError, TypeError):

return var(name)

SRM_parser = Parser(make_int = lambda x: SR(Integer(x)),

make_float = lambda x: SR(RealDoubleElement(x)),

make_var = _find_Mvar,

make_function = _find_func)

Coverage for local/lib/python2.7/site-packages/sage/calculus/calculus.py : 64%

416 statements 265 run 151 missing 0 excluded