1// Copyright 2017 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// Constants
8DATA ·log1pxlim<> + 0(SB)/4, $0xfff00000
9GLOBL ·log1pxlim<> + 0(SB), RODATA, $4
10DATA ·log1pxzero<> + 0(SB)/8, $0.0
11GLOBL ·log1pxzero<> + 0(SB), RODATA, $8
12DATA ·log1pxminf<> + 0(SB)/8, $0xfff0000000000000
13GLOBL ·log1pxminf<> + 0(SB), RODATA, $8
14DATA ·log1pxnan<> + 0(SB)/8, $0x7ff8000000000000
15GLOBL ·log1pxnan<> + 0(SB), RODATA, $8
16DATA ·log1pyout<> + 0(SB)/8, $0x40fce621e71da000
17GLOBL ·log1pyout<> + 0(SB), RODATA, $8
18DATA ·log1pxout<> + 0(SB)/8, $0x40f1000000000000
19GLOBL ·log1pxout<> + 0(SB), RODATA, $8
20DATA ·log1pxl2<> + 0(SB)/8, $0xbfda7aecbeba4e46
21GLOBL ·log1pxl2<> + 0(SB), RODATA, $8
22DATA ·log1pxl1<> + 0(SB)/8, $0x3ffacde700000000
23GLOBL ·log1pxl1<> + 0(SB), RODATA, $8
24DATA ·log1pxa<> + 0(SB)/8, $5.5
25GLOBL ·log1pxa<> + 0(SB), RODATA, $8
26DATA ·log1pxmone<> + 0(SB)/8, $-1.0
27GLOBL ·log1pxmone<> + 0(SB), RODATA, $8
28
29// Minimax polynomial approximations
30DATA ·log1pc8<> + 0(SB)/8, $0.212881813645679599E-07
31GLOBL ·log1pc8<> + 0(SB), RODATA, $8
32DATA ·log1pc7<> + 0(SB)/8, $-.148682720127920854E-06
33GLOBL ·log1pc7<> + 0(SB), RODATA, $8
34DATA ·log1pc6<> + 0(SB)/8, $0.938370938292558173E-06
35GLOBL ·log1pc6<> + 0(SB), RODATA, $8
36DATA ·log1pc5<> + 0(SB)/8, $-.602107458843052029E-05
37GLOBL ·log1pc5<> + 0(SB), RODATA, $8
38DATA ·log1pc4<> + 0(SB)/8, $0.397389654305194527E-04
39GLOBL ·log1pc4<> + 0(SB), RODATA, $8
40DATA ·log1pc3<> + 0(SB)/8, $-.273205381970859341E-03
41GLOBL ·log1pc3<> + 0(SB), RODATA, $8
42DATA ·log1pc2<> + 0(SB)/8, $0.200350613573012186E-02
43GLOBL ·log1pc2<> + 0(SB), RODATA, $8
44DATA ·log1pc1<> + 0(SB)/8, $-.165289256198351540E-01
45GLOBL ·log1pc1<> + 0(SB), RODATA, $8
46DATA ·log1pc0<> + 0(SB)/8, $0.181818181818181826E+00
47GLOBL ·log1pc0<> + 0(SB), RODATA, $8
48
49
50// Table of log10 correction terms
51DATA ·log1ptab<> + 0(SB)/8, $0.585235384085551248E-01
52DATA ·log1ptab<> + 8(SB)/8, $0.412206153771168640E-01
53DATA ·log1ptab<> + 16(SB)/8, $0.273839003221648339E-01
54DATA ·log1ptab<> + 24(SB)/8, $0.166383778368856480E-01
55DATA ·log1ptab<> + 32(SB)/8, $0.866678223433169637E-02
56DATA ·log1ptab<> + 40(SB)/8, $0.319831684989627514E-02
57DATA ·log1ptab<> + 48(SB)/8, $-.000000000000000000E+00
58DATA ·log1ptab<> + 56(SB)/8, $-.113006378583725549E-02
59DATA ·log1ptab<> + 64(SB)/8, $-.367979419636602491E-03
60DATA ·log1ptab<> + 72(SB)/8, $0.213172484510484979E-02
61DATA ·log1ptab<> + 80(SB)/8, $0.623271047682013536E-02
62DATA ·log1ptab<> + 88(SB)/8, $0.118140812789696885E-01
63DATA ·log1ptab<> + 96(SB)/8, $0.187681358930914206E-01
64DATA ·log1ptab<> + 104(SB)/8, $0.269985148668178992E-01
65DATA ·log1ptab<> + 112(SB)/8, $0.364186619761331328E-01
66DATA ·log1ptab<> + 120(SB)/8, $0.469505379381388441E-01
67GLOBL ·log1ptab<> + 0(SB), RODATA, $128
68
69// Log1p returns the natural logarithm of 1 plus its argument x.
70// It is more accurate than Log(1 + x) when x is near zero.
71//
72// Special cases are:
73//      Log1p(+Inf) = +Inf
74//      Log1p(±0) = ±0
75//      Log1p(-1) = -Inf
76//      Log1p(x < -1) = NaN
77//      Log1p(NaN) = NaN
78// The algorithm used is minimax polynomial approximation
79// with coefficients determined with a Remez exchange algorithm.
80
81TEXT	·log1pAsm(SB), NOSPLIT, $0-16
82	FMOVD	x+0(FP), F0
83	MOVD	$·log1pxmone<>+0(SB), R1
84	MOVD	·log1pxout<>+0(SB), R2
85	FMOVD	0(R1), F3
86	MOVD	$·log1pxa<>+0(SB), R1
87	MOVWZ	·log1pxlim<>+0(SB), R0
88	FMOVD	0(R1), F1
89	MOVD	$·log1pc8<>+0(SB), R1
90	FMOVD	0(R1), F5
91	MOVD	$·log1pc7<>+0(SB), R1
92	VLEG	$0, 0(R1), V20
93	MOVD	$·log1pc6<>+0(SB), R1
94	WFSDB	V0, V3, V4
95	VLEG	$0, 0(R1), V18
96	MOVD	$·log1pc5<>+0(SB), R1
97	VLEG	$0, 0(R1), V16
98	MOVD	R2, R5
99	LGDR	F4, R3
100	WORD	$0xC0190006	//iilf	%r1,425983
101	BYTE	$0x7F
102	BYTE	$0xFF
103	SRAD	$32, R3, R3
104	SUBW	R3, R1
105	SRW	$16, R1, R1
106	BYTE	$0x18	//lr	%r4,%r1
107	BYTE	$0x41
108	RISBGN	$0, $15, $48, R4, R2
109	RISBGN	$16, $31, $32, R4, R5
110	MOVW	R0, R6
111	MOVW	R3, R7
112	CMPBGT	R6, R7, L8
113	WFCEDBS	V4, V4, V6
114	MOVD	$·log1pxzero<>+0(SB), R1
115	FMOVD	0(R1), F2
116	BVS	LEXITTAGlog1p
117	WORD	$0xB3130044	// lcdbr %f4,%f4
118	WFCEDBS	V2, V4, V6
119	BEQ	L9
120	WFCHDBS	V4, V2, V2
121	BEQ	LEXITTAGlog1p
122	MOVD	$·log1pxnan<>+0(SB), R1
123	FMOVD	0(R1), F0
124	FMOVD	F0, ret+8(FP)
125	RET
126
127L8:
128	LDGR	R2, F2
129	FSUB	F4, F3
130	FMADD	F2, F4, F1
131	MOVD	$·log1pc4<>+0(SB), R2
132	WORD	$0xB3130041	// lcdbr %f4,%f1
133	FMOVD	0(R2), F7
134	FSUB	F3, F0
135	MOVD	$·log1pc3<>+0(SB), R2
136	FMOVD	0(R2), F3
137	MOVD	$·log1pc2<>+0(SB), R2
138	WFMDB	V1, V1, V6
139	FMADD	F7, F4, F3
140	WFMSDB	V0, V2, V1, V0
141	FMOVD	0(R2), F7
142	WFMADB	V4, V5, V20, V5
143	MOVD	$·log1pc1<>+0(SB), R2
144	FMOVD	0(R2), F2
145	FMADD	F7, F4, F2
146	WFMADB	V4, V18, V16, V4
147	FMADD	F3, F6, F2
148	WFMADB	V5, V6, V4, V5
149	FMUL	F6, F6
150	MOVD	$·log1pc0<>+0(SB), R2
151	WFMADB	V6, V5, V2, V6
152	FMOVD	0(R2), F4
153	WFMADB	V0, V6, V4, V6
154	RISBGZ	$57, $60, $3, R1, R1
155	MOVD	$·log1ptab<>+0(SB), R2
156	MOVD	$·log1pxl1<>+0(SB), R3
157	WORD	$0x68112000	//ld	%f1,0(%r1,%r2)
158	FMOVD	0(R3), F2
159	WFMADB	V0, V6, V1, V0
160	MOVD	$·log1pyout<>+0(SB), R1
161	LDGR	R5, F6
162	FMOVD	0(R1), F4
163	WFMSDB	V2, V6, V4, V2
164	MOVD	$·log1pxl2<>+0(SB), R1
165	FMOVD	0(R1), F4
166	FMADD	F4, F2, F0
167	FMOVD	F0, ret+8(FP)
168	RET
169
170L9:
171	MOVD	$·log1pxminf<>+0(SB), R1
172	FMOVD	0(R1), F0
173	FMOVD	F0, ret+8(FP)
174	RET
175
176
177LEXITTAGlog1p:
178	FMOVD	F0, ret+8(FP)
179	RET
180
181