xref: /illumos-gate/usr/src/lib/libmvec/common/vis/__vatan.S (revision 55fea89d)

/*
 * CDDL HEADER START
 *
 * The contents of this file are subject to the terms of the
 * Common Development and Distribution License (the "License").
 * You may not use this file except in compliance with the License.
 *
 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
 * or http://www.opensolaris.org/os/licensing.
 * See the License for the specific language governing permissions
 * and limitations under the License.
 *
 * When distributing Covered Code, include this CDDL HEADER in each
 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
 * If applicable, add the following below this CDDL HEADER, with the
 * fields enclosed by brackets "[]" replaced with your own identifying
 * information: Portions Copyright [yyyy] [name of copyright owner]
 *
 * CDDL HEADER END
 */
/*
 * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
 */
/*
 * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
 * Use is subject to license terms.
 */

	.file	"__vatan.S"

#include "libm.h"

	RO_DATA

! following is the C version of the ATAN algorithm
!  #include <math.h>
!  #include <stdio.h>
!  double jkatan(double *x)
!  {
!    double  f, z, ans, ansu, ansl, tmp, poly, conup, conlo, dummy;
!    int index, sign, intf, intz;
!    extern const double __vlibm_TBL_atan1[];
!    long *pf = (long *) &f, *pz = (long *) &z;
!
!  /*    Power series  atan(x) = x + p1*x**3 + p2*x**5 + p3*x**7
!   *    Error =  -3.08254E-18   On the interval  |x| < 1/64 */
!
!  /* define dummy names for readability.  Use parray to help compiler optimize loads */
!  #define p3    parray[0]
!  #define p2    parray[1]
!  #define p1    parray[2]
!  #define soffset      3
!
!    static const double parray[] = {
!     -1.428029046844299722E-01,		/* p[3]		*/
!      1.999999917247000615E-01, 		/* p[2]		*/
!     -3.333333333329292858E-01, 		/* p[1]		*/
!      1.0, 				/* not used for p[0], though		*/
!     -1.0,				/* used to flip sign of answer 		*/
!    };
!
!    f        = *x;				/* fetch argument		*/
!    intf     = pf[0];				/* grab upper half		*/
!    sign     = intf & 0x80000000;			/* sign of argument		*/
!    intf    ^= sign;				/* abs(upper argument)		*/
!    sign     = (unsigned) sign >> 31;		/* sign bit = 0 or 1		*/
!    pf[0]    = intf;
!
!    if( (intf > 0x43600000) || (intf < 0x3e300000) ) /* filter out special cases */
!    {
!      if(  (intf > 0x7ff00000) ||
!        ((intf == 0x7ff00000) &&  (pf[1] !=0)) ) return (*x-*x);/* return NaN if x=NaN*/
!      if( intf < 0x3e300000 ) 			/* avoid underflow for small arg */
!      {
!        dummy = 1.0e37 + f;
!        dummy = dummy;
!        return (*x);
!      }
!      if( intf > 0x43600000 )			/* avoid underflow for big arg  */
!      {
!        index = 2;
!        f     = __vlibm_TBL_atan1[index] + __vlibm_TBL_atan1[index+1];/* pi/2 up + pi/2 low   */
!        f     = parray[soffset + sign] * f;	/* put sign bit on ans		*/
!        return (f);
!      }
!    }
!
!    index    = 0;				/* points to 0,0 in table	*/
!    if (intf > 0x40500000)			/* if(|x| > 64               	*/
!    { f = -1.0/f;
!      index  = 2; 				/* point to pi/2 upper, lower	*/
!    }
!    else if( intf >= 0x3f900000 )		/* if |x| >= (1/64)... 		*/
!    {
!      intz   = (intf + 0x00008000) & 0x7fff0000;/* round arg, keep upper	*/
!      pz[0]  = intz;				/* store as a double (z)	*/
!      pz[1]  = 0;				/* ...lower			*/
!      f      = (f - z)/(1.0 + f*z); 		/* get reduced argument		*/
!      index  = (intz - 0x3f900000) >> 15;	/* (index >> 16) << 1)		*/
!      index += 4;				/* skip over 0,0,pi/2,pi/2	*/
!    }
!    conup    = __vlibm_TBL_atan1[index];	/* upper table 			*/
!    conlo    = __vlibm_TBL_atan1[index+1];	/* lower table 			*/
!    tmp      = f*f;
!    poly     = (f*tmp)*((p3*tmp + p2)*tmp + p1);
!    ansu     = conup + f;     			/* compute atan(f)  upper 	*/
!    ansl     = (((conup - ansu) + f) + poly) + conlo;
!    ans      = ansu + ansl;
!    ans      = parray[soffset + sign] * ans;
!    return ans;
!  }

