#ifndef jsmath_H
#define jsmath_H

#include <stdint.h>

/* J. Stolfi's miscellaneous general-purpose definitions */
/* Last edited on 2009-01-04 13:51:14 by stolfi */

#ifndef INF
#define INF INFINITY
#endif

int64_t ipow(int64_t x, unsigned int y);
  /* Returns {x^y}. Does not check for overflow. */

int64_t imod(int64_t x, int64_t y);
  /* Mathematical remainder of {x ÷ y}, always in {0..y-1}.
    Fails with divide-by-zero if {y} is zero or negative. */

int64_t ifloor(int64_t x, int64_t y);
  /* Rounds {x} down to the nearest multiple of {y}.
    Fails with divide-by-zero if {y} is zero or negative.
    Beware of overflow. */

int64_t iceil(int64_t x, int64_t y);
  /* Rounds {x} up to the nearest multiple of {y}.
    Fails with divide-by-zero if {y} is zero or negative.
    Beware of overflow. */

uint64_t gcd(uint64_t a, uint64_t b);
  /* Greatest common divisor of {a} and {b}.
    If both arguments are 0, returns 0. */
  
uint64_t lcm(uint64_t a, uint64_t b);
  /* Least common multiple of {a} and {b}.
    If either argument is 0, returns 0. */

int64_t imin(int64_t x, int64_t y);
int64_t imax(int64_t x, int64_t y);
  /* Returns maximum/minimum of {x,y}. */

uint64_t umin(uint64_t x, uint64_t y);
uint64_t umax(uint64_t x, uint64_t y);
  /* Returns maximum/minimum of {x,y}. */
  
void uint64_mul(uint64_t x, uint64_t y, uint64_t Z[]);
  /* Multiplies two 64-bit unsigned integers {x,y} to give a
    128-bit unsiged int.  The result is stored into {Z[0],Z[1]}
    and should be interpreted as the integer {Z[0] + Z[1]*2^64}. */

unsigned int digits(uint64_t x);
  /* Number of decimal digits needed to write {x}. 
    Namely, {1} if {x==0}, else {floor(log10(x)) + 1} */

double rel_diff(double x, double y);
  /* Returns the relative difference between {x} and {y}, namely
    {(x-y)/sqrt(x^2 + y^2)/2}. This number is zero when {x == y}, is
    0.5 when only one of them is zero, and has extremal values -1
    and +1 when {x} and {y} are equal and opposite).  It is {NaN}
    when both {x} and {y} are zero. */

double abs_rel_diff(double x, double y, double abs_tol, double rel_tol);
  /* Returns the difference {d = x - y} divided by 
    {D = hypot(abs_tol, rel_tol*hypot(x,y)/sqrt(2))}.
    
    If {fabs(x),fabs(y)} are both very small, the denominator {D} is
    approximately {abs_tol}. In that case, if the result {d/D} is 1.0 in
    absolute value, then {fabs(d)} is approximately {abs_tol}.
    
    If either of {fabs(x),fabs(y)} is large compared to {abs_tol}, {D}
    is between {0.7*M} and {M}, where {M = max{fabs(x),fabs(y)}};
    thus, if the result is less that 1.0 in absolute value, then
    {fabs(d) < rel_tol*fabs(x)} and {fabs(d) < rel_tol*fabs(y)}; if
    the result is greater than 1 in absolute value, then {fabs(d) >
    0.7*rel_tol*fabs(x)} and {fabs(d) > 0.7*rel_tol*fabs(y)}. */

#endif