1/* Copyright (C) 2004-2022 Free Software Foundation, Inc.
2   This file is part of the GNU C Library.
3
4   The GNU C Library is free software; you can redistribute it and/or
5   modify it under the terms of the GNU Lesser General Public
6   License as published by the Free Software Foundation; either
7   version 2.1 of the License, or (at your option) any later version.
8
9   The GNU C Library is distributed in the hope that it will be useful,
10   but WITHOUT ANY WARRANTY; without even the implied warranty of
11   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
12   Lesser General Public License for more details.
13
14   You should have received a copy of the GNU Lesser General Public
15   License along with the GNU C Library.  If not, see
16   <https://www.gnu.org/licenses/>.  */
17
18#include "div_libc.h"
19
20
21/* 64-bit unsigned long divide.  These are not normal C functions.  Argument
22   registers are t10 and t11, the result goes in t12.  Only t12 and AT may be
23   clobbered.
24
25   Theory of operation here is that we can use the FPU divider for virtually
26   all operands that we see: all dividend values between -2**53 and 2**53-1
27   can be computed directly.  Note that divisor values need not be checked
28   against that range because the rounded fp value will be close enough such
29   that the quotient is < 1, which will properly be truncated to zero when we
30   convert back to integer.
31
32   When the dividend is outside the range for which we can compute exact
33   results, we use the fp quotent as an estimate from which we begin refining
34   an exact integral value.  This reduces the number of iterations in the
35   shift-and-subtract loop significantly.
36
37   The FPCR save/restore is due to the fact that the EV6 _will_ set FPCR_INE
38   for cvttq/c even without /sui being set.  It will not, however, properly
39   raise the exception, so we don't have to worry about FPCR_INED being clear
40   and so dying by SIGFPE.  */
41
42	.text
43	.align	4
44	.globl	__divqu
45	.type	__divqu, @funcnoplt
46	.usepv	__divqu, no
47
48	cfi_startproc
49	cfi_return_column (RA)
50__divqu:
51	lda	sp, -FRAME(sp)
52	cfi_def_cfa_offset (FRAME)
53	CALL_MCOUNT
54
55	/* Get the fp divide insn issued as quickly as possible.  After
56	   that's done, we have at least 22 cycles until its results are
57	   ready -- all the time in the world to figure out how we're
58	   going to use the results.  */
59	beq	Y, DIVBYZERO
60
61	stt	$f0, 0(sp)
62	excb
63	stt	$f1, 8(sp)
64	stt	$f3, 48(sp)
65	cfi_rel_offset ($f0, 0)
66	cfi_rel_offset ($f1, 8)
67	cfi_rel_offset ($f3, 48)
68	mf_fpcr	$f3
69
70	_ITOFT2	X, $f0, 16, Y, $f1, 24
71	cvtqt	$f0, $f0
72	cvtqt	$f1, $f1
73
74	blt	X, $x_is_neg
75	divt/c	$f0, $f1, $f0
76
77	/* Check to see if Y was mis-converted as signed value.  */
78	ldt	$f1, 8(sp)
79	blt	Y, $y_is_neg
80
81	/* Check to see if X fit in the double as an exact value.  */
82	srl	X, 53, AT
83	bne	AT, $x_big
84
85	/* If we get here, we're expecting exact results from the division.
86	   Do nothing else besides convert and clean up.  */
87	cvttq/c	$f0, $f0
88	excb
89	mt_fpcr	$f3
90	_FTOIT	$f0, RV, 16
91
92	ldt	$f0, 0(sp)
93	ldt	$f3, 48(sp)
94	lda	sp, FRAME(sp)
95	cfi_remember_state
96	cfi_restore ($f0)
97	cfi_restore ($f1)
98	cfi_restore ($f3)
99	cfi_def_cfa_offset (0)
100	ret	$31, (RA), 1
101
102	.align	4
103	cfi_restore_state
104$x_is_neg:
105	/* If we get here, X is so big that bit 63 is set, which made the
106	   conversion come out negative.  Fix it up lest we not even get
107	   a good estimate.  */
108	ldah	AT, 0x5f80		/* 2**64 as float.  */
109	stt	$f2, 24(sp)
110	cfi_rel_offset ($f2, 24)
111	_ITOFS	AT, $f2, 16
112
113	.align	4
114	addt	$f0, $f2, $f0
115	unop
116	divt/c	$f0, $f1, $f0
117	unop
118
119	/* Ok, we've now the divide issued.  Continue with other checks.  */
120	ldt	$f1, 8(sp)
121	unop
122	ldt	$f2, 24(sp)
123	blt	Y, $y_is_neg
124	cfi_restore ($f1)
125	cfi_restore ($f2)
126	cfi_remember_state	/* for y_is_neg */
127
128	.align	4
129$x_big:
130	/* If we get here, X is large enough that we don't expect exact
131	   results, and neither X nor Y got mis-translated for the fp
132	   division.  Our task is to take the fp result, figure out how
133	   far it's off from the correct result and compute a fixup.  */
134	stq	t0, 16(sp)
135	stq	t1, 24(sp)
136	stq	t2, 32(sp)
137	stq	t3, 40(sp)
138	cfi_rel_offset (t0, 16)
139	cfi_rel_offset (t1, 24)
140	cfi_rel_offset (t2, 32)
141	cfi_rel_offset (t3, 40)
142
143#define Q	RV		/* quotient */
144#define R	t0		/* remainder */
145#define SY	t1		/* scaled Y */
146#define S	t2		/* scalar */
147#define QY	t3		/* Q*Y */
148
149	cvttq/c	$f0, $f0
150	_FTOIT	$f0, Q, 8
151	mulq	Q, Y, QY
152
153	.align	4
154	stq	t4, 8(sp)
155	excb
156	ldt	$f0, 0(sp)
157	mt_fpcr	$f3
158	cfi_rel_offset (t4, 8)
159	cfi_restore ($f0)
160
161	subq	QY, X, R
162	mov	Y, SY
163	mov	1, S
164	bgt	R, $q_high
165
166$q_high_ret:
167	subq	X, QY, R
168	mov	Y, SY
169	mov	1, S
170	bgt	R, $q_low
171
172$q_low_ret:
173	ldq	t4, 8(sp)
174	ldq	t0, 16(sp)
175	ldq	t1, 24(sp)
176	ldq	t2, 32(sp)
177
178	ldq	t3, 40(sp)
179	ldt	$f3, 48(sp)
180	lda	sp, FRAME(sp)
181	cfi_remember_state
182	cfi_restore (t0)
183	cfi_restore (t1)
184	cfi_restore (t2)
185	cfi_restore (t3)
186	cfi_restore (t4)
187	cfi_restore ($f3)
188	cfi_def_cfa_offset (0)
189	ret	$31, (RA), 1
190
191	.align	4
192	cfi_restore_state
193	/* The quotient that we computed was too large.  We need to reduce
194	   it by S such that Y*S >= R.  Obviously the closer we get to the
195	   correct value the better, but overshooting high is ok, as we'll
196	   fix that up later.  */
1970:
198	addq	SY, SY, SY
199	addq	S, S, S
200$q_high:
201	cmpult	SY, R, AT
202	bne	AT, 0b
203
204	subq	Q, S, Q
205	unop
206	subq	QY, SY, QY
207	br	$q_high_ret
208
209	.align	4
210	/* The quotient that we computed was too small.  Divide Y by the
211	   current remainder (R) and add that to the existing quotient (Q).
212	   The expectation, of course, is that R is much smaller than X.  */
213	/* Begin with a shift-up loop.  Compute S such that Y*S >= R.  We
214	   already have a copy of Y in SY and the value 1 in S.  */
2150:
216	addq	SY, SY, SY
217	addq	S, S, S
218$q_low:
219	cmpult	SY, R, AT
220	bne	AT, 0b
221
222	/* Shift-down and subtract loop.  Each iteration compares our scaled
223	   Y (SY) with the remainder (R); if SY <= R then X is divisible by
224	   Y's scalar (S) so add it to the quotient (Q).  */
2252:	addq	Q, S, t3
226	srl	S, 1, S
227	cmpule	SY, R, AT
228	subq	R, SY, t4
229
230	cmovne	AT, t3, Q
231	cmovne	AT, t4, R
232	srl	SY, 1, SY
233	bne	S, 2b
234
235	br	$q_low_ret
236
237	.align	4
238	cfi_restore_state
239$y_is_neg:
240	/* If we get here, Y is so big that bit 63 is set.  The results
241	   from the divide will be completely wrong.  Fortunately, the
242	   quotient must be either 0 or 1, so just compute it directly.  */
243	cmpule	Y, X, RV
244	excb
245	mt_fpcr	$f3
246	ldt	$f0, 0(sp)
247	ldt	$f3, 48(sp)
248	lda	sp, FRAME(sp)
249	cfi_restore ($f0)
250	cfi_restore ($f3)
251	cfi_def_cfa_offset (0)
252	ret	$31, (RA), 1
253
254	cfi_endproc
255	.size	__divqu, .-__divqu
256
257	DO_DIVBYZERO
258