/* 8 bytes = 1 double f.p. word */
#define WSIZE	8

	.align	32			!align with full D-cache line
.COEFFS:
	.double	0r-1.428029046844299722E-01 	!p[3]
        .double 0r1.999999917247000615E-01 	!p[2]
        .double 0r-3.333333333329292858E-01	!p[1]
        .double 0r-1.0, 			!constant -1.0
	.word   0x00008000,0x0			!for fp rounding of reduced arg
	.word   0x7fff0000,0x0			!for fp truncation
	.word	0x47900000,0			!a number close to 1.0E37
	.word   0x80000000,0x0			!mask for fp sign bit
	.word   0x3f800000,0x0			!1.0/128.0 dummy "safe" argument
	.type	.COEFFS,#object

	ENTRY(__vatan)
	save	%sp,-SA(MINFRAME)-16,%sp
	PIC_SETUP(g5)
	PIC_SET(g5,__vlibm_TBL_atan1,o4)
	PIC_SET(g5,.COEFFS,o0)
/*
   __vatan(int n, double *x, int stridex, double *y, stridey)
   computes y(i) = atan( x(i) ), for 1=1,n.  Stridex, stridey
   are the distance between x and y elements

	%i0    n
	%i1    address of x
	%i2    stride x
	%i3    address of y
	%i4    stride y
*/
	cmp	%i0,0			!if n <=0,
	ble,pn	%icc,.RETURN		!....then do nothing
	sll	%i2,3,%i2		!convert stride to byte count
	sll	%i4,3,%i4		!convert stride to byte count

/*  pre-load constants before beginning main loop */

	ldd	[%o0],%f58		!load p[3]
	mov	2,%i5			!argcount = 3

	ldd	[%o0+WSIZE],%f60	!load p[2]
	add	%fp,STACK_BIAS-8,%l1		!yaddr1 = &dummy
	fzero   %f18			!ansu1 = 0

	ldd	[%o0+2*WSIZE],%f62	!load p[1]
	add	%fp,STACK_BIAS-8,%l2		!yaddr2 = &dummy
	fzero   %f12			!(poly1) = 0

	ldd 	[%o0+3*WSIZE],%f56	!-1.0
	fzero   %f14			!tmp1 = 0

	ldd     [%o0+4*WSIZE],%f52	!load rounding mask
	fzero   %f16			!conup1 = 0

	ldd     [%o0+5*WSIZE],%f54	!load truncation mask
	fzero   %f36			!f1 = 0

	ldd	[%o0+6*WSIZE],%f50	!1.0e37
	fzero   %f38			!f2 = 0

	ldd	[%o0+7*WSIZE],%f32	!mask for sign bit

	ldd	[%o4+2*WSIZE],%f46	!pi/2 upper
	ldd	[%o4+(2*WSIZE+8)],%f48	!pi/2 lower
	sethi	%hi(0x40500000),%l6	!64.0
	sethi	%hi(0x3f900000),%l7	!1/64.0
	mov     0,%l4			!index1 = 0
	mov     0,%l5			!index2 = 0

.MAINLOOP:

    /*--------------------------------------------------------------------------*/
    /*--------------------------------------------------------------------------*/
    /*--------------------------------------------------------------------------*/

.LOOP0:
        deccc   %i0                     !--n
	bneg	1f
        mov	%i1,%o5			!xuse = x (delay slot)

	ba	2f
	nop				!delay slot
