Subversion Repositories planix.SVN

/*

 * hypot -- sqrt(p*p+q*q), but overflows only if the result does.

 * See Cleve Moler and Donald Morrison,

 * ``Replacing Square Roots by Pythagorean Sums,''

 * IBM Journal of Research and Development,

 * Vol. 27, Number 6, pp. 577-581, Nov. 1983

 */

#include <u.h>

#include <libc.h>

double

hypot(double p, double q)

{

	double r, s, pfac;

	if(p < 0)

		p = -p;

	if(q < 0)

		q = -q;

	if(p < q) {

		r = p;

		p = q;

		q = r;

	}

	if(p == 0)

		return 0;

	pfac = p;

	r = q = q/p;

	p = 1;

	for(;;) {

		r *= r;

		s = r+4;

		if(s == 4)

			return p*pfac;

		r /= s;

		p += 2*r*p;

		q *= r;

		r = q/p;

	}

}