CImg.h

      const unsigned int v = u|(1U<<(8*sizeof(unsigned int)-1)); // set sign bit to 1.
      // use memcpy instead of simple assignment to avoid undesired optimizations by C++-compiler.
      std::memcpy(&f,&v,sizeof(float));
      return f;
    }

    //! Return the value of a system timer, with a millisecond precision.
    /**
       \note The timer does not necessarily starts from \c 0.
    **/
    inline cimg_ulong time() {
#if cimg_OS==1
      struct timeval st_time;
      gettimeofday(&st_time,0);
      return (cimg_ulong)(st_time.tv_usec/1000 + st_time.tv_sec*1000);
#elif cimg_OS==2
      SYSTEMTIME st_time;
      GetLocalTime(&st_time);
      return (cimg_ulong)(st_time.wMilliseconds + 1000*(st_time.wSecond + 60*(st_time.wMinute + 60*st_time.wHour)));
#else
      return 0;
#endif
    }

    // Implement a tic/toc mechanism to display elapsed time of algorithms.
    inline cimg_ulong tictoc(const bool is_tic);

    //! Start tic/toc timer for time measurement between code instructions.
    /**
       \return Current value of the timer (same value as time()).
    **/
    inline cimg_ulong tic() {
      return cimg::tictoc(true);
    }

    //! End tic/toc timer and displays elapsed time from last call to tic().
    /**
       \return Time elapsed (in ms) since last call to tic().
    **/
    inline cimg_ulong toc() {
      return cimg::tictoc(false);
    }

    //! Sleep for a given numbers of milliseconds.
    /**
       \param milliseconds Number of milliseconds to wait for.
       \note This function frees the CPU ressources during the sleeping time.
       It can be used to temporize your program properly, without wasting CPU time.
    **/
    inline void sleep(const unsigned int milliseconds) {
#if cimg_OS==1
      struct timespec tv;
      tv.tv_sec = milliseconds/1000;
      tv.tv_nsec = (milliseconds%1000)*1000000;
      nanosleep(&tv,0);
#elif cimg_OS==2
      Sleep(milliseconds);
#else
      cimg::unused(milliseconds);
#endif
    }

    inline unsigned int _wait(const unsigned int milliseconds, cimg_ulong& timer) {
      if (!timer) timer = cimg::time();
      const cimg_ulong current_time = cimg::time();
      if (current_time>=timer + milliseconds) { timer = current_time; return 0; }
      const unsigned int time_diff = (unsigned int)(timer + milliseconds - current_time);
      timer = current_time + time_diff;
      cimg::sleep(time_diff);
      return time_diff;
    }

    //! Wait for a given number of milliseconds since the last call to wait().
    /**
       \param milliseconds Number of milliseconds to wait for.
       \return Number of milliseconds elapsed since the last call to wait().
       \note Same as sleep() with a waiting time computed with regard to the last call
       of wait(). It may be used to temporize your program properly, without wasting CPU time.
    **/
    inline cimg_long wait(const unsigned int milliseconds) {
      cimg::mutex(3);
      static cimg_ulong timer = 0;
      if (!timer) timer = cimg::time();
      cimg::mutex(3,0);
      return _wait(milliseconds,timer);
    }

    // Random number generators.
    // CImg may use its own Random Number Generator (RNG) if configuration macro 'cimg_use_rng' is set.
    // Use it for instance when you have to deal with concurrent threads trying to call std::srand()
    // at the same time!
#ifdef cimg_use_rng

#include <stdint.h>

    // Use a custom RNG.
    inline unsigned int _rand(const unsigned int seed=0, const bool set_seed=false) {
      static cimg_ulong next = 0xB16B00B5;
      cimg::mutex(4);
      if (set_seed) next = (cimg_ulong)seed;
      next = next*1103515245 + 12345U;
      cimg::mutex(4,0);
      return (unsigned int)(next&0xFFFFFFU);
    }

    inline void srand() {
      const unsigned int t = (unsigned int)cimg::time();
#if cimg_OS==1
      cimg::_rand(t + (unsigned int)getpid(),true);
#elif cimg_OS==2
      cimg::_rand(t + (unsigned int)_getpid(),true);
#else
      cimg::_rand(t,true);
#endif
    }

    inline void srand(const unsigned int seed) {
      _rand(seed,true);
    }

    inline double rand(const double val_min, const double val_max) {
      const double val = cimg::_rand()/16777215.;
      return val_min + (val_max - val_min)*val;
    }

#else

    // Use the system RNG.
    inline void srand() {
      const unsigned int t = (unsigned int)cimg::time();
#if cimg_OS==1 || defined(__BORLANDC__)
      std::srand(t + (unsigned int)getpid());
#elif cimg_OS==2
      std::srand(t + (unsigned int)_getpid());
#else
      std::srand(t);
#endif
    }

    inline void srand(const unsigned int seed) {
      std::srand(seed);
    }

    //! Return a random variable uniformely distributed between [val_min,val_max].
    /**
    **/
    inline double rand(const double val_min, const double val_max) {
      const double val = (double)std::rand()/RAND_MAX;
      return val_min + (val_max - val_min)*val;
    }
#endif

    //! Return a random variable uniformely distributed between [0,val_max].
    /**
     **/
    inline double rand(const double val_max=1) {
      return cimg::rand(0,val_max);
    }

    //! Return a random variable following a gaussian distribution and a standard deviation of 1.
    /**
    **/
    inline double grand() {
      double x1, w;
      do {
        const double x2 = cimg::rand(-1,1);
        x1 = cimg::rand(-1,1);
        w = x1*x1 + x2*x2;
      } while (w<=0 || w>=1.0);
      return x1*std::sqrt((-2*std::log(w))/w);
    }

    //! Return a random variable following a Poisson distribution of parameter z.
    /**
    **/
    inline unsigned int prand(const double z) {
      if (z<=1.0e-10) return 0;
      if (z>100) return (unsigned int)((std::sqrt(z) * cimg::grand()) + z);
      unsigned int k = 0;
      const double y = std::exp(-z);
      for (double s = 1.0; s>=y; ++k) s*=cimg::rand();
      return k - 1;
    }

    //! Cut (i.e. clamp) value in specified interval.
    template<typename T, typename t>
    inline T cut(const T& val, const t& val_min, const t& val_max) {
      return val<val_min?(T)val_min:val>val_max?(T)val_max:val;
    }

    //! Bitwise-rotate value on the left.
    template<typename T>
    inline T rol(const T& a, const unsigned int n=1) {
      return n?(T)((a<<n)|(a>>((sizeof(T)<<3) - n))):a;
    }

    inline float rol(const float a, const unsigned int n=1) {
      return (float)rol((int)a,n);
    }

    inline double rol(const double a, const unsigned int n=1) {
      return (double)rol((cimg_long)a,n);
    }

    inline double rol(const long double a, const unsigned int n=1) {
      return (double)rol((cimg_long)a,n);
    }

#ifdef cimg_use_half
    inline half rol(const half a, const unsigned int n=1) {
      return (half)rol((int)a,n);
    }
#endif

    //! Bitwise-rotate value on the right.
    template<typename T>
    inline T ror(const T& a, const unsigned int n=1) {
      return n?(T)((a>>n)|(a<<((sizeof(T)<<3) - n))):a;
    }

    inline float ror(const float a, const unsigned int n=1) {
      return (float)ror((int)a,n);
    }

    inline double ror(const double a, const unsigned int n=1) {
      return (double)ror((cimg_long)a,n);
    }

    inline double ror(const long double a, const unsigned int n=1) {
      return (double)ror((cimg_long)a,n);
    }

#ifdef cimg_use_half
    inline half ror(const half a, const unsigned int n=1) {
      return (half)ror((int)a,n);
    }