1:
	PIC_SET(g5,.COEFFS+8*WSIZE,o5)
	dec	%i5			!argcount--
2:
	sethi	%hi(0x80000000),%o7	!mask for sign bit
/*2 */  sethi	%hi(0x43600000),%o1	!big = 0x43600000,0
	ld	[%o5],%o0		!intf = pf[0] = f upper
	ldd     [%o4+%l5],%f26		!conup2 = __vlibm_TBL_atan1[index2]

	sethi	%hi(0x3e300000),%o2	!small = 0x3e300000,0
/*4 */  andn	%o0,%o7,%o0		!intf = fabs(intf)
	ldd 	[%o5],%f34		!f = *x into f34

	sub	%o1,%o0,%o1		!(-) if intf > big
/*6 */  sub	%o0,%o2,%o2		!(-) if intf < small
	fand 	%f34,%f32,%f40		!sign0 = sign bit
            fmuld   %f38,%f38,%f24      !tmp2= f2*f2

/*7 */  orcc	%o1,%o2,%g0		!(-) if either true
 	bneg,pn  %icc,.SPECIAL0		!if (-) goto special cases below
	fabsd	%f34,%f34		!abs(f)  (delay slot)
	!----------------------


	sethi   %hi(0x8000),%o7		!rounding bit
/*8 */  fpadd32 %f34,%f52,%f0		!intf + 0x00008000 (again)
            faddd   %f26,%f38,%f28      !ansu2 = conup2 + f2

 	add	%o0,%o7,%o0		!intf + 0x00008000  (delay slot)
/*9*/   fand    %f0,%f54,%f0		!pz[0] = intz = (intf + 0x00008000) & 0x7fff0000 (again)
            fmuld   %f58,%f24,%f22      !p[3]*tmp2

/*10 */ sethi   %hi(0x7fff0000),%o7     !mask for rounding argument
	fmuld	%f34,%f0,%f10		!f*z
	fsubd	%f34,%f0,%f20		!f - z
          add     %o4,%l4,%l4           !base addr + index1
          fmuld   %f14,%f12,%f12        !poly1 = (f1*tmp1)*((p3*tmp1 + p2)*tmp1 + p1)
          faddd   %f16,%f36,%f16        !(conup1 - ansu1) + f1

/*12 */ and	%o0,%o7,%o0		!intz = (intf + 0x00008000) & 0x7fff0000
            faddd   %f22,%f60,%f22      !p[3]*tmp2 + p[2]
          ldd     [%l4+WSIZE],%f14      !conlo1 = __vlibm_TBL_atan1[index+1]

/*13 */ sub	%o0,%l7,%o2		!intz - 0x3f900000
	fsubd	%f10,%f56,%f10		!(f*z - (-1.0))
          faddd   %f16,%f12,%f12        !((conup1 - ansu1) + f1) + poly1

	cmp	%o0,%l6			!(|f| > 64)
        ble	.ELSE0			!if(|f| > 64) then
/*15 */ sra	  %o2,15,%o3		!index  = (intz - 0x3f900000) >> 15
	  mov	  2,%o1			!index == 2, point to conup, conlo = pi/2 upper, lower
          ba	  .ENDIF0		!continue
/*16 */   fdivd	  %f56,%f34,%f34	!f = -1.0/f (delay slot)
        .ELSE0:				!else f( |x| >= (1/64))
	cmp	%o0,%l7			!if intf >= 1/64
	bl  	.ENDIF0			!if( |x| >= (1/64) ) then...
	mov	0,%o1			!index == 0 , point to conup,conlo = 0,0
          add	  %o3,4,%o1		!index = index + 4
/*16 */   fdivd	  %f20,%f10,%f34	!f = (f - z)/(1.0 + f*z), reduced argument
	.ENDIF0:

/*17*/  sll	%o1,3,%l3		!index0 = index
	mov	%i3,%l0			!yaddr0 = address of y
          faddd   %f12,%f14,%f12        !ansl1 = (((conup1 - ansu)1 + f1) + poly1) + conlo1
            fmuld   %f22,%f24,%f22      !(p3*tmp2 + p2)*tmp2
            fsubd   %f26,%f28,%f26      !conup2 - ansu2

