/*
  Copyright 2003 Yan Gerard, Isabelle Debled, Paul Zimmermann,
  LLAIC - LORIA/INRIA Lorraine.
  See http://www.loria.fr/~debled/plane/

  This program is free software; you can redistribute it and/or modify it
  under the terms of the GNU General Public License as published by the
  Free Software Foundation; either version 2 of the License, or (at your
  option) any later version.

  This program is distributed in the hope that it will be useful, but WITHOUT
  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for
  more details.

  You should have received a copy of the GNU General Public License along
  with this program; see the file COPYING.  If not, write to the Free
  Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
  02111-1307, USA.

  This program checks if a set of 3D points can be recognized as a discrete
  plan, i.e. if there exists integers (a, b, c, mu) such that all points 
  (x,y,z) satisfy:

                 mu <= a*x + b*y + c*z < mu + |c|

  To compile this program, you need first to install the GNU MP library.
  You can download it at <http://www.swox.se/gmp>, it it is not already
  installed on your computer. The type:

  gcc -O2 -g -I<gmpdir>/include -L<gmpdir>/lib reco.c -o reco -lgmp -lm

  where <gmpdir> is the directory when GNU is installed (/usr or /usr/local
  in most cases).

  To use the program, type "./reco nnn" where nnn is the number of vertices
  you want to generate at random. The program then prints the maximal and
  average time per try, and the maximal/average number of loops in the
  algorithm. For example:

  % ./reco 1000 
  Using double-precision floating-point numbers
  max. time =0 aver. time =0.000000
  max. loops=9 aver. loops=9.000000

  By default the program does only one random generation. Type 
  "./reco nnn kkk" if you want kkk random tries:

  % ./reco 1000 1000
  Using double-precision floating-point numbers
  max. time =10 aver. time =0.670000
  max. loops=15 aver. loops=10.855000

  By default the program uses double precision floating-point numbers.
  If you want to use GNU MP integers, compile the program with -DGMPZ:

  % ./reco 1000 100
  Using GNU MP integers
  max. time =20 aver. time =12.300000
  max. loops=15 aver. loops=10.800000
*/

// #define VERBOSE /* to print coordinates of points and found plane */

#include <stdio.h>
#include <stdlib.h>
#include <sys/types.h>
#include <sys/resource.h>
#include "gmp.h"

#ifdef GMPZ
#define TYPE mpz_t
#define print_coeff(x) mpz_out_str (stdout, 10, x)
#define coeff_init mpz_init
#define coeff_clear mpz_clear
#define coeff_sub mpz_sub
#define coeff_set mpz_set
#define coeff_set_ui mpz_set_ui
#define coeff_set_si mpz_set_si
#define coeff_swap mpz_swap
#define coeff_mul mpz_mul
#define coeff_mul_ui mpz_mul_ui
#define coeff_cmp mpz_cmp
#define coeff_neq(x,y) (mpz_cmp(x,y) != 0)
#define coeff_gt(x,y) (mpz_cmp(x,y) > 0)
#define coeff_ge(x,y) (mpz_cmp(x,y) >= 0)
#define coeff_lt(x,y) (mpz_cmp(x,y) < 0)
#define coeff_le(x,y) (mpz_cmp(x,y) <= 0)
#define coeff_submul mpz_submul
#define coeff_addmul mpz_addmul
#define coeff_sgn mpz_sgn
#define coeff_neg mpz_neg
#define coeff_set_z mpz_set
#else /* default is DOUBLE */
#define DOUBLE
#define TYPE double
#define print_coeff(x) printf ("%1.0f", x)
#define coeff_init(x) 
#define coeff_clear(x)
#define coeff_sub(x,y,z) (x) = (y) - (z)
#define coeff_set(x,y) (x) = (y)
#define coeff_set_ui(x,y) (x) = (double) (y)
#define coeff_set_si(x,y) (x) = (double) (y)
#define coeff_swap(x,y) { double t = x; x = y; y = t; }
#define coeff_mul(x,y,z) (x) = (y) * (z)
#define coeff_mul_ui(x,y,z) (x) = (y) * (double) (z)
#define coeff_cmp(x,y) (((x) > (y)) ? 1 : (((x) == (y)) ? 0 : -1))
#define coeff_neq(x,y) (x != y)
#define coeff_gt(x,y) (x > y)
#define coeff_ge(x,y) (x >= y)
#define coeff_lt(x,y) (x < y)
#define coeff_le(x,y) (x <= y)
#define coeff_submul(x,y,z) (x) -= (y) * (z)
#define coeff_addmul(x,y,z) (x) += (y) * (z)
#define coeff_sgn(x) (((x) > 0.0) ? 1 : (((x) == 0.0) ? 0 : -1))
#define coeff_neg(x,y) (x) = -(y)
#define coeff_set_z(x,y) x = mpz_get_d (y)
#endif

