/cln-1.3.2/src/complex/algebraic/cl_LF_hypot.cc

 1// cl_hypot().
 2
 3// General includes.
 4#include "base/cl_sysdep.h"
 5
 6// Specification.
 7#include "complex/cl_C.h"
 8
 9
10// Implementation.
11
12#include "cln/lfloat.h"
13#include "float/lfloat/cl_LF.h"
14#include "float/lfloat/cl_LF_impl.h"
15
16/* For inline version of minusp */
17#include "base/cl_inline.h"
18#include "float/lfloat/elem/cl_LF_minusp.cc"
19
20namespace cln {
21
22ALL_cl_LF_OPERATIONS_SAME_PRECISION()
23
24const cl_LF cl_hypot (const cl_LF& a, const cl_LF& b)
25{
26//  Zuerst a und b auf gleiche Länge bringen: den längeren runden.
27//  a=0.0 -> liefere abs(b).
28//  b=0.0 -> liefere abs(a).
29//  e:=max(exponent(a),exponent(b)).
30//  a':=a/2^e bzw. 0.0 bei Underflowmöglichkeit (beim Skalieren a':=a/2^e
31//      oder beim Quadrieren a'*a':  2*(e-exponent(a))>exp_mid-exp_low-1
32//      d.h. exponent(b)-exponent(a)>floor((exp_mid-exp_low-1)/2) ).
33//  b':=b/2^e bzw. 0.0 bei Underflowmöglichkeit (beim Skalieren b':=b/2^e
34//      oder beim Quadrieren b'*b':  2*(e-exponent(b))>exp_mid-exp_low-1
35//      d.h. exponent(a)-exponent(b)>floor((exp_mid-exp_low-1)/2) ).
36//  c':=a'*a'+b'*b', c':=sqrt(c'), liefere 2^e*c'.
37 {	Mutable(cl_LF,a);
38	Mutable(cl_LF,b);
39	{
40		var uintC a_len = TheLfloat(a)->len;
41		var uintC b_len = TheLfloat(b)->len;
42		if (!(a_len == b_len)) {
43			if (a_len < b_len)
44				b = shorten(b,a_len);
45			else
46				a = shorten(a,b_len);
47		}
48	}
49	var sintE a_exp;
50	var sintE b_exp;
51	{
52		// Exponenten von a holen:
53		var uintE uexp = TheLfloat(a)->expo;
54		if (uexp == 0)
55			// a=0.0 -> liefere (abs b) :
56			return (minusp_inline(b) ? -b : b);
57		a_exp = (sintE)(uexp - LF_exp_mid);
58	}
59	{
60		// Exponenten von b holen:
61		var uintE uexp = TheLfloat(b)->expo;
62		if (uexp == 0)
63			// b=0.0 -> liefere (abs a) :
64			return (minusp_inline(a) ? -a : a);
65		b_exp = (sintE)(uexp - LF_exp_mid);
66	}
67	// Nun a_exp = float_exponent(a), b_exp = float_exponent(b).
68	var sintE e = (a_exp > b_exp ? a_exp : b_exp); // Maximum der Exponenten
69	// a und b durch 2^e dividieren:
70	var cl_LF na = ((b_exp > a_exp) && ((uintE)(b_exp-a_exp) > (uintE)floor(LF_exp_mid-LF_exp_low-1,2)) ? encode_LF0(TheLfloat(a)->len) : scale_float(a,-e));
71	var cl_LF nb = ((a_exp > b_exp) && ((uintE)(a_exp-b_exp) > (uintE)floor(LF_exp_mid-LF_exp_low-1,2)) ? encode_LF0(TheLfloat(b)->len) : scale_float(b,-e));
72	// c' := a'*a'+b'*b' berechnen:
73	var cl_LF nc = square(na) + square(nb);
74	return scale_float(sqrt(nc),e); // c' := sqrt(c'), 2^e*c'
75}}
76
77}  // namespace cln