commons-numbers-arrays/src/main/java/org/apache/commons/numbers/arrays/SafeNorm.java

/*
 * Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The ASF licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
package org.apache.commons.numbers.arrays;

/**
 * Computes the Cartesian norm (2-norm), handling both overflow and underflow.
 * Translation of the <a href="http://www.netlib.org/minpack">minpack</a>
 * "enorm" subroutine.
 */
public final class SafeNorm {
    /** Constant. */
    private static final double R_DWARF = 3.834e-20;
    /** Constant. */
    private static final double R_GIANT = 1.304e+19;

    /** Private constructor. */
    private SafeNorm() {
        // intentionally empty.
    }

    /**
     * @param v Cartesian coordinates.
     * @return the 2-norm of the vector.
     */
    public static double value(double[] v) {
        double s1 = 0;
        double s2 = 0;
        double s3 = 0;
        double x1max = 0;
        double x3max = 0;
        final double floatn = v.length;
        final double agiant = R_GIANT / floatn;
        for (int i = 0; i < v.length; i++) {
            final double xabs = Math.abs(v[i]);
            if (xabs < R_DWARF || xabs > agiant) {
                if (xabs > R_DWARF) {
                    if (xabs > x1max) {
                        final double r = x1max / xabs;
                        s1 = 1 + s1 * r * r;
                        x1max = xabs;
                    } else {
                        final double r = xabs / x1max;
                        s1 += r * r;
                    }
                } else {
                    if (xabs > x3max) {
                        final double r = x3max / xabs;
                        s3 = 1 + s3 * r * r;
                        x3max = xabs;
                    } else {
                        if (xabs != 0) {
                            final double r = xabs / x3max;
                            s3 += r * r;
                        }
                    }
                }
            } else {
                s2 += xabs * xabs;
            }
        }
        double norm;
        if (s1 != 0) {
            norm = x1max * Math.sqrt(s1 + (s2 / x1max) / x1max);
        } else {
            if (s2 == 0) {
                norm = x3max * Math.sqrt(s3);
            } else {
                if (s2 >= x3max) {
                    norm = Math.sqrt(s2 * (1 + (x3max / s2) * (x3max * s3)));
                } else {
                    norm = Math.sqrt(x3max * ((s2 / x3max) + (x3max * s3)));
                }
            }
        }
        return norm;
    }
}