// gcc n.c -Wmost p.o -fnested-functions ../../../code/invert.o
#include "h.h"
typedef struct {R p1; R p2; R per;} Suv;
typedef struct {pt pos[N];} sv;
Suv S[] = {
// {0.347111, 0.532728, 6.324449},
{ 0.347116767811511,  0.532725000747374,   6.325913462237277},
//{0.306893, 0.125507, 6.235641},
{ 0.306893421600765,  0.125506565818490,   6.234674843059651},
//{0.392955, 0.097579, 7.003907},
{ 0.392955496410768,  0.097578967671380,   7.003709632880691},
{0.184279, 0.587188, 63.534541},
{0.080584, 0.588836, 21.270975}, {0.083300, 0.127889, 10.466818},
{0.464445, 0.396060, 14.893911}, {0.439166, 0.452968, 28.6700278},
{0.383444, 0.377364, 25.840618}, {0.405916, 0.230163, 13.865763},
//{0.350112, 0.079339, 79.475875},
{ 0.351587929760850,  0.078670498782693,  79.533325949045079},
{0.559064, 0.349192, 55.501762},
{0.513938, 0.304736, 17.328370}, {0.282699, 0.327209, 10.962563},
{0.416822, 0.330333, 55.789829}, {0.417343, 0.313100, 54.207599}
};
void pSuv(Suv p){printf("{%18.15f, %18.15f, %19.15f}\n", (d)p.p1, (d)p.p2, (d)p.per);}
void p3(char * m, d e[3]){int y; for(y=0; y<3; ++y) printf("%15e ", e[y]);
    printf(" %s\n", m);}
int solve(int n, double A[n][n], double B[n]);

int main(){
// printf("%d\n", (int)sizeof(R));
assert(N==3);
d ce(sv t, sv u){R e = 0; {int n=N; while(n--){int dd=D;
    while(dd--) e += s(t.pos[n].dv[dd].x  - u.pos[n].dv[dd].x )
                    +s(t.pos[n].dv[dd].xd - u.pos[n].dv[dd].xd);}}
    return sqrt(e);}
assert(N==3);
{int j; for(j=0; j<sizeof(S)/sizeof(S[0]); ++j) {
  R sc = 1., tc = 1/sqrtl(sc);
  void td(Suv T, d er[3], char *r){
    R p1 = tc*T.p1, p2 = tc*T.p2;
    sv t = {{{{{-sc, p1}, {0, p2}}, 1},
             {{{ sc, p1}, {0, p2}}, 1},
             {{{0, -2.*p1}, {0, -2.*p2}}, 1}}};
    sv u = t;
    spin(t.pos, sc/tc*T.per, 0, -1);
    if(r) {printf("%s error = %e\n", r, ce(u, t)); prx(r, t.pos);}
    er[0] = t.pos[0].dv[1].x;
    er[1] = t.pos[1].dv[1].x;
    er[2] = t.pos[2].dv[0].x;}
  {Suv s = S[j]; d J[3][3], e[3]; d const dx = 1.e-6; td(s, e, "initial"); pSuv(s);
    {Suv t = s; t.p1  += dx; td(t, J[0], 0);}
    {Suv t = s; t.p2  += dx; td(t, J[1], 0);}
    {Suv t = s; t.per += dx; td(t, J[2], 0);}
    {int y; for(y=0; y<3; ++y) p3("Ji", J[y]);}
    {int u, v; for(u=0; u<3; ++u) for(v=0; v<3; ++v)
       {J[u][v] -= e[v]; J[u][v] /= dx;}}
    {int i, j; for(i=0; i<3; ++i) for(j=i+1; j<3; ++j) // Transpose J.
       {d t=J[i][j]; J[i][j]=J[j][i]; J[j][i] = t;}} 
    {int y; for(y=0; y<3; ++y) p3("J", J[y]);}
    p3("e", e); solve(3, J, e); p3("dm", e);
    if(0) {d f[3]; int j; for(j=0; j<3; ++j) { // did I solve the equations?
      f[j]=0; {int k; for(k=0; k<3; ++k) f[j] += J[j][k]*e[k];}}
      p3("ec?", f);} // yes, ec = e
    {Suv t = s; t.p1 -= e[0]; t.p2 -= e[1]; t.per -= e[2]; pSuv(t);
    {d f[3]; td(t, f, "final"); printf("j = %d", j); p3("ne", f);}}}
  getchar();}} return 0;}
