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