[ublas] norm_2 bug when BOOST_UBLAS_SCALED_NORM is on

From: Vyacheslav E. Andrejev (mortituris_at_[hidden])
Date: 2007-03-09 04:04:15

Hello All,

norm_2 can perform zero by zero division if macro BOOST_UBLAS_SCALED_NORM 
is defined and the first component of a vector is zero.

Look at the implementation of vector_norm_2::apply (file $(BOOST)/boost/numeric/ublas/functional.hpp, 
line 437):

        template<class E>
        static BOOST_UBLAS_INLINE
        result_type apply (const vector_expression<E> &e) {
#ifndef BOOST_UBLAS_SCALED_NORM
            real_type t = real_type ();
            typedef typename E::size_type vector_size_type;
            vector_size_type size (e ().size ());
            for (vector_size_type i = 0; i < size; ++ i) {
                real_type u (type_traits<value_type>::norm_2 (e () (i)));
                t +=  u * u;
            }
            return type_traits<real_type>::type_sqrt (t);
#else
            real_type scale = real_type ();
            real_type sum_squares (1);
            size_type size (e ().size ());
            for (size_type i = 0; i < size; ++ i) {
                real_type u (type_traits<value_type>::norm_2 (e () (i)));
                if (scale < u) {
                    real_type v (scale / u);
                    sum_squares = sum_squares * v * v + real_type (1);
                    scale = u;
                } else {
                    real_type v (u / scale);
                    sum_squares += v * v;
                }
            }
            return scale * type_traits<real_type>::type_sqrt (sum_squares);
#endif

>From the code above it is obvious that when e()(0) is zero then on the first 
iteration u == 0 and scale == 0, and the result of line 458

                    real_type v (u / scale);

is undefined.

The result of the overall operation will be QNAN.

To fix it, it is enough to insert
                if (0 == u) continue;
to line 453:

        template<class E>
        static BOOST_UBLAS_INLINE
        result_type apply (const vector_expression<E> &e) {
#ifndef BOOST_UBLAS_SCALED_NORM
            real_type t = real_type ();
            typedef typename E::size_type vector_size_type;
            vector_size_type size (e ().size ());
            for (vector_size_type i = 0; i < size; ++ i) {
                real_type u (type_traits<value_type>::norm_2 (e () (i)));
                t +=  u * u;
            }
            return type_traits<real_type>::type_sqrt (t);
#else
            real_type scale = real_type ();
            real_type sum_squares (1);
            size_type size (e ().size ());
            for (size_type i = 0; i < size; ++ i) {
                real_type u (type_traits<value_type>::norm_2 (e () (i)));
                if (0 == u) continue;
                if (scale < u) {
                    real_type v (scale / u);
                    sum_squares = sum_squares * v * v + real_type (1);
                    scale = u;
                } else {
                    real_type v (u / scale);
                    sum_squares += v * v;
                }
            }
            return scale * type_traits<real_type>::type_sqrt (sum_squares);
#endif

Regards.
--
-- Vyacheslav E. Andrejev
-- System Architect, Optech International, Inc.
-- E-mail: mortituris_at_[hidden]