Fractional linear transformations in Asymptote

Question

I'm using Asymptote to create some figures illustrating constructions in the disc model of the hyperbolic plane, and I need to be able to transform paths by fractional linear transformations (FLTs), but haven't found any built-in functionality that does what I need.

There are two things I need to be able to do:

1. Given two pairs of points (a,b) and (c,d) on the unit circle, find the FLT that maps a to b, c to d, and preserves the unit circle.
2. Given an FLT and a path (in particular, an arc representing a geodesic), find the path that is the image of the original path under the FLT.

If there is a pre-existing package that does some or all of this, that would be ideal. If not, the following are the stumbling blocks to my writing an implementation of all of this myself.

• An FLT is determined by a 4-tuple, with scalar multiples of each other determining the same FLT; in particular, Item (1) above can be reduced to finding the kernel of a 4x4 matrix with rank 3. Asymptote includes a routine ("solve") to find the solution of Ax=b when A is non-singular, but I don't know how to use it to solve Ax=0 when A is singular. (This is the main difficulty I'm facing.)
• Once the coefficients determining the FLT are known, transforming points is just a matter of complex multiplication, division, and addition, so I can write a function F that takes as input the FLT's coefficients and the initial points, and gives as output the images of those points. What I really want in Item (2) is a similar function G that operates not on points but on paths. I can probably hack together a way of making the function operate on paths once the points version is available, but it'll take me some digging in the documentation, so if there's a natural (pre-defined?) way to define a function G that acts on paths by transforming all the control points and tangent vectors according to F (and DF), it would of course be better to use that.

Answer 1

Here's an ad-hoc function that returns a unit vector that, given a singular matrix A, returns a unit vector that is very close to being in the kernel of A. If A is nonsingular, the unit vector returned is nonsense. The basic idea is to invert A+e, where e is a small amount of random noise. If A is singular, then its kernel will be (very nearly) an eigenspace of (A+e)^{-1} with astronomically large eigenvalue. So for a random vector v, (A+e)^{-1}(v) will be huge and very nearly lie in the kernel of A.

Two caveats. First, I have tested this code on exactly one case. Second, in addition to the inherent randomness described above, it depends on the numerical stability of the built-in matrix solver.

import stats;

real length(real[] v) {
  real ans = 0;
  for (real a : v) ans += a*a;
  return sqrt(ans);
}

real[] kernel_vector(real[][] matrix4x4) {
  real[][] A = copy(matrix4x4);
  for (int i = 0; i < A.length; ++i) {
    for (int j = 0; j < A[i].length; ++j) {
      A[i][j] += 1e-5 * Gaussrand();
    }
  }
  real[] v = new real[A[0].length];
  for (int i = 0; i < v.length; ++i)
    v[i] = Gaussrand();
  real[] ans = solve(A, v);
  return ans / length(ans);
}

For your second point, applying nonlinear transformations to paths instead of just points: see https://tex.stackexchange.com/a/161146/484.