e_log2l.s

glibc 2.9,最新版的C语言库函数

{.mfb
  add r3=16,r2 
  // r=x*y-1
  fms.s1 f6=f41,f8,f1
  (p12) br.cond.spnt SPECIAL_log2l
}
{.mfi
  // load C_1
  ldfe f14=[r2],48
  // RN(x*y)
  fma.s1 f43=f41,f8,f0
  mov r23=0xff;;
}

{.mmi
  // load C_7, C_8
  ldfpd f10,f11=[r3],16
  // load C_3,C_4
  ldfpd f15,f42=[r2],16
  (p8) shr.u r28=r25,63-8;;
}


{.mfi
  // load C_5, C_6
  ldfpd f12,f13=[r3]
  // pseudo-zero ?
  fcmp.eq.s0 p7,p0=f7,f0
  // if first 9 bits after leading 1 are all zero, then p8=1
  cmp.ltu p8,p12=r25,r26
}
{.mfi
  // load C1l
  ldfe f34=[r2],16
  fmerge.se f7=f1,f7
  // get T_index
  and r28=r28,r23;;
}
{.mfi
  // r29=exponent-bias
  sub r29=r29,r27
  // if first 8 bits after leading bit are 0, use polynomial approx. only
  (p8) fms.s1 f6=f7,f1,f1
  // start address of T_low
  add r3=1024+16,r2
}
{.mfi
  // load C_2
  ldfe f35=[r2],16
  // x=1, return 0
  (p6) fma.s0 f8=f0,f0,f0
  // first 8 bits after leading 1 are all ones ?
  cmp.eq p10,p0=r23,r28;;
}

{.mfb
  // if first 8 bits after leading 1 are all ones, use polynomial approx. only
  // add 1 to the exponent additive term, and estimate log2(1-r)
  (p10) add r29=1,r29
  nop.f 0
  (p7) br.cond.spnt LOG2_PSEUDO_ZERO 
}
{.mfi
  // get T_low adress 
  shladd r3=r28,3,r3
  // if first 8 bits after leading 1 are all ones, use polynomial approx. only
  (p10) fms.s1 f6=f7,f36,f1
  // p10 --> p8=1, p12=0
  (p10) cmp.eq p8,p12=r0,r0;;
}

{.mfi
  // get T_high address
  shladd r2=r28,2,r2
  // L(x*y)=x*y-RN(x*y)
  fms.s1 f41=f41,f8,f43
  nop.i 0
}
{.mfi
  // p13=p12
  (p12) cmp.eq.unc p13,p0=r0,r0
  // RtH=RN(x*y)-1  (will eliminate rounding errors in r)
  fms.s1 f43=f43,f1,f1
  nop.i 0;;
}

.pred.rel "mutex",p8,p12
{.mfb
  // load T_high (unless first 9 bits after leading 1 are 0)
  (p12) ldfs f7=[r2]
  // set T_high=0 (if first 9 bits after leading 1 are 0)
  (p8) fma.s1 f7=f0,f0,f0
  // x=1, return
  (p6) br.ret.spnt b0
}
.pred.rel "mutex",p8,p12
{.mfi
  // p12: load T_low
  (p12) ldfd f36=[r3]
  // p8: set T_low=0
  (p8) fma.s1 f36=f0,f0,f0
  (p8) cmp.eq p8,p12=r29,r0;; //nop.i 0;;
}

.pred.rel "mutex",p8,p12
{.mfi
  // f8=expon - bias 
  setf.sig f8=r29
  // general case: 2^{16}+C1*r
  (p12) fma.s1 f33=f6,f14,f32
  nop.i 0
}
{.mfi
  // r26=1
  mov r26=1
  // p8 (mantissa is close to 1, or close to 2): 2^{-8}+C1*r
  (p8) fma.s1 f32=f6,f14,f33
  nop.i 0;;
}

