tangent lines to a circle
=========================

Given a fixed circle in the plane, compute the lines tangent to
the circle and passing through the origin.
In :numref:`figtouchcircle` we see a general line through the origin
and two lines touching the circle.

.. _figtouchcircle:

.. figure:: ./touchcircle.png
    :align: center

    A general line through a circle and two lines touching a circle.

The tangent lines are special: at the points where the lines touch
the circle, we have a double solution, a solution of multiplicity two.
One method is to consider the one parameter family of lines through
the origin and intersect this family with the polynomials which express
the singularity condition.

lines through the origin intersecting a circle
----------------------------------------------

The polynomials which express all lines through the origin
intersecting a fixed circle, fixed by its center and radius,
are returned by a function.

::

   def polynomials(a, b, r):
       """
       Returns string representations of two polynomials:
       1) a circle with radius r centered at (a, b);
       2) a line through the origin with slope s.
       """
       crc = '(x - %.15e)^2 + (y - %.15e)^2 - %.15e;' % (a, b, r**2)
       lin = 'y - s*x;'
       return [crc, lin]

There are two equations, one for the circle and one for the line.
The variables are two coordinates ``x``, ``y``, and the slope ``s``.
When given two equations in three variables we expect a one
dimensional solution set.  To represent this space curve,
we intersect the curve with a general hyperplane and compute
the points on the curve and on the hyperplane.

The code snippet below defines the problem for a circle
centered at ``(3, 2)`` with radius one.  The ``embed`` function
returns the original polynomials with one general hyperplane added
and also one slack variable.  
The blackbox solver computes the generic points.

::

   pols = polynomials(3, 2, 1)
   from phcpy.sets import embed
   from phcpy.solver import solve
   embpols = embed(3, 1, pols)
   embsols = solve(embpols)

The tuple ``(embpols, embsols)`` is a numerical representation for
the set of all lines through the origin intersecting a fixed circle.

As a sanity check, consider a point on the set of all lines as
in the left of :numref:`figtouchcircle`.
Such a point is for instance the line with slope one.
The coordinates for the intersection points,
as can be seen from :numref:`figtouchcircle` are ``(2, 2)`` and ``(3, 3)``.  
In the code below, the intersection point ``(2, 2)`` is joined with
the slope ``1`` in a solution string, called ``point``.

::

   from phcpy.solutions import make_solution
   point = make_solution(['x', 'y', 's'], [2, 2, 1])
   ismb = is_member(embpols, embsols, 1, point)

The call to ``is_member`` returns a boolean value,
so ``ismb`` should hold the value ``True`` for this point.

defining the equations for the singular locus
---------------------------------------------

The two tangent lines to the circle are two special solutions.
At any other line through the origin, the line intersects the
circle at two distinct complex solutions, but at the tangent lines,
the two intersection points collide into a double solution.
At a double solution, the Jacobian matrix of the system no longer
has full rank.  Instead of using the determinant of the matrix of
all first order partial derivatives, the equations we use express
that there is a nonzero combination of the columns of the Jacobian
matrix which yields the zero vector.

The equations for the singular locus are defined by the
function ``jacobian``.  For the circle centered at ``(3, 2)``,
the polynomial equations are obtained as follows:

::

    pols = jacobian(3, 2)
    for pol in pols:
        print(pol)

What is printed is

::

    2*(x-3.000000000000000e+00)*L1 + 2*(y-2.000000000000000e+00)*L2;
    -s*L1 + L2;
    (-0.94944388496-0.313936791907*i) +(0.253472461117-0.967342602936*i)*L1 \
    +(-0.209901989746-0.97772243234*i)*L2;

In the first equation we recognize the two partial derivatives
of :math:`(x-3)^2 + (y-2)^2`, multiplied with the multipliers
``L1`` and ``L2``.  The second equation is derived 
from :math:`y - s x = 0`.
If we have a nonzero combination of the columns of the Jacobian matrix
which yields the zero vector, then any nonzero multiple of the multipliers
also defines such a nonzero combination.
The last equation is a linear equation in the multipliers only,
requires that the multipliers are nonzero, and at the same time
fixing one combination among all nonzero multiples.

The definition of the function ``jacobian`` depends on another
function which returns a linear equation with random coefficients.

::

   def jacobian(a, b):
       """
       For the circle centered at (a, b),
       returns the equations which define the points
       where the Jacobian matrix is singular,
       as a random linear combination of the columns.
       Random complex coefficients are generated to
       scale the multiplier variables.
       """
       eq1 = '2*(x-%.15e)*L1 + 2*(y-%.15e)*L2;' % (a, b)
       eq2 = '-s*L1 + L2;'
       eq3 = random_hyperplane(['L1', 'L2'])
       return [eq1, eq2, eq3]

To avoid badly scaled coefficients, the complex numbers are generated
on the unit circle, but the function ``random_complex`` below.

