[Hi Peng, followup to topical group sci.math.num-analysis]

From: Peng Yu <Pe*******@gmail.com>

Date: Tue, 9 Sep 2008 21:24:10 -0700 (PDT)

Hi,

I have the following program which computes roots of a cubic

function.

The solution is sensitive to the type, which is due to the truncation

error. 'long double T' gives three solutions, and 'typedef double T'

gives one solutions. The correct number of solutions should be two, 1

and 2.

I know there is some trick to reduce the chance of under or overflow.

For example, std::abs(z) shall be implemented as

|x| * sqrt(1 + (y/x)*(y/x)) if |x|>|y|

and

|y| * sqrt(1 + (x/y)*(x/y)) if |y|>|x|,

where z is a complex number, and x and y are its real and complex

parts.

I'm wondering what trick can be played to reduce the truncation error

when solving the cubic polynomial equations.

BTW, I use the algorithm is shown at

http://www.hawaii.edu/suremath/jrootsCubic.html
Thanks,

Peng

#include <vector>

#include <complex>

#include <cmath>

#include <iostream>

template <typename T>

std::vector<Troots_of_cubic_function(const T &a2, const T &a1,

const

T &a0) {

const T one = static_cast<T>(1);

const T three = static_cast<T>(3);

T q = one / 3 * a1 - one / 9 * a2 * a2;

T r = one / 6 * (a1 * a2 - three * a0) - one / 27 * std::pow(a2,

3);

T Delta = std::pow(q, 3) + r * r;

std::cout << "Delta = " << Delta << std::endl;

std::vector<Tv;

if(Delta >= T()) {

T s1 = std::pow(r + sqrt(Delta), one/3);

T s2 = std::pow(r - sqrt(Delta), one/3);

v.push_back(s1 + s2 - a2 / 3);

if(Delta == T()) {

v.push_back(-.5 * (s1 + s2) - a2 / 3);

}

}

else {

std::complex<Ttemp = sqrt(std::complex<T>(Delta));

std::complex<Ts1 = std::pow(r + temp, one/3);

std::complex<Ts2 = std::pow(r - temp, one/3);

const T minus_half = - static_cast<T>(1)/2;

v.push_back((s1 + s2 - a2 / 3).real());

v.push_back((minus_half * (s1 + s2) - a2 / 3 + std::complex<T>(0,

sqrt(three)/2) * (s1 - s2)).real());

v.push_back((minus_half * (s1 + s2) - a2 / 3 - std::complex<T>(0,

sqrt(three)/2) * (s1 - s2)).real());

}

return v;

}

int main () {

//typedef long double T;

typedef double T;

const T a2 = -4;

const T a1 = 5;

const T a0 = -2;

std::vector<Tv = roots_of_cubic_function<T>(a2, a1, a0);

std::cout << "Solutions:" << std::endl;

for(std::vector<T>::const_iterator it = v.begin(); it != v.end();

++

it) {

T x = *it;

T f = ((x + a2) * x + a1) * x + a0;

std::cout << x << " " << f << std::endl;

}

}

--

Quidquid latine scriptum est, altum videtur.