Home | History | Annotate | Download | only in math
      1 // Copyright 2010 The Go Authors. All rights reserved.
      2 // Use of this source code is governed by a BSD-style
      3 // license that can be found in the LICENSE file.
      4 
      5 #include "textflag.h"
      6 
      7 #define HSqrt2 7.07106781186547524401e-01 // sqrt(2)/2
      8 #define Ln2Hi  6.93147180369123816490e-01 // 0x3fe62e42fee00000
      9 #define Ln2Lo  1.90821492927058770002e-10 // 0x3dea39ef35793c76
     10 #define L1     6.666666666666735130e-01   // 0x3FE5555555555593
     11 #define L2     3.999999999940941908e-01   // 0x3FD999999997FA04
     12 #define L3     2.857142874366239149e-01   // 0x3FD2492494229359
     13 #define L4     2.222219843214978396e-01   // 0x3FCC71C51D8E78AF
     14 #define L5     1.818357216161805012e-01   // 0x3FC7466496CB03DE
     15 #define L6     1.531383769920937332e-01   // 0x3FC39A09D078C69F
     16 #define L7     1.479819860511658591e-01   // 0x3FC2F112DF3E5244
     17 #define NaN    0x7FF8000000000001
     18 #define NegInf 0xFFF0000000000000
     19 #define PosInf 0x7FF0000000000000
     20 
     21 // func Log(x float64) float64
     22 TEXT Log(SB),NOSPLIT,$0
     23 	// test bits for special cases
     24 	MOVQ    x+0(FP), BX
     25 	MOVQ    $~(1<<63), AX // sign bit mask
     26 	ANDQ    BX, AX
     27 	JEQ     isZero
     28 	MOVQ    $0, AX
     29 	CMPQ    AX, BX
     30 	JGT     isNegative
     31 	MOVQ    $PosInf, AX
     32 	CMPQ    AX, BX
     33 	JLE     isInfOrNaN
     34 	// f1, ki := math.Frexp(x); k := float64(ki)
     35 	MOVQ    BX, X0
     36 	MOVQ    $0x000FFFFFFFFFFFFF, AX
     37 	MOVQ    AX, X2
     38 	ANDPD   X0, X2
     39 	MOVSD   $0.5, X0 // 0x3FE0000000000000
     40 	ORPD    X0, X2 // X2= f1
     41 	SHRQ    $52, BX
     42 	ANDL    $0x7FF, BX
     43 	SUBL    $0x3FE, BX
     44 	XORPS   X1, X1 // break dependency for CVTSL2SD
     45 	CVTSL2SD BX, X1 // x1= k, x2= f1
     46 	// if f1 < math.Sqrt2/2 { k -= 1; f1 *= 2 }
     47 	MOVSD   $HSqrt2, X0 // x0= 0.7071, x1= k, x2= f1
     48 	CMPSD   X2, X0, 5 // cmpnlt; x0= 0 or ^0, x1= k, x2 = f1
     49 	MOVSD   $1.0, X3 // x0= 0 or ^0, x1= k, x2 = f1, x3= 1
     50 	ANDPD   X0, X3 // x0= 0 or ^0, x1= k, x2 = f1, x3= 0 or 1
     51 	SUBSD   X3, X1 // x0= 0 or ^0, x1= k, x2 = f1, x3= 0 or 1
     52 	MOVSD   $1.0, X0 // x0= 1, x1= k, x2= f1, x3= 0 or 1
     53 	ADDSD   X0, X3 // x0= 1, x1= k, x2= f1, x3= 1 or 2
     54 	MULSD   X3, X2 // x0= 1, x1= k, x2= f1
     55 	// f := f1 - 1
     56 	SUBSD   X0, X2 // x1= k, x2= f
     57 	// s := f / (2 + f)
     58 	MOVSD   $2.0, X0
     59 	ADDSD   X2, X0
     60 	MOVAPD  X2, X3
     61 	DIVSD   X0, X3 // x1=k, x2= f, x3= s
     62 	// s2 := s * s
     63 	MOVAPD  X3, X4 // x1= k, x2= f, x3= s
     64 	MULSD   X4, X4 // x1= k, x2= f, x3= s, x4= s2
     65 	// s4 := s2 * s2
     66 	MOVAPD  X4, X5 // x1= k, x2= f, x3= s, x4= s2
     67 	MULSD   X5, X5 // x1= k, x2= f, x3= s, x4= s2, x5= s4
     68 	// t1 := s2 * (L1 + s4*(L3+s4*(L5+s4*L7)))
     69 	MOVSD   $L7, X6
     70 	MULSD   X5, X6
     71 	ADDSD   $L5, X6
     72 	MULSD   X5, X6
     73 	ADDSD   $L3, X6
     74 	MULSD   X5, X6
     75 	ADDSD   $L1, X6
     76 	MULSD   X6, X4 // x1= k, x2= f, x3= s, x4= t1, x5= s4
     77 	// t2 := s4 * (L2 + s4*(L4+s4*L6))
     78 	MOVSD   $L6, X6
     79 	MULSD   X5, X6
     80 	ADDSD   $L4, X6
     81 	MULSD   X5, X6
     82 	ADDSD   $L2, X6
     83 	MULSD   X6, X5 // x1= k, x2= f, x3= s, x4= t1, x5= t2
     84 	// R := t1 + t2
     85 	ADDSD   X5, X4 // x1= k, x2= f, x3= s, x4= R
     86 	// hfsq := 0.5 * f * f
     87 	MOVSD   $0.5, X0
     88 	MULSD   X2, X0
     89 	MULSD   X2, X0 // x0= hfsq, x1= k, x2= f, x3= s, x4= R
     90 	// return k*Ln2Hi - ((hfsq - (s*(hfsq+R) + k*Ln2Lo)) - f)
     91 	ADDSD   X0, X4 // x0= hfsq, x1= k, x2= f, x3= s, x4= hfsq+R
     92 	MULSD   X4, X3 // x0= hfsq, x1= k, x2= f, x3= s*(hfsq+R)
     93 	MOVSD   $Ln2Lo, X4
     94 	MULSD   X1, X4 // x4= k*Ln2Lo
     95 	ADDSD   X4, X3 // x0= hfsq, x1= k, x2= f, x3= s*(hfsq+R)+k*Ln2Lo
     96 	SUBSD   X3, X0 // x0= hfsq-(s*(hfsq+R)+k*Ln2Lo), x1= k, x2= f
     97 	SUBSD   X2, X0 // x0= (hfsq-(s*(hfsq+R)+k*Ln2Lo))-f, x1= k
     98 	MULSD   $Ln2Hi, X1 // x0= (hfsq-(s*(hfsq+R)+k*Ln2Lo))-f, x1= k*Ln2Hi
     99 	SUBSD   X0, X1 // x1= k*Ln2Hi-((hfsq-(s*(hfsq+R)+k*Ln2Lo))-f)
    100   	MOVSD   X1, ret+8(FP)
    101 	RET
    102 isInfOrNaN:
    103 	MOVQ    BX, ret+8(FP) // +Inf or NaN, return x
    104 	RET
    105 isNegative:
    106 	MOVQ    $NaN, AX
    107 	MOVQ    AX, ret+8(FP) // return NaN
    108 	RET
    109 isZero:
    110 	MOVQ    $NegInf, AX
    111 	MOVQ    AX, ret+8(FP) // return -Inf
    112 	RET
    113