/* See frb_proc.h */
/* Last edited on 2005-01-16 14:20:38 by stolfi */
/* Copyright © 2005 Universidade Estadual de Campinas. See note at end of file. */

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

#include <js.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;
      }
  }


/* COPYRIGHT AND AUTHORSHIP NOTICE

  Copyright © 2005 Universidade Estadual de Campinas (UNICAMP).
  Created by Helena C. G. Leitão and Jorge Stolfi in 1995--2005.
  
  This source file can be freely distributed, used, and modified,
  provided that this copyright and authorship notice is preserved in
  all copies, and that any modified versions of this file are clearly
  marked as such.
  
  This software has NO WARRANTY of correctness or applicability for
  any purpose. Neither the authors nor their employers shall be held
  responsible for any losses or damages that may result from its use.

  END OF NOTICE */