迭代与向量化正弦函数计算结果存在差异的问题排查

处理包含纬度值的NetCDF文件时，尝试计算纬度正弦值的倒数，分别用遍历逐个计算和数组向量化运算两种方式实现，但两者倒数的差值不为零，寻求差异原因。

代码如下：

import xarray as xr
import numpy as np

data = xr.open_dataset('ERA5_2000_01_01.nc')
lat = data.latitude.values
omega = 7.292e-5

f1 = np.array([2 * omega * np.sin(x * np.pi / 180.0) for x in lat])
f2 = 2 * omega * np.sin(lat * np.pi / 180.0)

diff = 1/f1 - 1/f2
print(diff)

得到的差值diff为：

[-0.00018776 -0.00050296 -0.00059619 -0.00035761  0.00011444 -0.00065717
 -0.00024923 -0.00060329  0.00018545 -0.00014504 -0.00064833 -0.00043263
  0.00061779  0.00013936 -0.00046115 -0.00041766 -0.0009715   0.0005565
 -0.00034907 -0.00123883 -0.00013027  0.00075469  0.00162418 -0.00066007
 -0.00037552 -0.00189229 -0.00101494  0.00089398  0.00063323  0.00030524
  0.00051748  0.00012772  0.00051868  0.00344283 -0.002541    0.00145573
         nan -0.00145573  0.002541   -0.00344283 -0.00051868 -0.00012772
 -0.00051748 -0.00030524 -0.00063323 -0.00089398  0.00101494  0.00189229
  0.00037552  0.00066007 -0.00162418 -0.00075469  0.00013027  0.00123883
  0.00034907 -0.0005565   0.0009715   0.00041766  0.00046115 -0.00013936
 -0.00061779  0.00043263  0.00064833  0.00014504 -0.00018545  0.00060329
  0.00024923  0.00065717 -0.00011444  0.00035761  0.00059619  0.00050296
  0.00018776]

• 浮点数精度累积差异：遍历计算时，每个纬度值单独完成弧度转换、正弦计算后再存入数组；向量化运算则是numpy对整个数组执行批量优化计算。两种方式的浮点数运算顺序、内部精度处理细节不同，会产生微小的数值偏差。
• 向量化运算的硬件优化：numpy向量化操作会利用CPU的SIMD指令集批量计算，这类优化可能调整运算顺序或舍入策略，和单个元素计算的精度表现略有不同，进而导致结果差异。
• 小数值的放大效应：omega是7.292e-5的小数值，计算出的f1、f2本身数值极小，此时两者间的微小绝对差异，在取倒数后会被放大为可观测的差值；尤其是赤道附近纬度为0时，sin(0)=0，倒数会出现NaN，对应diff中的NaN值。

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