diff --git a/reference/cmath/assoc_legendre.md b/reference/cmath/assoc_legendre.md
index dc96d1b09c..41ec51eb76 100644
--- a/reference/cmath/assoc_legendre.md
+++ b/reference/cmath/assoc_legendre.md
@@ -59,6 +59,30 @@ $$
 - `l >= 128` の場合、この関数の呼び出しの効果は実装定義である
 - (1) : C++23では、拡張浮動小数点数型を含む浮動小数点数型へのオーバーロードとして定義された
 
+### 負の m の対応
+この標準関数は $m$ が正の場合にしか対応していない。
+一方でルジャンドル陪関数はロドリゲスの公式を用いて負の $m$ に対して自然に拡張され、
+このことは球面調和関数を定義する上でも使われる。
+負の $m$ に対してもルジャンドル陪関数を計算する必要がある場合は、関係式
+$$ P_l^{-m}(x) = (-1)^m \frac{(l-m)!}{(l+m)!} P_l^m(x) $$
+を用いる必要がある。
+
+```cpp
+#include <cmath>
+
+// 負の m にも対応した実装例
+double assoc_legendre(unsigned l, int m, double x) {
+  if (m >= 0)
+    return std::assoc_legendre(l, (unsigned) m, x);
+  else
+    return std::pow(-1.0, m) * (std::tgamma(1.0 + l + m) / std::tgamma(1.0 + l - m)) * std::assoc_legendre(l, (unsigned) -m, x);
+}
+```
+* std::assoc_legendre[color ff0000]
+* std::tgamma[link tgamma.md]
+
+上記の例では簡単のために階乗をガンマ関数 $n! = \Gamma(n + 1)$ ([`tgamma`](tgamma.md)) で計算しているが、
+計算効率やオーバーフローなどを考えると、直接 $(l + |m|)\cdots(l - |m| + 1)$ で割り算したり、係数を事前計算しておくなど工夫すると良い。
 
 ## 例
 ```cpp example
diff --git a/reference/cmath/legendre.md b/reference/cmath/legendre.md
index 5e35b3d88c..9c9e54241c 100644
--- a/reference/cmath/legendre.md
+++ b/reference/cmath/legendre.md
@@ -103,6 +103,9 @@ legendre(3, 1) = 1
 - [ICC](/implementation.md#icc): ??
 - [Visual C++](/implementation.md#visual_cpp): ??
 
+### 備考
+- GCC 7, 8 には $x$ として [-1.0, 1.0] の範囲外の値を渡した時に、[`std::domain_error`](/reference/stdexcept.md) を投げるバグがあった (GCC 7.5, 8.2, 8.4 [mark verified])。
+  GCC 9 以降では直っている (GCC 9.3, 10.5, 11.4, 12.3, 13.2, 15.0 [mark verified])。
 
 ## 実装例
 ### 閉形式
diff --git a/reference/cmath/lgamma.md b/reference/cmath/lgamma.md
index bcb8dd5521..03027e27d0 100644
--- a/reference/cmath/lgamma.md
+++ b/reference/cmath/lgamma.md
@@ -35,7 +35,7 @@ namespace std {
 * Integral[italic]
 
 ## 概要
-ガンマ関数の絶対値の自然対数を求める。
+ガンマ関数 $\Gamma(x)$ ([`tgamma`](tgamma.md)) の絶対値の自然対数を求める。
 
 - (1) : `float`に対するオーバーロード
 - (2) : `double`に対するオーバーロード
@@ -60,7 +60,8 @@ namespace std {
     - `x = -∞` の場合、戻り値は `+∞` となる。
     - `x = +∞` の場合、戻り値は `+∞` となる。
 - C++23では、(1)、(2)、(3)が(4)に統合され、拡張浮動小数点数型を含む浮動小数点数型へのオーバーロードとして定義された
-
+- この関数はガンマ関数 ([`tgamma`](tgamma.md)) がオーバーフローするような場合に使う。
+  具体例については[ガンマ関数の備考](tgamma.md#remarks-lgamma)を参照のこと。
 
 ## 例
 ```cpp example
@@ -107,6 +108,8 @@ lgamma(+∞)  = inf
 
 - GCC 4.6.1 以上
 
