    package flare.util
    {
        /**
         * Utility methods for mathematics not found in the <code>Math</code> class.
         */
        public class Maths
        {
            private static var EPSILON:Number = 1e-30;

            /**
             * Constructor, throws an error if called, as this is an abstract class.
             */
            public function Maths() {
                throw new ArgumentError("This is an abstract class");
            }

            // -- Operators -------------------------------------------------------

            /**
             * Returns the logarithm with a given base value.
             * @param x the number for which to compute the logarithm
             * @param base the base of the logarithm
             * @return the logarithm value
             */
            public static function log(x:Number, base:Number):Number
            {
                return Math.log(x) / Math.log(base);
            }

            /**
             * Rounds an input value down according to its logarithm. The method
             *  takes the floor of the logarithm of the value and then uses the
             *  resulting value as an exponent for the base value. For example:
             * <code>
             * trace(Maths.logFloor(13,10)); // prints "10"
             * trace(Maths.logFloor(113,10)); // prints "100"
             * </code>
             * @param x the number for which to compute the logarithm floor
             * @param base the base of the logarithm
             * @return the rounded-by-logarithm value
             */
            public static function logFloor(x:Number, base:Number):Number {
                if (x > 0) {
                    return Math.pow(base, Math.floor(Maths.log(x, base)));
                } else {
                    return -Math.pow(base, -Math.floor(-Maths.log(-x, base)));
                }
            }

            /**
             * Rounds an input value up according to its logarithm. The method
             *  takes the ceiling of the logarithm of the value and then uses the
             *  resulting value as an exponent for the base value. For example:
             * <code>
             * trace(Maths.logCeil(13,10)); // prints "100"
             * trace(Maths.logCeil(113,10)); // prints "1000"
             * </code>
             * @param x the number for which to compute the logarithm ceiling
             * @param base the base of the logarithm
             * @return the rounded-by-logarithm value
             */
            public static function logCeil(x:Number, base:Number):Number {
                if (x > 0) {
                    return Math.pow(base, Math.ceil(Maths.log(x, base)));
                } else {
                    return -Math.pow(base, -Math.ceil(-Maths.log(-x, base)));
                }
            }

            /**
             * Computes a zero-symmetric logarithm. Computes the logarithm of the
             * absolute value of the input, and determines the sign of the output
             * according to the sign of the input value.
             * @param x the number for which to compute the logarithm
             * @param b the base of the logarithm
             * @return the symmetric log value.
             */
            public static function symLog(x:Number, b:Number) : Number
            {
                return x == 0 ? 0 : x > 0 ? log(x, b) : -log(-x, b);
            }

            /**
             * Computes a zero-symmetric logarithm, with adjustment to values
             * between zero and the logarithm base. This adjustment introduces
             * distortion for values less than the base number, but enables
             * simultaneous plotting of log-transformed data involving both
             * positive and negative numbers.
             * @param x the number for which to compute the logarithm
             * @param b the base of the logarithm
             * @return the adjusted, symmetric log value.
             */
            public static function adjLog(x:Number, b:Number) : Number
            {
                var neg:Boolean = (x < 0);
                if (neg) x = -x;
                if (x < b) x += (b-x)/b;
                return neg ? -log(x,b) : log(x,b);
            }

            /**
             * Computes a zero-symmetric square root. Computes the square root of
             * the absolute value of the input, and determines the sign of the
             * output according to the sign of the input value.
             * @param x the number for which to compute the square root
             * @return the symmetric square root value.
             */
            public static function symSqrt(x:Number) : Number
            {
                return (x > 0 ? Math.sqrt(x) : -Math.sqrt(-x));
            }

            /**
             * Computes a zero-symmetric Nth root. Computes the root of
             * the absolute value of the input, and determines the sign of the
             * output according to the sign of the input value.
             * @param x the number for which to compute the square root
             * @param p the root value (2 for square root, 3 for cubic root, etc)
             * @return the symmetric Nth root value.
             */
            public static function symRoot(x:Number, p:Number) : Number
            {
                return (x > 0 ? Math.pow(x, 1/p) : -Math.pow(-x, 1/p));
            }

            // -- Statistics ------------------------------------------------------

            /**
             * Computes the n-quantile boundaries for a set of values.
             * @param n the number of quantiles
             * @param values the values to break into quantiles
             * @param sort indicates if the values array needs to be sorted. If
             *  true, the array will be sorted prior to determining quantile
             *  boundaries. If false, the array must be in ascending sort order.
             * @return an n-length array of quantile boundaries
             */
            public static function quantile(n:uint, values:Array, sort:Boolean):Array
            {
                var len:uint = values.length-1;
                var qtls:Array = new Array(n);

                if (sort) {
                    values = Arrays.copy(values);
                    values.sort(Array.NUMERIC);
                   }

                for (var i:uint=1; i <= n; ++i)
                {
                    qtls[i-1] = values[uint((len*i)/n)];
                }
                return qtls;
            }

            // -- Forward Interpolation Routines ----------------------------------

            /**
             * Computes a linear interpolation between two values.
             * @param f the interpolation fraction (typically between 0 and 1)
             * @param min the minimum value (corresponds to f==0)
             * @param max the maximum value (corresponds to f==1)
             * @return the interpolated value
             */
            public static function linearInterp(f:Number, min:Number, max:Number):Number
            {
                return min + f * (max - min);
            }

            /**
             * Computes a logarithmic interpolation between two values. Uses a
             * zero-symmetric logarithm calculation (see <code>symLog</code>).
             * @param f the interpolation fraction (typically between 0 and 1)
             * @param min the minimum value (corresponds to f==0)
             * @param max the maximum value (corresponds to f==1)
             * @param b the base of the logarithm
             * @return the interpolated value
             */
            public static function logInterp(f:Number, min:Number, max:Number, b:Number):Number
            {
                min = symLog(min, b); max = symLog(max, b);
                f = min + f * (max - min);
                return f<0 ? -Math.pow(b, -f) : Math.pow(b, f);
            }

            /**
             * Computes a logarithmic interpolation between two values. Uses an
             * adjusted zero-symmetric logarithm calculation (see <code>adjLog</code>).
             * @param f the interpolation fraction (typically between 0 and 1)
             * @param min the minimum value (corresponds to f==0)
             * @param max the maximum value (corresponds to f==1)
             * @param b the base of the logarithm
             * @return the interpolated value
             */
            public static function adjLogInterp(f:Number, min:Number, max:Number, b:Number):Number
            {
                min = adjLog(min, b); max = adjLog(max, b);
                f = min + f * (max - min);
                var neg:Boolean = f < 0;
                f = neg ? Math.pow(b, -f) : Math.pow(b, f);
                f = b*(f-1) / (b-1);
                return neg ? -f : f;
            }

            /**
             * Computes a square root interpolation between two values. Uses a
             * zero-symmetric root calculation (see <code>symSqrt</code>).
             * @param f the interpolation fraction (typically between 0 and 1)
             * @param min the minimum value (corresponds to f==0)
             * @param max the maximum value (corresponds to f==1)
             * @return the interpolated value
             */
            public static function sqrtInterp(f:Number, min:Number, max:Number):Number
            {
                min = symSqrt(min); max = symSqrt(max);
                f = min + f * (max - min);
                return f<0 ? -(f*f) : f*f;
            }

            /**
             * Computes an Nth-root interpolation between two values. Uses a
             * zero-symmetric root calculation (see <code>symRoot</code>).
             * @param f the interpolation fraction (typically between 0 and 1)
             * @param min the minimum value (corresponds to f==0)
             * @param max the maximum value (corresponds to f==1)
             * @param p the root value (2 for square root, 3 for cubic root, etc)
             * @return the interpolated value
             */
            public static function rootInterp(f:Number, min:Number, max:Number, p:Number):Number
            {
                min = symRoot(min, p); max = symRoot(max, p);
                f = min + f*(max - min);
                var neg:Boolean = f < 0;
                return neg ? -Math.pow(-f, p) : Math.pow(f, p);
            }

            /**
             * Computes an interpolated value in a quantile scale.
             * @param f the interpolation fraction (typically between 0 and 1)
             * @param quantiles an array of quantile boundaries
             * @return the interpolated value
             */
            public static function quantileInterp(f:Number, quantiles:Array):Number
            {
                return quantiles[int(Math.round(f*(quantiles.length-1)))];
            }

            // -- Inverse Interpolation Routines ----------------------------------

            /**
             * Computes an inverse linear interpolation, returning an interpolation
             * fraction. Returns 0.5 if the min and max values are the same.
             * @param x the interpolated value
             * @param min the minimum value (corresponds to f==0)
             * @param min the maximum value (corresponds to f==1)
             * @return the inferred interpolation fraction
             */
            public static function invLinearInterp(x:Number, min:Number, max:Number):Number
            {
                var denom:Number = (max - min);
                return (denom < EPSILON && denom > -EPSILON ? 0 : (x - min) / denom);
            }

            /**
             * Computes an inverse logarithmic interpolation, returning an
             * interpolation fraction. Uses a zero-symmetric logarithm.
             * Returns 0.5 if the min and max values are the same.
             * @param x the interpolated value
             * @param min the minimum value (corresponds to f==0)
             * @param min the maximum value (corresponds to f==1)
             * @param b the base of the logarithm
             * @return the inferred interpolation fraction
             */
            public static function invLogInterp(x:Number, min:Number, max:Number, b:Number):Number
            {
                min = symLog(min, b);
                var denom:Number = symLog(max, b) - min;
                return (denom < EPSILON && denom > -EPSILON ? 0 : (symLog(x, b) - min) / denom);
            }

            /**
             * Computes an inverse logarithmic interpolation, returning an
             * interpolation fraction. Uses an adjusted zero-symmetric logarithm.
             * Returns 0.5 if the min and max values are the same.
             * @param x the interpolated value
             * @param min the minimum value (corresponds to f==0)
             * @param min the maximum value (corresponds to f==1)
             * @param b the base of the logarithm
             * @return the inferred interpolation fraction
             */
            public static function invAdjLogInterp(x:Number, min:Number, max:Number, b:Number):Number
            {
                min = adjLog(min, b);
                var denom:Number = adjLog(max, b) - min;
                return (denom < EPSILON && denom > -EPSILON ? 0 : (adjLog(x, b) - min) / denom);
            }

            /**
             * Computes an inverse square root interpolation, returning an
             * interpolation fraction. Uses a zero-symmetric square root.
             * Returns 0.5 if the min and max values are the same.
             * @param x the interpolated value
             * @param min the minimum value (corresponds to f==0)
             * @param min the maximum value (corresponds to f==1)
             * @return the inferred interpolation fraction
             */
            public static function invSqrtInterp(x:Number, min:Number, max:Number):Number
            {
                min = symSqrt(min);
                var denom:Number = symSqrt(max) - min;
                return (denom < EPSILON && denom > -EPSILON ? 0 : (symSqrt(x) - min) / denom);
            }

            /**
             * Computes an inverse Nth-root interpolation, returning an
             * interpolation fraction. Uses a zero-symmetric root.
             * Returns 0.5 if the min and max values are the same.
             * @param x the interpolated value
             * @param min the minimum value (corresponds to f==0)
             * @param min the maximum value (corresponds to f==1)
             * @param p the root value (2 for square root, 3 for cubic root, etc)
             * @return the inferred interpolation fraction
             */
            public static function invRootInterp(x:Number, min:Number, max:Number, p:Number):Number
            {
                min = symRoot(min,p);
                var denom:Number = symRoot(max,p) - min;
                return (denom < EPSILON  && denom > -EPSILON ? 0 : (symRoot(x,p) - min) / denom);
            }

            /**
             * Computes an inverse quantile scale interpolation, returning an
             * interpolation fraction.
             * @param x the interpolated value
             * @param quantiles an array of quantile boundaries
             * @return the inferred interpolation fraction
             */
            public static function invQuantileInterp(x:Number, quantiles:/*Number*/Array):Number
            {
                var a:uint = 0, b:uint = quantiles.length;
                var i:uint = b / 2;
                while (a < b)
                {   // binary search over the boundaries
                    if (quantiles[i] == x) {
                        for (; i>0 && quantiles[i-1] == x; --i);
                        break;
                    } else if (quantiles[i] < x) {
                        a = i+1;
                    } else {
                        b = i;
                    }
                    i = a + ((b - a) >> 1);
                }
                return Number(i) / (quantiles.length-1);
            }

            // -- Perlin Noise ----------------------------------------------------

            /**
             * Computes Perlin noise for a given (x, y, z) point.
             * @param x the x parameter
             * @param y the y parameter (default 0)
             * @param z the z parameter (default 0)
             * @return Perlin noise for the input co-ordinates
             */
            public static function noise(x:Number, y:Number=0, z:Number=0):Number
               {
                   var X:int,Y:int,Z:int,A:int,AA:int,AB:int,B:int,BA:int,BB:int;
                   var u:Number,v:Number,w:Number,a:Number,b:Number,c:Number,d:Number;
                   var xx:Number, yy:Number, zz:Number;

                   X = int(x); Y = int(y); Z = int(z);
                   x -= X;     y -= Y;     z -= Z;
                   X &= 255;   Y &= 255;   Z &= 255;
                   xx = x-1;   yy = y-1;   zz = z-1;

                   u = x * x * x * (x * (x*6 - 15) + 10);
                   v = y * y * y * (y * (y*6 - 15) + 10);
                   w = z * z * z * (z * (z*6 - 15) + 10);

                   A  = (_p[X  ]+Y);
                   AA = (_p[A  ]+Z);
                   AB = (_p[A+1]+Z);
                   B  = (_p[X+1]+Y);
                   BA = (_p[B  ]+Z);
                   BB = (_p[B+1]+Z);

                   // interpolate
                   a =       grad(_p[AA  ], x , y , z );
                   b =       grad(_p[AB  ], x , yy, z );
                   c =       grad(_p[AA+1], x , y , zz);
                   d =       grad(_p[AB+1], x , yy, zz);
                   a += u * (grad(_p[BA  ], xx, y , z ) - a);
                   b += u * (grad(_p[BB  ], xx, yy, z ) - b);
                   c += u * (grad(_p[BA+1], xx, y , zz) - c);
                   d += u * (grad(_p[BB+1], xx, yy, zz) - d);
                   a += v * (b - a);
                   c += v * (d - c);
                   return a + w * (c - a);
               }

            private static function grad(h:int, x:Number, y:Number, z:Number):Number
            {
                h &= 15;
                var u:Number = h<8 ? x : y;
                var v:Number = h<4 ? y : h==2||h==14 ? x : z;
                return ((h&1) == 0 ? u : -u) + ((h&2) == 0 ? v : -v);
            }

            private static const _p:Array = [151,160,137,91,90,15,131,13,201,95,96,
            53,194,233,7,225,140,36,103,30,69,142,8,99,37,240,21, 10,23,190, 6,148,
            247,120,234,75, 0,26,197,62,94,252,219,203,117,35, 11, 32,57,177,33,88,
            237,149,56,87,174, 20,125,136,171,168, 68,175,74,165, 71,134,139,48,27,
            166,77,146,158,231,83,111,229,122, 60,211,133,230,220,105,92, 41,55,46,
            245,40,244,102,143,54, 65,25,63,161, 1,216,80,73,209,76,132,187,208,89,
            18,169,200,196,135,130,116,188,159,86,164,100,109,198,173,186, 3,64,52,
            217,226,250,124,123, 5,202,38,147,118,126,255,82,85,212,207,206,59,227,
            47,16,58,17,182,189,28,42,223,183,170,213,119,248,152, 2,44,154,163,70,
            221,153,101,155,167, 43,172, 9,129,22,39,253, 19,98,108,110,79,113,224,
            232,178,185, 112,104,218,246,97,228,251,34,242,193,238,210,144, 12,191,
            179,162,241,81,51,145,235,249,14,239,107,49,192,214,31,181,199,106,157,
            184,84,204,176,115,121,50,45,127, 4,150,254,138,236,205,93,222,114, 67,
            29,24,72,243,141,128,195,78,66,215,61,156,180,// now we repeat the list
            151,160,137,91,90,15,131, 13,201,95,96,53,194,233, 7,225,140,36,103,30,
            69,142,8,99,37,240,21, 10,23,190, 6,148,247,120,234,75, 0,26,197,62,94,
            252,219,203,117,35,11,32,57,177,33,88,237,149,56,87,174,20,125,136,171,
            168, 68,175,74,165,71,134,139,48, 27,166,77,146,158,231,83,111,229,122,
            60,211,133,230,220,105,92,41,55,46,245,40,244,102,143,54, 65,25,63,161,
            1,216, 80,73,209,76,132,187,208, 89,18,169,200,196,135,130,116,188,159,
            86,164,100,109,198,173,186, 3,64,52,217,226,250,124,123, 5,202, 38,147,
            118,126,255,82,85,212,207,206,59,227,47,16,58,17,182,189,28,42,223,183,
            170,213,119,248,152, 2,44,154,163, 70,221,153,101,155,167,43,172,9,129,
            22,39,253, 19,98,108,110,79,113,224,232,178,185,112,104,218,246,97,228,
               251,34,242,193,238,210,144,12,191,179,162,241,81,51,145,235,249,14,239,
               107,49,192,214, 31,181,199,106,157,184,84,204,176,115,121,50,45,127, 4,
               150,254,138,236,205,93,222,114,67,29, 24,72,243,141,128,195, 78,66,215,
               61,156,180];

        } // end of class Maths
    }