--- /dev/null
+// 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