ieee-754/arc600-dsp/divdf3.S

/* Copyright (C) 2008-2017 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)