gh-35717: `sage.calculus`: Modularization fixes, doctest cosmetics, `# needs`
    
<!-- Please provide a concise, informative and self-explanatory title.
-->
<!-- Don't put issue numbers in the title. Put it in the Description
below. -->
<!-- For example, instead of "Fixes #12345", use "Add a new method to
multiply two integers" -->

### 📚 Description
Not all of `sage.calculus` is for symbolic calculus. We change some
import statements to `lazy_import` so that some modules can be loaded
even when the symbolic subsystem is not present, or move them into
methods.

We also add `# needs` tags to doctests to enable doctesting of
modularized distributions.


<!-- Describe your changes here in detail. -->
<!-- Why is this change required? What problem does it solve? -->
<!-- If this PR resolves an open issue, please link to it here. For
example "Fixes #12345". -->
- Part of: #29705
- Cherry-picked from: #35095
<!-- If your change requires a documentation PR, please link it
appropriately. -->

### 📝 Checklist

<!-- Put an `x` in all the boxes that apply. It should be `[x]` not `[x
]`. -->

- [x] The title is concise, informative, and self-explanatory.
- [x] The description explains in detail what this PR is about.
- [x] I have linked a relevant issue or discussion.
- [ ] I have created tests covering the changes.
- [ ] I have updated the documentation accordingly.

### ⌛ Dependencies

<!-- List all open PRs that this PR logically depends on
- #12345: short description why this is a dependency
- #34567: ...
-->

<!-- If you're unsure about any of these, don't hesitate to ask. We're
here to help! -->
    
URL: #35717
Reported by: Matthias Köppe
Reviewer(s): Kwankyu Lee, Matthias Köppe, Michael Orlitzky

Author: Release Manager

src/sage/calculus/calculus.py

@@ -874,8 +874,8 @@ def desolve_system(des, vars, ics=None, ivar=None, algorithm="maxima"):
 
     ::
 
-        sage: P1 = plot([solx,soly], (0,1))
-        sage: P2 = parametric_plot((solx,soly), (0,1))
+        sage: P1 = plot([solx,soly], (0,1))                                             # needs sage.plot
+        sage: P2 = parametric_plot((solx,soly), (0,1))                                  # needs sage.plot
 
     Now type ``show(P1)``, ``show(P2)`` to view these plots.
 
@@ -892,7 +892,8 @@ def desolve_system(des, vars, ics=None, ivar=None, algorithm="maxima"):
         sage: desolve_system([de1, de2], [x1, x2], ics=[1,1], ivar=t)
         Traceback (most recent call last):
         ...
-        ValueError: Initial conditions aren't complete: number of vars is different from number of dependent variables. Got ics = [1, 1], vars = [x1(t), x2(t)]
+        ValueError: Initial conditions aren't complete: number of vars is different
+        from number of dependent variables. Got ics = [1, 1], vars = [x1(t), x2(t)]
 
     Check that :trac:`9825` is fixed::
 
@@ -1013,9 +1014,9 @@ def eulers_method(f, x0, y0, h, x1, algorithm="table"):
     ::
 
         sage: pts = eulers_method(5*x+y-5,0,1,1/2,1,algorithm="none")
-        sage: P1 = list_plot(pts)
-        sage: P2 = line(pts)
-        sage: (P1+P2).show()
+        sage: P1 = list_plot(pts)                                                       # needs sage.plot
+        sage: P2 = line(pts)                                                            # needs sage.plot
+        sage: (P1 + P2).show()                                                          # needs sage.plot
 
     AUTHORS:
 
@@ -1163,7 +1164,7 @@ def eulers_method_2x2_plot(f, g, t0, x0, y0, h, t1):
 
         sage: from sage.calculus.desolvers import eulers_method_2x2_plot
         sage: f = lambda z : z[2]; g = lambda z : -sin(z[1])