/*20*/  add	%i1,%i2,%i1		!x     += stridex
        add	%i3,%i4,%i3		!y     += stridey
          faddd   %f18,%f12,%f36        !ans1 = ansu1 + ansl1
            fmuld   %f38,%f24,%f24      !f*tmp2
            faddd   %f22,%f62,%f22      !(p3*tmp2 + p2)*tmp2 + p1

/*23*/    for     %f36,%f42,%f36        !sign(ans1) = sign of argument
          std     %f36,[%l1]            !*yaddr1 = ans1
            add     %o4,%l5,%l5         !base addr + index2
            fmuld   %f24,%f22,%f22      !poly2 = (f2*tmp2)*((p3*tmp2 + p2)*tmp2 + p1)
            faddd   %f26,%f38,%f26      !(conup2 - ansu2) + f2
	cmp	%i5,0			!if argcount =0, we are done
	be	.RETURN
	  nop

    /*--------------------------------------------------------------------------*/
    /*--------------------------------------------------------------------------*/
    /*--------------------------------------------------------------------------*/

.LOOP1:
/*25*/  deccc   %i0                     !--n
        bneg    1f
        mov     %i1,%o5                 !xuse = x (delay slot)
        ba      2f
        nop                             !delay slot
1:
	PIC_SET(g5,.COEFFS+8*WSIZE,o5)
        dec     %i5                     !argcount--
2:

/*26*/  sethi	%hi(0x80000000),%o7	!mask for sign bit
	sethi	%hi(0x43600000),%o1	!big = 0x43600000,0
	ld	[%o5],%o0		!intf = pf[0] = f upper

/*28*/  sethi	%hi(0x3e300000),%o2	!small = 0x3e300000,0
	andn	%o0,%o7,%o0		!intf = fabs(intf)
	ldd 	[%o5],%f36		!f = *x into f36

/*30*/  sub	%o1,%o0,%o1		!(-) if intf > big
	sub	%o0,%o2,%o2		!(-) if intf < small
	fand 	%f36,%f32,%f42		!sign1 = sign bit

/*31*/  orcc	%o1,%o2,%g0		!(-) if either true
 	bneg,pn  %icc,.SPECIAL1		!if (-) goto special cases below
	fabsd	%f36,%f36		!abs(f)  (delay slot)
	!----------------------

/*32*/  fpadd32 %f36,%f52,%f0		!intf + 0x00008000 (again)
          ldd     [%l5+WSIZE],%f24      !conlo2 = __vlibm_TBL_atan1[index2+1]

/*33*/  fand    %f0,%f54,%f0		!pz[0] = intz = (intf + 0x00008000) & 0x7fff0000 (again)
        sethi   %hi(0x8000),%o7		!rounding bit
          faddd   %f26,%f22,%f22        !((conup2 - ansu2) + f2) + poly2

/*34*/  add	%o0,%o7,%o0		!intf + 0x00008000  (delay slot)
	sethi   %hi(0x7fff0000),%o7	!mask for rounding argument
	fmuld	%f36,%f0,%f10		!f*z
	fsubd	%f36,%f0,%f20		!f - z

/*35*/  and	%o0,%o7,%o0		!intz = (intf + 0x00008000) & 0x7fff0000
          faddd   %f22,%f24,%f22        !ansl2 = (((conup2 - ansu2) + f2) + poly2) + conlo2

/*37*/  sub	%o0,%l7,%o2		!intz - 0x3f900000
	fsubd	%f10,%f56,%f10		!(f*z - (-1.0))
      ldd     [%o4+%l3],%f6             !conup0 = __vlibm_TBL_atan1[index0]

        cmp	%o0,%l6			!(|f| > 64)
	ble	.ELSE1			!if(|f| > 64) then
/*38*/  sra	%o2,15,%o3		!index  = (intz - 0x3f900000) >> 15
          mov	  2,%o1			!index == 2, point to conup, conlo = pi/2 upper, lower
	  ba	  .ENDIF1		!continue
