| /* Copyright (C) 2008-2026 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. | |
| */ | |
| #include "arc-ieee-754.h" | |
| /* 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) | |
| __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 | |
| MPYHU r5,r4,r8 | |
| sub r7,r7,r11 | |
| asl r4,r4,12 | |
| b.d .Lpast_denorm_dbl1 | |
| asl r7,r7,20 | |
| ; wb stall | |
| .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_NaN ; denorm/large number -> 0 | |
| beq_s .Lret0_NaN | |
| 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 | |
| .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_NaN: | |
| xor_s DBL1H,DBL1H,DBL0H | |
| cmp_s r12,r9 | |
| mov_s DBL0L,0 | |
| bmsk DBL0H,DBL1H,30 | |
| xor_s DBL0H,DBL0H,DBL1H | |
| j_s.d [blink] | |
| sub.hi DBL0H,DBL0H,1 | |
| .Linf_nan_dbl1: ; Inf/Inf -> NaN x/Inf-> 0 x/NaN -> NaN | |
| not_s DBL0L,DBL1H | |
| cmp r6,r9 | |
| sub_s.ne DBL0L,DBL0L,DBL0L | |
| tst_s DBL0H,DBL0H | |
| add_s DBL0H,DBL1H,DBL0L | |
| j_s.d [blink] | |
| bxor.mi DBL0H,DBL0H,31 | |
| .Linf_nan_dbl0: | |
| tst_s DBL1H,DBL1H | |
| j_s.d [blink] | |
| bxor.mi DBL0H,DBL0H,31 | |
| .balign 4 | |
| .global __divdf3 | |
| /* N.B. the spacing between divtab and the add3 to get its address must | |
| be a multiple of 8. */ | |
| __divdf3: | |
| asl r8,DBL1H,12 | |
| lsr r12,DBL1L,20 | |
| lsr r4,r8,26 | |
| #if defined (__ARCHS__) || defined (__ARCEM__) | |
| add3 r10,pcl,60 ; (.Ldivtab-.) >> 3 | |
| #else | |
| add3 r10,pcl,59 ; (.Ldivtab-.) >> 3 | |
| #endif | |
| ld.as r4,[r10,r4] | |
| #if defined (__ARCHS__) || defined (__ARCEM__) | |
| ld.as r9,[pcl,182]; [pcl,(-((.-.L7ff00000) >> 2))] ; 0x7ff00000 | |
| #else | |
| ld.as r9,[pcl,180]; [pcl,(-((.-.L7ff00000) >> 2))] ; 0x7ff00000 | |
| #endif | |
| or r8,r8,r12 | |
| MPYHU r5,r4,r8 | |
| and.f r7,DBL1H,r9 | |
| asl r4,r4,12 ; having the asl here is a concession to the XMAC pipeline. | |
| beq.d .Ldenorm_dbl1 | |
| and r6,DBL0H,r9 | |
| .Lpast_denorm_dbl1: ; wb stall | |
| sub r4,r4,r5 | |
| MPYHU r5,r4,r4 | |
| breq.d r6,0,.Ldenorm_dbl0 | |
| lsr r8,r8,1 | |
| asl r12,DBL0H,11 | |
| lsr r10,DBL0L,21 | |
| .Lpast_denorm_dbl0: ; wb stall | |
| bset r8,r8,31 | |
| MPYHU r11,r5,r8 | |
| add_s r12,r12,r10 | |
| bset r5,r12,31 | |
| cmp r5,r8 | |
| cmp.eq DBL0L,DBL1L | |
| ; wb stall | |
| lsr.cc r5,r5,1 | |
| sub r4,r4,r11 ; u1.31 inverse, about 30 bit | |
| MPYHU r11,r5,r4 ; result fraction highpart | |
| breq r7,r9,.Linf_nan_dbl1 | |
| lsr r8,r8,2 ; u3.29 | |
| add r5,r6, /* wait for immediate / XMAC wb stall */ \ | |
| 0x3fe00000 | |
| ; wb stall (not for XMAC) | |
| breq r6,r9,.Linf_nan_dbl0 | |
| mpyu r12,r11,r8 ; u-28.31 | |
| asl_s DBL1L,DBL1L,9 ; u-29.23:9 | |
| sbc r6,r5,r7 | |
| ; resource conflict (not for XMAC) | |
| MPYHU r5,r11,DBL1L ; 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 | |
| ; wb stall (before 'and' for XMAC) | |
| lsr r7,r11,9 | |
| sub r5,DBL0L,r5 ; rest msw ; u-26.31:0 | |
| MPYH r12,r5,r4 ; result fraction lowpart | |
| xor.f 0,DBL0H,DBL1H | |
| and DBL0H,r6,r9 | |
| add_s DBL0H,DBL0H,r7 ; (XMAC wb stall) | |
| bxor.mi DBL0H,DBL0H,31 | |
| brhs r6, /* wb stall / wait for immediate */ \ | |
| 0x7fe00000,.Linf_denorm | |
| add.f r12,r12,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. */ | |
| sub.f DBL0L,DBL0L,1 | |
| asl r12,r9,2 ; u-22.30:2 | |
| mpyu r10,r11,DBL1L ; rest before considering r12 in r5 : -r10 | |
| sub.cs DBL0H,DBL0H,1 | |
| sub.f r12,r12,2 | |
| ; resource conflict (not for XMAC) | |
| MPYHU r7,r12,DBL1L ; u-51.32 | |
| asl r5,r5,25 ; s-51.7:25 | |
| lsr r10,r10,7 ; u-51.30:2 | |
| ; resource conflict (not for XMAC) | |
| ; resource conflict (not for XMAC) | |
| mpyu r9,r12,r8 ; 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 | |
| ; wb stall (one earlier for XMAC) | |
| sub r5,r5,r7 ; a highpart zero appears negative | |
| sub.f r5,r5,r9 ; rest msw | |
| add.pl.f DBL0L,DBL0L,1 | |
| j_s.d [blink] | |
| add.eq DBL0H,DBL0H,1 | |
| .balign 4 | |
| .Linf_denorm: | |
| brlo r6,0xc0000000,.Linf | |
| .Ldenorm: | |
| asr r6,r6,20 | |
| neg r9,r6 | |
| mov_s DBL0H,0 | |
| brhs.d r9,54,.Lret0 | |
| bxor.mi DBL0H,DBL0H,31 | |
| add_l r12,r12,1 | |
| and r12,r12,-4 | |
| rsub r7,r6,5 | |
| asr r10,r12,28 | |
| bmsk r4,r12,27 | |
| #if defined (__ARCHS__) || defined (__ARCEM__) | |
| min r7, r7, 31 | |
| asr DBL0L, r4, r7 | |
| #else | |
| asrs DBL0L,r4,r7 | |
| #endif | |
| add DBL1H,r11,r10 | |
| #if defined (__ARCHS__) || defined (__ARCEM__) | |
| abs.f r10, r4 | |
| sub.mi r10, r10, 1 | |
| #endif | |
| add.f r7,r6,32-5 | |
| #ifdef __ARC700__ | |
| abss r10,r4 | |
| #endif | |
| 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_l [blink] | |
| /* Calculation so far was not conclusive; calculate further rest. */ | |
| mpyu r11,r11,DBL1L ; rest before considering r12 in r5 : -r11 | |
| asr.f r12,r12,3 | |
| asl r5,r5,25 ; s-51.7:25 | |
| ; resource conflict (not for XMAC) | |
| mpyu DBL1H,r12,r8 ; u-51.31:1 | |
| and r9,DBL0L,1 ; tie-breaker: round to even | |
| lsr r11,r11,7 ; u-51.30:2 | |
| ; resource conflict (not for XMAC) | |
| MPYHU r8,r12,DBL1L ; u-51.32 | |
| sub.mi r11,r11,DBL1L ; signed multiply adjust for r12*DBL1L | |
| add_s DBL1H,DBL1H,r11 | |
| ; resource conflict (not for XMAC) | |
| ; resource conflict (not for XMAC) | |
| mpyu r12,r12,DBL1L ; u-83.30:2 | |
| sub DBL1H,DBL1H,r5 ; -rest msw | |
| add_s DBL1H,DBL1H,r8 ; -rest msw | |
| add.f 0,DBL1H,DBL1H ; can't ror.f by 32 :-( | |
| ; wb stall (XMAC: Before add.f) | |
| tst_s DBL1H,DBL1H | |
| cmp.eq r12,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 | |
| .balign 4 | |
| .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 | |
| .L7ff00000: | |
| .long 0x7ff00000 | |
| ENDFUNC(__divdf3) |