#include <cairo.h>
#include <math.h>
#include <float.h>
#include <stdio.h>

#undef MAX
#define MAX(a, b) ((a) > (b) ? (a) : (b))

static cairo_matrix_t m1 = {
  .xx = -38.835274957135283,
  .yx = -5.9002831204844393,
  .xy = 310.8899410274783,
  .yy = 43.796450025375023,
  .x0 = -266.44259062779747,
  .y0 = 242.52585686041357,
};
int x_offset1 = 0;
int y_offset1 = 0;

int *x_offset = &x_offset1, *y_offset = &y_offset1;
cairo_matrix_t *matrix = &m1;

int main (int argc, char *argv[])
{
  cairo_matrix_t m;

  m = *matrix;
  cairo_matrix_translate (&m, *x_offset, *y_offset);

  printf ("DBL_EPSILON: %g\n", DBL_EPSILON);

    if (m.x0 != 0.0 || m.y0 != 0.0) {
        double tx, ty, norm;
        int i, j;

        /* pixman also limits the [xy]_offset to 16 bits so evenly
         * spread the bits between the two.
         *
         * To do this, find the solutions of:
         *   |x| = |x*m.xx + y*m.xy + m.x0|
         *   |y| = |x*m.yx + y*m.yy + m.y0|
         *
         * and select the one whose maximum norm is smallest.
         */
        tx = m.x0;
        ty = m.y0;
        norm = MAX (fabs (tx), fabs (ty));

	for (i = -1; i < 2; i+=2) {
	    for (j = -1; j < 2; j+=2) {
		double x, y, den, new_norm;

		den = (m.xx - i) * (m.yy - j) - m.xy * m.yx;
		if (fabs (den) < DBL_EPSILON)
		    continue;

		x = m.y0 * m.xy - m.x0 * (m.yy - j);
		y = m.x0 * m.yx - m.y0 * (m.xx - i);

		den = 1 / den;
		x *= den;
		y *= den;

		new_norm = MAX (fabs (x), fabs (y));
		if (norm > new_norm) {
		    norm = new_norm;
		    tx = x;
		    ty = y;
		}
	    }
	}

	tx = floor (tx);
	ty = floor (ty);
	*x_offset = -tx;
	*y_offset = -ty;
	cairo_matrix_translate (&m, tx, ty);
    } else {
	*x_offset = 0;
	*y_offset = 0;
    }

    printf("matrix orig: xx/yx/xy/yy/x0/y0 %f %f %f %f %f %f\n",
           fabs (matrix->xx), fabs (matrix->yx), fabs (matrix->xy),
           fabs (matrix->yy), fabs (matrix->x0), fabs (matrix->y0));
    printf("matrix new : xx/yx/xy/yy/x0/y0 %f %f %f %f %f %f\n",
           fabs (m.xx), fabs (m.yx), fabs (m.xy),
           fabs (m.yy), fabs (m.x0), fabs (m.y0));
    printf("offset x/y: %d/%d\n", *x_offset, *y_offset);

  return 0;
}