typedef struct {
  TYPE x, y, z;
} 
__point_struct;

typedef __point_struct point_t[1];

void
PrintPoint (point_t M)
{
  putchar ('[');
  print_coeff (M->x);
  putchar (',');
  print_coeff (M->y);
  putchar (',');
  print_coeff (M->z);
  printf ("]");
}

void
InitPoint (point_t M)
{
  coeff_init (M->x);
  coeff_init (M->y);
  coeff_init (M->z);
}

void
ClearPoint (point_t M)
{
  coeff_clear (M->x);
  coeff_clear (M->y);
  coeff_clear (M->z);
}

void
SubPoint (point_t R, point_t M, point_t N)
{
  coeff_sub (R->x, M->x, N->x);
  coeff_sub (R->y, M->y, N->y);
  coeff_sub (R->z, M->z, N->z);
}

void
SetPoint (point_t R, point_t M)
{
  coeff_set (R->x, M->x);
  coeff_set (R->y, M->y);
  coeff_set (R->z, M->z);
}

void
SwapPoint (point_t R, point_t M)
{
  coeff_swap (R->x, M->x);
  coeff_swap (R->y, M->y);
  coeff_swap (R->z, M->z);
}

/* returns sign(Delta(A,B,C)) which is -1, 0, or 1 */
int
SignDelta (point_t A, point_t B, point_t C)
{
  TYPE t1, t2, t3;
  int res;

  coeff_init (t1);
  coeff_init (t2);
  coeff_init (t3);
  coeff_sub (t1, A->y, B->y);
  coeff_sub (t2, A->x, C->x);
  coeff_mul (t1, t1, t2);
  coeff_sub (t2, A->x, B->x);
  coeff_sub (t3, A->y, C->y);
  coeff_mul (t2, t2, t3);
  res = coeff_cmp (t1, t2);
  coeff_clear (t1);
  coeff_clear (t2);
  coeff_clear (t3);
  return res;
}

void
facet_to_normal (point_t w, point_t A, point_t B, point_t C)
{
  point_t u, v;
  TYPE t1;

  InitPoint (u);
  InitPoint (v);
  SubPoint (u, B, A);
  SubPoint (v, C, A);
  coeff_init (t1);
  coeff_mul (w->x, u->y, v->z);
  coeff_submul (w->x, u->z, v->y);
  coeff_mul (w->y, u->z, v->x);
  coeff_submul (w->y, u->x, v->z);
  coeff_mul (w->z, u->x, v->y);
  coeff_submul (w->z, u->y, v->x);
#ifdef GMPZZ
  mpz_gcd (t1, w->x, w->y);
  mpz_gcd (t1, t1, w->z);
  mpz_divexact (w->x, w->x, t1);
  mpz_divexact (w->y, w->y, t1);
  mpz_divexact (w->z, w->z, t1);
#endif
  ClearPoint (u);
  ClearPoint (v);
  coeff_clear (t1);
}

/* return 0 if ok, 1 if error */
int
find_starting_facet (point_t A, point_t B, point_t C, 
                     point_t *S, unsigned long n)
{
  unsigned long i;

  SetPoint (A, S[0]);
  for (i=1; i<n; i++)
    if (coeff_neq (A->x, S[i]->x) || coeff_neq (A->y, S[i]->y))
      break;
  if (i >= n)
    return 1; /* error */
  SetPoint (B, S[i]);
  for (; i<n; i++)
    if (SignDelta(A, B, S[i]) != 0)
      {
        point_t T;
        InitPoint(T);
        SetPoint (T, A);
        SubPoint (A, B, A);
        SubPoint (B, S[i], B);
        SubPoint (C, T, S[i]);
        ClearPoint(T);
        return 0;
      }
  return 1;
}

