Lecture 20: Computing with Functions
====================================

List comprehensions in Python take the form
``[f(x) for x in L]`` where ``L`` is some list and ``f`` some function.
The result is a new list that contains all function values of ``f(x)``
for all elements ``x`` in the list ``L``.  In SageMath we can apply list
comprehensions to build expressions of arbitrary size and shape. 

List Comprehensions
-------------------

We can use a :index:`list comprehension`
to create a sequence of 20 variables,
starting from ``x01``, ``x02``, .., ``x20``.

::

   v = [var('x' + '%02d' % k) for k in range(1,21)]
   v

Observe the :index:`string concatenation` with the proper formatting
of the integer index.  This leads to the list
``[x01, x02, x03, x04, x05, x06, x07, x08, x09, x10, x11, x12, 
x13, x14, x15, x16, x17, x18, x19, x20]``.
Now we raise every variable to the power three.

::

   p = [x^3 for x in v]
   p

and we computed
``[x01^3, x02^3, x03^3, x04^3, x05^3, x06^3, x07^3, x08^3, x09^3, x10^3, 
x11^3, x12^3, x13^3, x14^3, x15^3, x16^3, x17^3, x18^3, x19^3, x20^3]``.
We add up the powers into a sum to make an expression.

::

   e = sum(p)
   print(e, 'has type', type(e))

The sum leads to the expression 
``x01^3 + x02^3 + x03^3 + x04^3 + x05^3 + x06^3 + x07^3 + x08^3 + x09^3
+ x10^3 + x11^3 + x12^3 + x13^3 + x14^3 + x15^3 + x16^3 + x17^3 + x18^3
+ x19^3 + x20^3``.
We can return to the list representation via the ``operands()`` method.

::

   ope = e.operands()
   print(ope, 'has type', type(ope))

and this gives a list
``[x01^3, x02^3, x03^3, x04^3, x05^3, x06^3, x07^3, x08^3, x09^3, x10^3, 
x11^3, x12^3, x13^3, x14^3, x15^3, x16^3, x17^3, x18^3, x19^3, x20^3]``.
Raising the k-th operand to the power k goes as follows.

::

   r = [ope[k]^k for k in range(len(ope))]
   r

and then we see ``[1, x02^3, x03^6, x04^9, x05^12, x06^15, x07^18, x08^21,
x09^24, x10^27, x11^30, x12^33, x13^36, x14^39, x15^42, x16^45, x17^48,
x18^51, x19^54, x20^57]``.
Note that all operands are expressions.