/*40*/    fdivd	  %f56,%f36,%f36	!f = -1.0/f (delay slot)
        .ELSE1:				!else f( |x| >= (1/64))
	cmp	%o0,%l7			!if intf >= 1/64
	bl  	.ENDIF1			!if( |x| >= (1/64) ) then...
	mov	0,%o1			!index == 0 , point to conup,conlo = 0,0
          add	  %o3,4,%o1		!index = index + 4
/*40*/    fdivd	  %f20,%f10,%f36	!f = (f - z)/(1.0 + f*z), reduced argument
	.ENDIF1:

/*41*/sll     %o1,3,%l4		 	!index1 = index
      mov     %i3,%l1			!yaddr1 = address of y
      fmuld   %f34,%f34,%f4             !tmp0= f0*f0
          faddd   %f28,%f22,%f38        !ans2 = ansu2 + ansl2

/*44*/add	%i1,%i2,%i1		!x     += stridex
      add	%i3,%i4,%i3		!y     += stridey
    fmuld   %f58,%f4,%f2                !p[3]*tmp0
    faddd   %f6,%f34,%f8                !ansu0 = conup0 + f0
          for     %f38,%f44,%f38        !sign(ans2) = sign of argument
          std     %f38,[%l2]            !*yaddr2 = ans2
          cmp	%i5,0			!if argcount =0, we are done
	  be	.RETURN
	  nop

    /*--------------------------------------------------------------------------*/
    /*--------------------------------------------------------------------------*/
    /*--------------------------------------------------------------------------*/

.LOOP2:
/*46*/  deccc   %i0                     !--n
        bneg    1f
        mov     %i1,%o5                 !xuse = x (delay slot)
        ba      2f
        nop                             !delay slot
1:
	PIC_SET(g5,.COEFFS+8*WSIZE,o5)
        dec     %i5                     !argcount--
2:

/*47*/  sethi	%hi(0x80000000),%o7	!mask for sign bit
	sethi	%hi(0x43600000),%o1	!big = 0x43600000,0
	ld	[%o5],%o0		!intf = pf[0] = f upper

/*49*/  sethi	%hi(0x3e300000),%o2	!small = 0x3e300000,0
	andn	%o0,%o7,%o0		!intf = fabs(intf)
	ldd 	[%o5],%f38		!f = *x into f38

/*51*/  sub	%o1,%o0,%o1		!(-) if intf > big
	sub	%o0,%o2,%o2		!(-) if intf < small
	fand 	%f38,%f32,%f44		!sign2 = sign bit

/*52*/  orcc	%o1,%o2,%g0		!(-) if either true
 	bneg,pn  %icc,.SPECIAL2		!if (-) goto special cases below
	fabsd	%f38,%f38		!abs(f)  (delay slot)
	!----------------------

/*53*/  fpadd32 %f38,%f52,%f0		!intf + 0x00008000 (again)
      faddd   %f2,%f60,%f2              !p[3]*tmp0 + p[2]

/*54*/  sethi   %hi(0x8000),%o7		!rounding bit
	fand    %f0,%f54,%f0		!pz[0] = intz = (intf + 0x00008000) & 0x7fff0000 (again)

/*55*/  add	%o0,%o7,%o0		!intf + 0x00008000  (delay slot)
	sethi   %hi(0x7fff0000),%o7     !mask for rounding argument
	fmuld	%f38,%f0,%f10		!f*z
	fsubd	%f38,%f0,%f20		!f - z

/*56*/  and	%o0,%o7,%o0		!intz = (intf + 0x00008000) & 0x7fff0000
      fmuld   %f2,%f4,%f2               !(p3*tmp0 + p2)*tmp0
      fsubd   %f6,%f8,%f6               !conup0 - ansu0

/*58*/  sub	%o0,%l7,%o2		!intz - 0x3f900000
	fsubd	%f10,%f56,%f10		!(f*z - (-1.0))
            ldd     [%o4+%l4],%f16      !conup1 = __vlibm_TBL_atan1[index1]

	cmp	%o0,%l6			!(|f| > 64)
	ble	.ELSE2			!if(|f| > 64) then
