# Ublas :

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]
```