NFFT
3.3.1alpha
Main Page
Modules
Data Structures
Files
File List
Globals
kernel
util
float.c
1
/*
2
* Copyright (c) 2002, 2016 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
21
R Y(float_property)(
const
float_property p)
22
{
23
const
R base = FLT_RADIX;
24
const
R eps = EPSILON;
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 rmin */
38
{
39
const
INT n = 1 - MIN_EXP;
40
INT i;
41
for
(i = 0; i < n; i++)
42
rmin /= base;
43
}
44
45
/* Compute rmax */
46
{
47
INT i;
48
rmax -= eps;
49
for
(i = 0; i < emax; i++)
50
rmax *= base;
51
}
52
53
/* Compute sfmin */
54
{
55
R small = K(1.0) / rmax;
56
sfmin = rmin;
57
if
(small >= sfmin)
58
sfmin = small * (eps + K(1.0));
59
}
60
61
first = FALSE;
62
}
63
64
if
(p == NFFT_EPSILON)
65
return
eps;
66
else
if
(p == NFFT_SAFE__MIN)
67
return
sfmin;
68
else
if
(p == NFFT_BASE)
69
return
base;
70
else
if
(p == NFFT_PRECISION)
71
return
prec;
72
else
if
(p == NFFT_MANT_DIG)
73
return
t;
74
else
if
(p == NFFT_FLTROUND)
75
return
rnd;
76
else
if
(p == NFFT_E_MIN)
77
return
emin;
78
else
if
(p == NFFT_R_MIN)
79
return
rmin;
80
else
if
(p == NFFT_E_MAX)
81
return
emax;
82
else
if
(p == NFFT_R_MAX)
83
return
rmax;
84
else
85
CK(0
/* cannot happen */
);
86
87
return
K(-1.0);
88
}
/* dlamch_ */
89
91
R Y(prod_real)(R *vec, INT d)
92
{
93
INT t;
94
R prod;
95
96
prod = K(1.0);
97
for
(t = 0; t < d; t++)
98
prod *= vec[t];
99
100
return
prod;
101
}
Generated by
1.8.11