public static int log2(int x) throws ArithmeticException
   {
      if (x <= 0)
        throw new ArithmeticException("Cannot calculate log of a negative or null value");

      if (x < 256)
        return LOG2_4096[x] >> 2;
      
      if (x < 512)
        return 1024 + (LOG2_4096[x>>1] >> 2);
      
      int log = 8;

      for (long y=x+1; y>=512; y>>=1)
        log++;

      // x is a power of 2 ?
      if ((x & (x-1)) == 0)
         return log << 10;

      long base = 0;
      long z = (long) x - (1L << log);

      // Use the fact that log2(x) = log2(2^(log2(x)+1)*y) = log2(2^p)+ 1 + ln(1-z)/ln(2)
      // with z in ]0, 0.5[, it yields log2(x) = p + 1 - (z/1 + z^2/2 + z^3/3 ...)/ln(2)
      // To improve accuracy (keep z in ]0, 0.25[), one can choose either (1+z) or (1-z).
      // EG: log2(257) = log2(256) + log2(1+1/256) is better approximated with Taylor 
      // series expansion than log2(512) + log2(1-255/512)
      if (z > 1L << (log-1))
      {         
         // z in [0.5, 0.75[ => rescale x so that z in [0, 0.25[
         if (z < 1L << (log-1) + 1L << (log-2)) 
         {
            base = 497;
            x = (int) ((x * 5L) / 7L);   
         }
         
         // z in [0.75, 1[ => select 1 - x Taylor series expansion
         log++;      
      } 
      else 
      {
         // z in [0.25, 0.5[ => rescale x so that z in [0, 0.25[
         if (z >= 1L << (log-2))
         {
            base = 269;
            x = (int) ((x * 5L) / 6L);
         }            
         
         // select 1 + x Taylor series expansion
      }
     
      z = x - (1 << log);
      long z2 = (z * z) >> log;
      long taylor = z;
      taylor -= (z2 >> 1);
      taylor += (((z*z2) / 3) >> log);
      taylor -= ((z2*z2) >> (log+2));
      taylor = (taylor * 5909) >> (log+2); // rescale: 4096*1/log(2)
      return (int) (base + (log << 10) + taylor);
   }