// Copyright 2025 The Go Authors. All rights reserved.
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.

#include "textflag.h"

#define NearZero	0x3e30000000000000	// 2**-28
#define PosInf		0x7ff0000000000000
#define FracMask	0x000fffffffffffff
#define C1		0x3cb0000000000000	// 2**-52

DATA exprodata<>+0(SB)/8, $0.0
DATA exprodata<>+8(SB)/8, $0.5
DATA exprodata<>+16(SB)/8, $1.0
DATA exprodata<>+24(SB)/8, $2.0
DATA exprodata<>+32(SB)/8, $6.93147180369123816490e-01	// Ln2Hi
DATA exprodata<>+40(SB)/8, $1.90821492927058770002e-10	// Ln2Lo
DATA exprodata<>+48(SB)/8, $1.44269504088896338700e+00	// Log2e
DATA exprodata<>+56(SB)/8, $7.09782712893383973096e+02	// Overflow
DATA exprodata<>+64(SB)/8, $-7.45133219101941108420e+02	// Underflow
DATA exprodata<>+72(SB)/8, $1.0239999999999999e+03	// Overflow2
DATA exprodata<>+80(SB)/8, $-1.0740e+03			// Underflow2
DATA exprodata<>+88(SB)/8, $3.7252902984619141e-09	// NearZero
GLOBL exprodata<>+0(SB), NOPTR|RODATA, $96

DATA expmultirodata<>+0(SB)/8, $1.66666666666666657415e-01	// P1
DATA expmultirodata<>+8(SB)/8, $-2.77777777770155933842e-03	// P2
DATA expmultirodata<>+16(SB)/8, $6.61375632143793436117e-05	// P3
DATA expmultirodata<>+24(SB)/8, $-1.65339022054652515390e-06	// P4
DATA expmultirodata<>+32(SB)/8, $4.13813679705723846039e-08	// P5
GLOBL expmultirodata<>+0(SB), NOPTR|RODATA, $40

// Exp returns e**x, the base-e exponential of x.
// This is an assembly implementation of the method used for function Exp in file exp.go.
//
// func Exp(x float64) float64
TEXT ·archExp(SB),$0-16
	MOVD	x+0(FP), F0	// F0 = x

	MOVV	$exprodata<>+0(SB), R10
	MOVD	56(R10), F1	// Overflow
	MOVD	64(R10), F2	// Underflow
	MOVD	88(R10), F3	// NearZero
	MOVD	16(R10), F17	// 1.0

	CMPEQD	F0, F0, FCC0
	BFPF	isNaN		// x = NaN, return NaN

	CMPGTD	F0, F1, FCC0
	BFPT	overflow	// x > Overflow, return PosInf

	CMPGTD	F2, F0, FCC0
	BFPT	underflow	// x < Underflow, return 0

	ABSD	F0, F5
	CMPGTD	F3, F5, FCC0
	BFPT	nearzero	// fabs(x) < NearZero, return 1 + x

	// argument reduction, x = k*ln2 + r,  |r| <= 0.5*ln2
	// computed as r = hi - lo for extra precision.
	MOVD	0(R10), F5
	MOVD	8(R10), F3
	MOVD	48(R10), F2
	CMPGTD	F0, F5, FCC0
	BFPT	add		// x > 0
sub:
	FMSUBD	F3, F2, F0, F3	// Log2e*x - 0.5
	JMP	2(PC)
add:
	FMADDD	F3, F2, F0, F3	// Log2e*x + 0.5

	FTINTRZVD F3, F4	// float64 -> int64
	MOVV	F4, R5		// R5 = int(k)
	FFINTDV	F4, F3		// int64 -> float64

	MOVD	32(R10), F4
	MOVD	40(R10), F5
	FNMSUBD	F0, F3, F4, F4
	MULD	F3, F5, F5
	SUBD	F5, F4, F6
	MULD	F6, F6, F7

	// compute c
	MOVV	$expmultirodata<>+0(SB), R11
	MOVD	32(R11), F8
	MOVD	24(R11), F9
	FMADDD	F9, F8, F7, F13
	MOVD	16(R11), F10
	FMADDD	F10, F13, F7, F13
	MOVD	8(R11), F11
	FMADDD	F11, F13, F7, F13
	MOVD	0(R11), F12
	FMADDD	F12, F13, F7, F13
	FNMSUBD	F6, F13, F7, F13

	// compute y
	MOVD	24(R10), F14
	SUBD	F13, F14, F14
	MULD	F6, F13, F15
	DIVD	F14, F15, F15
	SUBD	F15, F5, F15
	SUBD	F4, F15, F15
	SUBD	F15, F17, F16

	// inline Ldexp(y, k), benefit:
	// 1, no parameter pass overhead.
	// 2, skip unnecessary checks for Inf/NaN/Zero
	MOVV	F16, R4
	MOVV	$FracMask, R9
	AND	R9, R4, R6	// fraction
	SRLV	$52, R4, R7	// exponent
	ADDV	R5, R7
	MOVV	$1, R12
	BGE	R7, R12, normal
	ADDV	$52, R7		// denormal
	MOVV	$C1, R8
	MOVV	R8, F17
