//===========================================================================//
// This file is part of the DIP Solver Framework.                            //
//                                                                           //
// DIP is distributed under the Eclipse Public License as part of the        //
// COIN-OR repository (http://www.coin-or.org).                              //
//                                                                           //
// Author: Matthew Galati, SAS Institute Inc. (matthew.galati@sas.com)       //
//                                                                           //
// Conceptual Design: Matthew Galati, SAS Institute Inc.                     //
//                    Ted Ralphs, Lehigh University                          //
//                                                                           //
// Copyright (C) 2002-2015, Lehigh University, Matthew Galati, Ted Ralphs    //
// All Rights Reserved.                                                      //
//===========================================================================//

//===========================================================================//
#include "UtilMacros.h"
#define SPACES " \t\r\n"
//===========================================================================//

// =========================================================================
// Graph Macros
// =========================================================================
// ------------------------------------------------------------------------- //
std::pair<int, int> UtilBothEndsU(const int index)
{
   int i = (int)(floor(sqrt(1 + 8.0 * index) / 2 + .500000001));
   return std::make_pair( i, index - (i * (i - 1) / 2) );
}

// =========================================================================
// Other Macros
// =========================================================================

// ------------------------------------------------------------------------- //
int UtilScaleDblToIntArr(const int      arrLen,
                         const double* arrDbl,
                         int*           arrInt,
                         const double   oneDbl,
                         int*           oneInt,
                         const double   epstol)
{
   //---
   //--- A very simple function to scale an array of doubles to integers.
   //---    Note: epstol denotes the preferred accuracy,
   //---    so, we will scale by 1.0/epstol, unless something smaller works.
   //--- It constructs the scaled array and returns the scale factor.
   //---
   //--- It can also scale oneDbl to oneInt wrt to the array (e.g..,
   //--- the rhs of row). If oneInt == NULL, then this part is skipped.
   //---
   int      i, scaleFactor = 1, n_aFrac = 0, factorTooBig = 0;
   double* arrFrac = NULL;
   double   fractionalPart;
   double   oneOverEps = 1.0 / epstol;
   arrFrac = new double[arrLen + 1];
   assert(arrFrac);

   for (i = 0; i < arrLen; i++) {
      fractionalPart = UtilFracPart(arrDbl[i]);

      if (!UtilIsZero(fractionalPart)) {
         fractionalPart     *= oneOverEps;
         arrFrac[n_aFrac++]  = static_cast<int>(round(fractionalPart))
                               * (double)epstol;
      }
   }

   if (oneInt) {
      fractionalPart = UtilFracPart(oneDbl);

      if (!UtilIsZero(fractionalPart)) {
         fractionalPart     *= oneOverEps;
         arrFrac[n_aFrac++]  = static_cast<int>(round(fractionalPart))
                               * (double)epstol;
      }
   }

   for (i = 0; i < n_aFrac; i++) {
      assert(arrFrac[i] < (INT_MAX / scaleFactor));
      arrFrac[i] *= scaleFactor;

      while (!UtilIsZero(UtilFracPart(arrFrac[i]))) {
         scaleFactor *= 10;

         if (scaleFactor >= oneOverEps) {
            factorTooBig = 1;
            break;
         }

         assert(arrFrac[i] < (INT_MAX / 10));
         arrFrac[i]    *= 10;
         assert(arrFrac[i] >= 0);
      }

      if (factorTooBig) {
         break;
      }
   }

   for (i = 0; i < arrLen; i++) {
      arrInt[i] = static_cast<int>(round(arrDbl[i] * scaleFactor));
   }

   if (oneInt) {
      *oneInt = static_cast<int>(round(oneDbl * scaleFactor));
   }

   UTIL_DELARR(arrFrac);
   return scaleFactor;
}