::

   def random_complex():
       """
       Returns a random complex number on the unit circle.
       """
       from math import cos, sin, pi
       from random import uniform
       theta = uniform(0, 2*pi)
       return complex(cos(theta), sin(theta))

The imaginary unit in Python is represented by ``j``
whereas for phcpy, the imaginary unit is represented by ``i`` and ``I``.
Therefore, the function ``random_hyperplane`` replaces the ``j`` by ``i``.

::

   def random_hyperplane(vars):
       """
       Returns a linear equation in the variables in
       the list vars, with random complex coefficients.
       """
       cf0 = str(random_complex())
       tf0 = cf0.replace('j', '*i')
       result = tf0
       for var in vars:
           cff = str(random_complex())
           tcf = cff.replace('j', '*i')
           result = result + '+' + tcf + '*' + var
       return result + ';'

The function ``jacobian(3, 2)`` returned three equations in the
two coordinates ``x``, ``y``, the slope ``s``, 
the multipliers ``L1``, and ``L2``; five variables in total.
In five dimensional space, three equations define a two dimensional set.

For a numerical representation of this two dimensional set,
two random linear equations are added with the ``embed`` function
and the generic points are computed with the blackbox solver
as done in the code snippet below.

::

    from phcpy.sets import embed
    from phcpy.solver import solve
    embpols = embed(5, 2, pols)
    embsols = solve(embpols)

The number of generic points equals three.

intersecting two algebraic sets
-------------------------------

We have two algebraic sets:

1. The set of all lines through the origin intersecting a fixed circle.
   The degree of this set is four.

2. The set of all intersection points of a line through the origin
   and a fixed circle where the Jacobian matrix is singular.
   The degree of this set is three.

Before we can intersect the two algebraic sets, we have to ensure that
their ambient space is the same.  The first set involves only the
variables ``x``, ``y``, and ``s``, but not the multiplier variables
``L1`` and ``L2`` which occur in the second algebraic set.
Therefore, to each generic point in the first one dimensional set
we add two values for ``L1`` and ``L2`` and two corresponding linear
equations.  So, the one dimensional set is upgraded to a three
dimensional sets in the same five dimensional space in where the
second two dimensional set lives.  Because we can choose any values
for ``L1`` and ``L2`` in this upgrade of the first set, the dimension
of the first set increase from one to three.

Add two variable names ``L1`` and ``L2``, both with values one
and two slack variables ``zz2`` and ``zz3`` with zero values
is done by the function ``extend_solutions``.

::

   def extend_solutions(sols):
       """
       To each solution in sols, adds L1 and L2 with values 1,
       and zz2 and zz3 with values zero.
       """
       from phcpy.solutions import make_solution, coordinates
       result = []
       for sol in sols:
           (vars, vals) = coordinates(sol)
           vars = vars + ['L1', 'L2', 'zz2', 'zz3']
           vals = vals + [1, 1, 0, 0]
           extsol = make_solution(vars, vals)
           result.append(extsol)
       return result

The function is called in the function ``extend`` which upgrades
the first set from a one dimensional to a three dimensional set,
raising its ambient space from a 3-space to the 5-space where the
second set lives.

::

   def extend(pols, sols):
       """
       Extends the witness set with two free variables
       L1 and L2, addition two linear equations,
       and two slack variables zz2 and zz3.
       """
       vars = ['zz2', 'zz3']
       eq1 = 'zz2;'
       eq2 = 'zz3;'
       eq3 = 'L1 - 1;'
       eq4 = 'L2 - 1;'
       extpols = pols[:-1] + [eq1, eq2, eq3, eq4, pols[-1]]
       extsols = extend_solutions(sols)
       return (extpols, extsols)

Note that the order of the equations is important.
The linear equations that cut down the positive dimensional solutions
to isolated points must occur at the end of the list of polynomials.

Also the order of the variables matters.
To ensure that the names of the variables line up in the same order
for both lists of polynomials, the first polynomial for both sets is
prepended with the string ``x-x+y-y+s-s+L1-L1+L2-L2``.

The relevant code snippet to intersect two sets with diagonal homotopies
is shown below.

::

    from phcpy.diagonal import diagonal_solver as diagsolve
    result = diagsolve(dim, w1d, w1eqs, w1sols, w2d, w2eqs, w2sols)
    (eqs, sols) = result

The polynomials and the corresponding generic points
for the first set are in ``w1eqs`` and ``w1sols`` respectively,
for the second set they are in ``w2eqs`` and ``w2sols``.
The dimensions of the two sets are in ``w1d`` and ``w2d``
(which respectively equal three and two)
and the ambient dimension (five) is given in ``dim``.

The number of solutions in list ``sols`` returned by the diagonal solver
equals two, defining the two tangent lines shown at the right 
of :numref:`figtouchcircle`.

The complete script which computes this use case in in the ``examples``
folder in the ``Python/PHCpy3`` directory of the source code.
