/* See frb_proc.h */
/* Last edited on 2005-01-01 11:41:46 by stolfi */

#include <frb_types.h>
#include <frb_proc.h>

#include <affirm.h>

#include <vec.h>
#include <stdlib.h>
#include <math.h>

double frb_dmax ( double A, double B)
  { return ((A) >= (B) ? (A) : (B)); }

double frb_dmin ( double A, double B )
  { return ((A) >= (B) ? (B) : (A)); }

int frb_imax ( int A, int B)
  { return ((A) >= (B) ? (A) : (B)); }

int frb_imin ( int A, int B )
  { return ((A) >= (B) ? (B) : (A)); }

int frb_round ( double x )
  { return (x < 0.0 ? (int)(x - 0.5) : (int)(x + 0.5)); }

int frb_imod ( int A, int B )
  { demand(B > 0, "nonpos divisor"); 
    int R = A % B;
    while (R < 0) { R += B; }
    return R;
  }
    
double  frb_adjust_unit
  ( double givenUnit, 
    double dev,
    double big 
  )
  {
    double eps, minUnit, newUnit;
    /* Compute desirable unit {eps}: */
    if (big == 0)
      { eps = 1; }
    else if (dev == 0)
      { eps = frb_dmax(big/1.0e-6, 1.0e-10); }
    else
      { eps = frb_dmax(dev/1.0e-6, 1.0e-10); }
    /* Compute minimum unit to avoid overflow: */
    minUnit = big/1.0e+9;
    /* Decide what to do with {givenUnit}: */
    if (givenUnit < minUnit)
      { /* Must increase {givenUnit} to avoid overflow: */
        newUnit = 1;
        while (newUnit/10 >= minUnit){ newUnit = newUnit / 10; }
        while (newUnit < minUnit){ newUnit = newUnit * 10; }
        fprintf(stderr, 
          "warning: unit increased from %g to %g to avoid overflow.\n",
          givenUnit, newUnit
        );
      }
    else if (givenUnit > eps)
      { /* Given unit is too big, try to reduce it. */
        /* Compute smallest power of 10 {p >= minUnit} and {p > eps/10}: */
        newUnit = 1;
        while ((newUnit > eps) && (newUnit/10 >= minUnit))
          { newUnit = newUnit / 10; }
        while ((newUnit <= eps/10 ) && (newUnit < minUnit))
          { newUnit = newUnit * 10; }
        if (newUnit < givenUnit)
          { fprintf(stderr, 
              "warning: reducing unit from %g to %g to reduce quantization errors.\n",
              givenUnit, newUnit
            );
          }
      }
    else
      { newUnit = givenUnit; }
    /* In any case, check whether result is OK: */
    if (newUnit > eps)
      { fprintf(stderr, 
          "warning: possible loss of information dev = %g unit = %g\n",
          dev, newUnit
        );
      }
    return newUnit;
  }

void frb_reduce_to_period
  ( double *t, /* IO */  
    double tLo,
    double tHi, 
    int *k  /* (OUT) */
  )
  { if ((*t >= tLo) && (*t < tHi))
      { *k = 0; }
    else
      { double T = tHi - tLo;
        affirm(T > 0, "neg period");
        double kk = floor(((*t) - tLo)/T);
        double tt = (*t) - ((double)kk)*T;
        while (tt < tLo){ tt += T; kk++; }
        while (tt >= tHi){ tt -= T; kk--; }
        *t = tt; *k = kk;
      }
  }