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
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
|
/* ieee-utils/fp-tru64.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 2 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.
*/
/*
* Under Compaq's Unix with the silly name, read the man pages for read_rnd,
* write_rnd, and ieee(3) for more information on the functions used here.
*
* Note that enabling control of dynamic rounding mode (via write_rnd) requires
* that you pass a special flag to your C compiler. For Compaq's C compiler
* the flag is `-fprm d', for gcc it's `-mfp-rounding-mode=d'.
*
* Enabling the trap control (via ieee_set_fp_control) also requires a
* flag be passed to the C compiler. The flag for Compaq's C compiler
* is `-ieee' and for gcc it's `-mieee'.
* We have not implemented the `inexact' case, since it is rarely used
* and requires the library being built with an additional compiler
* flag that can degrade performance for everything else. If you need
* to add support for `inexact' the relevant flag for Compaq's
* compiler is `-ieee_with_inexact', and the flag for gcc is
* `-mieee-with-inexact'.
*
* Problem have been reported with the "fixed" float.h installed with
* gcc-2.95 lacking some of the definitions in the system float.h (the
* symptoms are errors like: `FP_RND_RN' undeclared). To work around
* this we can include the system float.h before the gcc version, e.g.
*
* #include "/usr/include/float.h"
* #include <float.h>
*/
#include <float.h>
#ifndef FP_RND_RN
# undef _FLOAT_H_
# include "/usr/include/float.h"
# undef _FLOAT_H_
# include <float.h>
#endif
#include <machine/fpu.h>
#include <stdio.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_errno.h>
int
gsl_ieee_set_mode (int precision, int rounding, int exception_mask)
{
unsigned long int mode = 0 ;
unsigned int rnd = 0 ;
/* I'm actually not completely sure that the alpha only supports default
* precisions rounding, but I couldn't find any information regarding this, so
* it seems safe to assume this for now until it's proven otherwise.
*/
switch (precision)
{
case GSL_IEEE_SINGLE_PRECISION:
GSL_ERROR ("Tru64 Unix on the alpha only supports default precision rounding",
GSL_EUNSUP) ;
break ;
case GSL_IEEE_DOUBLE_PRECISION:
GSL_ERROR ("Tru64 Unix on the alpha only supports default precision rounding",
GSL_EUNSUP) ;
break ;
case GSL_IEEE_EXTENDED_PRECISION:
GSL_ERROR ("Tru64 Unix on the alpha only supports default precision rounding",
GSL_EUNSUP) ;
break ;
}
switch (rounding)
{
case GSL_IEEE_ROUND_TO_NEAREST:
rnd = FP_RND_RN ;
write_rnd (rnd) ;
break ;
case GSL_IEEE_ROUND_DOWN:
rnd = FP_RND_RM ;
write_rnd (rnd) ;
break ;
case GSL_IEEE_ROUND_UP:
rnd = FP_RND_RP ;
write_rnd (rnd) ;
break ;
case GSL_IEEE_ROUND_TO_ZERO:
rnd = FP_RND_RZ ;
write_rnd (rnd) ;
break ;
default:
rnd = FP_RND_RN ;
write_rnd (rnd) ;
}
/* Turn on all the exceptions apart from 'inexact' */
/* from the ieee(3) man page:
* IEEE_TRAP_ENABLE_INV -> Invalid operation
* IEEE_TRAP_ENABLE_DZE -> Divide by 0
* IEEE_TRAP_ENABLE_OVF -> Overflow
* IEEE_TRAP_ENABLE_UNF -> Underflow
* IEEE_TRAP_ENABLE_INE -> Inexact (requires special option to C compiler)
* IEEE_TRAP_ENABLE_DNO -> denormal operand
* Note: IEEE_TRAP_ENABLE_DNO is not supported on OSF 3.x or Digital Unix
* 4.0 - 4.0d(?).
* IEEE_TRAP_ENABLE_MASK -> mask of all the trap enables
* IEEE_MAP_DMZ -> map denormal inputs to zero
* IEEE_MAP_UMZ -> map underflow results to zero
*/
mode = IEEE_TRAP_ENABLE_INV | IEEE_TRAP_ENABLE_DZE | IEEE_TRAP_ENABLE_OVF
| IEEE_TRAP_ENABLE_UNF ;
if (exception_mask & GSL_IEEE_MASK_INVALID)
mode &= ~ IEEE_TRAP_ENABLE_INV ;
if (exception_mask & GSL_IEEE_MASK_DENORMALIZED)
{
#ifdef IEEE_TRAP_ENABLE_DNO
mode &= ~ IEEE_TRAP_ENABLE_DNO ;
#else
GSL_ERROR ("Sorry, this version of Digital Unix does not support denormalized operands", GSL_EUNSUP) ;
#endif
}
if (exception_mask & GSL_IEEE_MASK_DIVISION_BY_ZERO)
mode &= ~ IEEE_TRAP_ENABLE_DZE ;
if (exception_mask & GSL_IEEE_MASK_OVERFLOW)
mode &= ~ IEEE_TRAP_ENABLE_OVF ;
if (exception_mask & GSL_IEEE_MASK_UNDERFLOW)
mode &= ~ IEEE_TRAP_ENABLE_UNF ;
if (exception_mask & GSL_IEEE_TRAP_INEXACT)
{
/* To implement this would require a special flag to the C
compiler which can cause degraded performance */
GSL_ERROR ("Sorry, GSL does not implement trap-inexact for Tru64 Unix on the alpha - see fp-tru64.c for details", GSL_EUNSUP) ;
/* In case you need to add it, the appropriate line would be
*
* mode |= IEEE_TRAP_ENABLE_INE ;
*
*/
}
else
{
mode &= ~ IEEE_TRAP_ENABLE_INE ;
}
ieee_set_fp_control (mode) ;
return GSL_SUCCESS ;
}
|
