Fortran辛普森1/3算法real(kind=4)正常real(kind=16)段错误求解

段错误和real(kind=16)本身无关，是代码本身存在数组越界、参数设计不合理的问题，低精度下没触发崩溃只是因为越界访问的内存刚好没有被系统保护，属于未定义行为碰巧正常运行，核心问题如下：

• 数组越界：子例程中do k=2,s+1,2循环内访问B(k+2)，当k最大等于s时，k+2 = s+2，但B数组只分配了s+1个元素，直接访问越界内存触发段错误。
• 整数类型滥用：integer(kind=16)在多数常规平台支持度很差，完全不需要用128位整数存储区间数，用默认32位整数就可以满足需求。
• 参数设计冗余：主程序里声明的A、B、d、s1等变量都是子例程内部使用的临时变量，不需要作为参数传递，额外传递反而可能因为变量类型、分配状态不一致出问题。
• 浮点常量精度不匹配：计算d/3.0时3.0是单精度常量，和real16运算会有精度损失，应该写成对应精度的常量。
• 辛普森算法逻辑错误：当前的奇数项累加逻辑存在错误，就算不崩溃计算结果也不符合辛普森1/3积分法则的要求。

按照以下步骤修改即可适配real(kind=16)精度：

1. 删除冗余参数，子例程内部的临时变量都放在子例程内部声明，不要放在参数列表中。
2. 修正数组越界问题，重新编写奇偶项累加逻辑，不需要额外用B数组做移位累加，直接遍历节点函数值计算奇偶和即可。
3. 把integer(kind=16)改成默认整数类型，不要使用128位整数。
4. 所有浮点常量都加_16后缀保证精度匹配。
5. 增加辛普森1/3法则的前置校验：要求输入的区间数必须为偶数。

模块文件 #

module simpsonsonethirdmod
    implicit none
    public::si
contains
    subroutine si(f, low, high, n_intervals, result_)
        real(kind=16), external :: f
        real(kind=16), intent(in) :: low, high
        integer, intent(in) :: n_intervals
        real(kind=16), intent(out) :: result_
        real(kind=16), allocatable :: node_val(:)
        real(kind=16) :: step, s_even, s_odd
        integer :: k

        if (mod(n_intervals, 2) /= 0) error stop "辛普森1/3法则要求区间数必须为偶数"

        allocate(node_val(n_intervals + 1))
        step = (high - low) / real(n_intervals, kind=16)

        ! 计算所有节点的函数值
        do k = 1, n_intervals + 1
            node_val(k) = f(low + (k-1)*step)
        end do

        ! 计算奇数项、偶数项和（从第2个节点开始计数）
        s_odd = 0.0_16
        do k = 2, n_intervals, 2
            s_odd = s_odd + node_val(k)
        end do

        s_even = 0.0_16
        do k = 3, n_intervals-1, 2
            s_even = s_even + node_val(k)
        end do

        ! 计算最终积分结果
        result_ = (step / 3.0_16) * (node_val(1) + node_val(n_intervals+1) + 4.0_16 * s_odd + 2.0_16 * s_even)

        deallocate(node_val)
    end subroutine si
end module simpsonsonethirdmod

主程序文件 #

program simpsonsonethird
    use simpsonsonethirdmod
    implicit none
    real(kind=16), external :: f
    real(kind=16) :: i, e, result_, start_time, end_time
    integer :: s

    call cpu_time(start_time)

    print *, "enter the starting value of integral"
    read(*,*) i
    print *, "enter the final value of the integral"
    read(*,*) e
    print *, "enter the number of intervals you want (must be even)"
    read(*,*) s

    call si(f, i, e, s, result_)

    print *, "The integral by simpson's rule is", result_

    call cpu_time(end_time)
    print *, "Time taken in seconds =", end_time-start_time
end program simpsonsonethird

real(kind=16) function f(x)
    real(kind=16), intent(in) :: x
    f = x**2
end function f

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