-        sage: P = eulers_method_2x2_plot(f,g, 0.0, 0.75, 0.0, 0.1, 1.0)
+        sage: P = eulers_method_2x2_plot(f,g, 0.0, 0.75, 0.0, 0.1, 1.0)                 # needs sage.plot
     """
     from sage.plot.line import line
 
@@ -1430,13 +1431,14 @@ def desolve_system_rk4(des, vars, ics=None, ivar=None, end_points=None, step=0.1
     Lotka Volterra system::
 
         sage: from sage.calculus.desolvers import desolve_system_rk4
-        sage: x,y,t=var('x y t')
-        sage: P=desolve_system_rk4([x*(1-y),-y*(1-x)],[x,y],ics=[0,0.5,2],ivar=t,end_points=20)
-        sage: Q=[ [i,j] for i,j,k in P]
-        sage: LP=list_plot(Q)
+        sage: x,y,t = var('x y t')
+        sage: P = desolve_system_rk4([x*(1-y),-y*(1-x)], [x,y], ics=[0,0.5,2],
+        ....:                        ivar=t, end_points=20)
+        sage: Q = [[i,j] for i,j,k in P]
+        sage: LP = list_plot(Q)                                                         # needs sage.plot
 
-        sage: Q=[ [j,k] for i,j,k in P]
-        sage: LP=list_plot(Q)
+        sage: Q = [[j,k] for i,j,k in P]
+        sage: LP = list_plot(Q)                                                         # needs sage.plot
 
     ALGORITHM:
 
@@ -1570,49 +1572,52 @@ def desolve_odeint(des, ics, times, dvars, ivar=None, compute_jac=False, args=()
     Lotka Volterra Equations::
 
         sage: from sage.calculus.desolvers import desolve_odeint
-        sage: x,y=var('x,y')
-        sage: f=[x*(1-y),-y*(1-x)]
-        sage: sol=desolve_odeint(f,[0.5,2],srange(0,10,0.1),[x,y])
-        sage: p=line(zip(sol[:,0],sol[:,1]))
-        sage: p.show()
+        sage: x,y = var('x,y')
+        sage: f = [x*(1-y), -y*(1-x)]
+        sage: sol = desolve_odeint(f, [0.5,2], srange(0,10,0.1), [x,y])                 # needs scipy
+        sage: p = line(zip(sol[:,0],sol[:,1]))                                          # needs scipy sage.plot
+        sage: p.show()                                                                  # needs scipy sage.plot
 
     Lorenz Equations::
 
-        sage: x,y,z=var('x,y,z')
+        sage: x,y,z = var('x,y,z')
         sage: # Next we define the parameters
-        sage: sigma=10
-        sage: rho=28
-        sage: beta=8/3
+        sage: sigma = 10
+        sage: rho = 28
+        sage: beta = 8/3
         sage: # The Lorenz equations
-        sage: lorenz=[sigma*(y-x),x*(rho-z)-y,x*y-beta*z]
+        sage: lorenz = [sigma*(y-x),x*(rho-z)-y,x*y-beta*z]
         sage: # Time and initial conditions
-        sage: times=srange(0,50.05,0.05)
-        sage: ics=[0,1,1]
-        sage: sol=desolve_odeint(lorenz,ics,times,[x,y,z],rtol=1e-13,atol=1e-14)
+        sage: times = srange(0,50.05,0.05)
+        sage: ics = [0,1,1]
+        sage: sol = desolve_odeint(lorenz, ics, times, [x,y,z],                         # needs scipy
+        ....:                      rtol=1e-13, atol=1e-14)
 
     One-dimensional stiff system::
 
-        sage: y= var('y')
-        sage: epsilon=0.01
-        sage: f=y^2*(1-y)
-        sage: ic=epsilon
-        sage: t=srange(0,2/epsilon,1)
-        sage: sol=desolve_odeint(f,ic,t,y,rtol=1e-9,atol=1e-10,compute_jac=True)
-        sage: p=points(zip(t,sol[:,0]))
-        sage: p.show()
+        sage: y = var('y')
+        sage: epsilon = 0.01
+        sage: f = y^2*(1-y)
+        sage: ic = epsilon
+        sage: t = srange(0,2/epsilon,1)
+        sage: sol = desolve_odeint(f, ic, t, y,                                         # needs scipy
+        ....:                      rtol=1e-9, atol=1e-10, compute_jac=True)
+        sage: p = points(zip(t,sol[:,0]))                                               # needs scipy sage.plot
+        sage: p.show()                                                                  # needs scipy sage.plot
 
     Another stiff system with some optional parameters with no
     default value::
 
-        sage: y1,y2,y3=var('y1,y2,y3')
-        sage: f1=77.27*(y2+y1*(1-8.375*1e-6*y1-y2))
-        sage: f2=1/77.27*(y3-(1+y1)*y2)
-        sage: f3=0.16*(y1-y3)
-        sage: f=[f1,f2,f3]
-        sage: ci=[0.2,0.4,0.7]
-        sage: t=srange(0,10,0.01)
-        sage: v=[y1,y2,y3]
-        sage: sol=desolve_odeint(f,ci,t,v,rtol=1e-3,atol=1e-4,h0=0.1,hmax=1,hmin=1e-4,mxstep=1000,mxords=17)
+        sage: y1,y2,y3 = var('y1,y2,y3')
+        sage: f1 = 77.27*(y2+y1*(1-8.375*1e-6*y1-y2))
+        sage: f2 = 1/77.27*(y3-(1+y1)*y2)
+        sage: f3 = 0.16*(y1-y3)
+        sage: f = [f1,f2,f3]
+        sage: ci = [0.2,0.4,0.7]
+        sage: t = srange(0,10,0.01)
+        sage: v = [y1,y2,y3]
+        sage: sol = desolve_odeint(f, ci, t, v, rtol=1e-3, atol=1e-4,                   # needs scipy
+        ....:                      h0=0.1, hmax=1, hmin=1e-4, mxstep=1000, mxords=17)
 
     AUTHOR:
 

src/sage/calculus/desolvers.py

@@ -1,3 +1,4 @@
+# sage.doctest: needs sage.symbolic
 """
 Functional notation support for common calculus methods
 