#ifdef GMPZ
inline void
dotprod (TYPE s, point_t M, point_t v)
{
  coeff_mul (s, M->x, v->x);
  coeff_addmul (s, M->y, v->y);
  coeff_addmul (s, M->z, v->z);
}
#else
#define dotprod(a,b,c) a = b->x * c->x + b->y * c->y + b->z * c->z
#endif

/* puts in M = Mmax - Mmin, in m = maxi - mini */
void
find_minmax_scalprod (point_t M, TYPE *mini, TYPE *maxi,
                      point_t *S, unsigned long n, point_t v)
{
  TYPE s;
  unsigned Mmax, Mmin, i;

  coeff_init (s);
  dotprod (s, S[0], v);
  coeff_set (*mini, s);
  coeff_set (*maxi, s);
  Mmax = Mmin = 0; /* indices of extremal points */
  for (i=1; i<n; i++)
    {
      dotprod (s, S[i], v);
      if (coeff_gt (s, *maxi))
        {
          coeff_set (*maxi, s);
          Mmax = i;
        }
      else if (coeff_lt (s, *mini))
        {
          coeff_set (*mini, s);
          Mmin = i;
        }
    }
  SubPoint (M, S[Mmax], S[Mmin]);
  coeff_clear (s);
}

/* puts in (A,B,C) the new facet */
void
find_cutting_facet (point_t A, point_t B, point_t C, point_t M)
{
  int da, db, dc;
  point_t O;

  InitPoint (O);
  coeff_set_ui (O->x, 0);
  coeff_set_ui (O->y, 0);
  coeff_set_ui (O->z, 0);
  da = - SignDelta(A, O, M);
  db = - SignDelta(B, O, M);
  dc = - SignDelta(C, O, M);
  if (da >= 0 && db <= 0)
    SetPoint (C, M);
  else if (db >= 0 && dc <= 0)
    SetPoint (A, M);
  else if (dc >= 0 && da <= 0)
    SetPoint (B, M);
  else
    {
      fprintf (stderr, "cannot find cutting facet\n");
      exit (1);
    }
  ClearPoint (O);
}

/* return 0 if FAIL, otherwise number of loops */
int
get_separating_facet (point_t *S, unsigned long n)
{
  point_t A, B, C, v, M;
  unsigned long loop;
  TYPE sp, sp_ref, mini[1], maxi[1];

  InitPoint (A);
  InitPoint (B);
  InitPoint (C);
  InitPoint (M);
  InitPoint (v);
  coeff_init (sp_ref);
  coeff_init (sp);
  coeff_init (*mini);
  coeff_init (*maxi);
  if (find_starting_facet (A, B, C, S, n))
    {
      fprintf (stderr, "Error in find_starting_facet\n");
      exit (1);
    }
  for (loop=1; ; loop++)
    {
      facet_to_normal (v, A, B, C);
      if (coeff_sgn (v->z) < 0)
        {
          coeff_neg (v->x, v->x);
          coeff_neg (v->y, v->y);
          coeff_neg (v->z, v->z);
          SwapPoint (A, B);
        }
      dotprod (sp_ref, A, v);
      if (coeff_ge (sp_ref, v->z))
        {
          printf ("FAIL\n");
          return 0;
        }
      find_minmax_scalprod (M, mini, maxi, S, n, v);
      coeff_sub (sp, *maxi, *mini);
      if (coeff_le (sp, sp_ref))
        {
#ifdef VERBOSE
          printf ("mu:="); print_coeff (*mini);
          printf (": P:="); PrintPoint (v); printf (":\n");
          printf ("loop=%u\n", loop);
#endif
          return loop;
        }
      find_cutting_facet (A, B, C, M);
#ifdef DOUBLE
      if (loop > 100)
        return -1;
#endif
    }
  ClearPoint (A);
  ClearPoint (B);
  ClearPoint (C);
  ClearPoint (M);
  ClearPoint (v);
  coeff_clear (sp_ref);
  coeff_clear (sp);
  coeff_clear (*mini);
  coeff_clear (*maxi);
}

