Fortran中指数求和计算出现中间值溢出问题的技术咨询

你遇到的是统计物理玻尔兹曼权重计算的经典溢出问题，不需要复杂的对数转换，用指数项最大值平移法即可完全规避中间过程的大数溢出，特征值放大任意倍数都能稳定计算。

实现原理 #

直接计算exp(-E_i/T)时，低温度下指数自变量数值极大会直接超出浮点数表示上限触发溢出。利用比值计算的公共因子可约特性，做如下变换即可：

1. 先遍历所有项，找到所有指数自变量-E_i/T的最大值，记为max_arg
2. 把每一项的指数拆分为公共因子和余项：exp(-E_i/T) = exp(max_arg) * exp( (-E_i/T) - max_arg )
3. 配分函数Z可改写为Z = exp(max_arg) * sum( exp( (-E_i/T) - max_arg ) )，计算概率rho_i = exp(-E_i/T)/Z时，分子分母的exp(max_arg)会完全抵消，全程不需要计算这个超大的公共因子
4. 所有参与exp计算的自变量(-E_i/T) - max_arg都小于等于0，exp(非正数)的取值范围永远在(0,1]之间，从根源上杜绝上溢问题。

修改后的可运行代码 #

代码里已经默认把特征值放大300倍，运行时全程无溢出：

program testing
  implicit none
  integer, parameter :: q = SELECTED_REAL_KIND(10)
  integer, parameter :: qc = SELECTED_REAL_KIND(10)

  integer :: Ndim
  real (q)  :: Temperature
  real (q),allocatable  :: Eigenvalues (:),rho (:)
  real (q) :: Z, max_arg, arg_i
  integer :: i
  
  Ndim=4
  allocate (rho (Ndim),Eigenvalues(Ndim))
  ! 参数定义
  Temperature=0.00000158
  ! 特征值放大300倍，可根据需求继续调大
  Eigenvalues(1)=-0.000893 * 300
  Eigenvalues(2)=-0.000893 * 300
  Eigenvalues(3)=-0.000788 * 300
  Eigenvalues(4)=-0.000446 * 300

  Z = 0._q
  rho = 0._q

  ! 第一步：遍历找到所有指数项的最大自变量
  max_arg = -Eigenvalues(1)/Temperature
  do i = 2, Ndim
    arg_i = -Eigenvalues(i)/Temperature
    if (arg_i > max_arg) max_arg = arg_i
  end do

  ! 第二步：平移后累加配分函数，无溢出
  do i =1, Ndim
    arg_i = -Eigenvalues(i)/Temperature
    Z = Z + exp( arg_i - max_arg )
    write(*,*) "累加配分函数值：", Z
  end do

  ! 第三步：计算各能级占据概率，公共因子自动抵消
  do i = 1, Ndim
    arg_i = -Eigenvalues(i)/Temperature
    rho(i) = exp( arg_i - max_arg ) / Z
    write(*,*) "能级占据概率：", rho(i)
  end do

  deallocate(rho, Eigenvalues)
end

补充说明 #

• 该方法计算精度和无溢出场景下的直接计算结果完全一致，不需要升级浮点数精度
• 如果后续需要计算配分函数的对数值log(Z)，直接通过log(Z) = max_arg + log(平移后累加的Z值)计算即可，这就是工业界通用的log-sum-exp稳定实现，同样无溢出风险
• 哪怕温度继续降低、特征值放大多个数量级，只要最终概率结果是合法的0~1区间值，计算过程就不会触发溢出。

内容的提问来源于stack exchange，提问作者Román Picó