GNU Linux-libre 4.19.211-gnu1
[releases.git] / arch / parisc / math-emu / sfrem.c
1 /*
2  * Linux/PA-RISC Project (http://www.parisc-linux.org/)
3  *
4  * Floating-point emulation code
5  *  Copyright (C) 2001 Hewlett-Packard (Paul Bame) <bame@debian.org>
6  *
7  *    This program is free software; you can redistribute it and/or modify
8  *    it under the terms of the GNU General Public License as published by
9  *    the Free Software Foundation; either version 2, or (at your option)
10  *    any later version.
11  *
12  *    This program is distributed in the hope that it will be useful,
13  *    but WITHOUT ANY WARRANTY; without even the implied warranty of
14  *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15  *    GNU General Public License for more details.
16  *
17  *    You should have received a copy of the GNU General Public License
18  *    along with this program; if not, write to the Free Software
19  *    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
20  */
21 /*
22  * BEGIN_DESC
23  *
24  *  File:
25  *      @(#)    pa/spmath/sfrem.c               $Revision: 1.1 $
26  *
27  *  Purpose:
28  *      Single Precision Floating-point Remainder
29  *
30  *  External Interfaces:
31  *      sgl_frem(srcptr1,srcptr2,dstptr,status)
32  *
33  *  Internal Interfaces:
34  *
35  *  Theory:
36  *      <<please update with a overview of the operation of this file>>
37  *
38  * END_DESC
39 */
40
41
42
43 #include "float.h"
44 #include "sgl_float.h"
45
46 /*
47  *  Single Precision Floating-point Remainder
48  */
49
50 int
51 sgl_frem (sgl_floating_point * srcptr1, sgl_floating_point * srcptr2,
52           sgl_floating_point * dstptr, unsigned int *status)
53 {
54         register unsigned int opnd1, opnd2, result;
55         register int opnd1_exponent, opnd2_exponent, dest_exponent, stepcount;
56         register boolean roundup = FALSE;
57
58         opnd1 = *srcptr1;
59         opnd2 = *srcptr2;
60         /*
61          * check first operand for NaN's or infinity
62          */
63         if ((opnd1_exponent = Sgl_exponent(opnd1)) == SGL_INFINITY_EXPONENT) {
64                 if (Sgl_iszero_mantissa(opnd1)) {
65                         if (Sgl_isnotnan(opnd2)) {
66                                 /* invalid since first operand is infinity */
67                                 if (Is_invalidtrap_enabled()) 
68                                         return(INVALIDEXCEPTION);
69                                 Set_invalidflag();
70                                 Sgl_makequietnan(result);
71                                 *dstptr = result;
72                                 return(NOEXCEPTION);
73                         }
74                 }
75                 else {
76                         /*
77                          * is NaN; signaling or quiet?
78                          */
79                         if (Sgl_isone_signaling(opnd1)) {
80                                 /* trap if INVALIDTRAP enabled */
81                                 if (Is_invalidtrap_enabled()) 
82                                         return(INVALIDEXCEPTION);
83                                 /* make NaN quiet */
84                                 Set_invalidflag();
85                                 Sgl_set_quiet(opnd1);
86                         }
87                         /* 
88                          * is second operand a signaling NaN? 
89                          */
90                         else if (Sgl_is_signalingnan(opnd2)) {
91                                 /* trap if INVALIDTRAP enabled */
92                                 if (Is_invalidtrap_enabled()) 
93                                         return(INVALIDEXCEPTION);
94                                 /* make NaN quiet */
95                                 Set_invalidflag();
96                                 Sgl_set_quiet(opnd2);
97                                 *dstptr = opnd2;
98                                 return(NOEXCEPTION);
99                         }
100                         /*
101                          * return quiet NaN
102                          */
103                         *dstptr = opnd1;
104                         return(NOEXCEPTION);
105                 }
106         } 
107         /*
108          * check second operand for NaN's or infinity
109          */
110         if ((opnd2_exponent = Sgl_exponent(opnd2)) == SGL_INFINITY_EXPONENT) {
111                 if (Sgl_iszero_mantissa(opnd2)) {
112                         /*
113                          * return first operand
114                          */
115                         *dstptr = opnd1;
116                         return(NOEXCEPTION);
117                 }
118                 /*
119                  * is NaN; signaling or quiet?
120                  */
121                 if (Sgl_isone_signaling(opnd2)) {
122                         /* trap if INVALIDTRAP enabled */
123                         if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
124                         /* make NaN quiet */
125                         Set_invalidflag();
126                         Sgl_set_quiet(opnd2);
127                 }
128                 /*
129                  * return quiet NaN
130                  */
131                 *dstptr = opnd2;
132                 return(NOEXCEPTION);
133         }
134         /*
135          * check second operand for zero
136          */
137         if (Sgl_iszero_exponentmantissa(opnd2)) {
138                 /* invalid since second operand is zero */
139                 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
140                 Set_invalidflag();
141                 Sgl_makequietnan(result);
142                 *dstptr = result;
143                 return(NOEXCEPTION);
144         }
145
146         /* 
147          * get sign of result
148          */
149         result = opnd1;  
150
151         /* 
152          * check for denormalized operands
153          */
154         if (opnd1_exponent == 0) {
155                 /* check for zero */
156                 if (Sgl_iszero_mantissa(opnd1)) {
157                         *dstptr = opnd1;
158                         return(NOEXCEPTION);
159                 }
160                 /* normalize, then continue */
161                 opnd1_exponent = 1;
162                 Sgl_normalize(opnd1,opnd1_exponent);
163         }
164         else {
165                 Sgl_clear_signexponent_set_hidden(opnd1);
166         }
167         if (opnd2_exponent == 0) {
168                 /* normalize, then continue */
169                 opnd2_exponent = 1;
170                 Sgl_normalize(opnd2,opnd2_exponent);
171         }
172         else {
173                 Sgl_clear_signexponent_set_hidden(opnd2);
174         }
175
176         /* find result exponent and divide step loop count */
177         dest_exponent = opnd2_exponent - 1;
178         stepcount = opnd1_exponent - opnd2_exponent;
179
180         /*
181          * check for opnd1/opnd2 < 1
182          */
183         if (stepcount < 0) {
184                 /*
185                  * check for opnd1/opnd2 > 1/2
186                  *
187                  * In this case n will round to 1, so 
188                  *    r = opnd1 - opnd2 
189                  */
190                 if (stepcount == -1 && Sgl_isgreaterthan(opnd1,opnd2)) {
191                         Sgl_all(result) = ~Sgl_all(result);   /* set sign */
192                         /* align opnd2 with opnd1 */
193                         Sgl_leftshiftby1(opnd2); 
194                         Sgl_subtract(opnd2,opnd1,opnd2);
195                         /* now normalize */
196                         while (Sgl_iszero_hidden(opnd2)) {
197                                 Sgl_leftshiftby1(opnd2);
198                                 dest_exponent--;
199                         }
200                         Sgl_set_exponentmantissa(result,opnd2);
201                         goto testforunderflow;
202                 }
203                 /*
204                  * opnd1/opnd2 <= 1/2
205                  *
206                  * In this case n will round to zero, so 
207                  *    r = opnd1
208                  */
209                 Sgl_set_exponentmantissa(result,opnd1);
210                 dest_exponent = opnd1_exponent;
211                 goto testforunderflow;
212         }
213
214         /*
215          * Generate result
216          *
217          * Do iterative subtract until remainder is less than operand 2.
218          */
219         while (stepcount-- > 0 && Sgl_all(opnd1)) {
220                 if (Sgl_isnotlessthan(opnd1,opnd2))
221                         Sgl_subtract(opnd1,opnd2,opnd1);
222                 Sgl_leftshiftby1(opnd1);
223         }
224         /*
225          * Do last subtract, then determine which way to round if remainder 
226          * is exactly 1/2 of opnd2 
227          */
228         if (Sgl_isnotlessthan(opnd1,opnd2)) {
229                 Sgl_subtract(opnd1,opnd2,opnd1);
230                 roundup = TRUE;
231         }
232         if (stepcount > 0 || Sgl_iszero(opnd1)) {
233                 /* division is exact, remainder is zero */
234                 Sgl_setzero_exponentmantissa(result);
235                 *dstptr = result;
236                 return(NOEXCEPTION);
237         }
238
239         /* 
240          * Check for cases where opnd1/opnd2 < n 
241          *
242          * In this case the result's sign will be opposite that of
243          * opnd1.  The mantissa also needs some correction.
244          */
245         Sgl_leftshiftby1(opnd1);
246         if (Sgl_isgreaterthan(opnd1,opnd2)) {
247                 Sgl_invert_sign(result);
248                 Sgl_subtract((opnd2<<1),opnd1,opnd1);
249         }
250         /* check for remainder being exactly 1/2 of opnd2 */
251         else if (Sgl_isequal(opnd1,opnd2) && roundup) { 
252                 Sgl_invert_sign(result);
253         }
254
255         /* normalize result's mantissa */
256         while (Sgl_iszero_hidden(opnd1)) {
257                 dest_exponent--;
258                 Sgl_leftshiftby1(opnd1);
259         }
260         Sgl_set_exponentmantissa(result,opnd1);
261
262         /* 
263          * Test for underflow
264          */
265     testforunderflow:
266         if (dest_exponent <= 0) {
267                 /* trap if UNDERFLOWTRAP enabled */
268                 if (Is_underflowtrap_enabled()) {
269                         /*
270                          * Adjust bias of result
271                          */
272                         Sgl_setwrapped_exponent(result,dest_exponent,unfl);
273                         *dstptr = result;
274                         /* frem is always exact */
275                         return(UNDERFLOWEXCEPTION);
276                 }
277                 /*
278                  * denormalize result or set to signed zero
279                  */
280                 if (dest_exponent >= (1 - SGL_P)) {
281                         Sgl_rightshift_exponentmantissa(result,1-dest_exponent);
282                 }
283                 else {
284                         Sgl_setzero_exponentmantissa(result);
285                 }
286         }
287         else Sgl_set_exponent(result,dest_exponent);
288         *dstptr = result;
289         return(NOEXCEPTION);
290 }