::

   print(r[2], 'has type', type(r[2]))
   print('its degree in' , v[2], ':', r[2].degree(v[2])

With the proper selection of the variable in the argument
of ``degree()`` we get the right degree, which is ``6`` for ``r[2]``.
We can compute the degrees in their variables.

::

   [(r[k]).degree(v[k]) for k in range(len(ope))]

which leads to ``[0, 3, 6, 9, 12, 15, 18, 21, 24, 27, 30,
33, 36, 39, 42, 45, 48, 51, 54, 57]``.

Suppose we want to select those operands of degree less than 13.
Recall the unit_step() function.

::

   print(r[2], ':', unit_step(13 - r[2].degree(v[2])))
   print(r[8], ':', unit_step(13 - r[8].degree(v[8])))

and we see ``x03^6 : 1`` and ``x09^24 : 0``.
Now we can remove all operands with degree higher than 13.

.. index:: filter, lambda

::

   filter = lambda x, k: x*unit_step(13 - x.degree(v[k]))
   s = [filter(r[k],k) for k in range(len(r))]
   s

and we see ``[1, x02^3, x03^6, x04^9, x05^12, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]``
and then build the expression again with ``sum()``.

::

   t = sum(s)
   t

and we obtain ``x05^12 + x04^9 + x03^6 + x02^3 + 1``.

Composing Functions
-------------------

Consider the trajectory of a projectile in space modeled by a parabola,
subject to the following constraints.  At time t = 0 it is launched from
the ground and at time t = 45 it hits the ground 120 miles further.
Assuming constant horizontal speed we create a function f(t) to give
the altitude of the projectile in function of time.
First we model the shape of the trajectory.

::

   y(x) = x*(120 - x)
   plot(y(x), (x, 0, 120), aspect_ratio=1/100, thickness=3, color='red')

.. index:: aspect_ratio, thickness, color

The parameters ``thickness`` and ``color`` embellish the plot,
adjusting the default value for the ``aspect_ratio`` to ``1/100``
alters the display of the shape of the trajectory,
so we do not have to worry about units.
The plot is displayed in :numref:`figparabola`.

.. _figparabola:

.. figure:: ./figparabola.png
    :align: center

    A parabolic trajectory.

Now we want a function in time t.
The assumption of constant horizontal speed implies that x is just a 
rescaling of t.
This means that when t = 45, x must be 120.

::

   x(t) = 120/45*t
   print('x at 0 :', x(0))
   print('x at 45 :', x(45))

and indeed, we see ``x at 0 : 0``
and ``x at 45 : 120`` printed.
Then the altitude is given as the composition of y after x.

::

   f(t) = y(x(t))
   print('altitude at 0 :', f(0))
   print('altitude at 22.5 :', f(22.5))
   print('altitude at 45 :', f(45))

as altitudes at 0, 22.5, and 45, we respectively see ``0``, 
``3600.00000000000`` and ``0``.
With the composition of functions we have separated the geometric 
shape of the trajectory and its evolution in time.
Suppose we want to halve the speed.

::

   h(t) = t/2
   hf(t) = f(h(t))
   print('altitude at 45 :', hf(45))
   print('altitude at 90 :', hf(90))

and the prints confirm that the projectile reaches its peak at 45
and falls to the ground at 90.

We can also make a 3d-plot of the space curve.

::

   parametric_plot3d([x(t), 0, f(t)], (t, 0, 45), thickness=5, color='red')

and this shows the following image:

.. _figparabolain3dims:

.. figure:: ./figparabolain3dims.png
    :width: 50%
    :align: center

    The parabolic trajectory plotted in three dimensions.

Functions Returning Functions
-----------------------------

We can define functions that return functions.
Suppose we are interested in solving equations with Newton's method.

::

   def newton_step(f, x):
       """
       Returns the function that will perform one Newton step
       to solve the equation f(x) = 0.
       """
       n(x) = x - f(x)/diff(f,x)
       return n

Let us try it out on the square function,
so we can compute the square root of two.

::

   sqr(x) = x^2 - 2
   our_sqrt = newton_step(sqr, x)
   our_sqrt

Then we see ``x |--> x - 1/2*(x^2 - 2)/x``.
Observe that if we start at an integer value,
we get rational approximations for the square root of two.

::

   our_sqrt(our_sqrt(our_sqrt(2)))

Starting at ``2``, three iterations of Newton's method yields ``577/408``.
For a repeated application of a function, we multiply strings.

::

   s = 'our_sqrt('*5 + '2' + ')'*5
   s

We will then execute
``our_sqrt(our_sqrt(our_sqrt(our_sqrt(our_sqrt(2)))))``
and we evaluate the expression.

::

   v = eval(s)
   v

This gives ``886731088897/627013566048``.
To compare with an actual approximation for sqrt(2).

::

   print(v.n(digits=30))
   print(sqrt(2).n(digits=30))

and we see

::

   1.41421356237309504880168962350
   1.41421356237309504880168872421

Assignments
-----------

1. Type ``L = [randint(0,99) for i in range(10)]`` to generate a list of 10
   random numbers between 0 and 99.
   Give the SageMath commands for the following operations.

   (a) divide every element in the list by 100;
   (b) convert the list to a list of floating-point numbers of
       a precision of 3 decimal places;
   (c) select all elements in the list that are larger than 0.5;
   (d) compute the sum of the elements in the list.

2. The command ``is_prime()`` applied to a number returns ``True``
   if the number is prime and ``False`` otherwise.
   Use ``is_prime()`` to make a list of all primes less than 1000.
   How many 3-digit primes are there?

3. Do ``P.<x> = PolynomialRing(RR, sparse=False)`` and then

   ``q = P.random_element(degree=10)``.

   Select from ``q`` all terms with positive coefficients
   and make a new polynomial with those terms.

4. For a positive integer :math:`n`, and a list :math:`X` of :math:`n`
   variables :math:`X` = ``[`` ``x1``, ``x2``, .. , ``xn`` ``]``,
   define the function

   .. math::

      H(n,X,k,c) = 2 n X[k] - c X[k] 
      \left( 1 + \sum_{i=0}^{n-1} \frac{i+1}{i+3} X[i] \right)
      - 2 n,

   where :math:`c` is some variable parameter.
   Show the definition of your function is correct for a list of 9 variables.

5. Let \ :math:`n` and \ :math:`k` be positive natural numbers 
   with \ :math:`k < n`
   and consider the polynomial \ :math:`{\displaystyle \sum_{i=0}^{n-1}
   \prod_{j=0}^{k-1} z_{(i+j) \mbox{ mod } n}}.`

   Define a function that takes on input a list of variables
   (where \ :math:`n` is the length of the list) and a value
   for \ :math:`k`.  On output is the expression as defined
   by the polynomial above.  Hint: use ``prod`` to make the product
   of a list of expressions.

   Show the result of your function
   on a sequence of 10 variables and for \ :math:`k = 3`.

6. Execute the ``newton_step`` function on the ``our_sqr`` function,
   defined again as ``x^2 - 2``.  Newton's method has the property
   of converging quadratically, that is: in each step the correct
   number of decimal places doubles.  Starting at ``2``, perform 10
   steps with Newton's method and verify the progess of the accuracy
   of the rational approximations for the square root of 2?
   Do you observe quadratic convergence?

7. Define a Python function ``makeHorner`` which takes on input
   the coefficients of a polynomial and the symbol of the variable
   in the polynomial.  For example, if ``C`` is ``[c0, c1, c2, c3, c4]``
   and the variable is ``x``, then ``makeHorner(C, x)`` returns

   .. code::

      (((c4*x + c3)*x + c2)*x + c1)*x + c0

   Your function ``makeHorner`` should directly build the polynomial
   and not apply the ``horner()`` to an expanded polynomial.