+## 関連項目
+- ガンマ関数 [`tgamma`](tgamma.md)
 
 ## 参照
 - [P1467R9 Extended floating-point types and standard names](https://www.open-std.org/jtc1/sc22/wg21/docs/papers/2022/p1467r9.html)
diff --git a/reference/cmath/sph_legendre.md b/reference/cmath/sph_legendre.md
index 1f32cba3dd..daa2b669b5 100644
--- a/reference/cmath/sph_legendre.md
+++ b/reference/cmath/sph_legendre.md
@@ -47,7 +47,7 @@ namespace std {
 
 
 ## 戻り値
-引数 `l`, `m`, `theta` について $Y_l^m(\theta, 0)$ を返す。
+引数 `l`, `m`, `theta` について $Y_l^m(\theta, 0)$ (ただし $0 \le m \le l$) を返す。
 $Y_l^m(\theta, \phi)$ は球面調和関数
 $$
 Y_l^m(\theta, \phi) = (-1)^m \sqrt{\frac{2l + 1}{4\pi} \frac{(l - m)!}{(l + m)!}} P_l^m(\cos \theta) \exp(i m \phi)
@@ -56,23 +56,85 @@ $$
 である。
 $P_l^m$ はルジャンドル陪関数 ([`assoc_legendre`](assoc_legendre.md)) である。
 
-
 ## 備考
 - `l >= 128` の場合、この関数の呼び出しの効果は実装定義である
 - (1) : C++23では、拡張浮動小数点数型を含む浮動小数点数型へのオーバーロードとして定義された
 
+### 球面調和関数
+球面調和関数として使う場合には $\phi$ 依存性は自分で与える必要がある。
+また、この標準関数は $m\ge0$ にしか対応していないので、$m < 0$ の時の球面調和関数を計算したければ自分で $Y_l^{|m|}(\theta,0)$ の値を調節して使う必要がある。
+ルジャンドル陪関数の性質 $P_l^{-m}(x) = (-1)^m [(l - m)!/(l + m)!] P_l^m(x)$ より、一般の $m$ の球面調和関数は
+$$
+\begin{align*}
+Y_l^m(\theta, \phi)
+  &= (-1)^{(m+|m|)/2} \sqrt{\frac{2l + 1}{4\pi} \frac{(l - |m|)!}{(l + |m|)!}} P_l^{|m|}(\cos \theta) e^{i m \phi} \\
+  &= \begin{cases}
+    Y_l^{m}(\theta, 0) e^{im\phi}, & (0 \le m \le l), \\
+    (-1)^{|m|} Y_l^{|m|}(\theta, 0) e^{im\phi}, & (-l \le m < 0).
+  \end{cases}
+\end{align*}
+$$
+で与えられる。
+
+```cpp
+#include <cmath>
+#include <complex>
+#include <numbers>
+
+// 球面調和関数の実装例
+std::complex<double> sph_harmonics(unsigned l, int m, double theta, double phi) {
+  if (m >= 0)
+    return std::sph_legendre(l, (unsigned) m, theta) * std::polar(1.0, m * phi);
+  else
+    return std::sph_legendre(l, (unsigned) -m, theta) * std::polar(1.0, m * (phi - std::numbers::pi));
+}
+```
+* std::sph_legendre[color ff0000]
+* std::polar[link /reference/complex/complex/polar.md]
+* std::numbers::pi[link /reference/numbers/pi.md]
+
+また線形結合を取り直して実数にした、実数球面調和関数 $Y_{lm}(\theta,\phi)$ を計算することもできる。
+
+$$
+Y_{lm}(\theta,\phi)
+= \begin{cases}
+  \frac{(-1)^m Y_l^{|m|}(\theta,\phi) - Y_l^{-|m|}(\theta,\phi)}{\sqrt{2} i} = \sqrt{2} (-1)^{|m|} Y_l^{|m|}(\theta, 0) \sin(|m|\phi), & (-l \le m < 0), \\
+  Y_l^0 = Y_l^0(\theta, 0), & (m = 0), \\
+  \frac{(-1)^m Y_l^{|m|}(\theta,\phi) + Y_l^{-|m|}(\theta,\phi)}{\sqrt{2}} = \sqrt{2} (-1)^{|m|} Y_l^{|m|}(\theta, 0) \cos(|m|\phi), & (0 < m \le l).
+\end{cases}
+$$
+
+```cpp
+#include <cmath>
+#include <numbers>
+
+// 実数球面調和関数の実装例
+double real_sph_harmonics(unsigned l, int m, double theta, double phi) {
+  if (m == 0)
+    return std::sph_legendre(l, 0u, theta);
+  else if (m > 0)
+    return std::numbers::sqrt2 * std::sph_legendre(l, (unsigned) m, theta) * std::cos(m * (phi - std::numbers::pi));
+  else
+    return std::numbers::sqrt2 * std::sph_legendre(l, (unsigned) -m, theta) * std::sin(-m * (phi - std::numbers::pi));
+}
+```
+* std::sph_legendre[color ff0000]
+* std::numbers::pi[link /reference/numbers/pi.md]
+* std::numbers::sqrt2[link /reference/numbers/sqrt2.md]
 
 ## 例
 ```cpp example
 #include <cmath>
 #include <complex>
+#include <numbers>
 #include <iostream>
 
-constexpr double pi = 3.141592653589793;
-
 // 球面調和関数
-std::complex<double> sph_harmonics(unsigned l, unsigned m, double theta, double phi) {
-  return std::sph_legendre(l, m, theta) * std::polar(1.0, m * phi);
+std::complex<double> sph_harmonics(unsigned l, int m, double theta, double phi) {
+  if (m >= 0)
+    return std::sph_legendre(l, (unsigned) m, theta) * std::polar(1.0, m * phi);
+  else
+    return std::sph_legendre(l, (unsigned) -m, theta) * std::polar(1.0, m * (phi - std::numbers::pi));
 }
 
 int main() {
@@ -82,9 +144,9 @@ int main() {
 
   std::cout << "#θ / π\tφ / π\tY_" << l << "^" << m << "(θ, φ)\n";
   for (double t : {0., 0.25, 0.5, 0.75, 1.}) {
-    double theta = t * pi;
+    double theta = t * std::numbers::pi;
     for (double p : {0., 0.25, 0.5, 0.75, 1., 1.25, 1.5, 1.75, 2.}) {
-      double phi = p * pi / 4;
+      double phi = p * std::numbers::pi / 4;
       std::cout << t << "\t" << p << "\t" << sph_harmonics(l, m, theta, phi) << "\n";
       if (t == 0 || t == 1) break;
     }
@@ -93,6 +155,7 @@ int main() {
 ```
 * std::sph_legendre[color ff0000]
 * std::polar[link /reference/complex/complex/polar.md]
+* std::numbers::pi[link /reference/numbers/pi.md]
 
 ### 出力例
 ```
diff --git a/reference/cmath/tgamma.md b/reference/cmath/tgamma.md
index 4f3e397eb7..615ff77179 100644
--- a/reference/cmath/tgamma.md
+++ b/reference/cmath/tgamma.md
@@ -47,12 +47,12 @@ namespace std {
 
 
 ## 戻り値
-引数 `x` のガンマ関数を返す。
+引数 `x` のガンマ関数 $\Gamma(x)$ を返す。
 
+$$ \Gamma (x) = \int_0^\infty t^{x-1} e^{-t} dt $$
 
 ## 備考
-- $$ f(x) = \Gamma (x) $$
-- ガンマ関数は階乗の一般化である。
+- ガンマ関数は階乗の一般化である。階乗はガンマ関数を用いて $n! = \Gamma(n + 1)$ で計算できる。
 - C++11 以降では、処理系が IEC 60559 に準拠している場合（[`std::numeric_limits`](../limits/numeric_limits.md)`<T>::`[`is_iec559`](../limits/numeric_limits/is_iec559.md)`() != false`）、以下の規定が追加される。
     - `x = ±0` の場合、戻り値は `±∞` となり、
     [`FE_DIVBYZERO`](../cfenv/fe_divbyzero.md)（ゼロ除算浮動小数点例外）が発生する。
@@ -64,6 +64,32 @@ namespace std {
 - `gamma` という関数は既にあったが処理系によって定義が違ったため、本当の (true) ガンマ関数 `tgamma` と名付けられた。
 - C++23では、(1)、(2)、(3)が(4)に統合され、拡張浮動小数点数型を含む浮動小数点数型へのオーバーロードとして定義された
 
+### <a id="remarks-lgamma" href="#remarks-lgamma">lgamma との使い分け</a>
+ガンマ関数は急激に増加し容易にオーバーフローするので、代わりにガンマ関数の結果を自然対数で返す関数 [`lgamma`](lgamma.md) を用いた方が良いことが多くある。
+例えばガンマ関数の比を計算する場合には、 ガンマ関数の対数の差を取ってから指数関数 [`std::exp`](exp.md) を適用するのが賢明である。
+
+$$ \frac{\Gamma(2026)}{\Gamma(2025)} = \exp[\ln\Gamma(2026) - \ln\Gamma(2025)] $$
+
+```cpp example
+#include <cmath>
+#include <iostream>
+int main() {
+  std::cout << std::tgamma(2026.0) / std::tgamma(2025.0) << std::endl;
+  std::cout << std::exp(std::lgamma(2026.0) - std::lgamma(2025.0)) << std::endl;
+}
+```
+* std::tgamma[color ff0000]
+* std::lgamma[color 0000ff][link lgamma.md]
+
+出力例
+```
+-nan
+2025
+```
+
+上の結果では、直接ガンマ関数を計算した場合はオーバーフローによって inf / inf となり最終結果が -nan になっているが、`lgamma` を使った場合には正しい値が計算できている。
+ただし、`lgamma` は飽くまでガンマ関数の「絶対値」の対数であることに注意する。
+ガンマ関数の引数が負になる場合はガンマ関数が負の値を取りうるので符号は別に求める必要がある。
 
 ## 例
 ```cpp example
@@ -110,6 +136,8 @@ tgamma(+∞)  = inf
 
 - GCC 4.6.1 以上
 
+## 関連項目
+- ガンマ関数の絶対値の自然対数 [`lgamma`](lgamma.md)
 
 ## 参照
 - [P1467R9 Extended floating-point types and standard names](https://www.open-std.org/jtc1/sc22/wg21/docs/papers/2022/p1467r9.html)