1 /* Copyright (C) 2008-2020 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 64 bit multiply results and/or mmed .
61 #include "../arc-ieee-754.h"
63 /* N.B. fp-bit.c does double rounding on denormal numbers. */
89 bic.f 0,0x7ff80000,r12 ; both NaN -> OK
93 #define __divdf3 __divdf3_asm
166 __divdf3_support: /* This label makes debugger output saner. */
170 0x43500000,.Linf_NaN ; large number / denorm -> Inf
175 norm.f r11,r12 ; flag for x/0 -> Inf check
188 bxor.mi DBL1H,DBL1H,31
194 b.d .Lpast_denorm_dbl1
203 norm.f r11,r12 ; flag for 0/x -> 0 check
205 0x43500000, .Lret0_2 ; denorm/large number -> 0
219 b.d .Lpast_denorm_dbl0
223 tst_s DBL0L,DBL0L ; 0/0 -> NaN
224 xor_s DBL1H,DBL1H,DBL0H
225 bclr.eq.f DBL0H,DBL0H,31
227 xor_s DBL0H,DBL0H,DBL1H
234 xor_s DBL1H,DBL1H,DBL0H
238 xor_s DBL0H,DBL0H,DBL1H
241 /* N.B. the spacing between divtab and the sub3 to get its address must
242 be a multiple of 8. */
246 sub3 r10,pcl,61; (.-.Ldivtab) >> 3
247 ld.as r9,[pcl,-124]; [pcl,(-((.-.L7ff00000) >> 2))] ; 0x7ff00000
256 breq.d r7,r9,.Linf_nan_dbl1
262 breq.d r6,r9,.Linf_nan_dbl0
273 sub r4,r4,mhi ; u1.31 inverse, about 30 bit
274 mulu64 r5,r4 ; result fraction highpart
276 add r5,r6, /* wait for immediate */ \
278 mov r11,mhi ; result fraction highpart
279 mulu64 r11,r8 ; u-28.31
280 asl_s DBL1L,DBL1L,9 ; u-29.23:9
282 mov r12,mlo ; u-28.31
283 mulu64 r11,DBL1L ; mhi: u-28.23:9
284 add.cs DBL0L,DBL0L,DBL0L
285 asl_s DBL0L,DBL0L,6 ; u-26.25:7
287 sub_l DBL0L,DBL0L,r12
289 sub r5,DBL0L,mhi ; rest msw ; u-26.31:0
290 mul64 r5,r4 ; mhi: result fraction lowpart
294 bclr r12,r9,20 ; 0x7fe00000
295 brhs.d r6,r12,.Linf_denorm
296 bxor.mi DBL0H,DBL0H,31
304 /* work out exact rounding if we fall through here. */
305 /* We know that the exact result cannot be represented in double
306 precision. Find the mid-point between the two nearest
307 representable values, multiply with the divisor, and check if
308 the result is larger than the dividend. Since we want to know
309 only the sign bit, it is sufficient to calculate only the
310 highpart of the lower 64 bits. */
311 mulu64 r11,DBL1L ; rest before considering r12 in r5 : -mlo
313 asl r12,r9,2 ; u-22.30:2
316 mov r10,mlo ; rest before considering r12 in r5 : -r10
317 mulu64 r12,DBL1L ; mhi: u-51.32
318 asl r5,r5,25 ; s-51.7:25
319 lsr r10,r10,7 ; u-51.30:2
321 mulu64 r12,r8 ; mlo: u-51.31:1
323 add.mi r5,r5,DBL1L ; signed multiply adjust for r12*DBL1L
324 bset r7,r7,0 ; make sure that the result is not zero, and that
325 sub r5,r5,r7 ; a highpart zero appears negative
326 sub.f r5,r5,mlo ; rest msw
327 add.pl.f DBL0L,DBL0L,1
331 .Linf_nan_dbl1: ; 0/Inf -> NaN Inf/Inf -> NaN x/Inf-> 0 x/NaN -> NaN
335 sub_s.ne DBL0L,DBL0L,DBL0L
337 add_s DBL0H,DBL1H,DBL0L
339 bxor.mi DBL0H,DBL0H,31
343 bxor.mi DBL0H,DBL0H,31
353 bxor.mi DBL0H,DBL0H,31
371 or.pnz DBL0H,DBL0H,r7
374 add.f DBL0L,r10,DBL0L
375 add.cs.f DBL0H,DBL0H,1 ; carry clear after this point
377 add.pnz.f DBL0L,DBL0L,1
378 add.cs.f DBL0H,DBL0H,1
380 /* Calculation so far was not conclusive; calculate further rest. */
381 mulu64 r11,DBL1L ; rest before considering r12 in r5 : -mlo
383 asl r5,r5,25 ; s-51.7:25
384 mov r11,mlo ; rest before considering r12 in r5 : -r11
385 mulu64 r12,r8 ; u-51.31:1
386 and r9,DBL0L,1 ; tie-breaker: round to even
387 lsr r11,r11,7 ; u-51.30:2
388 mov DBL1H,mlo ; u-51.31:1
389 mulu64 r12,DBL1L ; u-51.62:2
390 sub.mi r11,r11,DBL1L ; signed multiply adjust for r12*DBL1L
391 add_s DBL1H,DBL1H,r11
392 sub DBL1H,DBL1H,r5 ; -rest msw
393 add_s DBL1H,DBL1H,mhi ; -rest msw
394 add.f 0,DBL1H,DBL1H ; can't ror.f by 32 :-(
397 add.cs.f DBL0L,DBL0L,1
409 bxor.mi DBL0H,DBL0H,31