[GR-24066] Fix ComplexBuiltin.__abs__ performance

PullRequest: graalpython/1019

Author: qunaibit

graalpython/com.oracle.graal.python/src/com/oracle/graal/python/builtins/objects/complex/ComplexBuiltins.java

@@ -63,7 +63,6 @@
 import static com.oracle.graal.python.nodes.SpecialMethodNames.__STR__;
 import static com.oracle.graal.python.nodes.SpecialMethodNames.__SUB__;
 import static com.oracle.graal.python.nodes.SpecialMethodNames.__TRUEDIV__;
-import static com.oracle.graal.python.runtime.exception.PythonErrorType.OverflowError;
 
 import java.util.List;
 
@@ -104,27 +103,143 @@ public abstract static class AbsNode extends PythonUnaryBuiltinNode {
 
         public abstract double executeDouble(Object arg);
 
-        public static AbsNode create() {
-            return ComplexBuiltinsFactory.AbsNodeFactory.create();
-        }
-
         @Specialization
         double abs(PComplex c) {
-            if (!Double.isFinite(c.getReal()) || !Double.isFinite(c.getImag())) {
-                if (Double.isInfinite(c.getReal())) {
-                    return Math.abs(c.getReal());
-                } else if (Double.isInfinite(c.getImag())) {
-                    return Math.abs(c.getImag());
+            double x = c.getReal();
+            double y = c.getImag();
+            if (Double.isInfinite(x) || Double.isInfinite(y)) {
+                return Double.POSITIVE_INFINITY;
+            } else if (Double.isNaN(x) || Double.isNaN(y)) {
+                return Double.NaN;
+            } else {
+
+                final int expX = getExponent(x);
+                final int expY = getExponent(y);
+                if (expX > expY + 27) {
+                    // y is neglectible with respect to x
+                    return abs(x);
+                } else if (expY > expX + 27) {
+                    // x is neglectible with respect to y
+                    return abs(y);
                 } else {
-                    return Double.NaN;
+
+                    // find an intermediate scale to avoid both overflow and
+                    // underflow
+                    final int middleExp = (expX + expY) / 2;
+
+                    // scale parameters without losing precision
+                    final double scaledX = scalb(x, -middleExp);
+                    final double scaledY = scalb(y, -middleExp);
+
+                    // compute scaled hypotenuse
+                    final double scaledH = Math.sqrt(scaledX * scaledX + scaledY * scaledY);
+
+                    // remove scaling
+                    double r = scalb(scaledH, middleExp);
+                    if (Double.isInfinite(r)) {
+                        throw raise(PythonErrorType.OverflowError, ErrorMessages.ABSOLUTE_VALUE_TOO_LARGE);
+                    }
+                    return r;
                 }
             }
-            double r = Math.hypot(c.getReal(), c.getImag());
-            if (Double.isInfinite(r)) {
-                throw raise(OverflowError, ErrorMessages.ABSOLUTE_VALUE_TOO_LARGE);
+        }
+
+        private static final long MASK_NON_SIGN_LONG = 0x7fffffffffffffffL;
+
+        static double abs(double x) {
+            return Double.longBitsToDouble(MASK_NON_SIGN_LONG & Double.doubleToRawLongBits(x));
+        }
+
+        static double scalb(final double d, final int n) {
+
+            // first simple and fast handling when 2^n can be represented using
+            // normal numbers
+            if ((n > -1023) && (n < 1024)) {
+                return d * Double.longBitsToDouble(((long) (n + 1023)) << 52);
+            }
+
+            // handle special cases
+            if (Double.isNaN(d) || Double.isInfinite(d) || (d == 0)) {
+                return d;
             }
-            return r;
+            if (n < -2098) {
+                return (d > 0) ? 0.0 : -0.0;
+            }
+            if (n > 2097) {
+                return (d > 0) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
+            }
+
+            // decompose d
+            final long bits = Double.doubleToRawLongBits(d);
+            final long sign = bits & 0x8000000000000000L;
+            int exponent = ((int) (bits >>> 52)) & 0x7ff;
+            long mantissa = bits & 0x000fffffffffffffL;
+
+            // compute scaled exponent
+            int scaledExponent = exponent + n;
+
+            if (n < 0) {
+                // we are really in the case n <= -1023
+                if (scaledExponent > 0) {
+                    // both the input and the result are normal numbers, we only
+                    // adjust the exponent
+                    return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa);
+                } else if (scaledExponent > -53) {
+                    // the input is a normal number and the result is a subnormal
+                    // number
+
+                    // recover the hidden mantissa bit
+                    mantissa = mantissa | (1L << 52);
+
+                    // scales down complete mantissa, hence losing least significant
+                    // bits
+                    final long mostSignificantLostBit = mantissa & (1L << (-scaledExponent));
+                    mantissa = mantissa >>> (1 - scaledExponent);
+                    if (mostSignificantLostBit != 0) {
+                        // we need to add 1 bit to round up the result
+                        mantissa++;
+                    }
+                    return Double.longBitsToDouble(sign | mantissa);
+
+                } else {
+                    // no need to compute the mantissa, the number scales down to 0
+                    return (sign == 0L) ? 0.0 : -0.0;
+                }
+            } else {
+                // we are really in the case n >= 1024
+                if (exponent == 0) {
+
+                    // the input number is subnormal, normalize it
+                    while ((mantissa >>> 52) != 1) {
+                        mantissa = mantissa << 1;
+                        --scaledExponent;
+                    }
+                    ++scaledExponent;
+                    mantissa = mantissa & 0x000fffffffffffffL;
+
+                    if (scaledExponent < 2047) {
+                        return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa);
+                    } else {
+                        return (sign == 0L) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
+                    }
+
+                } else if (scaledExponent < 2047) {
+                    return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa);
+                } else {
+                    return (sign == 0L) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
+                }
+            }
+        }
+
+        static int getExponent(final double d) {
+            // NaN and Infinite will return 1024 anywho so can use raw bits
+            return (int) ((Double.doubleToRawLongBits(d) >>> 52) & 0x7ff) - 1023;
+        }
+
+        public static AbsNode create() {
+            return ComplexBuiltinsFactory.AbsNodeFactory.create();
         }
+
     }
 
     @Builtin(name = __ADD__, minNumOfPositionalArgs = 2)