1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
|
//go:build !tinygo && !noasm
// Copyright 2014 Xuanyi Chew. All rights reserved.
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.
//
// The original code is lifted from the Go standard library which is governed by
// a BSD-style licence which can be found here: https://golang.org/LICENSE
#include "textflag.h"
// The method is based on a paper by Naoki Shibata: "Efficient evaluation
// methods of elementary functions suitable for SIMD computation", Proc.
// of International Supercomputing Conference 2010 (ISC'10), pp. 25 -- 32
// (May 2010). The paper is available at
// http://www.springerlink.com/content/340228x165742104/
//
// The original code and the constants below are from the author's
// implementation available at http://freshmeat.net/projects/sleef.
// The README file says, "The software is in public domain.
// You can use the software without any obligation."
//
// This code is a simplified version of the original.
// The magic numbers for the float32 are lifted from the same project
#define LN2 0.693147182464599609375 // log_e(2)
#define LOG2E 1.44269502162933349609375 // 1/LN2
#define LN2U 0.693145751953125 // upper half LN2
#define LN2L 1.428606765330187045e-06 // lower half LN2
#define T0 1.0
#define T1 0.5
#define T2 0.166665524244308471679688
#define T3 0.0416710823774337768554688
#define T4 0.00836596917361021041870117
#define PosInf 0x7F800000
#define NegInf 0xFF800000
// func archExp(x float32) float32
TEXT ·archExp(SB),NOSPLIT,$0
// test bits for not-finite
MOVL x+0(FP), BX
MOVQ $~(1<<31), AX // sign bit mask
MOVL BX, DX
ANDL AX, DX
MOVL $PosInf, AX
CMPL AX, DX
JLE notFinite
MOVL BX, X0
MOVSS $LOG2E, X1
MULSS X0, X1
CVTSS2SL X1, BX // BX = exponent
CVTSL2SS BX, X1
MOVSS $LN2U, X2
MULSS X1, X2
SUBSS X2, X0
MOVSS $LN2L, X2
MULSS X1, X2
SUBSS X2, X0
// reduce argument
MULSS $0.0625, X0
// Taylor series evaluation
ADDSS $T4, X1
MULSS X0, X1
ADDSS $T3, X1
MULSS X0, X1
ADDSS $T2, X1
MULSS X0, X1
ADDSS $T1, X1
MULSS X0, X1
ADDSS $T0, X1
MULSS X1, X0
MOVSS $2.0, X1
ADDSS X0, X1
MULSS X1, X0
MOVSS $2.0, X1
ADDSS X0, X1
MULSS X1, X0
MOVSS $2.0, X1
ADDSS X0, X1
MULSS X1, X0
MOVSS $2.0, X1
ADDSS X0, X1
MULSS X1, X0
ADDSS $1.0, X0
// return fr * 2**exponent
MOVL $0x7F, AX // bias
ADDL AX, BX
JLE underflow
CMPL BX, $0xFF
JGE overflow
MOVL $23, CX
SHLQ CX, BX
MOVL BX, X1
MULSS X1, X0
MOVSS X0, ret+8(FP)
RET
notFinite:
// test bits for -Inf
MOVL $NegInf, AX
CMPQ AX, BX
JNE notNegInf
// -Inf, return 0
underflow: // return 0
MOVL $0, AX
MOVL AX, ret+8(FP)
RET
overflow: // return +Inf
MOVL $PosInf, BX
notNegInf: // NaN or +Inf, return x
MOVL BX, ret+8(FP)
RET
TEXT ·archExp2(SB),NOSPLIT,$0
JMP ·exp2(SB)
|
