1 /* Copyright (C) 2008-2016 Free Software Foundation, Inc.
2 Contributor: Joern Rennecke <joern.rennecke@embecosm.com>
3 on behalf of Synopsys Inc.
5 This file is part of GCC.
7 GCC is free software; you can redistribute it and/or modify it under
8 the terms of the GNU General Public License as published by the Free
9 Software Foundation; either version 3, or (at your option) any later
12 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
13 WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
17 Under Section 7 of GPL version 3, you are granted additional
18 permissions described in the GCC Runtime Library Exception, version
19 3.1, as published by the Free Software Foundation.
21 You should have received a copy of the GNU General Public License and
22 a copy of the GCC Runtime Library Exception along with this program;
23 see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
24 <http://www.gnu.org/licenses/>. */
27 - calculate 15..18 bit inverse using a table of approximating polynoms.
28 precision is higher for polynoms used to evaluate input with larger
30 - do one newton-raphson iteration step to double the precision,
31 then multiply this with the divisor
32 -> more time to decide if dividend is subnormal
33 - the worst error propagation is on the side of the value range
34 with the least initial defect, thus giving us about 30 bits precision.
36 #include "../arc-ieee-754.h"
39 #define mul64(b,c) mullw 0,b,c` machlw 0,b,c
40 #define mulu64(b,c) mululw 0,b,c` machulw 0,b,c
67 #define __divsf3 __divsf3_asm
140 .global __divsf3_support
143 bclr.f 0,r0,31 ; 0/0 -> NaN
155 /* N.B. the spacing between divtab and the sub3 to get its address must
156 be a multiple of 8. */
158 ld.as r9,[pcl,-9]; [pcl,(-((.-.L7f800000) >> 2))] ; 0x7f800000
159 sub3 r3,pcl,37;(.-.Ldivtab) >> 3
169 breq.d r11,r9,.Linf_nan_fp1
173 breq.d r2,r9,.Linf_nan_fp0
185 add_s r2,r2, /* wait for immediate */ \
187 sub r7,r7,mhi ; u1.31 inverse, about 30 bit
192 bclr r3,r9,23 ; 0x7f000000
193 brhs.d r2,r3,.Linf_denorm
196 add r3,mhi,0x22 ; round to nearest or higher
197 tst r3,0x3c ; check if rounding was unsafe
199 jne.d [blink] ; return if rounding was safe.
201 /* work out exact rounding if we fall through here. */
202 /* We know that the exact result cannot be represented in single
203 precision. Find the mid-point between the two nearest
204 representable values, multiply with the divisor, and check if
205 the result is larger than the dividend. */
209 asr.f 0,r0,1 ; for round-to-even in case this is a denorm
219 bne_s 0f ; inf/inf -> nan
220 brne r2,r9,.Lsigned0 ; x/inf -> 0, but x/nan -> nan
230 /* For denormal results, it is possible that an exact result needs
231 rounding, and thus the round-to-even rule has to come into play. */
233 brlo r2,0xc0000000,.Linf
238 brlo.d r9,25,.Lpast_denorm
240 /* Fall through: return +- 0 */
247 norm.f r12,r6 ; flag for x/0 -> Inf check
253 0x60000000,r0 ; large number / denorm -> Inf
260 bne.d .Lpast_denorm_fp1
265 norm.f r3,r12 ; flag for 0/x -> 0 check
266 bic.ne.f 0,0x60000000,r1 ; denorm/large number -> 0
272 b.d .Lpast_denorm_fp0