Fortran数组与非数组存储差异致分箱计数偏差问题排查

问题背景 #

在分箱统计数据的Fortran程序中，改用数组循环计算分箱上下界后，部分分箱的计数结果和之前直接赋值变量的版本偏差1-5个。测试发现：

• 直接赋值的分箱边界（如第5箱的-1.8）输出为-1.79999995
• 数组迭代计算的lower(5)输出为-1.80000019
  这种微小的浮点数差异，导致部分数据被错误划分到相邻分箱，最终计数偏差。

原分析程序 #

program analysis
        implicit none
        
        integer i, j, k, l
        double precision a, b, c, d, e
        integer binb(1:20),binc(1:20),bind1(1:20),bine(1:20)
        real lower(1:20),upper(1:20), next
    
        character(100) event
        
        upper(1)=-2.7
        lower(1)=-3.0
        binb(1:20)=0
        binc(1:20)=0
        bind1(1:20)=0
        bine(1:20)=0
        
        next=.3
        
        do l=2,20
            lower(l)=lower(l-1)+next
            upper(l)=upper(l-1)+next
        end do

        open(unit = 7, file="zpc_initial_momenta.dat")
            
        do k=1, 10
            read(7,'(A)') event
            do j=1,4000
                read(7,*) a, b, c, d, e
                
                do i=1, 20
                    if(b>=lower(i) .and. b<upper(i)) then
                        binb(i)=binb(i)+1
                    end if
                    if(c>=lower(i) .and. c<upper(i)) then 
                        binc(i)=binc(i)+1
                    end if
                    if(d>=lower(i) .and. d<upper(i)) then
                        bind1(i)=bind1(i)+1
                    end if
                    if(e>=lower(i) .and. e<upper(i)) then 
                        bine(i)=bine(i)+1
                    end if
                end do
            end do
        end do
            
        close(7)
        
        open(unit = 8, file="outputanalysis2.dat")
            Write(8,*) 'The bins in each column are as follows:'
            Write(8,*) 'FIRST COLUMN (MOMENTUM IN X DIRECTION)'
            write(8,*) binb(1:20)
            write(8,*) 'THE TOTAL COUNT IS', SUM(binb)
            Write(8,*) 'SECOND COLUMN (MOMENTUM IN Y DIRECTION)'
            write(8,*) binc(1:20)
            write(8,*) 'THE TOTAL COUNT IS', SUM(binc)
            Write(8,*) 'THIRD COLUMN (MOMENTUM IN Z DIRECTION)'
            write(8,*) bind1(1:20)
            write(8,*) 'THE TOTAL COUNT IS', SUM(bind1)
            Write(8,*) 'FOURTH COLUMN (KINETIC ENERGY)'
            write(8,*) bine(1:20)
            write(8,*) 'THE TOTAL COUNT IS', SUM(bine)
        close(8)
        
end program

复现测试程序 #

program help
    implicit none
    
    integer i
    real lower(1:5),next,nonarraylower2
    
    lower(1)=-2.1
    nonarraylower2=-0.9
    next=0.3
    
    do i=2,5
        lower(i)=lower(i-1)+next
    end do
    
    write(*,*) lower(5),nonarraylower2
    
end program

补充数据示例 #

21  0.1314E+00 -0.1232E+01  0.6258E+00  0.1388E+01
21 -0.3478E+00 -0.1268E+01 -0.3834E+00  0.1370E+01
21  0.2138E+00  0.1995E+01 -0.9115E-01  0.2009E+01
21  0.8936E+00 -0.5758E-01  0.7773E+00  0.1186E+01
21 -0.5949E+00 -0.1999E+00  0.3787E+00  0.7330E+00

原因分析 #

问题根源是二进制浮点数的精度限制：

1. 十进制的0.3无法用二进制浮点数精确表示，存储为单精度real时是一个近似值。
2. 循环迭代累加next=0.3时，每次都会引入微小的误差，多次迭代后误差累积，导致数组计算的边界值和直接赋值的理论值（编译器转换的近似值）出现差异。
3. 原程序中lower/upper是单精度real，而统计的b/c/d/e是双精度double precision，类型不匹配导致隐式转换，进一步放大了误差影响。

解决方法 #

1. 用公式直接计算边界，避免迭代累加 #

不要通过前一个值累加得到当前边界，而是用初始值加上步长的倍数，只做一次运算，彻底避免累积误差：

! 修改为双精度，统一类型
double precision lower(1:20), upper(1:20), next
lower(1) = -3.0d0
upper(1) = -2.7d0
next = 0.3d0

do l=2,20
    lower(l) = lower(1) + (l-1)*next  ! 直接计算第l个边界
    upper(l) = upper(1) + (l-1)*next
end do

2. 统一数据类型 #

把所有和边界、统计相关的变量都改成双精度double precision，避免单/双精度转换带来的额外误差：

! 原程序中real类型改为double precision
double precision lower(1:20), upper(1:20), next

3. 调整分箱判断逻辑 #

对浮点数的边界判断，避免严格的>=和<组合，可以引入极小的epsilon值调整边界，或者统一判断规则（比如所有分箱用>左边界、<=右边界，避免边界值被重复统计或漏统计）：

! 示例：用epsilon调整左边界，避免因精度问题漏判
double precision, parameter :: eps = epsilon(1.0d0)
if(b > (lower(i)-eps) .and. b <= upper(i)) then
    binb(i) = binb(i) + 1
end if

4. 手动预定义边界值 #

如果分箱数量不多，直接手动赋值所有lower和upper的元素，完全避免计算误差：

double precision lower(1:20), upper(1:20)
lower = [-3.0d0, -2.7d0, -2.4d0, ...]  ! 按顺序手动列出所有边界
upper = [-2.7d0, -2.4d0, -2.1d0, ...]

内容的提问来源于stack exchange，提问作者Garrett Leigh