/*60*/  sra	%o2,15,%o3		!index  = (intz - 0x3f900000) >> 15
	  mov	  2,%o1			!index == 2, point to conup, conlo = pi/2 upper, lower
	  ba	  .ENDIF2		!continue
/*61*/    fdivd	  %f56,%f38,%f38	!f = -1.0/f (delay slot)
        .ELSE2:				!else f( |x| >= (1/64))
	cmp	%o0,%l7			!if intf >= 1/64
	bl  	.ENDIF2			!if( |x| >= (1/64) ) then...
	mov	0,%o1			!index == 0 , point to conup,conlo = 0,0
          add	  %o3,4,%o1		!index = index + 4
/*61*/    fdivd	  %f20,%f10,%f38	!f = (f - z)/(1.0 + f*z), reduced argument
	.ENDIF2:


/*62*/  sll	%o1,3,%l5		!index2 = index
	mov	%i3,%l2			!yaddr2 = address of y
      fmuld   %f34,%f4,%f4              !f0*tmp0
      faddd   %f2,%f62,%f2              !(p3*tmp0 + p2)*tmp0 + p1
            fmuld   %f36,%f36,%f14      !tmp1= f1*f1

/*65*/add     %o4,%l3,%l3               !base addr + index0
      fmuld   %f4,%f2,%f2               !poly0 = (f0*tmp0)*((p3*tmp0 + p2)*tmp0 + p1)
      faddd   %f6,%f34,%f6              !(conup0 - ansu0) + f0
            fmuld   %f58,%f14,%f12      !p[3]*tmp1
            faddd   %f16,%f36,%f18      !ansu1 = conup1 + f1
      ldd     [%l3+WSIZE],%f4           !conlo0 = __vlibm_TBL_atan1[index0+1]

/*68*/  add	%i1,%i2,%i1		!x     += stridex
	add	%i3,%i4,%i3		!y     += stridey
      faddd   %f6,%f2,%f2               !((conup0 - ansu0) + f0) + poly0
            faddd   %f12,%f60,%f12      !p[3]*tmp1 + p[2]

/*71*/faddd   %f2,%f4,%f2               !ansl0 = (((conup0 - ansu)0 + f0) + poly0) + conlo0
            fmuld   %f12,%f14,%f12      !(p3*tmp1 + p2)*tmp1
            fsubd   %f16,%f18,%f16      !conup1 - ansu1

/*74*/faddd   %f8,%f2,%f34              !ans0 = ansu0 + ansl0
          fmuld   %f36,%f14,%f14        !f1*tmp1
          faddd   %f12,%f62,%f12        !(p3*tmp1 + p2)*tmp1 + p1

/*77*/  for     %f34,%f40,%f34          !sign(ans0) = sign of argument
        std     %f34,[%l0]              !*yaddr0 = ans, always gets stored (delay slot)
        cmp	%i5,0			!if argcount =0, we are done
	bg	.MAINLOOP
	  nop

    /*--------------------------------------------------------------------------*/
    /*--------------------------------------------------------------------------*/
    /*--------------------------------------------------------------------------*/

.RETURN:
	ret
	restore	%g0,%g0,%g0

    /*--------------------------------------------------------------------------*/
    /*------------SPECIAL CASE HANDLING FOR LOOP0 ------------------------------*/
    /*--------------------------------------------------------------------------*/

/* at this point
   %i1     x address
   %o0     intf
   %o2     intf - 0x3e300000
   %f34,36,38     f0,f1,f2
   %f40,42,44     sign0,sign1,sign2
*/

 	.align	32				!align on I-cache boundary
.SPECIAL0:
	orcc	%o2,%g0,%g0			!(-) if intf < 0x3e300000
	bpos	1f  				!if >=...continue
	sethi	%hi(0x7ff00000),%g1	!upper word of Inf (we use 64-bit wide int for this)
	  ba	  3f
	  faddd	  %f34,%f50,%f30		!dummy op just to generate exception (delay slot)
