// clang a.c -O3 ...
// gcc a.c -Wmost -std=c99 -fnested-functions
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <unistd.h>
#define C const

int main(){
FILE * rp = fopen("rprt", "w");
int C X=9; int C Y=9;
typedef double R;
typedef struct{R x; R y; R xd; R yd;} pt;
pt a[X][Y];
short fx(double x){union {double d; long int i;} q;
   q.d = (x/Y)*750 + 0x3.p51; return q.i;}
{short w[] = {X, Y}; write(1, w, 2*2);}
void plt(){for(int x=0; x<X; ++x) {short w[2*Y];
   for(int y=0; y<Y; ++y) {w[2*y] = fx(a[x][y].x); w[2*y+1] = fx(a[x][y].y);}
   write(1, w, 4*Y);}}
{int i=X; while(i--){int j=Y; while(j--) 
   a[i][j] = (pt){i + 0.5*j, (sqrt(3)/2)*j, 0, 0};}}
 plt();
 double t = 0, dt = 0.1, tp = -dt/2, Xm = a[X-1][Y-1].x, Ym = a[X-1][Y-1].y;
    while(0<1) {
      if(t > tp) {tp += 1; fprintf(rp, "t = %11.7f\n", t);
        for(int i=0; i<X; ++i) for(int j=0; j<Y; ++j) {pt p = a[i][j];
          fprintf(rp, "%2d %2d %11.7f %11.7f %11.7f %11.7f\n",
              i, j, p.x, p.y, p.xd, p.yd);}}
// Symmetries??
if(1) {double c6 = sqrt(3)/2; R eps=0000001;
if(X!=Y) exit(fprintf(rp, "Sym!!\n"));
for(int i=0; i<X; ++i) for(int j=0; j<X; ++j) {
  pt p=a[i][j], q=a[j][i], r=a[X-1-i][Y-1-j];
  // Verify symmetry about long diagonal:
  if(fabs(p.x - ( 0.5*q.x + c6*q.y))> eps) exit(fprintf(rp, "Fox %d %d\n", i, j));
  if(fabs(p.y - (-0.5*q.y + c6*q.x))> eps) exit(fprintf(rp, "Foy %d %d\n", i, j));
  // Verify central symmetry.
  if(fabs(p.x + r.x - Xm) > eps) exit(fprintf(rp, "Fux %d %d\n", i, j));
  if(fabs(p.y + r.y - Ym) > eps) exit(fprintf(rp, "Fuy %d %d\n", i, j));
  }}
       if(t > 10.55) break;
       R dx[Y-1], dy[Y-1];
       for(int j=0; j<Y-1; ++j) {
          dx[j] = a[0][j+1].x - a[0][j].x;
          dy[j] = a[0][j+1].y - a[0][j].y;}
       for(int i=0; i<X-1; ++i) {double
          dex = a[i+1][0].x - a[i][0].x,
          dey = a[i+1][0].y - a[i][0].y;
       for(int j=0; j<Y-1; ++j) {double
          p0 = dt/(dy[j]*dex - dx[j]*dey),
          ddx = a[i][j+1].x - a[i+1][j].x,
          ddy = a[i][j+1].y - a[i+1][j].y;
          a[i][j].xd -= p0*ddy; a[i][j].yd += p0*ddx;
          a[i][j+1].xd -= p0*dey; a[i][j+1].yd += p0*dex;
          a[i+1][j].xd += p0*dy[j]; a[i+1][j].yd -= p0*dx[j];
          dx[j] = a[i+1][j+1].x - a[i+1][j].x;
          dy[j] = a[i+1][j+1].y - a[i+1][j].y;
          dex = a[i+1][j+1].x - a[i][j+1].x;
          dey = a[i+1][j+1].y - a[i][j+1].y;
        R p1 = dt/(dy[j]*dex - dx[j]*dey);
          a[i][j+1].xd -= p1*dy[j]; a[i][j+1].yd += p1*dx[j];
          a[i+1][j].xd += p1*dey;  a[i+1][j].yd -= p1*dex;
          a[i+1][j+1].xd += p1*ddy; a[i+1][j+1].yd -= p1*ddx;
          a[i][j].x += dt*a[i][j].xd; a[i][j].y += dt*a[i][j].yd;}
       a[i][Y-1].x += dt*a[i][Y-1].xd; a[i][Y-1].y += dt*a[i][Y-1].yd;}
       for(int j=0; j<Y; ++j) {
          a[X-1][j].x += dt*a[X-1][j].xd; a[X-1][j].y += dt*a[X-1][j].yd;}
  plt(); t += dt;}}