NFFT 3.5.3alpha
float.c
1/*
2 * Copyright (c) 2002, 2017 Jens Keiner, Stefan Kunis, Daniel Potts
3 *
4 * This program is free software; you can redistribute it and/or modify it under
5 * the terms of the GNU General Public License as published by the Free Software
6 * Foundation; either version 2 of the License, or (at your option) any later
7 * version.
8 *
9 * This program is distributed in the hope that it will be useful, but WITHOUT
10 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
12 * details.
13 *
14 * You should have received a copy of the GNU General Public License along with
15 * this program; if not, write to the Free Software Foundation, Inc., 51
16 * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17 */
18
19#include "infft.h"
20
21R Y(float_property)(const float_property p)
22{
23 const R base = FLT_RADIX;
24 static R eps = K(1.0);
25 const R t = MANT_DIG;
26 const R emin = MIN_EXP;
27 const R emax = MAX_EXP;
28 const R prec = eps * base;
29 static R rmin = K(1.0);
30 static R rmax = K(1.0);
31 const R rnd = FLTROUND;
32 static R sfmin = K(-1.0);
33 static short first = TRUE;
34
35 if (first)
36 {
37 /* Compute eps = 2^(1-MANT_DIG).
38 * The usual definition of EPSILON is too small for double-double arithmetic on PowerPC. */
39 for (INT i=0; i<MANT_DIG-1; i++)
40 eps /= K(2.0);
41
42 /* Compute rmin */
43 {
44 const INT n = 1 - MIN_EXP;
45 INT i;
46 for (i = 0; i < n; i++)
47 rmin /= base;
48 }
49
50 /* Compute rmax */
51 {
52 INT i;
53 rmax -= eps;
54 for (i = 0; i < emax; i++)
55 rmax *= base;
56 }
57
58 /* Compute sfmin */
59 {
60 R small = K(1.0) / rmax;
61 sfmin = rmin;
62 if (small >= sfmin)
63 sfmin = small * (eps + K(1.0));
64 }
65
66 first = FALSE;
67 }
68
69 if (p == NFFT_EPSILON)
70 return eps;
71 else if (p == NFFT_SAFE__MIN)
72 return sfmin;
73 else if (p == NFFT_BASE)
74 return base;
75 else if (p == NFFT_PRECISION)
76 return prec;
77 else if (p == NFFT_MANT_DIG)
78 return t;
79 else if (p == NFFT_FLTROUND)
80 return rnd;
81 else if (p == NFFT_E_MIN)
82 return emin;
83 else if (p == NFFT_R_MIN)
84 return rmin;
85 else if (p == NFFT_E_MAX)
86 return emax;
87 else if (p == NFFT_R_MAX)
88 return rmax;
89 else
90 CK(0 /* cannot happen */);
91
92 return K(-1.0);
93} /* dlamch_ */
94
96R Y(prod_real)(R *vec, INT d)
97{
98 INT t;
99 R prod;
100
101 prod = K(1.0);
102 for (t = 0; t < d; t++)
103 prod *= vec[t];
104
105 return prod;
106}
Internal header file for auxiliary definitions and functions.