1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
|
/*
-- MAGMA (version 2.9.0) --
Univ. of Tennessee, Knoxville
Univ. of California, Berkeley
Univ. of Colorado, Denver
@date January 2025
*/
#include "magma_internal.h"
// ========================================
// complex support
// --------------------
// propagates inf and nan correctly:
// abs( [xxx, nan] ) == nan, for any xxx
// abs( [nan, xxx] ) == nan, for any xxx
// abs( [inf, xxx] ) == inf, for any xxx except nan
// abs( [xxx, inf] ) == inf, for any xxx except nan
extern "C" double
magma_cabs(magmaDoubleComplex z)
{
double x = fabs( real( z ));
double y = fabs( imag( z ));
double big, little; // "small" reserved in Windows. Ugh.
if ( x > y ) {
big = x;
little = y;
}
else {
big = y;
little = x;
}
if ( big == 0 || isinf(big) ) {
return big + little; // add to propagate nan
}
little /= big;
return big * sqrt( 1 + little*little );
}
// --------------------
extern "C" float
magma_cabsf(magmaFloatComplex z)
{
float x = fabsf( real( z ));
float y = fabsf( imag( z ));
float big, little;
if ( x > y ) {
big = x;
little = y;
}
else {
big = y;
little = x;
}
if ( big == 0 || isinf(big) ) {
return big + little; // add to propagate nan
}
little /= big;
return big * sqrt( 1 + little*little );
}
|
