/* See {hr2_pmap_opt_affine.h}. */
/* Last edited on 2024-09-05 19:20:56 by stolfi */

#define _GNU_SOURCE
#include <stdint.h>
#include <assert.h>
#include <math.h>
 
#include <bool.h>
#include <r2.h>
#include <affirm.h>

#include <hr2_pmap_opt_affine.h>

void hr2_pmap_encode_affine(hr2_pmap_t *M, double y[])
  { double A00 = M->dir.c[0][0];
  
    y[0] = M->dir.c[0][1] / A00;
    y[1] = M->dir.c[0][2] / A00;
    y[2] = M->dir.c[1][1] / A00;
    y[3] = M->dir.c[1][2] / A00;
    y[4] = M->dir.c[2][1] / A00;
    y[5] = M->dir.c[2][2] / A00;
    
    if (y[2]*y[5]-y[3]*y[4] < 0)
      { for (int32_t j = 0; j <= 1; j++)
          { double tmp = y[2+j]; y[2+j] = y[4+j]; y[4+j] = tmp; }
      }
  }

void hr2_pmap_decode_affine(double y[], sign_t sgn, hr2_pmap_t *M)
  { 
    r3x3_ident(&(M->dir));
    M->dir.c[0][1] = y[0];
    M->dir.c[0][2] = y[1];
    
    if (sgn*(y[2]*y[5]-y[3]*y[4]) > 0)
      { 
        M->dir.c[1][1] = y[2];
        M->dir.c[1][2] = y[3];

        M->dir.c[2][1] = y[4];
        M->dir.c[2][2] = y[5];
      }
    else
      { M->dir.c[1][1] = y[4];
        M->dir.c[1][2] = y[5];

        M->dir.c[2][1] = y[2];
        M->dir.c[2][2] = y[3];
      }
    
    r3x3_inv(&(M->dir), &(M->inv));

    /* Just in case: */
    M->inv.c[1][0] = 0.0;
    M->inv.c[2][0] = 0.0;
  }