为什么Fortran中的一元Horner比NumPy的更快，而二元Horner是n

Question

我想用Python执行多项式演算。numpy中的polynomial包对我来说不够快。因此，我决定重写Fortran中的几个函数，并使用f2py创建可轻松导入Python的共享库。目前，我正在对照numpy对应的单变量和二元多项式求值例程。在

在单变量例程中，我使用Horner's method，就像numpy.polynomial.polynomial.polyval一样。我已经观察到Fortran例程比numpy对应程序快的因素随着多项式阶数的增加而增加。在

在双变量程序中，我使用霍纳方法两次。首先在y中，然后在x中。不幸的是，我注意到对于增加多项式阶，numpy对应的函数会赶上并最终超过我的Fortran例程。由于numpy.polynomial.polynomial.polyval2d使用了与我类似的方法，我认为第二个观察结果很奇怪。在

我希望这个结果源于我对Fortran和f2py缺乏经验。有人知道为什么一元程序总是显得优越，而二元程序只对低阶多项式优越？在

编辑 以下是我最新更新的代码、基准脚本和性能图：

多项式

subroutine polyval(p, x, pval, nx)

    implicit none

    real(8), dimension(nx), intent(in) :: p
    real(8), intent(in) :: x
    real(8), intent(out) :: pval
    integer, intent(in) :: nx
    integer :: i

    pval = 0.0d0
    do i = nx, 1, -1
        pval = pval*x + p(i)
    end do

end subroutine polyval

subroutine polyval2(p, x, y, pval, nx, ny)

    implicit none

    real(8), dimension(nx, ny), intent(in) :: p
    real(8), intent(in) :: x, y
    real(8), intent(out) :: pval
    integer, intent(in) :: nx, ny
    real(8) :: tmp
    integer :: i, j

    pval = 0.0d0
    do j = ny, 1, -1
        tmp = 0.0d0
        do i = nx, 1, -1
            tmp = tmp*x + p(i, j)
        end do
        pval = pval*y + tmp
    end do

end subroutine polyval2

subroutine polyval3(p, x, y, z, pval, nx, ny, nz)

    implicit none

    real(8), dimension(nx, ny, nz), intent(in) :: p
    real(8), intent(in) :: x, y, z
    real(8), intent(out) :: pval
    integer, intent(in) :: nx, ny, nz
    real(8) :: tmp, tmp2
    integer :: i, j, k

    pval = 0.0d0
    do k = nz, 1, -1
        tmp2 = 0.0d0
        do j = ny, 1, -1
            tmp = 0.0d0
            do i = nx, 1, -1
                tmp = tmp*x + p(i, j, k)
            end do
            tmp2 = tmp2*y + tmp
        end do
        pval = pval*z + tmp2
    end do

end subroutine polyval3

基准.py（使用此脚本生成绘图）

^{pr2}$

结果

编辑对steabert提案的修正

subroutine polyval(p, x, pval, nx)

    implicit none

    real*8, dimension(nx), intent(in) :: p
    real*8, intent(in) :: x
    real*8, intent(out) :: pval
    integer, intent(in) :: nx

    integer, parameter :: simd = 8
    real*8 :: tmp(simd), xpower(simd), maxpower
    integer :: i, j, k

    xpower(1) = x
    do i = 2, simd
        xpower(i) = xpower(i-1)*x
    end do
    maxpower = xpower(simd)

    tmp = 0.0d0
    do i = nx+1, simd+2, -simd
        do j = 1, simd
            tmp(j) = tmp(j)*maxpower + p(i-j)*xpower(simd-j+1)
        end do
    end do

    k = mod(nx-1, simd)
    if (k == 0) then
        pval = sum(tmp) + p(1)
    else
        pval = sum(tmp) + p(k+1)
        do i = k, 1, -1
            pval = pval*x + p(i)
        end do
    end if

end subroutine polyval

编辑测试代码，以验证上面的代码对x>；1的结果不佳

import polynomial as P
import numpy.polynomial.polynomial as PP

import numpy as np

for n in xrange(2,100):
    poly1n = np.random.rand(n)
    poly1f = np.asfortranarray(poly1n)

    x = 2

    print np.linalg.norm(P.polyval(poly1f, x) - PP.polyval(x, poly1n)), '\n'

Answer 1

我猜，你的tmp数组太大了，以至于它需要二级、三级甚至主存访问，而不是缓存。最好把这些循环分开，一次只处理其中的大块（露天开采）。在