#endif

    //! Return absolute value of a value.
    template<typename T>
    inline T abs(const T& a) {
      return a>=0?a:-a;
    }
    inline bool abs(const bool a) {
      return a;
    }
    inline int abs(const unsigned char a) {
      return (int)a;
    }
    inline int abs(const unsigned short a) {
      return (int)a;
    }
    inline int abs(const unsigned int a) {
      return (int)a;
    }
    inline int abs(const int a) {
      return std::abs(a);
    }
    inline cimg_int64 abs(const cimg_uint64 a) {
      return (cimg_int64)a;
    }
    inline double abs(const double a) {
      return std::fabs(a);
    }
    inline float abs(const float a) {
      return (float)std::fabs((double)a);
    }

    //! Return square of a value.
    template<typename T>
    inline T sqr(const T& val) {
      return val*val;
    }

    //! Return <tt>1 + log_10(x)</tt> of a value \c x.
    inline int xln(const int x) {
      return x>0?(int)(1 + std::log10((double)x)):1;
    }

    //! Return the minimum between three values.
    template<typename t>
    inline t min(const t& a, const t& b, const t& c) {
      return std::min(std::min(a,b),c);
    }

    //! Return the minimum between four values.
    template<typename t>
    inline t min(const t& a, const t& b, const t& c, const t& d) {
      return std::min(std::min(a,b),std::min(c,d));
    }

    //! Return the maximum between three values.
    template<typename t>
    inline t max(const t& a, const t& b, const t& c) {
      return std::max(std::max(a,b),c);
    }

    //! Return the maximum between four values.
    template<typename t>
    inline t max(const t& a, const t& b, const t& c, const t& d) {
      return std::max(std::max(a,b),std::max(c,d));
    }

    //! Return the sign of a value.
    template<typename T>
    inline T sign(const T& x) {
      return (T)(x<0?-1:x>0);
    }

    //! Return the nearest power of 2 higher than given value.
    template<typename T>
    inline cimg_ulong nearest_pow2(const T& x) {
      cimg_ulong i = 1;
      while (x>i) i<<=1;
      return i;
    }

    //! Return the sinc of a given value.
    inline double sinc(const double x) {
      return x?std::sin(x)/x:1;
    }

    //! Return the modulo of a value.
    /**
       \param x Input value.
       \param m Modulo value.
       \note This modulo function accepts negative and floating-points modulo numbers, as well as variables of any type.
    **/
    template<typename T>
    inline T mod(const T& x, const T& m) {
      const double dx = (double)x, dm = (double)m;
      return (T)(dx - dm * std::floor(dx / dm));
    }
    inline int mod(const bool x, const bool m) {
      return m?(x?1:0):0;
    }
    inline int mod(const unsigned char x, const unsigned char m) {
      return x%m;
    }
    inline int mod(const char x, const char m) {
#if defined(CHAR_MAX) && CHAR_MAX==255
      return x%m;
#else
      return x>=0?x%m:(x%m?m + x%m:0);
#endif
    }
    inline int mod(const unsigned short x, const unsigned short m) {
      return x%m;
    }
    inline int mod(const short x, const short m) {
      return x>=0?x%m:(x%m?m + x%m:0);
    }
    inline int mod(const unsigned int x, const unsigned int m) {
      return (int)(x%m);
    }
    inline int mod(const int x, const int m) {
      return x>=0?x%m:(x%m?m + x%m:0);
    }
    inline cimg_int64 mod(const cimg_uint64 x, const cimg_uint64 m) {
      return x%m;
    }
    inline cimg_int64 mod(const cimg_int64 x, const cimg_int64 m) {
      return x>=0?x%m:(x%m?m + x%m:0);
    }

    //! Return the min-mod of two values.
    /**
       \note <i>minmod(\p a,\p b)</i> is defined to be:
       - <i>minmod(\p a,\p b) = min(\p a,\p b)</i>, if \p a and \p b have the same sign.
       - <i>minmod(\p a,\p b) = 0</i>, if \p a and \p b have different signs.
    **/
    template<typename T>
    inline T minmod(const T& a, const T& b) {
      return a*b<=0?0:(a>0?(a<b?a:b):(a<b?b:a));
    }

    //! Return base-2 logarithm of a value.
    inline double log2(const double x) {
      const double base = std::log(2.0);
      return std::log(x)/base;
    }

    template<typename T>
    inline T round(const T& x) {
      return (T)std::floor((_cimg_Tfloat)x + 0.5f);
    }

    //! Return rounded value.
    /**
       \param x Value to be rounded.
       \param y Rounding precision.
       \param rounding_type Type of rounding operation (\c 0 = nearest, \c -1 = backward, \c 1 = forward).
       \return Rounded value, having the same type as input value \c x.
    **/
    template<typename T>
    inline T round(const T& x, const double y, const int rounding_type=0) {
      if (y<=0) return x;
      if (y==1) switch (rounding_type) {
        case 0 : return round(x);
        case 1 : return (T)std::ceil((_cimg_Tfloat)x);
        default : return (T)std::floor((_cimg_Tfloat)x);
        }
      const double sx = (double)x/y, floor = std::floor(sx), delta =  sx - floor;
      return (T)(y*(rounding_type<0?floor:rounding_type>0?std::ceil(sx):delta<0.5?floor:std::ceil(sx)));
    }

    //! Return x^(1/3).
    template<typename T>
    inline double cbrt(const T& x) {
#if cimg_use_cpp11==1
      return std::cbrt(x);
#else
      return x>=0?std::pow((double)x,1.0/3):-std::pow(-(double)x,1.0/3);
#endif
    }

    // Code to compute fast median from 2,3,5,7,9,13,25 and 49 values.
    // (contribution by RawTherapee: http://rawtherapee.com/).
    template<typename T>
    inline T median(T val0, T val1) {
      return (val0 + val1)/2;
    }

    template<typename T>
    inline T median(T val0, T val1, T val2) {
      return std::max(std::min(val0,val1),std::min(val2,std::max(val0,val1)));
    }

    template<typename T>
    inline T median(T val0, T val1, T val2, T val3, T val4) {
      T tmp = std::min(val0,val1);
      val1 = std::max(val0,val1); val0 = tmp; tmp = std::min(val3,val4); val4 = std::max(val3,val4);
      val3 = std::max(val0,tmp);  val1 = std::min(val1,val4); tmp = std::min(val1,val2); val2 = std::max(val1,val2);
      val1 = tmp; tmp = std::min(val2,val3);
      return std::max(val1,tmp);
    }

    template<typename T>
    inline T median(T val0, T val1, T val2, T val3, T val4, T val5, T val6) {
      T tmp = std::min(val0,val5);
      val5 = std::max(val0,val5); val0 = tmp; tmp = std::min(val0,val3); val3 = std::max(val0,val3); val0 = tmp;
      tmp = std::min(val1,val6); val6 = std::max(val1,val6); val1 = tmp; tmp = std::min(val2,val4);
      val4 = std::max(val2,val4); val2 = tmp; val1 = std::max(val0,val1); tmp = std::min(val3,val5);
      val5 = std::max(val3,val5); val3 = tmp; tmp = std::min(val2,val6); val6 = std::max(val2,val6);
      val3 = std::max(tmp,val3); val3 = std::min(val3,val6); tmp = std::min(val4,val5); val4 = std::max(val1,tmp);
      tmp = std::min(val1,tmp); val3 = std::max(tmp,val3);
      return std::min(val3,val4);
    }

    template<typename T>
    inline T median(T val0, T val1, T val2, T val3, T val4, T val5, T val6, T val7, T val8) {
      T tmp = std::min(val1,val2);
      val2 = std::max(val1,val2); val1 = tmp; tmp = std::min(val4,val5);
      val5 = std::max(val4,val5); val4 = tmp; tmp = std::min(val7,val8);
      val8 = std::max(val7,val8); val7 = tmp; tmp = std::min(val0,val1);
      val1 = std::max(val0,val1); val0 = tmp; tmp = std::min(val3,val4);
      val4 = std::max(val3,val4); val3 = tmp; tmp = std::min(val6,val7);
      val7 = std::max(val6,val7); val6 = tmp; tmp = std::min(val1,val2);
      val2 = std::max(val1,val2); val1 = tmp; tmp = std::min(val4,val5);
      val5 = std::max(val4,val5); val4 = tmp; tmp = std::min(val7,val8);
      val8 = std::max(val7,val8); val3 = std::max(val0,val3); val5 = std::min(val5,val8);
      val7 = std::max(val4,tmp); tmp = std::min(val4,tmp); val6 = std::max(val3,val6);
      val4 = std::max(val1,tmp); val2 = std::min(val2,val5); val4 = std::min(val4,val7);
      tmp = std::min(val4,val2); val2 = std::max(val4,val2); val4 = std::max(val6,tmp);
      return std::min(val4,val2);
    }

    template<typename T>
    inline T median(T val0, T val1, T val2, T val3, T val4, T val5, T val6, T val7, T val8, T val9, T val10, T val11,
                    T val12) {
      T tmp = std::min(val1,val7);
      val7 = std::max(val1,val7); val1 = tmp; tmp = std::min(val9,val11); val11 = std::max(val9,val11); val9 = tmp;
      tmp = std::min(val3,val4);  val4 = std::max(val3,val4); val3 = tmp; tmp = std::min(val5,val8);
      val8 = std::max(val5,val8); val5 = tmp; tmp = std::min(val0,val12); val12 = std::max(val0,val12);
      val0 = tmp; tmp = std::min(val2,val6); val6 = std::max(val2,val6); val2 = tmp; tmp = std::min(val0,val1);
      val1 = std::max(val0,val1); val0 = tmp; tmp = std::min(val2,val3); val3 = std::max(val2,val3); val2 = tmp;
      tmp = std::min(val4,val6);  val6 = std::max(val4,val6); val4 = tmp; tmp = std::min(val8,val11);
      val11 = std::max(val8,val11); val8 = tmp; tmp = std::min(val7,val12); val12 = std::max(val7,val12); val7 = tmp;
      tmp = std::min(val5,val9); val9 = std::max(val5,val9); val5 = tmp; tmp = std::min(val0,val2);
      val2 = std::max(val0,val2); val0 = tmp; tmp = std::min(val3,val7); val7 = std::max(val3,val7); val3 = tmp;
      tmp = std::min(val10,val11); val11 = std::max(val10,val11); val10 = tmp; tmp = std::min(val1,val4);
      val4 = std::max(val1,val4); val1 = tmp; tmp = std::min(val6,val12); val12 = std::max(val6,val12); val6 = tmp;
      tmp = std::min(val7,val8); val8 = std::max(val7,val8); val7 = tmp; val11 = std::min(val11,val12);
      tmp = std::min(val4,val9); val9 = std::max(val4,val9); val4 = tmp; tmp = std::min(val6,val10);
      val10 = std::max(val6,val10); val6 = tmp; tmp = std::min(val3,val4); val4 = std::max(val3,val4); val3 = tmp;
      tmp = std::min(val5,val6); val6 = std::max(val5,val6); val5 = tmp; val8 = std::min(val8,val9);
      val10 = std::min(val10,val11); tmp = std::min(val1,val7); val7 = std::max(val1,val7); val1 = tmp;
      tmp = std::min(val2,val6); val6 = std::max(val2,val6); val2 = tmp; val3 = std::max(val1,val3);
      tmp = std::min(val4,val7); val7 = std::max(val4,val7); val4 = tmp; val8 = std::min(val8,val10);
      val5 = std::max(val0,val5); val5 = std::max(val2,val5); tmp = std::min(val6,val8); val8 = std::max(val6,val8);
      val5 = std::max(val3,val5); val7 = std::min(val7,val8); val6 = std::max(val4,tmp); tmp = std::min(val4,tmp);
      val5 = std::max(tmp,val5); val6 = std::min(val6,val7);
      return std::max(val5,val6);
    }

    template<typename T>
    inline T median(T val0, T val1, T val2, T val3, T val4,
                    T val5, T val6, T val7, T val8, T val9,
                    T val10, T val11, T val12, T val13, T val14,
                    T val15, T val16, T val17, T val18, T val19,
                    T val20, T val21, T val22, T val23, T val24) {
      T tmp = std::min(val0,val1);
      val1 = std::max(val0,val1); val0 = tmp; tmp = std::min(val3,val4); val4 = std::max(val3,val4);
      val3 = tmp; tmp = std::min(val2,val4); val4 = std::max(val2,val4); val2 = std::min(tmp,val3);
      val3 = std::max(tmp,val3); tmp = std::min(val6,val7); val7 = std::max(val6,val7); val6 = tmp;
      tmp = std::min(val5,val7); val7 = std::max(val5,val7); val5 = std::min(tmp,val6); val6 = std::max(tmp,val6);
      tmp = std::min(val9,val10); val10 = std::max(val9,val10); val9 = tmp; tmp = std::min(val8,val10);
      val10 = std::max(val8,val10); val8 = std::min(tmp,val9); val9 = std::max(tmp,val9);
      tmp = std::min(val12,val13); val13 = std::max(val12,val13); val12 = tmp; tmp = std::min(val11,val13);
      val13 = std::max(val11,val13); val11 = std::min(tmp,val12); val12 = std::max(tmp,val12);
      tmp = std::min(val15,val16); val16 = std::max(val15,val16); val15 = tmp; tmp = std::min(val14,val16);
      val16 = std::max(val14,val16); val14 = std::min(tmp,val15); val15 = std::max(tmp,val15);
      tmp = std::min(val18,val19); val19 = std::max(val18,val19); val18 = tmp; tmp = std::min(val17,val19);
      val19 = std::max(val17,val19); val17 = std::min(tmp,val18); val18 = std::max(tmp,val18);
      tmp = std::min(val21,val22); val22 = std::max(val21,val22); val21 = tmp; tmp = std::min(val20,val22);
      val22 = std::max(val20,val22); val20 = std::min(tmp,val21); val21 = std::max(tmp,val21);
      tmp = std::min(val23,val24); val24 = std::max(val23,val24); val23 = tmp; tmp = std::min(val2,val5);
      val5 = std::max(val2,val5); val2 = tmp; tmp = std::min(val3,val6); val6 = std::max(val3,val6); val3 = tmp;
      tmp = std::min(val0,val6); val6 = std::max(val0,val6); val0 = std::min(tmp,val3); val3 = std::max(tmp,val3);
      tmp = std::min(val4,val7); val7 = std::max(val4,val7); val4 = tmp; tmp = std::min(val1,val7);
      val7 = std::max(val1,val7); val1 = std::min(tmp,val4); val4 = std::max(tmp,val4); tmp = std::min(val11,val14);
      val14 = std::max(val11,val14); val11 = tmp; tmp = std::min(val8,val14); val14 = std::max(val8,val14);
      val8 = std::min(tmp,val11); val11 = std::max(tmp,val11); tmp = std::min(val12,val15);
      val15 = std::max(val12,val15); val12 = tmp; tmp = std::min(val9,val15); val15 = std::max(val9,val15);
      val9 = std::min(tmp,val12); val12 = std::max(tmp,val12); tmp = std::min(val13,val16);
      val16 = std::max(val13,val16); val13 = tmp; tmp = std::min(val10,val16); val16 = std::max(val10,val16);
      val10 = std::min(tmp,val13); val13 = std::max(tmp,val13); tmp = std::min(val20,val23);
      val23 = std::max(val20,val23); val20 = tmp; tmp = std::min(val17,val23); val23 = std::max(val17,val23);
      val17 = std::min(tmp,val20); val20 = std::max(tmp,val20); tmp = std::min(val21,val24);
      val24 = std::max(val21,val24); val21 = tmp; tmp = std::min(val18,val24); val24 = std::max(val18,val24);
      val18 = std::min(tmp,val21); val21 = std::max(tmp,val21); tmp = std::min(val19,val22);
      val22 = std::max(val19,val22); val19 = tmp; val17 = std::max(val8,val17); tmp = std::min(val9,val18);
      val18 = std::max(val9,val18); val9 = tmp; tmp = std::min(val0,val18); val18 = std::max(val0,val18);
      val9 = std::max(tmp,val9); tmp = std::min(val10,val19); val19 = std::max(val10,val19); val10 = tmp;
      tmp = std::min(val1,val19); val19 = std::max(val1,val19); val1 = std::min(tmp,val10);
      val10 = std::max(tmp,val10); tmp = std::min(val11,val20); val20 = std::max(val11,val20); val11 = tmp;
      tmp = std::min(val2,val20); val20 = std::max(val2,val20); val11 = std::max(tmp,val11);
      tmp = std::min(val12,val21); val21 = std::max(val12,val21); val12 = tmp; tmp = std::min(val3,val21);
      val21 = std::max(val3,val21); val3 = std::min(tmp,val12); val12 = std::max(tmp,val12);
      tmp = std::min(val13,val22); val22 = std::max(val13,val22); val4 = std::min(val4,val22);
      val13 = std::max(val4,tmp); tmp = std::min(val4,tmp); val4 = tmp; tmp = std::min(val14,val23);
      val23 = std::max(val14,val23); val14 = tmp; tmp = std::min(val5,val23); val23 = std::max(val5,val23);
      val5 = std::min(tmp,val14); val14 = std::max(tmp,val14); tmp = std::min(val15,val24);
      val24 = std::max(val15,val24); val15 = tmp; val6 = std::min(val6,val24); tmp = std::min(val6,val15);
      val15 = std::max(val6,val15); val6 = tmp; tmp = std::min(val7,val16); val7 = std::min(tmp,val19);
      tmp = std::min(val13,val21); val15 = std::min(val15,val23); tmp = std::min(val7,tmp);
      val7 = std::min(tmp,val15); val9 = std::max(val1,val9); val11 = std::max(val3,val11);
      val17 = std::max(val5,val17); val17 = std::max(val11,val17); val17 = std::max(val9,val17);
      tmp = std::min(val4,val10); val10 = std::max(val4,val10); val4 = tmp; tmp = std::min(val6,val12);
      val12 = std::max(val6,val12); val6 = tmp; tmp = std::min(val7,val14); val14 = std::max(val7,val14);
      val7 = tmp; tmp = std::min(val4,val6); val6 = std::max(val4,val6); val7 = std::max(tmp,val7);
      tmp = std::min(val12,val14); val14 = std::max(val12,val14); val12 = tmp; val10 = std::min(val10,val14);
      tmp = std::min(val6,val7); val7 = std::max(val6,val7); val6 = tmp; tmp = std::min(val10,val12);
      val12 = std::max(val10,val12); val10 = std::max(val6,tmp); tmp = std::min(val6,tmp);
      val17 = std::max(tmp,val17); tmp = std::min(val12,val17); val17 = std::max(val12,val17); val12 = tmp;
      val7 = std::min(val7,val17); tmp = std::min(val7,val10); val10 = std::max(val7,val10); val7 = tmp;
      tmp = std::min(val12,val18); val18 = std::max(val12,val18); val12 = std::max(val7,tmp);
      val10 = std::min(val10,val18); tmp = std::min(val12,val20); val20 = std::max(val12,val20); val12 = tmp;
      tmp = std::min(val10,val20);
      return std::max(tmp,val12);
    }

    template<typename T>
    inline T median(T val0, T val1, T val2, T val3, T val4, T val5, T val6,
                    T val7, T val8, T val9, T val10, T val11, T val12, T val13,
                    T val14, T val15, T val16, T val17, T val18, T val19, T val20,
                    T val21, T val22, T val23, T val24, T val25, T val26, T val27,
                    T val28, T val29, T val30, T val31, T val32, T val33, T val34,
                    T val35, T val36, T val37, T val38, T val39, T val40, T val41,
                    T val42, T val43, T val44, T val45, T val46, T val47, T val48) {
      T tmp = std::min(val0,val32);
      val32 = std::max(val0,val32); val0 = tmp; tmp = std::min(val1,val33); val33 = std::max(val1,val33); val1 = tmp;
      tmp = std::min(val2,val34); val34 = std::max(val2,val34); val2 = tmp; tmp = std::min(val3,val35);
      val35 = std::max(val3,val35); val3 = tmp; tmp = std::min(val4,val36); val36 = std::max(val4,val36); val4 = tmp;
      tmp = std::min(val5,val37); val37 = std::max(val5,val37); val5 = tmp; tmp = std::min(val6,val38);
      val38 = std::max(val6,val38); val6 = tmp; tmp = std::min(val7,val39); val39 = std::max(val7,val39); val7 = tmp;
      tmp = std::min(val8,val40); val40 = std::max(val8,val40); val8 = tmp; tmp = std::min(val9,val41);
      val41 = std::max(val9,val41); val9 = tmp; tmp = std::min(val10,val42); val42 = std::max(val10,val42);
      val10 = tmp; tmp = std::min(val11,val43); val43 = std::max(val11,val43); val11 = tmp;
      tmp = std::min(val12,val44); val44 = std::max(val12,val44); val12 = tmp; tmp = std::min(val13,val45);
      val45 = std::max(val13,val45); val13 = tmp; tmp = std::min(val14,val46); val46 = std::max(val14,val46);
      val14 = tmp; tmp = std::min(val15,val47); val47 = std::max(val15,val47); val15 = tmp;
      tmp = std::min(val16,val48); val48 = std::max(val16,val48); val16 = tmp; tmp = std::min(val0,val16);
      val16 = std::max(val0,val16); val0 = tmp; tmp = std::min(val1,val17); val17 = std::max(val1,val17);
      val1 = tmp; tmp = std::min(val2,val18); val18 = std::max(val2,val18); val2 = tmp; tmp = std::min(val3,val19);
      val19 = std::max(val3,val19); val3 = tmp; tmp = std::min(val4,val20); val20 = std::max(val4,val20); val4 = tmp;
      tmp = std::min(val5,val21); val21 = std::max(val5,val21); val5 = tmp; tmp = std::min(val6,val22);
      val22 = std::max(val6,val22); val6 = tmp; tmp = std::min(val7,val23); val23 = std::max(val7,val23); val7 = tmp;
      tmp = std::min(val8,val24); val24 = std::max(val8,val24); val8 = tmp; tmp = std::min(val9,val25);
      val25 = std::max(val9,val25); val9 = tmp; tmp = std::min(val10,val26); val26 = std::max(val10,val26);
      val10 = tmp; tmp = std::min(val11,val27); val27 = std::max(val11,val27); val11 = tmp;
      tmp = std::min(val12,val28); val28 = std::max(val12,val28); val12 = tmp; tmp = std::min(val13,val29);
      val29 = std::max(val13,val29); val13 = tmp; tmp = std::min(val14,val30); val30 = std::max(val14,val30);
      val14 = tmp; tmp = std::min(val15,val31); val31 = std::max(val15,val31); val15 = tmp;
      tmp = std::min(val32,val48); val48 = std::max(val32,val48); val32 = tmp; tmp = std::min(val16,val32);
      val32 = std::max(val16,val32); val16 = tmp; tmp = std::min(val17,val33); val33 = std::max(val17,val33);
      val17 = tmp; tmp = std::min(val18,val34); val34 = std::max(val18,val34); val18 = tmp;
      tmp = std::min(val19,val35); val35 = std::max(val19,val35); val19 = tmp; tmp = std::min(val20,val36);
      val36 = std::max(val20,val36); val20 = tmp; tmp = std::min(val21,val37); val37 = std::max(val21,val37);
      val21 = tmp; tmp = std::min(val22,val38); val38 = std::max(val22,val38); val22 = tmp;
      tmp = std::min(val23,val39); val39 = std::max(val23,val39); val23 = tmp; tmp = std::min(val24,val40);
      val40 = std::max(val24,val40); val24 = tmp; tmp = std::min(val25,val41); val41 = std::max(val25,val41);
      val25 = tmp; tmp = std::min(val26,val42); val42 = std::max(val26,val42); val26 = tmp;
      tmp = std::min(val27,val43); val43 = std::max(val27,val43); val27 = tmp; tmp = std::min(val28,val44);
      val44 = std::max(val28,val44); val28 = tmp; tmp = std::min(val29,val45); val45 = std::max(val29,val45);
      val29 = tmp; tmp = std::min(val30,val46); val46 = std::max(val30,val46); val30 = tmp;
      tmp = std::min(val31,val47); val47 = std::max(val31,val47); val31 = tmp; tmp = std::min(val0,val8);
      val8 = std::max(val0,val8); val0 = tmp; tmp = std::min(val1,val9); val9 = std::max(val1,val9); val1 = tmp;
      tmp = std::min(val2,val10); val10 = std::max(val2,val10); val2 = tmp; tmp = std::min(val3,val11);
      val11 = std::max(val3,val11); val3 = tmp; tmp = std::min(val4,val12); val12 = std::max(val4,val12); val4 = tmp;
      tmp = std::min(val5,val13); val13 = std::max(val5,val13); val5 = tmp; tmp = std::min(val6,val14);
      val14 = std::max(val6,val14); val6 = tmp; tmp = std::min(val7,val15); val15 = std::max(val7,val15); val7 = tmp;
      tmp = std::min(val16,val24); val24 = std::max(val16,val24); val16 = tmp; tmp = std::min(val17,val25);
      val25 = std::max(val17,val25); val17 = tmp; tmp = std::min(val18,val26); val26 = std::max(val18,val26);
      val18 = tmp; tmp = std::min(val19,val27); val27 = std::max(val19,val27); val19 = tmp;
      tmp = std::min(val20,val28); val28 = std::max(val20,val28); val20 = tmp; tmp = std::min(val21,val29);
      val29 = std::max(val21,val29); val21 = tmp; tmp = std::min(val22,val30); val30 = std::max(val22,val30);
      val22 = tmp; tmp = std::min(val23,val31); val31 = std::max(val23,val31); val23 = tmp;
      tmp = std::min(val32,val40); val40 = std::max(val32,val40); val32 = tmp; tmp = std::min(val33,val41);
      val41 = std::max(val33,val41); val33 = tmp; tmp = std::min(val34,val42); val42 = std::max(val34,val42);
      val34 = tmp; tmp = std::min(val35,val43); val43 = std::max(val35,val43); val35 = tmp;
      tmp = std::min(val36,val44); val44 = std::max(val36,val44); val36 = tmp; tmp = std::min(val37,val45);
      val45 = std::max(val37,val45); val37 = tmp; tmp = std::min(val38,val46); val46 = std::max(val38,val46);
      val38 = tmp; tmp = std::min(val39,val47); val47 = std::max(val39,val47); val39 = tmp;
      tmp = std::min(val8,val32); val32 = std::max(val8,val32); val8 = tmp; tmp = std::min(val9,val33);
      val33 = std::max(val9,val33); val9 = tmp; tmp = std::min(val10,val34); val34 = std::max(val10,val34);
      val10 = tmp; tmp = std::min(val11,val35); val35 = std::max(val11,val35); val11 = tmp;
      tmp = std::min(val12,val36); val36 = std::max(val12,val36); val12 = tmp; tmp = std::min(val13,val37);
      val37 = std::max(val13,val37); val13 = tmp; tmp = std::min(val14,val38); val38 = std::max(val14,val38);
      val14 = tmp; tmp = std::min(val15,val39); val39 = std::max(val15,val39); val15 = tmp;
      tmp = std::min(val24,val48); val48 = std::max(val24,val48); val24 = tmp; tmp = std::min(val8,val16);
      val16 = std::max(val8,val16); val8 = tmp; tmp = std::min(val9,val17); val17 = std::max(val9,val17);
      val9 = tmp; tmp = std::min(val10,val18); val18 = std::max(val10,val18); val10 = tmp;
      tmp = std::min(val11,val19); val19 = std::max(val11,val19); val11 = tmp; tmp = std::min(val12,val20);
      val20 = std::max(val12,val20); val12 = tmp; tmp = std::min(val13,val21); val21 = std::max(val13,val21);
      val13 = tmp; tmp = std::min(val14,val22); val22 = std::max(val14,val22); val14 = tmp;
      tmp = std::min(val15,val23); val23 = std::max(val15,val23); val15 = tmp; tmp = std::min(val24,val32);
      val32 = std::max(val24,val32); val24 = tmp; tmp = std::min(val25,val33); val33 = std::max(val25,val33);
      val25 = tmp; tmp = std::min(val26,val34); val34 = std::max(val26,val34); val26 = tmp;
      tmp = std::min(val27,val35); val35 = std::max(val27,val35); val27 = tmp; tmp = std::min(val28,val36);
      val36 = std::max(val28,val36); val28 = tmp; tmp = std::min(val29,val37); val37 = std::max(val29,val37);
      val29 = tmp; tmp = std::min(val30,val38); val38 = std::max(val30,val38); val30 = tmp;
      tmp = std::min(val31,val39); val39 = std::max(val31,val39); val31 = tmp; tmp = std::min(val40,val48);
      val48 = std::max(val40,val48); val40 = tmp; tmp = std::min(val0,val4); val4 = std::max(val0,val4);
      val0 = tmp; tmp = std::min(val1,val5); val5 = std::max(val1,val5); val1 = tmp; tmp = std::min(val2,val6);
      val6 = std::max(val2,val6); val2 = tmp; tmp = std::min(val3,val7); val7 = std::max(val3,val7); val3 = tmp;
      tmp = std::min(val8,val12); val12 = std::max(val8,val12); val8 = tmp; tmp = std::min(val9,val13);
      val13 = std::max(val9,val13); val9 = tmp; tmp = std::min(val10,val14); val14 = std::max(val10,val14);
      val10 = tmp; tmp = std::min(val11,val15); val15 = std::max(val11,val15); val11 = tmp;
      tmp = std::min(val16,val20); val20 = std::max(val16,val20); val16 = tmp; tmp = std::min(val17,val21);
      val21 = std::max(val17,val21); val17 = tmp; tmp = std::min(val18,val22); val22 = std::max(val18,val22);
      val18 = tmp; tmp = std::min(val19,val23); val23 = std::max(val19,val23); val19 = tmp;
      tmp = std::min(val24,val28); val28 = std::max(val24,val28); val24 = tmp; tmp = std::min(val25,val29);
      val29 = std::max(val25,val29); val25 = tmp; tmp = std::min(val26,val30); val30 = std::max(val26,val30);
      val26 = tmp; tmp = std::min(val27,val31); val31 = std::max(val27,val31); val27 = tmp;
      tmp = std::min(val32,val36); val36 = std::max(val32,val36); val32 = tmp; tmp = std::min(val33,val37);
      val37 = std::max(val33,val37); val33 = tmp; tmp = std::min(val34,val38); val38 = std::max(val34,val38);
      val34 = tmp; tmp = std::min(val35,val39); val39 = std::max(val35,val39); val35 = tmp;
      tmp = std::min(val40,val44); val44 = std::max(val40,val44); val40 = tmp; tmp = std::min(val41,val45);
      val45 = std::max(val41,val45); val41 = tmp; tmp = std::min(val42,val46); val46 = std::max(val42,val46);
      val42 = tmp; tmp = std::min(val43,val47); val47 = std::max(val43,val47); val43 = tmp;
      tmp = std::min(val4,val32); val32 = std::max(val4,val32); val4 = tmp; tmp = std::min(val5,val33);
      val33 = std::max(val5,val33); val5 = tmp; tmp = std::min(val6,val34); val34 = std::max(val6,val34);
      val6 = tmp; tmp = std::min(val7,val35); val35 = std::max(val7,val35); val7 = tmp;
      tmp = std::min(val12,val40); val40 = std::max(val12,val40); val12 = tmp; tmp = std::min(val13,val41);
      val41 = std::max(val13,val41); val13 = tmp; tmp = std::min(val14,val42); val42 = std::max(val14,val42);
      val14 = tmp; tmp = std::min(val15,val43); val43 = std::max(val15,val43); val15 = tmp;
      tmp = std::min(val20,val48); val48 = std::max(val20,val48); val20 = tmp; tmp = std::min(val4,val16);
      val16 = std::max(val4,val16); val4 = tmp; tmp = std::min(val5,val17); val17 = std::max(val5,val17);
      val5 = tmp; tmp = std::min(val6,val18); val18 = std::max(val6,val18); val6 = tmp;
      tmp = std::min(val7,val19); val19 = std::max(val7,val19); val7 = tmp; tmp = std::min(val12,val24);
      val24 = std::max(val12,val24); val12 = tmp; tmp = std::min(val13,val25); val25 = std::max(val13,val25);
      val13 = tmp; tmp = std::min(val14,val26); val26 = std::max(val14,val26); val14 = tmp;
      tmp = std::min(val15,val27); val27 = std::max(val15,val27); val15 = tmp; tmp = std::min(val20,val32);
      val32 = std::max(val20,val32); val20 = tmp; tmp = std::min(val21,val33); val33 = std::max(val21,val33);
      val21 = tmp; tmp = std::min(val22,val34); val34 = std::max(val22,val34); val22 = tmp;
      tmp = std::min(val23,val35); val35 = std::max(val23,val35); val23 = tmp; tmp = std::min(val28,val40);
      val40 = std::max(val28,val40); val28 = tmp; tmp = std::min(val29,val41); val41 = std::max(val29,val41);
      val29 = tmp; tmp = std::min(val30,val42); val42 = std::max(val30,val42); val30 = tmp;
      tmp = std::min(val31,val43); val43 = std::max(val31,val43); val31 = tmp; tmp = std::min(val36,val48);
      val48 = std::max(val36,val48); val36 = tmp; tmp = std::min(val4,val8); val8 = std::max(val4,val8);
      val4 = tmp; tmp = std::min(val5,val9); val9 = std::max(val5,val9); val5 = tmp; tmp = std::min(val6,val10);
      val10 = std::max(val6,val10); val6 = tmp; tmp = std::min(val7,val11); val11 = std::max(val7,val11); val7 = tmp;
      tmp = std::min(val12,val16); val16 = std::max(val12,val16); val12 = tmp; tmp = std::min(val13,val17);
      val17 = std::max(val13,val17); val13 = tmp; tmp = std::min(val14,val18); val18 = std::max(val14,val18);
      val14 = tmp; tmp = std::min(val15,val19); val19 = std::max(val15,val19); val15 = tmp;
      tmp = std::min(val20,val24); val24 = std::max(val20,val24); val20 = tmp; tmp = std::min(val21,val25);
      val25 = std::max(val21,val25); val21 = tmp; tmp = std::min(val22,val26); val26 = std::max(val22,val26);
      val22 = tmp; tmp = std::min(val23,val27); val27 = std::max(val23,val27); val23 = tmp;
      tmp = std::min(val28,val32); val32 = std::max(val28,val32); val28 = tmp; tmp = std::min(val29,val33);
      val33 = std::max(val29,val33); val29 = tmp; tmp = std::min(val30,val34); val34 = std::max(val30,val34);
      val30 = tmp; tmp = std::min(val31,val35); val35 = std::max(val31,val35); val31 = tmp;
      tmp = std::min(val36,val40); val40 = std::max(val36,val40); val36 = tmp; tmp = std::min(val37,val41);
      val41 = std::max(val37,val41); val37 = tmp; tmp = std::min(val38,val42); val42 = std::max(val38,val42);
      val38 = tmp; tmp = std::min(val39,val43); val43 = std::max(val39,val43); val39 = tmp;
      tmp = std::min(val44,val48); val48 = std::max(val44,val48); val44 = tmp; tmp = std::min(val0,val2);
      val2 = std::max(val0,val2); val0 = tmp; tmp = std::min(val1,val3); val3 = std::max(val1,val3); val1 = tmp;
      tmp = std::min(val4,val6); val6 = std::max(val4,val6); val4 = tmp; tmp = std::min(val5,val7);
      val7 = std::max(val5,val7); val5 = tmp; tmp = std::min(val8,val10); val10 = std::max(val8,val10); val8 = tmp;
      tmp = std::min(val9,val11); val11 = std::max(val9,val11); val9 = tmp; tmp = std::min(val12,val14);
      val14 = std::max(val12,val14); val12 = tmp; tmp = std::min(val13,val15); val15 = std::max(val13,val15);
      val13 = tmp; tmp = std::min(val16,val18); val18 = std::max(val16,val18); val16 = tmp;
      tmp = std::min(val17,val19); val19 = std::max(val17,val19); val17 = tmp; tmp = std::min(val20,val22);
      val22 = std::max(val20,val22); val20 = tmp; tmp = std::min(val21,val23); val23 = std::max(val21,val23);
      val21 = tmp; tmp = std::min(val24,val26); val26 = std::max(val24,val26); val24 = tmp;
      tmp = std::min(val25,val27); val27 = std::max(val25,val27); val25 = tmp; tmp = std::min(val28,val30);
      val30 = std::max(val28,val30); val28 = tmp; tmp = std::min(val29,val31); val31 = std::max(val29,val31);
      val29 = tmp; tmp = std::min(val32,val34); val34 = std::max(val32,val34); val32 = tmp;
      tmp = std::min(val33,val35); val35 = std::max(val33,val35); val33 = tmp; tmp = std::min(val36,val38);
      val38 = std::max(val36,val38); val36 = tmp; tmp = std::min(val37,val39); val39 = std::max(val37,val39);
      val37 = tmp; tmp = std::min(val40,val42); val42 = std::max(val40,val42); val40 = tmp;
      tmp = std::min(val41,val43); val43 = std::max(val41,val43); val41 = tmp; tmp = std::min(val44,val46);
      val46 = std::max(val44,val46); val44 = tmp; tmp = std::min(val45,val47); val47 = std::max(val45,val47);
      val45 = tmp; tmp = std::min(val2,val32); val32 = std::max(val2,val32); val2 = tmp; tmp = std::min(val3,val33);
      val33 = std::max(val3,val33); val3 = tmp; tmp = std::min(val6,val36); val36 = std::max(val6,val36); val6 = tmp;
      tmp = std::min(val7,val37); val37 = std::max(val7,val37); val7 = tmp; tmp = std::min(val10,val40);
      val40 = std::max(val10,val40); val10 = tmp; tmp = std::min(val11,val41); val41 = std::max(val11,val41);
      val11 = tmp; tmp = std::min(val14,val44); val44 = std::max(val14,val44); val14 = tmp;
      tmp = std::min(val15,val45); val45 = std::max(val15,val45); val15 = tmp; tmp = std::min(val18,val48);
      val48 = std::max(val18,val48); val18 = tmp; tmp = std::min(val2,val16); val16 = std::max(val2,val16);
      val2 = tmp; tmp = std::min(val3,val17); val17 = std::max(val3,val17); val3 = tmp;
      tmp = std::min(val6,val20); val20 = std::max(val6,val20); val6 = tmp; tmp = std::min(val7,val21);
      val21 = std::max(val7,val21); val7 = tmp; tmp = std::min(val10,val24); val24 = std::max(val10,val24);
      val10 = tmp; tmp = std::min(val11,val25); val25 = std::max(val11,val25); val11 = tmp;
      tmp = std::min(val14,val28); val28 = std::max(val14,val28); val14 = tmp; tmp = std::min(val15,val29);
      val29 = std::max(val15,val29); val15 = tmp; tmp = std::min(val18,val32); val32 = std::max(val18,val32);
      val18 = tmp; tmp = std::min(val19,val33); val33 = std::max(val19,val33); val19 = tmp;
      tmp = std::min(val22,val36); val36 = std::max(val22,val36); val22 = tmp; tmp = std::min(val23,val37);
      val37 = std::max(val23,val37); val23 = tmp; tmp = std::min(val26,val40); val40 = std::max(val26,val40);
      val26 = tmp; tmp = std::min(val27,val41); val41 = std::max(val27,val41); val27 = tmp;
      tmp = std::min(val30,val44); val44 = std::max(val30,val44); val30 = tmp; tmp = std::min(val31,val45);
      val45 = std::max(val31,val45); val31 = tmp; tmp = std::min(val34,val48); val48 = std::max(val34,val48);
      val34 = tmp; tmp = std::min(val2,val8); val8 = std::max(val2,val8); val2 = tmp; tmp = std::min(val3,val9);
      val9 = std::max(val3,val9); val3 = tmp; tmp = std::min(val6,val12); val12 = std::max(val6,val12); val6 = tmp;
      tmp = std::min(val7,val13); val13 = std::max(val7,val13); val7 = tmp; tmp = std::min(val10,val16);
      val16 = std::max(val10,val16); val10 = tmp; tmp = std::min(val11,val17); val17 = std::max(val11,val17);
      val11 = tmp; tmp = std::min(val14,val20); val20 = std::max(val14,val20); val14 = tmp;
      tmp = std::min(val15,val21); val21 = std::max(val15,val21); val15 = tmp; tmp = std::min(val18,val24);
      val24 = std::max(val18,val24); val18 = tmp; tmp = std::min(val19,val25); val25 = std::max(val19,val25);
      val19 = tmp; tmp = std::min(val22,val28); val28 = std::max(val22,val28); val22 = tmp;
      tmp = std::min(val23,val29); val29 = std::max(val23,val29); val23 = tmp; tmp = std::min(val26,val32);
      val32 = std::max(val26,val32); val26 = tmp; tmp = std::min(val27,val33); val33 = std::max(val27,val33);
      val27 = tmp; tmp = std::min(val30,val36); val36 = std::max(val30,val36); val30 = tmp;
      tmp = std::min(val31,val37); val37 = std::max(val31,val37); val31 = tmp; tmp = std::min(val34,val40);
      val40 = std::max(val34,val40); val34 = tmp; tmp = std::min(val35,val41); val41 = std::max(val35,val41);
      val35 = tmp; tmp = std::min(val38,val44); val44 = std::max(val38,val44); val38 = tmp;
      tmp = std::min(val39,val45); val45 = std::max(val39,val45); val39 = tmp; tmp = std::min(val42,val48);
      val48 = std::max(val42,val48); val42 = tmp; tmp = std::min(val2,val4); val4 = std::max(val2,val4);
      val2 = tmp; tmp = std::min(val3,val5); val5 = std::max(val3,val5); val3 = tmp; tmp = std::min(val6,val8);
      val8 = std::max(val6,val8); val6 = tmp; tmp = std::min(val7,val9); val9 = std::max(val7,val9); val7 = tmp;
      tmp = std::min(val10,val12); val12 = std::max(val10,val12); val10 = tmp; tmp = std::min(val11,val13);
      val13 = std::max(val11,val13); val11 = tmp; tmp = std::min(val14,val16); val16 = std::max(val14,val16);
      val14 = tmp; tmp = std::min(val15,val17); val17 = std::max(val15,val17); val15 = tmp;
      tmp = std::min(val18,val20); val20 = std::max(val18,val20); val18 = tmp; tmp = std::min(val19,val21);
      val21 = std::max(val19,val21); val19 = tmp; tmp = std::min(val22,val24); val24 = std::max(val22,val24);
      val22 = tmp; tmp = std::min(val23,val25); val25 = std::max(val23,val25); val23 = tmp;
      tmp = std::min(val26,val28); val28 = std::max(val26,val28); val26 = tmp; tmp = std::min(val27,val29);
      val29 = std::max(val27,val29); val27 = tmp; tmp = std::min(val30,val32); val32 = std::max(val30,val32);
      val30 = tmp; tmp = std::min(val31,val33); val33 = std::max(val31,val33); val31 = tmp;
      tmp = std::min(val34,val36); val36 = std::max(val34,val36); val34 = tmp; tmp = std::min(val35,val37);
      val37 = std::max(val35,val37); val35 = tmp; tmp = std::min(val38,val40); val40 = std::max(val38,val40);
      val38 = tmp; tmp = std::min(val39,val41); val41 = std::max(val39,val41); val39 = tmp;
      tmp = std::min(val42,val44); val44 = std::max(val42,val44); val42 = tmp; tmp = std::min(val43,val45);
      val45 = std::max(val43,val45); val43 = tmp; tmp = std::min(val46,val48); val48 = std::max(val46,val48);
      val46 = tmp; val1 = std::max(val0,val1); val3 = std::max(val2,val3); val5 = std::max(val4,val5);
      val7 = std::max(val6,val7); val9 = std::max(val8,val9); val11 = std::max(val10,val11);
      val13 = std::max(val12,val13); val15 = std::max(val14,val15); val17 = std::max(val16,val17);
      val19 = std::max(val18,val19); val21 = std::max(val20,val21); val23 = std::max(val22,val23);
      val24 = std::min(val24,val25); val26 = std::min(val26,val27); val28 = std::min(val28,val29);
      val30 = std::min(val30,val31); val32 = std::min(val32,val33); val34 = std::min(val34,val35);
      val36 = std::min(val36,val37); val38 = std::min(val38,val39); val40 = std::min(val40,val41);
      val42 = std::min(val42,val43); val44 = std::min(val44,val45); val46 = std::min(val46,val47);
      val32 = std::max(val1,val32); val34 = std::max(val3,val34); val36 = std::max(val5,val36);
      val38 = std::max(val7,val38); val9 = std::min(val9,val40); val11 = std::min(val11,val42);
      val13 = std::min(val13,val44); val15 = std::min(val15,val46); val17 = std::min(val17,val48);
      val24 = std::max(val9,val24); val26 = std::max(val11,val26); val28 = std::max(val13,val28);
      val30 = std::max(val15,val30); val17 = std::min(val17,val32); val19 = std::min(val19,val34);
      val21 = std::min(val21,val36); val23 = std::min(val23,val38); val24 = std::max(val17,val24);
      val26 = std::max(val19,val26); val21 = std::min(val21,val28); val23 = std::min(val23,val30);
      val24 = std::max(val21,val24); val23 = std::min(val23,val26);
      return std::max(val23,val24);
    }

    //! Return sqrt(x^2 + y^2).
    template<typename T>
    inline T hypot(const T x, const T y) {
      return std::sqrt(x*x + y*y);
    }

    template<typename T>
    inline T hypot(const T x, const T y, const T z) {
      return std::sqrt(x*x + y*y + z*z);
    }

    template<typename T>
    inline T _hypot(const T x, const T y) { // Slower but more precise version
      T nx = cimg::abs(x), ny = cimg::abs(y), t;
      if (nx<ny) { t = nx; nx = ny; } else t = ny;
      if (nx>0) { t/=nx; return nx*std::sqrt(1 + t*t); }
      return 0;
    }

    //! Return the factorial of n
    inline double factorial(const int n) {
      if (n<0) return cimg::type<double>::nan();
      if (n<2) return 1;
      double res = 2;
      for (int i = 3; i<=n; ++i) res*=i;
      return res;
    }

    //! Return the number of permutations of k objects in a set of n objects.
    inline double permutations(const int k, const int n, const bool with_order) {
      if (n<0 || k<0) return cimg::type<double>::nan();
      if (k>n) return 0;
      double res = 1;
      for (int i = n; i>=n - k + 1; --i) res*=i;
      return with_order?res:res/cimg::factorial(k);
    }

    inline double _fibonacci(int exp) {
      double
        base = (1 + std::sqrt(5.0))/2,
        result = 1/std::sqrt(5.0);
      while (exp) {
        if (exp&1) result*=base;
        exp>>=1;
        base*=base;
      }
      return result;
    }

    //! Calculate fibonacci number.
    // (Precise up to n = 78, less precise for n>78).
    inline double fibonacci(const int n) {
      if (n<0) return cimg::type<double>::nan();
      if (n<3) return 1;
      if (n<11) {
        cimg_uint64 fn1 = 1, fn2 = 1, fn = 0;
        for (int i = 3; i<=n; ++i) { fn = fn1 + fn2; fn2 = fn1; fn1 = fn; }
        return (double)fn;
      }
      if (n<75) // precise up to n = 74, faster than the integer calculation above for n>10
        return (double)((cimg_uint64)(_fibonacci(n) + 0.5));

      if (n<94) { // precise up to n = 78, less precise for n>78 up to n = 93, overflows for n>93
        cimg_uint64
          fn1 = (cimg_uint64)1304969544928657U,
          fn2 = (cimg_uint64)806515533049393U,
          fn = 0;
        for (int i = 75; i<=n; ++i) { fn = fn1 + fn2; fn2 = fn1; fn1 = fn; }
        return (double)fn;
      }
      return _fibonacci(n); // Not precise, but better than the wrong overflowing calculation
    }

    //! Convert ascii character to lower case.
    inline char lowercase(const char x) {
      return (char)((x<'A'||x>'Z')?x:x - 'A' + 'a');
    }
    inline double lowercase(const double x) {
      return (double)((x<'A'||x>'Z')?x:x - 'A' + 'a');
    }

    //! Convert C-string to lower case.
    inline void lowercase(char *const str) {
      if (str) for (char *ptr = str; *ptr; ++ptr) *ptr = lowercase(*ptr);
    }

    //! Convert ascii character to upper case.
    inline char uppercase(const char x) {
      return (char)((x<'a'||x>'z')?x:x - 'a' + 'A');
    }

    inline double uppercase(const double x) {
      return (double)((x<'a'||x>'z')?x:x - 'a' + 'A');
    }

    //! Convert C-string to upper case.
    inline void uppercase(char *const str) {
      if (str) for (char *ptr = str; *ptr; ++ptr) *ptr = uppercase(*ptr);
    }

    //! Read value in a C-string.
    /**
       \param str C-string containing the float value to read.
       \return Read value.
       \note Same as <tt>std::atof()</tt> extended to manage the retrieval of fractions from C-strings,
       as in <em>"1/2"</em>.
    **/
    inline double atof(const char *const str) {
      double x = 0, y = 1;
      return str && cimg_sscanf(str,"%lf/%lf",&x,&y)>0?x/y:0;
    }

    //! Compare the first \p l characters of two C-strings, ignoring the case.
    /**
       \param str1 C-string.
       \param str2 C-string.
       \param l Number of characters to compare.
       \return \c 0 if the two strings are equal, something else otherwise.
       \note This function has to be defined since it is not provided by all C++-compilers (not ANSI).
    **/
    inline int strncasecmp(const char *const str1, const char *const str2, const int l) {
      if (!l) return 0;
      if (!str1) return str2?-1:0;
      const char *nstr1 = str1, *nstr2 = str2;
      int k, diff = 0; for (k = 0; k<l && !(diff = lowercase(*nstr1) - lowercase(*nstr2)); ++k) { ++nstr1; ++nstr2; }
      return k!=l?diff:0;
    }

    //! Compare two C-strings, ignoring the case.
    /**
       \param str1 C-string.
       \param str2 C-string.
       \return \c 0 if the two strings are equal, something else otherwise.
       \note This function has to be defined since it is not provided by all C++-compilers (not ANSI).
    **/
    inline int strcasecmp(const char *const str1, const char *const str2) {
      if (!str1) return str2?-1:0;
      const int
        l1 = (int)std::strlen(str1),
        l2 = (int)std::strlen(str2);
      return cimg::strncasecmp(str1,str2,1 + (l1<l2?l1:l2));
    }

    //! Ellipsize a string.
    /**
       \param str C-string.
       \param l Max number of characters.
       \param is_ending Tell if the dots are placed at the end or at the center of the ellipsized string.
    **/
    inline char *strellipsize(char *const str, const unsigned int l=64,
                              const bool is_ending=true) {
      if (!str) return str;
      const unsigned int nl = l<5?5:l, ls = (unsigned int)std::strlen(str);
      if (ls<=nl) return str;
      if (is_ending) std::strcpy(str + nl - 5,"(...)");
      else {
        const unsigned int ll = (nl - 5)/2 + 1 - (nl%2), lr = nl - ll - 5;
        std::strcpy(str + ll,"(...)");
        std::memmove(str + ll + 5,str + ls - lr,lr);
      }
      str[nl] = 0;
      return str;
    }

    //! Ellipsize a string.
    /**
       \param str C-string.
       \param res output C-string.
       \param l Max number of characters.
       \param is_ending Tell if the dots are placed at the end or at the center of the ellipsized string.
    **/
    inline char *strellipsize(const char *const str, char *const res, const unsigned int l=64,
                              const bool is_ending=true) {
      const unsigned int nl = l<5?5:l, ls = (unsigned int)std::strlen(str);
      if (ls<=nl) { std::strcpy(res,str); return res; }
      if (is_ending) {
        std::strncpy(res,str,nl - 5);
        std::strcpy(res + nl -5,"(...)");
      } else {
        const unsigned int ll = (nl - 5)/2 + 1 - (nl%2), lr = nl - ll - 5;
        std::strncpy(res,str,ll);
        std::strcpy(res + ll,"(...)");
        std::strncpy(res + ll + 5,str + ls - lr,lr);
      }
      res[nl] = 0;
      return res;
    }

    //! Remove delimiters on the start and/or end of a C-string.
    /**
       \param[in,out] str C-string to work with (modified at output).
       \param delimiter Delimiter character code to remove.
       \param is_symmetric Tells if the removal is done only if delimiters are symmetric
       (both at the beginning and the end of \c s).
       \param is_iterative Tells if the removal is done if several iterations are possible.
       \return \c true if delimiters have been removed, \c false otherwise.
   **/
    inline bool strpare(char *const str, const char delimiter,
                        const bool is_symmetric, const bool is_iterative) {
      if (!str) return false;
      const int l = (int)std::strlen(str);
      int p, q;
      if (is_symmetric) for (p = 0, q = l - 1; p<q && str[p]==delimiter && str[q]==delimiter; ) {
          --q; ++p; if (!is_iterative) break;
        } else {
        for (p = 0; p<l && str[p]==delimiter; ) { ++p; if (!is_iterative) break; }
        for (q = l - 1; q>p && str[q]==delimiter; ) { --q; if (!is_iterative) break; }
      }
      const int n = q - p + 1;
      if (n!=l) { std::memmove(str,str + p,(unsigned int)n); str[n] = 0; return true; }
      return false;
    }

    //! Remove white spaces on the start and/or end of a C-string.
    inline bool strpare(char *const str, const bool is_symmetric, const bool is_iterative) {
      if (!str) return false;
      const int l = (int)std::strlen(str);