diff --git a/doc/doxygen/Doxyfile.in b/doc/doxygen/Doxyfile.in
index d454898f4e..e905d5a64d 100644
--- a/doc/doxygen/Doxyfile.in
+++ b/doc/doxygen/Doxyfile.in
@@ -2,7 +2,7 @@
 
 DOXYFILE_ENCODING      = UTF-8
 PROJECT_NAME           = "LAMMPS Programmer's Guide"
-PROJECT_NUMBER         = "24 August 2020"
+PROJECT_NUMBER         = "4 May 2022"
 PROJECT_BRIEF          = "Documentation of the LAMMPS library interface and Python wrapper"
 PROJECT_LOGO           = lammps-logo.png
 CREATE_SUBDIRS         = NO
@@ -437,6 +437,8 @@ INPUT                  = @LAMMPS_SOURCE_DIR@/utils.cpp                 \
                          @LAMMPS_SOURCE_DIR@/math_eigen.h              \
                          @LAMMPS_SOURCE_DIR@/platform.h                \
                          @LAMMPS_SOURCE_DIR@/platform.cpp              \
+                         @LAMMPS_SOURCE_DIR@/math_special.h            \
+                         @LAMMPS_SOURCE_DIR@/math_special.cpp          \
 
 # The EXCLUDE_SYMLINKS tag can be used to select whether or not files or
 # directories that are symbolic links (a Unix file system feature) are excluded
diff --git a/doc/src/Developer_utils.rst b/doc/src/Developer_utils.rst
index d7b1f077dd..0883d6b3c8 100644
--- a/doc/src/Developer_utils.rst
+++ b/doc/src/Developer_utils.rst
@@ -246,6 +246,23 @@ Customized standard functions
 
 ---------------------------
 
+Special Math functions
+----------------------
+
+The ``MathSpecial`` namespace implements a selection of custom and optimized
+mathematical functions for a variety of applications.
+
+.. doxygenfunction:: factorial
+   :project: progguide
+
+.. doxygenfunction:: exp2_x86
+   :project: progguide
+
+.. doxygenfunction:: fm_exp
+   :project: progguide
+
+---------------------------
+
 Tokenizer classes
 -----------------
 
diff --git a/src/math_special.cpp b/src/math_special.cpp
index da99c5850a..fd47aec9e9 100644
--- a/src/math_special.cpp
+++ b/src/math_special.cpp
@@ -706,6 +706,7 @@ static const double fm_exp2_p[] = {
 
 double MathSpecial::exp2_x86(double x)
 {
+#if defined(__BYTE_ORDER__) && (__BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__)
     double   ipart, fpart, px, qx;
     udi_t    epart;
 
@@ -726,6 +727,9 @@ double MathSpecial::exp2_x86(double x)
 
     x = 1.0 + 2.0*(px/(qx-px));
     return epart.f*x;
+#else
+    return pow(2.0, x);
+#endif
 }
 
 double MathSpecial::fm_exp(double x)
diff --git a/src/math_special.h b/src/math_special.h
index fa39d46321..db780a240c 100644
--- a/src/math_special.h
+++ b/src/math_special.h
@@ -20,20 +20,50 @@ namespace LAMMPS_NS {
 
 namespace MathSpecial {
 
-  // tabulated factorial function
+  /*! Fast tabulated factorial function
+   *
+   *  This function looks up precomputed factorial values for arguments of n = 0
+   *  to a maximum of 167, which is the maximal value representable by a double
+   *  precision floating point number.  For other values of n a NaN value is returned.
+   *
+   *  \param   n  argument (valid: 0 <= n <= 167)
+   *  \return  value of n! as double precision number or NaN */
 
-  extern double factorial(const int);
+  extern double factorial(const int n);
+
+  /*! Fast implementation of 2^x without argument checks for little endian CPUs
+   *
+   *  This function implements an optimized version of pow(2.0, x) that does not
+   *  check for valid arguments and thus may only be used where arguments are well
+   *  behaved.  The implementation makes assumptions about the layout of double
+   *  precision floating point numbers in memory and thus will only work on little
+   *  endian CPUs.  If little endian cannot be safely detected, the result of
+   *  calling pow(2.0, x) will be returned.  This function also is the basis for
+   *  the fast exponential fm_exp(x).
+   *
+   *  \param   x argument
+   *  \return  value of 2^x as double precision number */
+
+  extern double exp2_x86(double x);
+
+  /*! Fast implementation of exp(x) for little endian CPUs
+   *
+   *  This function implements an optimized version of exp(x) for little endian CPUs.
+   *  It calls the exp2_x86(x) function with a suitable prefactor to x to return exp(x).
+   *  The implementation makes assumptions about the layout of double
+   *  precision floating point numbers in memory and thus will only work on little
+   *  endian CPUs.  If little endian cannot be safely detected, the result of
+   *  calling the exp(x) implementation in the libm math library will be returned.
+   *
+   *  \param   x argument
+   *  \return  value of e^x as double precision number */
+
+  extern double fm_exp(double x);
 
   // support function for scaled error function complement
 
   extern double erfcx_y100(const double y100);
 
-  // fast 2**x function without argument checks for little endian CPUs
-  extern double exp2_x86(double x);
-
-  // fast e**x function for little endian CPUs, falls back to libc on other platforms
-  extern double fm_exp(double x);
-
   // scaled error function complement exp(x*x)*erfc(x) for coul/long styles
 
   static inline double my_erfcx(const double x)