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
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
|
/* ieee-utils/fp-aix.c
*
* Copyright (C) 1996, 1997, 1998, 1999, 2000 Tim Mooney
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 3 of the License, or (at
* your option) any later version.
*
* This program is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
#include <math.h>
#include <fptrap.h>
#include <float.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_errno.h>
/* GCC uses a common float.h for all platforms, and ignores all vendor
specific code in float.h. That causes problems with fprnd_t,
FP_RND_RZ, FP_RND_RN and FP_RND_RP on AIX.
Personally I consider that a bug, and was advised by a GCC developer
to report this as a bug, which I did.
http://gcc.gnu.org/bugzilla/show_bug.cgi?id=46155
However, there is not universal agreement whether this is a gcc bug, so
the issue needs to be worked around on AIX.
David Kirkby, 26th October 2010. */
#if !HAVE_DECL_FPRND_T
typedef unsigned short fprnd_t;
#endif
#ifndef FP_RND_RZ
#define FP_RND_RZ 0
#endif
#ifndef FP_RND_RN
#define FP_RND_RN 1
#endif
#ifndef FP_RND_RP
#define FP_RND_RP 2
#endif
#ifndef FP_RND_RM
#define FP_RND_RM 3
#endif
int
gsl_ieee_set_mode (int precision, int rounding, int exception_mask)
{
fptrap_t mode = 0 ;
fprnd_t rnd = 0 ;
switch (precision)
{
/* I'm not positive about AIX only supporting default precision rounding,
* but this is the best assumption until it's proven otherwise. */
case GSL_IEEE_SINGLE_PRECISION:
GSL_ERROR ("AIX only supports default precision rounding",
GSL_EUNSUP) ;
break ;
case GSL_IEEE_DOUBLE_PRECISION:
GSL_ERROR ("AIX only supports default precision rounding",
GSL_EUNSUP) ;
break ;
case GSL_IEEE_EXTENDED_PRECISION:
GSL_ERROR ("AIX only supports default precision rounding",
GSL_EUNSUP) ;
break ;
}
switch (rounding)
{
case GSL_IEEE_ROUND_TO_NEAREST:
rnd = FP_RND_RN ;
fp_swap_rnd (rnd) ;
break ;
case GSL_IEEE_ROUND_DOWN:
rnd = FP_RND_RM ;
fp_swap_rnd (rnd) ;
break ;
case GSL_IEEE_ROUND_UP:
rnd = FP_RND_RP ;
fp_swap_rnd (rnd) ;
break ;
case GSL_IEEE_ROUND_TO_ZERO:
rnd = FP_RND_RZ ;
fp_swap_rnd (rnd) ;
break ;
default:
rnd = FP_RND_RN ;
fp_swap_rnd (rnd) ;
}
/* Turn on all the exceptions apart from 'inexact' */
mode = TRP_INVALID | TRP_DIV_BY_ZERO | TRP_OVERFLOW | TRP_UNDERFLOW ;
if (exception_mask & GSL_IEEE_MASK_INVALID)
mode &= ~ TRP_INVALID ;
if (exception_mask & GSL_IEEE_MASK_DENORMALIZED)
{
/* do nothing */
}
else
{
GSL_ERROR ("AIX does not support the denormalized operand exception. "
"Use 'mask-denormalized' to work around this.",
GSL_EUNSUP) ;
}
if (exception_mask & GSL_IEEE_MASK_DIVISION_BY_ZERO)
mode &= ~ TRP_DIV_BY_ZERO ;
if (exception_mask & GSL_IEEE_MASK_OVERFLOW)
mode &= ~ TRP_OVERFLOW ;
if (exception_mask & GSL_IEEE_MASK_UNDERFLOW)
mode &= ~ TRP_UNDERFLOW ;
if (exception_mask & GSL_IEEE_TRAP_INEXACT)
{
mode |= TRP_INEXACT ;
}
else
{
mode &= ~ TRP_INEXACT ;
}
/* AIX appears to require two steps -- first enable floating point traps
* in general... */
fp_trap(FP_TRAP_SYNC);
/* next, enable the traps we're interested in */
fp_enable(mode);
return GSL_SUCCESS ;
}
|