src/sage/calculus/functional.py

@@ -1,9 +1,12 @@
+# sage.doctest: needs sage.symbolic
 r"""
 Calculus functions
 """
-from sage.matrix.constructor import matrix
+from sage.misc.lazy_import import lazy_import
 from sage.structure.element import is_Matrix, is_Vector, Expression
 
+lazy_import('sage.matrix.constructor', 'matrix')
+
 from .functional import diff
 
 

src/sage/calculus/functions.py

@@ -1,3 +1,4 @@
+# sage.doctest: needs sage.symbolic
 """
 Numerical Integration
 

src/sage/calculus/integration.pyx

@@ -55,9 +55,9 @@ cdef class Spline:
 
     This example is in the GSL documentation::
 
-        sage: v = [(i + sin(i)/2, i+cos(i^2)) for i in range(10)]
+        sage: v = [(i + RDF(i).sin()/2, i + RDF(i^2).cos()) for i in range(10)]
         sage: s = spline(v)
-        sage: show(point(v) + plot(s,0,9, hue=.8))
+        sage: show(point(v) + plot(s,0,9, hue=.8))                                      # needs sage.plot
 
     We compute the area underneath the spline::
 
@@ -77,13 +77,13 @@ cdef class Spline:
     We compute the first and second-order derivatives at a few points::
 
         sage: s.derivative(5)
-        -0.16230085261803...
+        -0.1623008526180...
         sage: s.derivative(6)
-        0.20997986285714...
+        0.2099798628571...
         sage: s.derivative(5, order=2)
-        -3.08747074561380...
+        -3.0874707456138...
         sage: s.derivative(6, order=2)
-        2.61876848274853...
+        2.6187684827485...
 
     Only the first two derivatives are supported::
 
@@ -118,7 +118,7 @@ cdef class Spline:
 
         Replace `0`-th point, which changes the spline::
 
-            sage: S[0]=(0,1); S
+            sage: S[0] = (0,1); S
             [(0, 1), (2, 3), (4, 5)]
             sage: S(1.5)
             2.5

src/sage/calculus/interpolation.pyx

@@ -1,3 +1,4 @@
+# sage.doctest: needs numpy
 """
 Complex Interpolation
 
