1// SPDX-License-Identifier: GPL-2.0-or-later
2/*
3 * Linux/PA-RISC Project (http://www.parisc-linux.org/)
4 *
5 * Floating-point emulation code
6 *  Copyright (C) 2001 Hewlett-Packard (Paul Bame) <bame@debian.org>
7 */
8/*
9 * BEGIN_DESC
10 *
11 *  File:
12 *	@(#)	pa/spmath/dfrem.c		$Revision: 1.1 $
13 *
14 *  Purpose:
15 *	Double Precision Floating-point Remainder
16 *
17 *  External Interfaces:
18 *	dbl_frem(srcptr1,srcptr2,dstptr,status)
19 *
20 *  Internal Interfaces:
21 *
22 *  Theory:
23 *	<<please update with a overview of the operation of this file>>
24 *
25 * END_DESC
26*/
27
28
29
30#include "float.h"
31#include "dbl_float.h"
32
33/*
34 *  Double Precision Floating-point Remainder
35 */
36
37int
38dbl_frem (dbl_floating_point * srcptr1, dbl_floating_point * srcptr2,
39	  dbl_floating_point * dstptr, unsigned int *status)
40{
41	register unsigned int opnd1p1, opnd1p2, opnd2p1, opnd2p2;
42	register unsigned int resultp1, resultp2;
43	register int opnd1_exponent, opnd2_exponent, dest_exponent, stepcount;
44	register boolean roundup = FALSE;
45
46	Dbl_copyfromptr(srcptr1,opnd1p1,opnd1p2);
47	Dbl_copyfromptr(srcptr2,opnd2p1,opnd2p2);
48	/*
49	 * check first operand for NaN's or infinity
50	 */
51	if ((opnd1_exponent = Dbl_exponent(opnd1p1)) == DBL_INFINITY_EXPONENT) {
52		if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
53			if (Dbl_isnotnan(opnd2p1,opnd2p2)) {
54				/* invalid since first operand is infinity */
55				if (Is_invalidtrap_enabled())
56                                	return(INVALIDEXCEPTION);
57                                Set_invalidflag();
58                                Dbl_makequietnan(resultp1,resultp2);
59				Dbl_copytoptr(resultp1,resultp2,dstptr);
60				return(NOEXCEPTION);
61			}
62		}
63		else {
64                	/*
65                 	 * is NaN; signaling or quiet?
66                 	 */
67                	if (Dbl_isone_signaling(opnd1p1)) {
68                        	/* trap if INVALIDTRAP enabled */
69                        	if (Is_invalidtrap_enabled())
70                            		return(INVALIDEXCEPTION);
71                        	/* make NaN quiet */
72                        	Set_invalidflag();
73                        	Dbl_set_quiet(opnd1p1);
74                	}
75			/*
76			 * is second operand a signaling NaN?
77			 */
78			else if (Dbl_is_signalingnan(opnd2p1)) {
79                        	/* trap if INVALIDTRAP enabled */
80                        	if (Is_invalidtrap_enabled())
81                            		return(INVALIDEXCEPTION);
82                        	/* make NaN quiet */
83                        	Set_invalidflag();
84                        	Dbl_set_quiet(opnd2p1);
85				Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
86                		return(NOEXCEPTION);
87			}
88                	/*
89                 	 * return quiet NaN
90                 	 */
91			Dbl_copytoptr(opnd1p1,opnd1p2,dstptr);
92                	return(NOEXCEPTION);
93		}
94	}
95	/*
96	 * check second operand for NaN's or infinity
97	 */
98	if ((opnd2_exponent = Dbl_exponent(opnd2p1)) == DBL_INFINITY_EXPONENT) {
99		if (Dbl_iszero_mantissa(opnd2p1,opnd2p2)) {
100			/*
101			 * return first operand
102			 */
103			Dbl_copytoptr(opnd1p1,opnd1p2,dstptr);
104			return(NOEXCEPTION);
105		}
106                /*
107                 * is NaN; signaling or quiet?
108                 */
109                if (Dbl_isone_signaling(opnd2p1)) {
110                        /* trap if INVALIDTRAP enabled */
111                        if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
112                        /* make NaN quiet */
113                        Set_invalidflag();
114                        Dbl_set_quiet(opnd2p1);
115                }
116                /*
117                 * return quiet NaN
118                 */
119		Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
120                return(NOEXCEPTION);
121	}
122	/*
123	 * check second operand for zero
124	 */
125	if (Dbl_iszero_exponentmantissa(opnd2p1,opnd2p2)) {
126		/* invalid since second operand is zero */
127		if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
128                Set_invalidflag();
129                Dbl_makequietnan(resultp1,resultp2);
130		Dbl_copytoptr(resultp1,resultp2,dstptr);
131		return(NOEXCEPTION);
132	}
133
134	/*
135	 * get sign of result
136	 */
137	resultp1 = opnd1p1;
138
139	/*
140	 * check for denormalized operands
141	 */
142	if (opnd1_exponent == 0) {
143		/* check for zero */
144		if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
145			Dbl_copytoptr(opnd1p1,opnd1p2,dstptr);
146			return(NOEXCEPTION);
147		}
148		/* normalize, then continue */
149		opnd1_exponent = 1;
150		Dbl_normalize(opnd1p1,opnd1p2,opnd1_exponent);
151	}
152	else {
153		Dbl_clear_signexponent_set_hidden(opnd1p1);
154	}
155	if (opnd2_exponent == 0) {
156		/* normalize, then continue */
157		opnd2_exponent = 1;
158		Dbl_normalize(opnd2p1,opnd2p2,opnd2_exponent);
159	}
160	else {
161		Dbl_clear_signexponent_set_hidden(opnd2p1);
162	}
163
164	/* find result exponent and divide step loop count */
165	dest_exponent = opnd2_exponent - 1;
166	stepcount = opnd1_exponent - opnd2_exponent;
167
168	/*
169	 * check for opnd1/opnd2 < 1
170	 */
171	if (stepcount < 0) {
172		/*
173		 * check for opnd1/opnd2 > 1/2
174		 *
175		 * In this case n will round to 1, so
176		 *    r = opnd1 - opnd2
177		 */
178		if (stepcount == -1 &&
179		    Dbl_isgreaterthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) {
180			/* set sign */
181			Dbl_allp1(resultp1) = ~Dbl_allp1(resultp1);
182			/* align opnd2 with opnd1 */
183			Dbl_leftshiftby1(opnd2p1,opnd2p2);
184			Dbl_subtract(opnd2p1,opnd2p2,opnd1p1,opnd1p2,
185			 opnd2p1,opnd2p2);
186			/* now normalize */
187                	while (Dbl_iszero_hidden(opnd2p1)) {
188                        	Dbl_leftshiftby1(opnd2p1,opnd2p2);
189                        	dest_exponent--;
190			}
191			Dbl_set_exponentmantissa(resultp1,resultp2,opnd2p1,opnd2p2);
192			goto testforunderflow;
193		}
194		/*
195		 * opnd1/opnd2 <= 1/2
196		 *
197		 * In this case n will round to zero, so
198		 *    r = opnd1
199		 */
200		Dbl_set_exponentmantissa(resultp1,resultp2,opnd1p1,opnd1p2);
201		dest_exponent = opnd1_exponent;
202		goto testforunderflow;
203	}
204
205	/*
206	 * Generate result
207	 *
208	 * Do iterative subtract until remainder is less than operand 2.
209	 */
210	while (stepcount-- > 0 && (Dbl_allp1(opnd1p1) || Dbl_allp2(opnd1p2))) {
211		if (Dbl_isnotlessthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) {
212			Dbl_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2,opnd1p1,opnd1p2);
213		}
214		Dbl_leftshiftby1(opnd1p1,opnd1p2);
215	}
216	/*
217	 * Do last subtract, then determine which way to round if remainder
218	 * is exactly 1/2 of opnd2
219	 */
220	if (Dbl_isnotlessthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) {
221		Dbl_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2,opnd1p1,opnd1p2);
222		roundup = TRUE;
223	}
224	if (stepcount > 0 || Dbl_iszero(opnd1p1,opnd1p2)) {
225		/* division is exact, remainder is zero */
226		Dbl_setzero_exponentmantissa(resultp1,resultp2);
227		Dbl_copytoptr(resultp1,resultp2,dstptr);
228		return(NOEXCEPTION);
229	}
230
231	/*
232	 * Check for cases where opnd1/opnd2 < n
233	 *
234	 * In this case the result's sign will be opposite that of
235	 * opnd1.  The mantissa also needs some correction.
236	 */
237	Dbl_leftshiftby1(opnd1p1,opnd1p2);
238	if (Dbl_isgreaterthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) {
239		Dbl_invert_sign(resultp1);
240		Dbl_leftshiftby1(opnd2p1,opnd2p2);
241		Dbl_subtract(opnd2p1,opnd2p2,opnd1p1,opnd1p2,opnd1p1,opnd1p2);
242	}
243	/* check for remainder being exactly 1/2 of opnd2 */
244	else if (Dbl_isequal(opnd1p1,opnd1p2,opnd2p1,opnd2p2) && roundup) {
245		Dbl_invert_sign(resultp1);
246	}
247
248	/* normalize result's mantissa */
249        while (Dbl_iszero_hidden(opnd1p1)) {
250                dest_exponent--;
251                Dbl_leftshiftby1(opnd1p1,opnd1p2);
252        }
253	Dbl_set_exponentmantissa(resultp1,resultp2,opnd1p1,opnd1p2);
254
255        /*
256         * Test for underflow
257         */
258    testforunderflow:
259	if (dest_exponent <= 0) {
260                /* trap if UNDERFLOWTRAP enabled */
261                if (Is_underflowtrap_enabled()) {
262                        /*
263                         * Adjust bias of result
264                         */
265                        Dbl_setwrapped_exponent(resultp1,dest_exponent,unfl);
266			/* frem is always exact */
267			Dbl_copytoptr(resultp1,resultp2,dstptr);
268			return(UNDERFLOWEXCEPTION);
269                }
270                /*
271                 * denormalize result or set to signed zero
272                 */
273                if (dest_exponent >= (1 - DBL_P)) {
274			Dbl_rightshift_exponentmantissa(resultp1,resultp2,
275			 1-dest_exponent);
276                }
277                else {
278			Dbl_setzero_exponentmantissa(resultp1,resultp2);
279		}
280	}
281	else Dbl_set_exponent(resultp1,dest_exponent);
282	Dbl_copytoptr(resultp1,resultp2,dstptr);
283	return(NOEXCEPTION);
284}
285