normal:
	SLLV	$52, R7
	OR	R7, R6, R4
	MOVV	R4, F0
	MULD	F17, F0		// return m * x
	MOVD	F0, ret+8(FP)
	RET
nearzero:
	ADDD	F17, F0, F0
isNaN:
	MOVD	F0, ret+8(FP)
	RET
underflow:
	MOVV	R0, ret+8(FP)
	RET
overflow:
	MOVV	$PosInf, R4
	MOVV	R4, ret+8(FP)
	RET


// Exp2 returns 2**x, the base-2 exponential of x.
// This is an assembly implementation of the method used for function Exp2 in file exp.go.
//
// func Exp2(x float64) float64
TEXT ·archExp2(SB),$0-16
	MOVD	x+0(FP), F0	// F0 = x

	MOVV	$exprodata<>+0(SB), R10
	MOVD	72(R10), F1	// Overflow2
	MOVD	80(R10), F2	// Underflow2
	MOVD	88(R10), F3	// NearZero

	CMPEQD	F0, F0, FCC0
	BFPF	isNaN		// x = NaN, return NaN

	CMPGTD	F0, F1, FCC0
	BFPT	overflow	// x > Overflow, return PosInf

	CMPGTD	F2, F0, FCC0
	BFPT	underflow	// x < Underflow, return 0

	// argument reduction; x = r*lg(e) + k with |r| <= ln(2)/2
	// computed as r = hi - lo for extra precision.
	MOVD	0(R10), F10
	MOVD	8(R10), F2
	CMPGTD	F0, F10, FCC0
	BFPT	add
sub:
	SUBD	F2, F0, F3	// x - 0.5
	JMP	2(PC)
add:
	ADDD	F2, F0, F3	// x + 0.5

	FTINTRZVD F3, F4
	MOVV	F4, R5
	FFINTDV	F4, F3

	MOVD	32(R10), F4
	MOVD	40(R10), F5
	SUBD	F3, F0, F3
	MULD	F3, F4
	FNMSUBD	F10, F3, F5, F5
	SUBD	F5, F4, F6
	MULD	F6, F6, F7

	// compute c
	MOVV	$expmultirodata<>+0(SB), R11
	MOVD	32(R11), F8
	MOVD	24(R11), F9
	FMADDD	F9, F8, F7, F13
	MOVD	16(R11), F10
	FMADDD	F10, F13, F7, F13
	MOVD	8(R11), F11
	FMADDD	F11, F13, F7, F13
	MOVD	0(R11), F12
	FMADDD	F12, F13, F7, F13
	FNMSUBD	F6, F13, F7, F13

	// compute y
	MOVD	24(R10), F14
	SUBD	F13, F14, F14
	MULD	F6, F13, F15
	DIVD	F14, F15

	MOVD	16(R10), F17
	SUBD	F15, F5, F15
	SUBD	F4, F15, F15
	SUBD	F15, F17, F16

	// inline Ldexp(y, k), benefit:
	// 1, no parameter pass overhead.
	// 2, skip unnecessary checks for Inf/NaN/Zero
	MOVV	F16, R4
	MOVV	$FracMask, R9
	SRLV	$52, R4, R7	// exponent
	AND	R9, R4, R6	// fraction
	ADDV	R5, R7
	MOVV	$1, R12
	BGE	R7, R12, normal

	ADDV	$52, R7		// denormal
	MOVV	$C1, R8
	MOVV	R8, F17
normal:
	SLLV	$52, R7
	OR	R7, R6, R4
	MOVV	R4, F0
	MULD	F17, F0
isNaN:
	MOVD	F0, ret+8(FP)
	RET
underflow:
	MOVV	R0, ret+8(FP)
	RET
overflow:
	MOVV	$PosInf, R4
	MOVV	R4, ret+8(FP)
	RET