@@ -49,9 +50,10 @@ def polygon_spline(pts):
         sage: ps = polygon_spline(pts)
         sage: fx = lambda x: ps.value(x).real
         sage: fy = lambda x: ps.value(x).imag
-        sage: show(parametric_plot((fx, fy), (0, 2*pi)))
-        sage: m = Riemann_Map([lambda x: ps.value(real(x))], [lambda x: ps.derivative(real(x))],0)
-        sage: show(m.plot_colored() + m.plot_spiderweb())
+        sage: show(parametric_plot((fx, fy), (0, 2*pi)))                                # needs sage.plot
+        sage: m = Riemann_Map([lambda x: ps.value(real(x))],
+        ....:                 [lambda x: ps.derivative(real(x))], 0)
+        sage: show(m.plot_colored() + m.plot_spiderweb())                               # needs sage.plot
 
     Polygon approximation of an circle::
 
@@ -119,7 +121,7 @@ cdef class PSpline:
             sage: ps = polygon_spline(pts)
             sage: ps.value(.5)
             (-0.363380227632...-1j)
-            sage: ps.value(0) - ps.value(2*pi)
+            sage: ps.value(0) - ps.value(2*RDF.pi())
             0j
             sage: ps.value(10)
             (0.26760455264...+1j)
@@ -152,6 +154,7 @@ cdef class PSpline:
             sage: ps = polygon_spline(pts)
             sage: ps.derivative(1 / 3)
             (1.27323954473...+0j)
+            sage: from math import pi
             sage: ps.derivative(0) - ps.derivative(2*pi)
             0j
             sage: ps.derivative(10)
@@ -171,7 +174,7 @@ def complex_cubic_spline(pts):
 
     INPUT:
 
-    - ``pts`` A list or array of complex numbers, or tuples of the form
+    - ``pts`` -- A list or array of complex numbers, or tuples of the form
       `(x,y)`.
 
     EXAMPLES:
@@ -182,13 +185,16 @@ def complex_cubic_spline(pts):
         sage: cs = complex_cubic_spline(pts)
         sage: fx = lambda x: cs.value(x).real
         sage: fy = lambda x: cs.value(x).imag
-        sage: show(parametric_plot((fx, fy), (0, 2*pi)))
-        sage: m = Riemann_Map([lambda x: cs.value(real(x))], [lambda x: cs.derivative(real(x))], 0)
-        sage: show(m.plot_colored() + m.plot_spiderweb())
+        sage: from math import pi
+        sage: show(parametric_plot((fx, fy), (0, 2*pi)))                                # needs sage.plot
+        sage: m = Riemann_Map([lambda x: cs.value(real(x))],
+        ....:                 [lambda x: cs.derivative(real(x))], 0)
+        sage: show(m.plot_colored() + m.plot_spiderweb())                               # needs sage.plot
 
     Polygon approximation of a circle::
 
-        sage: pts = [e^(I*t / 25) for t in range(25)]
+        sage: from cmath import exp
+        sage: pts = [exp(1j * t / 25) for t in range(25)]
         sage: cs = complex_cubic_spline(pts)
         sage: cs.derivative(2)
         (-0.0497765406583...+0.151095006434...j)
@@ -273,6 +279,7 @@ cdef class CCSpline:
             sage: cs = complex_cubic_spline(pts)
             sage: cs.value(4 / 7)
             (-0.303961332787...-1.34716728183...j)
+            sage: from math import pi
             sage: cs.value(0) - cs.value(2*pi)
             0j
             sage: cs.value(-2.73452)
@@ -304,6 +311,7 @@ cdef class CCSpline:
             sage: cs = complex_cubic_spline(pts)
             sage: cs.derivative(3 / 5)
             (1.40578892327...-0.225417136326...j)
+            sage: from math import pi
             sage: cs.derivative(0) - cs.derivative(2 * pi)
             0j
             sage: cs.derivative(-6)