/*==========================================================================*/
int UtilScaleDblToIntArr(const int      arrLen,
                         const double* arrDbl,
                         int*           arrInt,
                         const double   epstol)
{
   //---
   //--- A very simple function to scale an array of doubles to integers.
   //---    Note: epstol denotes the preferred accuracy,
   //---    so, we will scale by 1.0/epstol, unless something smaller works.
   //--- It constructs the scaled array and returns the scale factor.
   //---
   //--- It can also scale oneDbl to oneInt wrt to the array (e.g..,
   //--- the rhs of row). If oneInt == NULL, then this part is skipped.
   //---
   int      i, scaleFactor = 1, n_aFrac = 0, factorTooBig = 0;
   double* arrFrac = NULL;
   double   fractionalPart;
   double   oneOverEps = 1.0 / epstol;
   //TODO: pass in arrFrac?
   arrFrac = new double[arrLen];
   assert(arrFrac);

   for (i = 0; i < arrLen; i++) {
      fractionalPart = UtilFracPart(arrDbl[i]);

      //printf("arrDbl[%d]=%10.5f fracPart=%6.5f\n",
      //     i, arrDbl[i], fractionalPart);
      if (!UtilIsZero(fractionalPart)) {
         fractionalPart     *= oneOverEps;
         //printf("fracPart is not zero oneOverEps= %10.5f fracPart= %10.5f\n",
         //oneOverEps, fractionalPart);
         arrFrac[n_aFrac++]  = static_cast<int>(round(fractionalPart))
                               * (double)epstol;
         //printf("arrFrac[%d] = %10.5f\n", (n_aFrac-1), arrFrac[n_aFrac-1]);
      }
   }

   for (i = 0; i < n_aFrac; i++) {
      assert(arrFrac[i] < (INT_MAX / scaleFactor));
      arrFrac[i] *= scaleFactor;

      while (!UtilIsZero(UtilFracPart(arrFrac[i]))) {
         scaleFactor *= 10;

         if (scaleFactor >= oneOverEps) {
            factorTooBig = 1;
            break;
         }

         assert(arrFrac[i] < (INT_MAX / 10));
         arrFrac[i]    *= 10;
         assert(arrFrac[i] >= 0);
      }

      if (factorTooBig) {
         break;
      }
   }

   //---
   //--- must be careful not to trunc here
   //---   so, we want to round
   //---
   for (i = 0; i < arrLen; i++) {
      arrInt[i] = static_cast<int>(round(arrDbl[i] * scaleFactor));
   }

   UTIL_DELARR(arrFrac);
   return scaleFactor;
}

#if 0
/*==========================================================================*/
//---
//--- taken from Concorde's integerize_vector in LOCALCUT/first.c
//---
/*==========================================================================*/
#define INTEGERIZE_MUL (16*9*5*7*11*13*17)
#define CC_SWAP(a,b,t) (((t)=(a)),((a)=(b)),((b)=(t)))

/*==========================================================================*/
static int CCutil_our_gcd (int a, int b)
{
   int c;

   if (a < 0) {
      a = -a;
   }

   if (b < 0) {
      b = -b;
   }

   if (a > b) {
      CC_SWAP (a, b, c);
   }

   while (a) {
      c = b % a;
      b = a;
      a = c;
   }

   return b;
}
/*==========================================================================*/
static void integerize_check (int count, double* dvec, int* ivec)
{
   double error;
   double scale;
   int i;
   scale = 1.0;

   for (i = 0; i < count; i++) {
      if (ivec[i]) {
         scale = fabs(dvec[i] / ivec[i]);
         break;
      }
   }

   error = 0.0;

   for (i = 0; i < count; i++) {
      error += fabs(dvec[i] / scale - ivec[i]);
   }

   if (error > 0.001) {
      printf ("bad integerize error %f\n", error);
      printf ("from:");

      for (i = 0; i < count; i++) {
         printf (" %f", dvec[i]);
      }

      printf ("\n");
      printf ("to:");

      for (i = 0; i < count; i++) {
         printf (" %d", ivec[i]);
      }

      printf ("\n");
      fflush (stdout);
   }
}

