NFFT  3.5.3
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 
21 R 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 
96 R 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.