void
mpz_randoms (mpz_t x, mpz_t n, int neg)
{
  mpz_random (x, mpz_size (n));
  mpz_mod (x, x, n);
  if (neg)
    mpz_neg (x, x);
}

void
randplan (point_t *S, unsigned long m)
{
  unsigned long i;
  mpz_t n, a, b, c, mu, x, y, z, r;
  mpz_init (n);
  mpz_init (a);
  mpz_init (b);
  mpz_init (c);
  mpz_init (mu);
  mpz_init (x);
  mpz_init (y);
  mpz_init (z);
  mpz_init (r);

#if 0
  mpz_set_ui (n, m);
  mpz_mul_ui (n, n, 100);
#else
  mpz_set_ui (n, 57250);
#endif
  do
    {
      mpz_randoms (a, n, rand () % 2);
      mpz_randoms (b, n, rand () % 2);
      mpz_randoms (c, n, 0);
    }
  while (mpz_cmp_ui (c, 0) == 0);
  mpz_randoms (mu, n, rand () % 2);
  // printf("%d %d %d %d\n", mpz_get_si (a), mpz_get_si (b), mpz_get_si (c), mpz_get_si (mu));

  /* print number of points */
  /* generate m points on (a,b,c,mu) */
  for (i=0; i<m; i++)
    {
      mpz_randoms (x, n, rand () % 2);
      mpz_randoms (y, n, rand () % 2);
      /* we have mu, a, x, b, y <= n in absolute value
         thus |z| <= 2n2+n */
      mpz_set (z, mu);
      mpz_submul (z, a, x);
      mpz_submul (z, b, y);
      mpz_fdiv_q (z, z, c);
      coeff_set_z (S[i]->x, x);
      coeff_set_z (S[i]->y, y);
      coeff_set_z (S[i]->z, z);
    }
  mpz_clear (n);
  mpz_clear (a);
  mpz_clear (b);
  mpz_clear (c);
  mpz_clear (mu);
  mpz_clear (x);
  mpz_clear (y);
  mpz_clear (z);
  mpz_clear (r);
}

int
cputime ()
{
  struct rusage rus;

  getrusage (0, &rus);
  return rus.ru_utime.tv_sec * 1000 + rus.ru_utime.tv_usec / 1000;
}

/* reco3d N [k] where n is the number of points, k is the number of tests.
*/
int
main (int argc, char *argv[])
{
  unsigned long N, i, k, j;
  long a, b, c;
  point_t *S;
  unsigned long loops, maxloops=0, totloops=0, errors=0;
  int st, maxtime=0, tottime=0;

#ifdef GMPZ
  printf ("Using GNU MP integers\n");
#else
  printf ("Using double-precision floating-point numbers\n");
#endif  

  N = atoi (argv[1]);
  k = (argc > 2) ? atoi (argv[2]) : 1;

  S = (point_t*) malloc (N * sizeof (point_t));
  for (i=0; i<N; i++)
    {
      coeff_init (S[i]->x);
      coeff_init (S[i]->y);
      coeff_init (S[i]->z);
    }

  for (i=0; i<k; i++)
    {
      randplan (S, N);
#ifdef VERBOSE
      printf ("S:=[\n");
      for (j=0; j<N; j++)
        {
          PrintPoint (S[j]);
          if (j + 1 < N)
            printf (",");
          printf ("\n");
        }
      printf ("]:\n");
#endif
      st = cputime ();
      loops = get_separating_facet (S, N);
      if (loops == -1)
        errors ++;
      else if (loops == 0)
        {
          printf ("i=%u\n", i);
          for (j=0; j<N; j++)
            PrintPoint (S[j]);
        }
      else
        {
          st = cputime () - st;
          if (st > maxtime)
            maxtime = st;
          tottime += st;
          if (loops > maxloops)
            maxloops = loops;
          totloops += loops;
        }
    }

  printf ("max. time =%d aver. time =%f\n", maxtime , (double) tottime  / k);
  printf ("max. loops=%u aver. loops=%f\n", maxloops, (double) totloops / k);
  if (errors)
    printf ("%u error(s)\n", errors);

  for (i=0; i<N; i++)
    {
      coeff_clear (S[i]->x);
      coeff_clear (S[i]->y);
      coeff_clear (S[i]->z);
    }

  free (S);
  return 0;
}