Overflow/underflow in quaternion and octanion division operator

[tk-yoshimura]

Overflow or underflow occurs when the divisor is a giant or minute number in quaternion and octanion division.

The following code, which is expected to yield 1, unexpectedly yields NaN.

#include <iostream>
#include <boost/math/quaternion.hpp>

int main(){
    boost::math::quaternion<double> q1(1e-200, 2e-200, 3e-200, 4e-200);
    boost::math::quaternion<double> q2 = q1 / q1;

    std::cout << q2 << std::endl;
}

I consider exponential normalization to be necessary in the following code.
denominator underflows or overflows when rhs is a minute or huge number.

Quaternion:

math/include/boost/math/quaternion.hpp

Lines 398 to 420 in c3afa49

	BOOST_MATH_CXX14_CONSTEXPR quaternion<T> & operator /= (::std::complex<T> const & rhs)
	{
	T ar = rhs.real();
	T br = rhs.imag();
	T denominator = ar*ar+br*br;
	quaternion<T> result((+a*ar + b*br) / denominator, (-a*br + b*ar) / denominator, (+c*ar - d*br) / denominator, (+c*br + d*ar) / denominator);
	swap(result);
	return(*this);
	}
	template<typename X>
	BOOST_MATH_CXX14_CONSTEXPR quaternion<T> & operator /= (quaternion<X> const & rhs)
	{
	T ar = static_cast<T>(rhs.R_component_1());
	T br = static_cast<T>(rhs.R_component_2());
	T cr = static_cast<T>(rhs.R_component_3());
	T dr = static_cast<T>(rhs.R_component_4());
	T denominator = ar*ar+br*br+cr*cr+dr*dr;
	quaternion<T> result((+a*ar+b*br+c*cr+d*dr)/denominator, (-a*br+b*ar-c*dr+d*cr)/denominator, (-a*cr+b*dr+c*ar-d*br)/denominator, (-a*dr-b*cr+c*br+d*ar)/denominator);
	swap(result);
	return(*this);
	}

Octanion:

math/include/boost/math/octonion.hpp

Lines 561 to 654 in c3afa49

	octonion<T> & operator /= (::std::complex<T> const & rhs)
	{
	T ar = rhs.real();
	T br = rhs.imag();
	T denominator = ar*ar+br*br;
	T at = (+a*ar-b*br)/denominator;
	T bt = (-a*br+b*ar)/denominator;
	T ct = (+c*ar-d*br)/denominator;
	T dt = (+c*br+d*ar)/denominator;
	T et = (+e*ar-f*br)/denominator;
	T ft = (+e*br+f*ar)/denominator;
	T gt = (+g*ar+h*br)/denominator;
	T ht = (+g*br+h*ar)/denominator;
	a = at;
	b = bt;
	c = ct;
	d = dt;
	e = et;
	f = ft;
	g = gt;
	h = ht;
	return(*this);
	}
	octonion<T> & operator /= (::boost::math::quaternion<T> const & rhs)
	{
	T ar = rhs.R_component_1();
	T br = rhs.R_component_2();
	T cr = rhs.R_component_2();
	T dr = rhs.R_component_2();
	T denominator = ar*ar+br*br+cr*cr+dr*dr;
	T at = (+a*ar+b*br+c*cr+d*dr)/denominator;
	T bt = (-a*br+b*ar-c*dr+d*cr)/denominator;
	T ct = (-a*cr+b*dr+c*ar-d*br)/denominator;
	T dt = (-a*dr-b*cr+c*br+d*ar)/denominator;
	T et = (+e*ar-f*br-g*cr-h*dr)/denominator;
	T ft = (+e*br+f*ar+g*dr-h*cr)/denominator;
	T gt = (+e*cr-f*dr+g*ar+h*br)/denominator;
	T ht = (+e*dr+f*cr-g*br+h*ar)/denominator;
	a = at;
	b = bt;
	c = ct;
	d = dt;
	e = et;
	f = ft;
	g = gt;
	h = ht;
	return(*this);
	}
	template<typename X>
	octonion<T> & operator /= (octonion<X> const & rhs)
	{
	T ar = static_cast<T>(rhs.R_component_1());
	T br = static_cast<T>(rhs.R_component_2());
	T cr = static_cast<T>(rhs.R_component_3());
	T dr = static_cast<T>(rhs.R_component_4());
	T er = static_cast<T>(rhs.R_component_5());
	T fr = static_cast<T>(rhs.R_component_6());
	T gr = static_cast<T>(rhs.R_component_7());
	T hr = static_cast<T>(rhs.R_component_8());
	T denominator = ar*ar+br*br+cr*cr+dr*dr+er*er+fr*fr+gr*gr+hr*hr;
	T at = (+a*ar+b*br+c*cr+d*dr+e*er+f*fr+g*gr+h*hr)/denominator;
	T bt = (-a*br+b*ar-c*dr+d*cr-e*fr+f*er+g*hr-h*gr)/denominator;
	T ct = (-a*cr+b*dr+c*ar-d*br-e*gr-f*hr+g*er+h*fr)/denominator;
	T dt = (-a*dr-b*cr+c*br+d*ar-e*hr+f*gr-g*fr+h*er)/denominator;
	T et = (-a*er+b*fr+c*gr+d*hr+e*ar-f*br-g*cr-h*dr)/denominator;
	T ft = (-a*fr-b*er+c*hr-d*gr+e*br+f*ar+g*dr-h*cr)/denominator;
	T gt = (-a*gr-b*hr-c*er+d*fr+e*cr-f*dr+g*ar+h*br)/denominator;
	T ht = (-a*hr+b*gr-c*fr-d*er+e*dr+f*cr-g*br+h*ar)/denominator;
	a = at;
	b = bt;
	c = ct;
	d = dt;
	e = et;
	f = ft;
	g = gt;
	h = ht;
	return(*this);
	}

[mborland]

While this should get fixed have you tried Boost.QVM (https://www.boost.org/doc/libs/1_86_0/libs/qvm/doc/html/index.html)? I believe it has a much richer set of support than our quaternions (which is the oldest part of the library)

[tk-yoshimura]

Thank you for sharing about QVM.

Here is how I discovered this problem.
While implementing complex Bessel functions, I encountered a similar problem with the division of the micro-norm.
I verified that boost does not have a similar flaw and reported it because it seems to have the same mistake in quaternion and octonion.