/*==========================================================================*/
static void cfrac_approx (double d, int* p_num, int* p_den)
{
   int sgn = 1;
   int a;
   int g2 = 0;
   int g1 = 1;
   int g = 0;
   int h2 = 1;
   int h1 = 0;
   int h = 1;
   int bestg = 0;
   int besth = 1;
   double besterr = 1.0;
   int lim;

   if (d < 0.0) {
      d = -d;
      sgn = -1;
   }

   for (;;) {
      a = (int) d;

      if (a) {
         lim = INT_MAX / a;
      } else {
         lim = INT_MAX;
      }

      if (g1 > lim || h1 > lim ||
            g2 > INT_MAX - g1 * a || h2 > INT_MAX - h1 * a) {
         break;
      }

      g = a * g1 + g2;
      h = a * h1 + h2;
      g2 = g1;
      h2 = h1;
      g1 = g;
      h1 = h;

      if (d - (double) a < besterr) {
         bestg = g;
         besth = h;
         besterr = d - (double) a;

         if (besterr < 0.0001) {
            break;
         }
      }

      d = 1 / (d - a);
   }

   *p_num = bestg * sgn;
   *p_den = besth;
}

/*==========================================================================*/
static void integerize_vector2 (int count, double* dvec, int* ivec)
{
   int i;
   double scale;
   int divv;
#if 0
   printf ("integerize(2):");

   for (i = 0; i < count; i++) {
      printf (" %f", dvec[i]);
   }

   printf ("\n");
#endif
   scale = fabs(dvec[0]);

   for (i = 1; i < count; i++) {
      if (fabs(dvec[i]) > scale) {
         scale = fabs(dvec[i]);
      }
   }

   if (scale == 0.0) {
      for (i = 0; i < count; i++) {
         ivec[i] = 0;
      }

      return;
   }

   for (i = 0; i < count; i++) {
      if (dvec[i] < 0.0) {
         ivec[i] = -((int) (-dvec[i] * INTEGERIZE_MUL / scale + 0.5));
      } else {
         ivec[i] = ((int) (dvec[i] * INTEGERIZE_MUL / scale + 0.5));
      }
   }

   divv = ivec[0];

   for (i = 1; i < count; i++) {
      divv = CCutil_our_gcd (divv, ivec[i]);
   }

   if (divv > 1) {
      for (i = 0; i < count; i++) {
         ivec[i] /= divv;
      }
   }

#if 0
   printf ("integerized");

   for (i = 0; i < count; i++) {
      printf (" %d", ivec[i]);
   }

   printf ("\n");
   integerize_check (count, dvec, ivec);
#endif
}

/*==========================================================================*/
static void integerize_vector (int count, double* dvec, int* ivec)
{
   int i;
   int num;
   int den;
   int g;
   int scale = 1;
   double dmax;
#if 0
   printf ("integerize:");

   for (i = 0; i < count; i++) {
      printf (" %f", dvec[i]);
   }

   printf ("\n");
#endif
   dmax = 0.0;

   for (i = 0; i < count; i++) {
      if (fabs(dvec[i]) > dmax) {
         dmax = fabs(dvec[i]);
      }
   }

   dmax = 1.0 / dmax;

   for (i = 0; i < count; i++) {
      cfrac_approx (dvec[i]*dmax, &num, &den);
      g = CCutil_our_gcd (den, scale);
      scale /= g;

      if (INT_MAX / abs(scale) < abs(den)) {
         goto FAILURE;
      }

      scale *= den;
   }

   for (i = 0; i < count; i++) {
      cfrac_approx (dvec[i]*dmax, &num, &den);
      den = scale / den;

      if (INTEGERIZE_MUL / abs(den) < abs(num)) {
         goto FAILURE;
      }

      ivec[i] = num * den;
   }

#if 0
   printf ("integerized");

   for (i = 0; i < count; i++) {
      printf (" %d", ivec[i]);
   }

   printf ("\n");
   integerize_check (count, dvec, ivec);
#endif
   return;
FAILURE:
   printf ("Overflow in integerize_vector.  Calling integerize_vector2\n");
   integerize_vector2 (count, dvec, ivec);
}



/*==========================================================================*/
//---
//--- taken from concorde integerize_vector in LOCALCUT/first.c
//---
void UtilScaleDblToIntArrCC(const int      arrLen,
                            double* arrDbl,
                            int*           arrInt)
{
   integerize_vector(arrLen, arrDbl, arrInt);
}
#endif
