dfrem.c

  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}