1:
	ld	[%o5+4],%o5			!load x lower word
	sllx	%o0,32,%o0			!left justify intf
	sllx	%g1,32,%g1		!left justify Inf
	or	%o0,%o5,%o0			!merge in lower intf
	cmp	%o0,%g1				!if intf > 0x7ff00000 00000000
	ble,pt	%xcc,2f				!pass thru if NaN
	  nop
	  fmuld	  %f34,%f34,%f34		!...... (x*x) trigger invalid exception
	  ba      3f
	  nop
2:
	faddd	%f46,%f48,%f34			!ans = pi/2 upper + pi/2 lower
3:
	add	%i1,%i2,%i1			!x += stridex
	for     %f34,%f40,%f34			!sign(ans) = sign of argument
	std	%f34,[%i3]			!*y = ans
	ba	.LOOP0				!keep looping
	add	%i3,%i4,%i3			!y += stridey (delay slot)

    /*--------------------------------------------------------------------------*/
    /*-----------SPECIAL CASE HANDLING FOR LOOP1 -------------------------------*/
    /*--------------------------------------------------------------------------*/

 	.align	32				!align on I-cache boundary
.SPECIAL1:
	orcc	%o2,%g0,%g0			!(-) if intf < 0x3e300000
	bpos	1f  				!if >=...continue
	sethi	%hi(0x7ff00000),%g1	!upper word of Inf (we use 64-bit wide int for this)
	  ba	  3f
	  faddd	  %f36,%f50,%f30		!dummy op just to generate exception (delay slot)
1:
	ld	[%o5+4],%o5			!load x lower word
	sllx	%o0,32,%o0			!left justify intf
	sllx	%g1,32,%g1		!left justify Inf
	or	%o0,%o5,%o0			!merge in lower intf
	cmp	%o0,%g1				!if intf > 0x7ff00000 00000000
	ble,pt	%xcc,2f				!pass thru if NaN
	  nop
	  fmuld	  %f36,%f36,%f36		!...... (x*x) trigger invalid exception
	  ba      3f
	  nop
2:
	faddd	%f46,%f48,%f36			!ans = pi/2 upper + pi/2 lower
3:
	add	%i1,%i2,%i1			!x += stridex
	for     %f36,%f42,%f36			!sign(ans) = sign of argument
	std	%f36,[%i3]			!*y = ans
	ba      .LOOP1 				!keep looping
	add	%i3,%i4,%i3			!y += stridey (delay slot)

    /*--------------------------------------------------------------------------*/
    /*------------SPECIAL CASE HANDLING FOR LOOP2 ------------------------------*/
    /*--------------------------------------------------------------------------*/

 	.align	32				!align on I-cache boundary
.SPECIAL2:
	orcc	%o2,%g0,%g0			!(-) if intf < 0x3e300000
	bpos	1f  				!if >=...continue
	sethi	%hi(0x7ff00000),%g1	!upper word of Inf (we use 64-bit wide int for this)
	  ba	  3f
	  faddd	  %f38,%f50,%f30		!dummy op just to generate exception (delay slot)
1:
	ld	[%o5+4],%o5			!load x lower word
	sllx	%o0,32,%o0			!left justify intf
	sllx	%g1,32,%g1		!left justify Inf
	or	%o0,%o5,%o0			!merge in lower intf
	cmp	%o0,%g1				!if intf > 0x7ff00000 00000000
	ble,pt	%xcc,2f				!pass thru if NaN
	  nop
	  fmuld	  %f38,%f38,%f38		!...... (x*x) trigger invalid exception
	  ba      3f
	  nop
2:
	faddd	%f46,%f48,%f38			!ans = pi/2 upper + pi/2 lower
3:
	add	%i1,%i2,%i1			!x += stridex
	for     %f38,%f44,%f38			!sign(ans) = sign of argument
	std	%f38,[%i3]			!*y = ans
	ba      .LOOP2 				!keep looping
	add	%i3,%i4,%i3			!y += stridey

    /*--------------------------------------------------------------------------*/
    /*--------------------------------------------------------------------------*/
    /*--------------------------------------------------------------------------*/

	SET_SIZE(__vatan)

!	.ident	"03-20-96 Sparc V9 3-way-unrolled version"