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 to calculate a := b/x as b*y, with y := 1/x:
28 - x is in the range [1..2)
29 - calculate 15..18 bit inverse y0 using a table of approximating polynoms.
30 Precision is higher for polynoms used to evaluate input with larger
32 - Do one newton-raphson iteration step to double the precision,
33 then multiply this with the divisor
34 -> more time to decide if dividend is subnormal
35 - the worst error propagation is on the side of the value range
36 with the least initial defect, thus giving us about 30 bits precision.
37 The truncation error for the either is less than 1 + x/2 ulp.
38 A 31 bit inverse can be simply calculated by using x with implicit 1
39 and chaining the multiplies. For a 32 bit inverse, we multiply y0^2
40 with the bare fraction part of x, then add in y0^2 for the implicit
42 - If calculating a 31 bit inverse, the systematic error is less than
43 -1 ulp; likewise, for 32 bit, it is less than -2 ulp.
44 - If we calculate our seed with a 32 bit fraction, we can archive a
45 tentative result strictly better than -2 / +2.5 (1) ulp/128, i.e. we
46 only need to take the step to calculate the 2nd stage rest and
47 rounding adjust 1/32th of the time. However, if we use a 20 bit
48 fraction for the seed, the negative error can exceed -2 ulp/128, (2)
49 thus for a simple add / tst check, we need to do the 2nd stage
50 rest calculation/ rounding adjust 1/16th of the time.
51 (1): The inexactness of the 32 bit inverse contributes an error in the
52 range of (-1 .. +(1+x/2) ) ulp/128. Leaving out the low word of the
53 rest contributes an error < +1/x ulp/128 . In the interval [1,2),
55 (2): Unless proven otherwise. I have not actually looked for an
56 example where -2 ulp/128 is exceeded, and my calculations indicate
57 that the excess, if existent, is less than -1/512 ulp.
58 ??? The algorithm is still based on the ARC700 optimized code.
59 Maybe we could make better use of 32x16 bit multiply, or 64 bit multiply
62 #include "../arc-ieee-754.h"
65 #define mul64(b,c) mullw 0,b,c` machlw 0,b,c
66 #define mulu64(b,c) mululw 0,b,c` machulw 0,b,c
68 /* N.B. fp-bit.c does double rounding on denormal numbers. */
94 bic.f 0,0x7ff80000,r12 ; both NaN -> OK
98 #define __divdf3 __divdf3_asm
171 __divdf3_support: /* This label makes debugger output saner. */
175 0x43500000,.Linf_NaN ; large number / denorm -> Inf
180 norm.f r11,r12 ; flag for x/0 -> Inf check
193 bxor.mi DBL1H,DBL1H,31
199 b.d .Lpast_denorm_dbl1
203 tst_s DBL0L,DBL0L ; 0/0 -> NaN
204 xor_s DBL1H,DBL1H,DBL0H
205 bclr.eq.f DBL0H,DBL0H,31
207 xor_s DBL0H,DBL0H,DBL1H
214 xor_s DBL1H,DBL1H,DBL0H
218 xor_s DBL0H,DBL0H,DBL1H
221 /* N.B. the spacing between divtab and the sub3 to get its address must
222 be a multiple of 8. */
226 sub3 r10,pcl,51;(.-.Ldivtab) >> 3
227 ld.as r9,[pcl,-104]; [pcl,(-((.-.L7ff00000) >> 2))] ; 0x7ff00000
236 breq.d r7,r9,.Linf_nan_dbl1
250 norm.f r11,r12 ; flag for 0/x -> 0 check
252 0x43500000, .Lret0_2 ; denorm/large number -> 0
266 b.d .Lpast_denorm_dbl0
270 .Linf_nan_dbl1: ; 0/Inf -> NaN Inf/Inf -> NaN x/Inf-> 0 x/NaN -> NaN
274 sub_s.ne DBL0L,DBL0L,DBL0L
276 add_s DBL0H,DBL1H,DBL0L
278 bxor.mi DBL0H,DBL0H,31
282 breq.d r6,r9,.Linf_nan_dbl0
293 sub r4,r4,mhi ; u1.31 inverse, about 30 bit
295 machulw r11,r5,r4 ; result fraction highpart
297 add r5,r6, /* wait for immediate */ \
299 mulu64 (r11,r8) ; u-28.31
300 asl_s DBL1L,DBL1L,9 ; u-29.23:9
302 mov r12,mlo ; u-28.31
303 mulu64 (r11,DBL1L) ; mhi: u-28.23:9
304 add.cs DBL0L,DBL0L,DBL0L
305 asl_s DBL0L,DBL0L,6 ; u-26.25:7
307 sub_l DBL0L,DBL0L,r12
309 sub r5,DBL0L,mhi ; rest msw ; u-26.31:0
310 mul64 (r5,r4) ; mhi: result fraction lowpart
314 bclr r12,r9,20 ; 0x7fe00000
315 brhs.d r6,r12,.Linf_denorm
316 bxor.mi DBL0H,DBL0H,31
324 /* work out exact rounding if we fall through here. */
325 /* We know that the exact result cannot be represented in double
326 precision. Find the mid-point between the two nearest
327 representable values, multiply with the divisor, and check if
328 the result is larger than the dividend. Since we want to know
329 only the sign bit, it is sufficient to calculate only the
330 highpart of the lower 64 bits. */
331 mulu64 (r11,DBL1L) ; rest before considering r12 in r5 : -mlo
333 asl r12,r9,2 ; u-22.30:2
336 mov r10,mlo ; rest before considering r12 in r5 : -r10
338 machulw r7,r12,DBL1L ; mhi: u-51.32
339 asl r5,r5,25 ; s-51.7:25
340 lsr r10,r10,7 ; u-51.30:2
341 mulu64 (r12,r8) ; mlo: u-51.31:1
343 add.mi r5,r5,DBL1L ; signed multiply adjust for r12*DBL1L
344 bset r7,r7,0 ; make sure that the result is not zero, and that
345 sub r5,r5,r7 ; a highpart zero appears negative
346 sub.f r5,r5,mlo ; rest msw
347 add.pl.f DBL0L,DBL0L,1
354 bxor.mi DBL0H,DBL0H,31
364 bxor.mi DBL0H,DBL0H,31
382 or.pnz DBL0H,DBL0H,r7
385 add.f DBL0L,r10,DBL0L
386 add.cs.f DBL0H,DBL0H,1 ; carry clear after this point
388 add.pnz.f DBL0L,DBL0L,1
389 add.cs.f DBL0H,DBL0H,1
391 /* Calculation so far was not conclusive; calculate further rest. */
392 mulu64 (r11,DBL1L) ; rest before considering r12 in r5 : -mlo
394 asl r5,r5,25 ; s-51.7:25
395 mov r11,mlo ; rest before considering r12 in r5 : -r11
396 mulu64 (r12,r8) ; u-51.31:1
397 and r9,DBL0L,1 ; tie-breaker: round to even
398 lsr r11,r11,7 ; u-51.30:2
399 mov DBL1H,mlo ; u-51.31:1
400 mulu64 (r12,DBL1L) ; u-51.62:2
401 sub.mi r11,r11,DBL1L ; signed multiply adjust for r12*DBL1L
402 add_s DBL1H,DBL1H,r11
403 sub DBL1H,DBL1H,r5 ; -rest msw
404 add_s DBL1H,DBL1H,mhi ; -rest msw
405 add.f 0,DBL1H,DBL1H ; can't ror.f by 32 :-(
408 add.cs.f DBL0L,DBL0L,1
420 bxor.mi DBL0H,DBL0H,31