001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017package org.apache.commons.numbers.arrays;
018
019/**
020 * Computes the Cartesian norm (2-norm), handling both overflow and underflow.
021 * Translation of the <a href="http://www.netlib.org/minpack">minpack</a>
022 * "enorm" subroutine.
023 */
024public final class SafeNorm {
025    /** Constant. */
026    private static final double R_DWARF = 3.834e-20;
027    /** Constant. */
028    private static final double R_GIANT = 1.304e+19;
029
030    /** Private constructor. */
031    private SafeNorm() {
032        // intentionally empty.
033    }
034
035    /**
036     * @param v Cartesian coordinates.
037     * @return the 2-norm of the vector.
038     */
039    public static double value(double[] v) {
040        double s1 = 0;
041        double s2 = 0;
042        double s3 = 0;
043        double x1max = 0;
044        double x3max = 0;
045        final double floatn = v.length;
046        final double agiant = R_GIANT / floatn;
047        for (int i = 0; i < v.length; i++) {
048            final double xabs = Math.abs(v[i]);
049            if (xabs < R_DWARF || xabs > agiant) {
050                if (xabs > R_DWARF) {
051                    if (xabs > x1max) {
052                        final double r = x1max / xabs;
053                        s1 = 1 + s1 * r * r;
054                        x1max = xabs;
055                    } else {
056                        final double r = xabs / x1max;
057                        s1 += r * r;
058                    }
059                } else {
060                    if (xabs > x3max) {
061                        final double r = x3max / xabs;
062                        s3 = 1 + s3 * r * r;
063                        x3max = xabs;
064                    } else {
065                        if (xabs != 0) {
066                            final double r = xabs / x3max;
067                            s3 += r * r;
068                        }
069                    }
070                }
071            } else {
072                s2 += xabs * xabs;
073            }
074        }
075        double norm;
076        if (s1 != 0) {
077            norm = x1max * Math.sqrt(s1 + (s2 / x1max) / x1max);
078        } else {
079            if (s2 == 0) {
080                norm = x3max * Math.sqrt(s3);
081            } else {
082                if (s2 >= x3max) {
083                    norm = Math.sqrt(s2 * (1 + (x3max / s2) * (x3max * s3)));
084                } else {
085                    norm = Math.sqrt(x3max * ((s2 / x3max) + (x3max * s3)));
086                }
087            }
088        }
089        return norm;
090    }
091}