| /* Copyright (C) 2008-2023 Free Software Foundation, Inc. |
| Contributor: Joern Rennecke <joern.rennecke@embecosm.com> |
| on behalf of Synopsys Inc. |
| |
| This file is part of GCC. |
| |
| GCC is free software; you can redistribute it and/or modify it under |
| the terms of the GNU General Public License as published by the Free |
| Software Foundation; either version 3, or (at your option) any later |
| version. |
| |
| GCC is distributed in the hope that it will be useful, but WITHOUT ANY |
| WARRANTY; without even the implied warranty of MERCHANTABILITY or |
| FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
| for more details. |
| |
| Under Section 7 of GPL version 3, you are granted additional |
| permissions described in the GCC Runtime Library Exception, version |
| 3.1, as published by the Free Software Foundation. |
| |
| You should have received a copy of the GNU General Public License and |
| a copy of the GCC Runtime Library Exception along with this program; |
| see the files COPYING3 and COPYING.RUNTIME respectively. If not, see |
| <http://www.gnu.org/licenses/>. */ |
| |
| /* |
| to calculate a := b/x as b*y, with y := 1/x: |
| - x is in the range [1..2) |
| - calculate 15..18 bit inverse y0 using a table of approximating polynoms. |
| Precision is higher for polynoms used to evaluate input with larger |
| value. |
| - Do one newton-raphson iteration step to double the precision, |
| then multiply this with the divisor |
| -> more time to decide if dividend is subnormal |
| - the worst error propagation is on the side of the value range |
| with the least initial defect, thus giving us about 30 bits precision. |
| The truncation error for the either is less than 1 + x/2 ulp. |
| A 31 bit inverse can be simply calculated by using x with implicit 1 |
| and chaining the multiplies. For a 32 bit inverse, we multiply y0^2 |
| with the bare fraction part of x, then add in y0^2 for the implicit |
| 1 of x. |
| - If calculating a 31 bit inverse, the systematic error is less than |
| -1 ulp; likewise, for 32 bit, it is less than -2 ulp. |
| - If we calculate our seed with a 32 bit fraction, we can archive a |
| tentative result strictly better than -2 / +2.5 (1) ulp/128, i.e. we |
| only need to take the step to calculate the 2nd stage rest and |
| rounding adjust 1/32th of the time. However, if we use a 20 bit |
| fraction for the seed, the negative error can exceed -2 ulp/128, (2) |
| thus for a simple add / tst check, we need to do the 2nd stage |
| rest calculation/ rounding adjust 1/16th of the time. |
| (1): The inexactness of the 32 bit inverse contributes an error in the |
| range of (-1 .. +(1+x/2) ) ulp/128. Leaving out the low word of the |
| rest contributes an error < +1/x ulp/128 . In the interval [1,2), |
| x/2 + 1/x <= 1.5 . |
| (2): Unless proven otherwise. I have not actually looked for an |
| example where -2 ulp/128 is exceeded, and my calculations indicate |
| that the excess, if existent, is less than -1/512 ulp. |
| ??? The algorithm is still based on the ARC700 optimized code. |
| Maybe we could make better use of 32x16 bit multiply, or 64 bit multiply |
| results. |
| */ |
| #include "../arc-ieee-754.h" |
| #define mlo acc2 |
| #define mhi acc1 |
| #define mul64(b,c) mullw 0,b,c` machlw 0,b,c |
| #define mulu64(b,c) mululw 0,b,c` machulw 0,b,c |
| |
| /* N.B. fp-bit.c does double rounding on denormal numbers. */ |
| #if 0 /* DEBUG */ |
| .global __divdf3 |
| FUNC(__divdf3) |
| .balign 4 |
| __divdf3: |
| push_s blink |
| push_s r2 |
| push_s r3 |
| push_s r0 |
| bl.d __divdf3_c |
| push_s r1 |
| ld_s r2,[sp,12] |
| ld_s r3,[sp,8] |
| st_s r0,[sp,12] |
| st_s r1,[sp,8] |
| pop_s r1 |
| bl.d __divdf3_asm |
| pop_s r0 |
| pop_s r3 |
| pop_s r2 |
| pop_s blink |
| cmp r0,r2 |
| cmp.eq r1,r3 |
| jeq_s [blink] |
| and r12,DBL0H,DBL1H |
| bic.f 0,0x7ff80000,r12 ; both NaN -> OK |
| jeq_s [blink] |
| bl abort |
| ENDFUNC(__divdf3) |
| #define __divdf3 __divdf3_asm |
| #endif /* DEBUG */ |
| |
| FUNC(__divdf3) |
| .balign 4 |
| .L7ff00000: |
| .long 0x7ff00000 |
| .Ldivtab: |
| .long 0xfc0fffe1 |
| .long 0xf46ffdfb |
| .long 0xed1ffa54 |
| .long 0xe61ff515 |
| .long 0xdf7fee75 |
| .long 0xd91fe680 |
| .long 0xd2ffdd52 |
| .long 0xcd1fd30c |
| .long 0xc77fc7cd |
| .long 0xc21fbbb6 |
| .long 0xbcefaec0 |
| .long 0xb7efa100 |
| .long 0xb32f92bf |
| .long 0xae8f83b7 |
| .long 0xaa2f7467 |
| .long 0xa5ef6479 |
| .long 0xa1cf53fa |
| .long 0x9ddf433e |
| .long 0x9a0f3216 |
| .long 0x965f2091 |
| .long 0x92df0f11 |
| .long 0x8f6efd05 |
| .long 0x8c1eeacc |
| .long 0x88eed876 |
| .long 0x85dec615 |
| .long 0x82eeb3b9 |
| .long 0x800ea10b |
| .long 0x7d3e8e0f |
| .long 0x7a8e7b3f |
| .long 0x77ee6836 |
| .long 0x756e5576 |
| .long 0x72fe4293 |
| .long 0x709e2f93 |
| .long 0x6e4e1c7f |
| .long 0x6c0e095e |
| .long 0x69edf6c5 |
| .long 0x67cde3a5 |
| .long 0x65cdd125 |
| .long 0x63cdbe25 |
| .long 0x61ddab3f |
| .long 0x600d991f |
| .long 0x5e3d868c |
| .long 0x5c6d7384 |
| .long 0x5abd615f |
| .long 0x590d4ecd |
| .long 0x576d3c83 |
| .long 0x55dd2a89 |
| .long 0x545d18e9 |
| .long 0x52dd06e9 |
| .long 0x516cf54e |
| .long 0x4ffce356 |
| .long 0x4e9cd1ce |
| .long 0x4d3cbfec |
| .long 0x4becae86 |
| .long 0x4aac9da4 |
| .long 0x496c8c73 |
| .long 0x483c7bd3 |
| .long 0x470c6ae8 |
| .long 0x45dc59af |
| .long 0x44bc4915 |
| .long 0x43ac3924 |
| .long 0x428c27fb |
| .long 0x418c187a |
| .long 0x407c07bd |
| |
| __divdf3_support: /* This label makes debugger output saner. */ |
| .balign 4 |
| .Ldenorm_dbl1: |
| brge r6, \ |
| 0x43500000,.Linf_NaN ; large number / denorm -> Inf |
| bmsk.f r12,DBL1H,19 |
| mov.eq r12,DBL1L |
| mov.eq DBL1L,0 |
| sub.eq r7,r7,32 |
| norm.f r11,r12 ; flag for x/0 -> Inf check |
| beq_s .Linf_NaN |
| mov.mi r11,0 |
| add.pl r11,r11,1 |
| add_s r12,r12,r12 |
| asl r8,r12,r11 |
| rsub r12,r11,31 |
| lsr r12,DBL1L,r12 |
| tst_s DBL1H,DBL1H |
| or r8,r8,r12 |
| lsr r4,r8,26 |
| lsr DBL1H,r8,12 |
| ld.as r4,[r10,r4] |
| bxor.mi DBL1H,DBL1H,31 |
| sub r11,r11,11 |
| asl DBL1L,DBL1L,r11 |
| sub r11,r11,1 |
| mulu64 (r4,r8) |
| sub r7,r7,r11 |
| b.d .Lpast_denorm_dbl1 |
| asl r7,r7,20 |
| |
| .Linf_NaN: |
| tst_s DBL0L,DBL0L ; 0/0 -> NaN |
| xor_s DBL1H,DBL1H,DBL0H |
| bclr.eq.f DBL0H,DBL0H,31 |
| bmsk DBL0H,DBL1H,30 |
| xor_s DBL0H,DBL0H,DBL1H |
| sub.eq DBL0H,DBL0H,1 |
| mov_s DBL0L,0 |
| j_s.d [blink] |
| or DBL0H,DBL0H,r9 |
| .balign 4 |
| .Lret0_2: |
| xor_s DBL1H,DBL1H,DBL0H |
| mov_s DBL0L,0 |
| bmsk DBL0H,DBL1H,30 |
| j_s.d [blink] |
| xor_s DBL0H,DBL0H,DBL1H |
| .balign 4 |
| .global __divdf3 |
| /* N.B. the spacing between divtab and the sub3 to get its address must |
| be a multiple of 8. */ |
| __divdf3: |
| asl r8,DBL1H,12 |
| lsr r4,r8,26 |
| sub3 r10,pcl,51;(.-.Ldivtab) >> 3 |
| ld.as r9,[pcl,-104]; [pcl,(-((.-.L7ff00000) >> 2))] ; 0x7ff00000 |
| ld.as r4,[r10,r4] |
| lsr r12,DBL1L,20 |
| and.f r7,DBL1H,r9 |
| or r8,r8,r12 |
| mulu64 (r4,r8) |
| beq.d .Ldenorm_dbl1 |
| .Lpast_denorm_dbl1: |
| and.f r6,DBL0H,r9 |
| breq.d r7,r9,.Linf_nan_dbl1 |
| asl r4,r4,12 |
| sub r4,r4,mhi |
| mululw 0,r4,r4 |
| machulw r5,r4,r4 |
| bne.d .Lnormal_dbl0 |
| lsr r8,r8,1 |
| |
| .balign 4 |
| .Ldenorm_dbl0: |
| bmsk.f r12,DBL0H,19 |
| ; wb stall |
| mov.eq r12,DBL0L |
| sub.eq r6,r6,32 |
| norm.f r11,r12 ; flag for 0/x -> 0 check |
| brge r7, \ |
| 0x43500000, .Lret0_2 ; denorm/large number -> 0 |
| beq_s .Lret0_2 |
| mov.mi r11,0 |
| add.pl r11,r11,1 |
| asl r12,r12,r11 |
| sub r6,r6,r11 |
| add.f 0,r6,31 |
| lsr r10,DBL0L,r6 |
| mov.mi r10,0 |
| add r6,r6,11+32 |
| neg.f r11,r6 |
| asl DBL0L,DBL0L,r11 |
| mov.pl DBL0L,0 |
| sub r6,r6,32-1 |
| b.d .Lpast_denorm_dbl0 |
| asl r6,r6,20 |
| |
| .balign 4 |
| .Linf_nan_dbl1: ; 0/Inf -> NaN Inf/Inf -> NaN x/Inf-> 0 x/NaN -> NaN |
| or.f 0,r6,DBL0L |
| cmp.ne r6,r9 |
| not_s DBL0L,DBL1H |
| sub_s.ne DBL0L,DBL0L,DBL0L |
| tst_s DBL0H,DBL0H |
| add_s DBL0H,DBL1H,DBL0L |
| j_s.d [blink] |
| bxor.mi DBL0H,DBL0H,31 |
| |
| .balign 4 |
| .Lnormal_dbl0: |
| breq.d r6,r9,.Linf_nan_dbl0 |
| asl r12,DBL0H,11 |
| lsr r10,DBL0L,21 |
| .Lpast_denorm_dbl0: |
| bset r8,r8,31 |
| mulu64 (r5,r8) |
| add_s r12,r12,r10 |
| bset r5,r12,31 |
| cmp r5,r8 |
| cmp.eq DBL0L,DBL1L |
| lsr.cc r5,r5,1 |
| sub r4,r4,mhi ; u1.31 inverse, about 30 bit |
| mululw 0,r5,r4 |
| machulw r11,r5,r4 ; result fraction highpart |
| lsr r8,r8,2 ; u3.29 |
| add r5,r6, /* wait for immediate */ \ |
| 0x3fe00000 |
| mulu64 (r11,r8) ; u-28.31 |
| asl_s DBL1L,DBL1L,9 ; u-29.23:9 |
| sbc r6,r5,r7 |
| mov r12,mlo ; u-28.31 |
| mulu64 (r11,DBL1L) ; mhi: u-28.23:9 |
| add.cs DBL0L,DBL0L,DBL0L |
| asl_s DBL0L,DBL0L,6 ; u-26.25:7 |
| asl r10,r11,23 |
| sub_l DBL0L,DBL0L,r12 |
| lsr r7,r11,9 |
| sub r5,DBL0L,mhi ; rest msw ; u-26.31:0 |
| mul64 (r5,r4) ; mhi: result fraction lowpart |
| xor.f 0,DBL0H,DBL1H |
| and DBL0H,r6,r9 |
| add_s DBL0H,DBL0H,r7 |
| bclr r12,r9,20 ; 0x7fe00000 |
| brhs.d r6,r12,.Linf_denorm |
| bxor.mi DBL0H,DBL0H,31 |
| add.f r12,mhi,0x11 |
| asr r9,r12,5 |
| sub.mi DBL0H,DBL0H,1 |
| add.f DBL0L,r9,r10 |
| tst r12,0x1c |
| jne.d [blink] |
| add.cs DBL0H,DBL0H,1 |
| /* work out exact rounding if we fall through here. */ |
| /* We know that the exact result cannot be represented in double |
| precision. Find the mid-point between the two nearest |
| representable values, multiply with the divisor, and check if |
| the result is larger than the dividend. Since we want to know |
| only the sign bit, it is sufficient to calculate only the |
| highpart of the lower 64 bits. */ |
| mulu64 (r11,DBL1L) ; rest before considering r12 in r5 : -mlo |
| sub.f DBL0L,DBL0L,1 |
| asl r12,r9,2 ; u-22.30:2 |
| sub.cs DBL0H,DBL0H,1 |
| sub.f r12,r12,2 |
| mov r10,mlo ; rest before considering r12 in r5 : -r10 |
| mululw 0,r12,DBL1L |
| machulw r7,r12,DBL1L ; mhi: u-51.32 |
| asl r5,r5,25 ; s-51.7:25 |
| lsr r10,r10,7 ; u-51.30:2 |
| mulu64 (r12,r8) ; mlo: u-51.31:1 |
| sub r5,r5,r10 |
| add.mi r5,r5,DBL1L ; signed multiply adjust for r12*DBL1L |
| bset r7,r7,0 ; make sure that the result is not zero, and that |
| sub r5,r5,r7 ; a highpart zero appears negative |
| sub.f r5,r5,mlo ; rest msw |
| add.pl.f DBL0L,DBL0L,1 |
| j_s.d [blink] |
| add.eq DBL0H,DBL0H,1 |
| |
| .Linf_nan_dbl0: |
| tst_s DBL1H,DBL1H |
| j_s.d [blink] |
| bxor.mi DBL0H,DBL0H,31 |
| .balign 4 |
| .Linf_denorm: |
| lsr r12,r6,28 |
| brlo.d r12,0xc,.Linf |
| .Ldenorm: |
| asr r6,r6,20 |
| neg r9,r6 |
| mov_s DBL0H,0 |
| brhs.d r9,54,.Lret0 |
| bxor.mi DBL0H,DBL0H,31 |
| add r12,mhi,1 |
| and r12,r12,-4 |
| rsub r7,r6,5 |
| asr r10,r12,28 |
| bmsk r4,r12,27 |
| min r7,r7,31 |
| asr DBL0L,r4,r7 |
| add DBL1H,r11,r10 |
| abs.f r10,r4 |
| sub.mi r10,r10,1 |
| add.f r7,r6,32-5 |
| asl r4,r4,r7 |
| mov.mi r4,r10 |
| add.f r10,r6,23 |
| rsub r7,r6,9 |
| lsr r7,DBL1H,r7 |
| asl r10,DBL1H,r10 |
| or.pnz DBL0H,DBL0H,r7 |
| or.mi r4,r4,r10 |
| mov.mi r10,r7 |
| add.f DBL0L,r10,DBL0L |
| add.cs.f DBL0H,DBL0H,1 ; carry clear after this point |
| bxor.f 0,r4,31 |
| add.pnz.f DBL0L,DBL0L,1 |
| add.cs.f DBL0H,DBL0H,1 |
| jne_s [blink] |
| /* Calculation so far was not conclusive; calculate further rest. */ |
| mulu64 (r11,DBL1L) ; rest before considering r12 in r5 : -mlo |
| asr.f r12,r12,3 |
| asl r5,r5,25 ; s-51.7:25 |
| mov r11,mlo ; rest before considering r12 in r5 : -r11 |
| mulu64 (r12,r8) ; u-51.31:1 |
| and r9,DBL0L,1 ; tie-breaker: round to even |
| lsr r11,r11,7 ; u-51.30:2 |
| mov DBL1H,mlo ; u-51.31:1 |
| mulu64 (r12,DBL1L) ; u-51.62:2 |
| sub.mi r11,r11,DBL1L ; signed multiply adjust for r12*DBL1L |
| add_s DBL1H,DBL1H,r11 |
| sub DBL1H,DBL1H,r5 ; -rest msw |
| add_s DBL1H,DBL1H,mhi ; -rest msw |
| add.f 0,DBL1H,DBL1H ; can't ror.f by 32 :-( |
| tst_s DBL1H,DBL1H |
| cmp.eq mlo,r9 |
| add.cs.f DBL0L,DBL0L,1 |
| j_s.d [blink] |
| add.cs DBL0H,DBL0H,1 |
| |
| .Lret0: |
| /* return +- 0 */ |
| j_s.d [blink] |
| mov_s DBL0L,0 |
| .Linf: |
| mov_s DBL0H,r9 |
| mov_s DBL0L,0 |
| j_s.d [blink] |
| bxor.mi DBL0H,DBL0H,31 |
| ENDFUNC(__divdf3) |