{.mfi
  nop.m 0
  // P78=C_7+C_8*r
  fma.s1 f10=f11,f6,f10
  // r26=2^{63}
  shl r26=r26,63
}
{.mfi
  nop.m 0
  // P34=C_3+r*C_4
  fma.s1 f15=f42,f6,f15
  nop.i 0;;
}
{.mfi
  nop.m 0
  // r2=r*r
  fma.s1 f11=f6,f6,f0
  nop.i 0
}
{.mfi
  nop.m 0
  // P56=C_5+C_6*r
  fma.s1 f13=f13,f6,f12
  nop.i 0;;
}

{.mfi
  nop.m 0
  // Rth-r
  (p13) fms.s1 f43=f43,f1,f6
  nop.i 0
}
{.mfi
  // significand(x)=1 ?
  cmp.eq p0,p6=r25,r26
  // P12=C1l+C_2*r
  fma.s1 f34=f35,f6,f34
  nop.i 0;;
}

.pred.rel "mutex",p8,p12
{.mfi
  nop.m 0
  // p12: C1r=(2^{16}+C1*r)-2^{16}
  (p12) fms.s1 f32=f33,f1,f32
  nop.i 0
}
{.mfi
  nop.m 0
  // p8: C1r=C1*r (double extended)
  (p8) fms.s1 f32=f32,f1,f33
  nop.i 0;;
}

{.mfi
  nop.m 0
  // L(x*y)*C_1+T_low
  (p13) fma.s1 f36=f41,f14,f36
  nop.i 0
}
{.mfi
  nop.m 0
  // P58=P56+r2*P78
  fma.s1 f13=f11,f10,f13
  nop.i 0;;
}
{.mfi
  nop.m 0
  // P14=P12+r2*P34
  fma.s1 f15=f15,f11,f34
  nop.i 0
}
{.mfi
  nop.m 0
  // r4=r2*r2
  fma.s1 f11=f11,f11,f0
  nop.i 0;;
}

{.mfi
  nop.m 0
  // normalize additive term (l=exponent of x)
  fcvt.xf f8=f8
  nop.i 0;;
}


{.mfi
  nop.m 0
  // D=C1*r-C1r
  (p6) fms.s1 f12=f14,f6,f32
  nop.i 0;;
}

{.mfi
  nop.m 0
  // T_low'=(Rth-r)*C1+(L(x*y)*C1+T_low)
  (p13) fma.s1 f36=f43,f14,f36
  nop.i 0;;
}
{.mfi
  nop.m 0
  // P18=P14+r4*P58
  (p6) fma.s1 f13=f11,f13,f15
  nop.i 0;;
}

{.mfi
  nop.m 0
  // add T_high+l
  (p6) fma.s1 f8=f8,f1,f7
  nop.i 0;;
}


{.mfi
  nop.m 0
  // D+T_low
  (p6) fma.s1 f12=f12,f1,f36
  nop.i 0;;
}


{.mfi
  nop.m 0
  // (T_high+l)+C1r
  (p6) fma.s1 f8=f8,f1,f32
  nop.i 0
}
{.mfi
  nop.m 0
  // (D+T_low)+r*P18
  (p6) fma.s1 f13=f13,f6,f12
  nop.i 0;;
}

//{.mfb
//nop.m 0
//mov f8=f36
//fma.s0 f8=f13,f6,f0
//br.ret.sptk b0;;
//}


{.mfb
  nop.m 0
  // result=((T_high+l)+C1r)+((D+T_low)+r*P18)
  (p6) fma.s0 f8=f13,f1,f8
  // return
  br.ret.sptk b0;;
}