/*
initial error = 3.405540e-03
{ 0.347116767811511,  0.532725000747374,   6.325913462237277}
final error = 3.521189e-08
j = 0  -2.405090e-09    2.962021e-10   -1.450535e-08  ne
{ 0.347116888201945,  0.532724945343319,   6.325913983190598}
final error = 8.363796e-11
j = 0  -9.260606e-12    2.181972e-11   -5.650207e-12  ne

initial error = 1.883057e-03
{ 0.306893421600765,  0.125506565818490,   6.234674843059651}
final error = 1.454109e-09
j = 1   3.924495e-10    2.194644e-10   -1.820373e-10  ne
{ 0.306893420461400,  0.125506567025518,   6.234674838735611}
final error = 1.110243e-11
j = 1   2.376870e-15    1.497913e-14    3.399990e-14  ne


initial error = 4.125304e-04
{ 0.392955496410768,  0.097578967671380,   7.003709632880691}
final error = 1.975702e-09
j = 2  -1.611148e-12    8.466832e-12    1.535382e-10  ne
{ 0.392955495765450,  0.097578967947484,   7.003709625621720}
final error = 1.796268e-09
j = 2   2.185806e-15    6.841302e-15   -2.638560e-14  ne

initial error = 4.420275e-04
{ 0.184278498839971,  0.587188172147161,  63.534353305621387}
final error = 2.468039e-09
j = 3  -1.112306e-09    1.124441e-09   -3.398075e-10  ne

initial error = 3.101838e-03
{ 0.080584226622695,  0.588836089776186,  21.272337407346157}
final error = 5.054873e-09
j = 4   1.251614e-09    1.296648e-09   -1.515461e-09  ne

initial error = 3.558149e-03
{ 0.083300072234192,  0.127889254772242,  10.464849503339641}
final error = 3.630933e-08
j = 5  -2.045542e-09   -1.324287e-09    4.754263e-09  ne

initial error = 8.893301e-04
{ 0.464445173012218,  0.396060014316101,  14.894305165792927}
final error = 6.364526e-09
j = 6   1.166769e-09    2.007906e-10    1.544340e-09  ne

initial error = 1.674701e-03
{ 0.439165917738381,  0.452967643050678,  28.669270887120557}
final error = 9.210015e-09
j = 7   4.673497e-09   -2.302974e-09    1.718763e-09  ne

initial error = 2.953727e-03
{ 0.383443518308644,  0.377363695824407,  25.839236327190715}
final error = 2.567724e-08
j = 8  -1.539128e-09   -5.440134e-09   -1.100489e-08  ne

initial error = 2.884289e-03
{ 0.405915567153430,  0.230163125722042,  13.867123429011103}
final error = 5.948509e-09
j = 9  -3.627975e-11    2.283602e-10    2.190852e-09  ne

initial error = 3.068110e-01
{ 0.351587929760850,  0.078670498782693,  79.533325949045079}
final error = 1.101187e-01
j = 10   3.765722e-03    2.910091e-03   -4.047374e-02  ne
{ 0.348787753732397,  0.079904483693584,  79.263017178155820}
final error = 5.892616e-01
j = 10   1.779153e-02    1.371193e-02   -1.851782e-01  ne

initial error = 1.598727e-03
{ 0.559064255492877,  0.349191564810128,  55.502423599029491}
final error = 3.613604e-08
j = 11   5.980765e-09   -1.374165e-08   -7.928083e-09  ne

initial error = 1.003007e-03
{ 0.513938537496037,  0.304735919326180,  17.328834020206632}
final error = 1.077812e-09
j = 12   4.630994e-10   -2.959527e-10   -3.871984e-10  ne

initial error = 1.437527e-03
{ 0.282702100444260,  0.327208972140476,  10.963303244900368}
final error = 1.386489e-08
j = 13  -1.107879e-09   -6.303651e-10   -7.343089e-09  ne

initial error = 1.461054e-03
{ 0.416821663676514,  0.330333291922373,  55.789149974024340}
final error = 9.120777e-06
j = 14   2.604301e-10   -8.320222e-10    5.385103e-11  ne

initial error = 6.213411e-04
{ 0.417342634073892,  0.313099936304987,  54.207834463918431}
final error = 8.656209e-06
j = 15   2.311253e-11   -4.760028e-10   -9.786899e-11  ne
*/