From: ivan@siamics.netREMOVE.invalid   
      
   >>>>> On 2025-09-19, I wrote:   
      
    > m : matrix ([ a, b, 0 ],   
    >             [ c, d, 0 ],   
    >             [ 0, 0, 1 ])$   
    > s : linsolve (flatten ([  args (m . [ u, v, 1 ] - [ 1, 0, 1 ]),   
    >                           args (m . [ w, t, 1 ] - [ 0, 1, 1 ]) ]),   
    >               [ a, b, c, d ])$   
      
    > [Upon re-reading this article before submission I've found that   
    > I've managed to break my promise above to use the point p[k] so   
    > that p[k] - p[1] is the longest one.  Hence, the resulting ellipse   
    > is defined by /a/ diameter p[1], p[2] - and a "side" point p[3].   
    > I /think/ that by using (p[1] + p[3]) / 2 for the center point   
    > here, I'd get an ellipse with p[1], p[3] for its major axis.]   
      
   	More important are symmetry considerations, though.  The   
   	source figure we're transforming, a circle, is symmetric with   
   	respect to any line going through the origin.  When we take two   
   	/orthogonal/ unit vectors from the origin, scale along their   
   	directions, then rotate about origin, still then translate to   
   	some new origin, we at each point preserve these reflectional   
   	symmetries of the figure.  Hence, we get an ellipse with those   
   	two vectors, now transformed, as its axes of symmetry.   
      
   	Above, however, we also apply shearing, which, in general,   
   	doesn't preserve the symmetries.   
      
   	I believe the correct way to achieve the desired result is   
   	as follows.   
      
   C : matrix ([ 1, 0,  0 ],   
               [ 0, 1,  0 ],   
               [ 0, 0, -1 ])$   
   m : matrix ([ 'a, 'b, 0 ],   
               [ 'c, 'd, 0 ],   
               [  0,  0, 1 ])$   
   s : solve (flatten ([  args (m . [   'u, 'v, 1 ] - [ 1,  0, 1 ]),   
                          args (m . [ - 'v, 'u, 1 ] - [ 0, 'q, 1 ]),   
                          block ([ z ], z : m . [ 'w, 't, 1 ],   
                                 z . C . transpose (z)) ]),   
              [ 'a, 'b, 'c, 'd, 'q ])$   
   S : subst (s[1], m)$   
   optimize (S);   
   /**   
                                          2        2          1   
       block ([%1, %2, %3, %4, %5], %1 : u , %2 : v , %3 : -------,   
                                                           %2 + %1   
                      1   
   %4 : -----------------------------,   
          3                         3   
        (v  + %1 v) w - t u %2 - t u   
                     2                 4            2        4   
   %5 : sqrt ((- %1 w ) - 2 t u v w + v  + (2 %1 - t ) %2 + u ),   
   [  u %3      v %3     0 ]   
   [                       ]   
   [ v %4 %5  - u %4 %5  0 ])   
   [                       ]   
   [    0         0      1 ]   
    **/   
      
   	The resulting S matrix isn't particularly elegant, but isn't   
   	too complex, either.   
      
   	Note that the argument of sqrt () above might happen to be   
   	negative.  So far as I can tell, it happens when [ w, t ] ^^ 2   
   	is greater than [ u, v ] ^^ 2, which is something I intend to   
   	make sure never happens in my code (as in, I'm going to swap   
   	vectors in such a case before computing the matrix.)   
      
   	Let's apply it to our points:   
      
   p : matrix ([  0, 6, 2 ],   
               [ -5, 0, 1 ],   
               [  1, 1, 1 ])$   
   q : flatten (args (col (p, 1) + col (p, 2))) / 2;   
   /**       5   
       [3, - -, 1]   
             2     **/   
      
   R :   at (S, [ u = p[1, 2] - q[1], v = p[2, 2] - q[2],   
                  w = p[1, 3] - q[1], t = p[2, 3] - q[2] ])   
       . at (L, [ x = q[1], y = q[2] ])$   
   10309 * R;   
   /** [          2028                 1690               - 1859       ]   
       [                                                               ]   
       [                             3/2                 3/2           ]   
       [ - 65 sqrt(2) sqrt(399)  39 2    sqrt(399)  195 2    sqrt(399) ]   
       [                                                               ]   
       [           0                     0                10309        ] **/   
      
   Q : transpose (R) . C . R$   
   10309 * Q;   
   /** [  726    - 60  - 2328 ]   
       [                      ]   
       [  - 60   748    2050  ]   
       [                      ]   
       [ - 2328  2050   1800  ] **/   
      
   maplist (lambda ([ i ], transpose (col (p, i)) . Q . col (p, i)), [ 1, 2, 3 ]);   
   /** [0, 0, 0] */   
      
   plot2d (0 = [ x, y, 1 ] . Q . transpose ([ x, y, 1 ]),   
           [ x, -2, 7 ], [ y, -7, 2 ])$   
   /**           2 +------------------------------------------------+   
                   |    +     |    +     +    +     +    +     +    |   
                 1 |-+        |       ********************        +-|   
                   |          |   *****                   ****      |   
                 0 |----------|-------------------------------------|   
                   |         *|*                                **  |   
                -1 |-+     ***|                                  *+-|   
                   |       *  |                                   **|   
                -2 |-+    *   |                                   +-|   
              y    |      *   |                                    *|   
                -3 |-+    *   |                                   +-|   
                   |      *   |                                  ** |   
                -4 |-+    **  |                                 **+-|   
                   |       ***|                               **    |   
                -5 |-+       *|*                           ***    +-|   
                   |          |*****                   *****        |   
                -6 |-+        |    ********************           +-|   
                   |    +     |    +     +    +     +    +     +    |   
                -7 +------------------------------------------------+   
                  -2   -1     0    1     2    3     4    5     6    7   
                                           x   
    **/   
      
   	Aside of the actual code, the problem seems to be solved.   
      
   	Thoughts?   
      
   --- SoupGate-Win32 v1.05   
    * Origin: you cannot sedate... all the things you hate (1:229/2)