SPECIAL_log2l:
{.mfi
  nop.m 0
  mov FR_X=f8
  nop.i 0
}
{.mfi 
  nop.m 0
  // x=+Infinity ?
  fclass.m p7,p0=f8,0x21
  nop.i 0;;
}
{.mfi
  nop.m 0
  // x=+/-Zero ?
  fclass.m p8,p0=f7,0x7
  nop.i 0;;
}
{.mfi
  nop.m 0
  // x=-Infinity, -normal, -denormal ?
  fclass.m p6,p0=f8,0x3a
  nop.i 0;;
}
{.mfb
  nop.m 0
  // log2l(+Infinity)=+Infinity
  nop.f 0
  (p7) br.ret.spnt b0;;
}
{.mfi
  (p8) mov GR_Parameter_TAG = 168                          
  // log2l(+/-0)=-infinity, raises Divide by Zero
  // set f8=-0
  (p8) fmerge.ns f8=f0,f8
  nop.i 0;;
}
{.mfb
  nop.m 0
  (p8) frcpa.s0 f8,p0=f1,f8
  (p8) br.cond.sptk __libm_error_region;;
}
{.mfb
  (p6) mov GR_Parameter_TAG = 169 
  // x<0: return NaN, raise Invalid
  (p6) frcpa.s0 f8,p0=f0,f0
  (p6) br.cond.sptk __libm_error_region;;
}                          
  

{.mfb
  nop.m 0
  // Remaining cases: NaNs
  fma.s0 f8=f8,f1,f0
  br.ret.sptk b0;;
}

LOG2_PSEUDO_ZERO:

{.mfi
  nop.m 0
  mov FR_X=f8
  nop.i 0
}
{.mfi
  mov GR_Parameter_TAG = 168                          
  // log2l(+/-0)=-infinity, raises Divide by Zero
  // set f8=-0
  fmerge.ns f8=f0,f8
  nop.i 0;;
}
{.mfb
  nop.m 0
  frcpa.s0 f8,p0=f1,f8
  br.cond.sptk __libm_error_region;;
}


GLOBAL_IEEE754_END(log2l)


LOCAL_LIBM_ENTRY(__libm_error_region)
.prologue
{ .mfi
        add   GR_Parameter_Y=-32,sp             // Parameter 2 value
        nop.f 0
.save   ar.pfs,GR_SAVE_PFS
        mov  GR_SAVE_PFS=ar.pfs                 // Save ar.pfs 
}
{ .mfi
.fframe 64 
        add sp=-64,sp                           // Create new stack
        nop.f 0
        mov GR_SAVE_GP=gp                       // Save gp
};;
{ .mmi
        stfe [GR_Parameter_Y] = FR_Y,16         // STORE Parameter 2 on stack
        add GR_Parameter_X = 16,sp              // Parameter 1 address
.save   b0, GR_SAVE_B0                      
        mov GR_SAVE_B0=b0                       // Save b0 
};;
.body
{ .mib
        stfe [GR_Parameter_X] = FR_X                  // STORE Parameter 1 on stack 
        add   GR_Parameter_RESULT = 0,GR_Parameter_Y  // Parameter 3 address 
	nop.b 0                                      
}
{ .mib
        stfe [GR_Parameter_Y] = FR_RESULT             // STORE Parameter 3 on stack
        add   GR_Parameter_Y = -16,GR_Parameter_Y  
        br.call.sptk b0=__libm_error_support#         // Call error handling function
};;
{ .mmi
        nop.m 0
        nop.m 0
        add   GR_Parameter_RESULT = 48,sp
};;
{ .mmi
        ldfe  f8 = [GR_Parameter_RESULT]       // Get return result off stack
.restore sp
        add   sp = 64,sp                       // Restore stack pointer
        mov   b0 = GR_SAVE_B0                  // Restore return address
};;
{ .mib
        mov   gp = GR_SAVE_GP                  // Restore gp 
        mov   ar.pfs = GR_SAVE_PFS             // Restore ar.pfs
        br.ret.sptk     b0                     // Return
};; 

LOCAL_LIBM_END(__libm_error_region)
.type   __libm_error_support#,@function
.global __libm_error_support#