数据累加误差

Kieran · 发表于 2024-4-30 08:11:44

大家好，

我想向大家请一个问题。

我用双精度的数组从一个源文件中读取浮点数，方式如下。

[Fortran] 纯文本查看 复制代码

INTEGER，PARAMETER :: dp = SELECTED_REAL_KIND(15,14)

REAL (KIND=dp) :: arr(100,3)

INTEGER :: i

OPEN (UNIT=3, FILE='filename.dat', STATUS='OLD')

DO i = 1, 100, 1

      READ (UNIT=3, FMT=*) arr(i,:)

END DO

读取这个数组后，我想对这个数组做如下处理，将arr数组里的数据复制10000次并写入brr数组，每次复制时，都保持arr(i,1:2)不变，但会改变arr(i,3)增加一个浮点数的变量inc。最终将这个新数组brr写入一个新的文件中。方式如下。

[Fortran] 纯文本查看 复制代码

REAL (KIND=dp) :: inc = 31.494732290178

REAL (KIND=dp) :: brr(1000000,3)

DO i = 1, 10000, 1

     brr(1+(i-1)*100:i*100,1:2) = arr(1:100,1:2)

     brr(1+(i-1)*100:i*100,3) = arr(1:100,3)+DBLE(i)*inc

END DO

OPEN (UNIT=4， FILE=‘new.dat’, STATUS='UNKNOWN')

DO i = 1, 1000000, 1

     WRITE (UNIT=4, FMT=*) brr(i,:)

END DO

CLOSE (UNIT=4)

用这样的方式完成了文件new.dat的书写。new.dat文件中每相邻的两个100行数据中，第3列上的元素数值就会相差一个数值，但这个差值不是固定不变的，随着间隔的次数增多，差值会逐渐增大，由最初的31.494732290178，增大到31.794732290178，甚至更大。

而我想要做的事情是，让这个差值永远保持不变，文件new.dat中，每相邻的两个100行数据中，第3列元素的差值永远都是变量inc的数值31.494732290178。

请问要如何调整代码，实现这个一点呢？能麻烦大家给些建议吗？

谢谢。

kyra · 发表于 2024-4-30 09:43:23

本帖最后由 kyra 于 2024-4-30 09:44 编辑

REAL (KIND=dp) :: inc = 31.494732290178_dp
这是 Fortran 新手很容易弄错的一点，Fortran 中不仅变量有 kind 值，常量也有。尤其对双精度的浮点数需要特别注意。

Kieran · 发表于 2024-4-30 10:39:38

kyra 发表于 2024-4-30 09:43
REAL (KIND=dp) :: inc = 31.494732290178_dp
这是 Fortran 新手很容易弄错的一点，Fortran 中不仅变量有 k ...

非常感谢你的回复和指点。

我想问下，如果我要是对数组arr也做了一些数值上的累加。也同样不希望在累加的过程中，数组产生任何累加误差，是不是也可以在声明数组变量时也做这样限制呢？比如，arr(100,3)_dp，或者brr(1000000,3)_dp

kyra · 发表于 2024-4-30 15:05:08

变量申明的时候 real(kind=dp) :: arr 就行了。
常量才用 _dp 的形式。

Kieran · 发表于 2024-5-4 20:44:08

kyra 发表于 2024-4-30 15:05
变量申明的时候 real(kind=dp) :: arr 就行了。
常量才用 _dp 的形式。

谢谢你的回复。

我还是想再详细解释下我遇到的问题。

这是一个数据文件，文件名叫fil.dat，内容如下。

[Fortran] 纯文本查看 复制代码

     0.305999987         0.000000000         0.250000000
     0.693999988         0.000000000         0.749999982
     0.000000000         0.305999988         0.250000000
     0.000000016         0.694000020         0.749999982
     0.694000004         0.694000020         0.250000000
     0.305999988         0.305999988         0.749999982
     0.972666681         0.333333336         0.583333376
     0.360666706         0.333333336         0.083333339
     0.666666685         0.639333382         0.583333376
     0.666666680         0.027333348         0.083333339
     0.360666693         0.027333348         0.583333376
     0.972666697         0.639333382         0.083333339
     0.639333349         0.666666672         0.916666661
     0.027333349         0.666666672         0.416666697
     0.333333339         0.972666690         0.916666661
     0.333333323         0.360666684         0.416666697
     0.027333336         0.360666684         0.916666661
     0.639333351         0.972666690         0.416666697
     0.000000000         0.000000000         0.352000004
     0.000000000         0.000000000         0.647999996
    -0.000000000        -0.000000000         0.148000006
     0.000000000         0.000000000         0.851999967
     0.666666668         0.333333336         0.685333361
     0.666666668         0.333333336         0.981333354
     0.666666668         0.333333336         0.481333318
     0.666666668         0.333333336         0.185333251
     0.333333336         0.666666672         0.018666744
     0.333333336         0.666666672         0.314666675
     0.333333336         0.666666672         0.814666675
     0.333333336         0.666666672         0.518666682

我想用下面的代码读取这个文件里的数据，并对第3列的数据内容做处理，加上一个常数0.49870538。

[Fortran] 纯文本查看 复制代码

PROGRAM PROCESS
IMPLICIT NONE
INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15,14)
INTEGER :: i, j
REAL (KIND=dp) :: coor(30,3)
REAL (KIND=dp) :: inte

OPEN (UNIT=3, FILE='fil.dat', STATUS='OLD')
OPEN (UNIT=4, FILE='processed.dat', STATUS='UNKNOWN')

inte = 0.49870538_dp

DO i = 1, 30, 1
   READ (UNIT=3, FMT=*) coor(i,:)
END DO

DO i = 1, 30 ,1
   WRITE (UNIT=4, FMT=*) coor(i,1:2), coor(i,3)+inte
END DO

CLOSE (UNIT=3)
CLOSE (UNIT=4)

STOP
END PROGRAM PROCESS

比如，对第一行的第3各数字（0.250000000）加上常数0.49870538，应该得到的数值是0.748705380。但当我运行程序，并将数据输出后，这个数字却是0.74870537999999998。全部的输出数据内容如下。

[Fortran] 纯文本查看 复制代码

  0.30599998699999997        0.0000000000000000       0.74870537999999998     
  0.69399998799999996        0.0000000000000000        1.2487053619999999     
   0.0000000000000000       0.30599998800000000       0.74870537999999998     
   1.6000000000000001E-008  0.69400002000000005        1.2487053619999999     
  0.69400000399999995       0.69400002000000005       0.74870537999999998     
  0.30599998800000000       0.30599998800000000        1.2487053619999999     
  0.97266668099999998       0.33333333599999998        1.0820387560000000     
  0.36066670600000000       0.33333333599999998       0.58203871900000004     
  0.66666668500000004       0.63933338200000001        1.0820387560000000     
  0.66666667999999996        2.7333348000000000E-002  0.58203871900000004     
  0.36066669299999998        2.7333348000000000E-002   1.0820387560000000     
  0.97266669699999997       0.63933338200000001       0.58203871900000004     
  0.63933334900000005       0.66666667199999996        1.4153720409999999     
   2.7333349000000000E-002  0.66666667199999996       0.91537207700000001     
  0.33333333900000001       0.97266668999999994        1.4153720409999999     
  0.33333332300000001       0.36066668400000002       0.91537207700000001     
   2.7333336000000000E-002  0.36066668400000002        1.4153720409999999     
  0.63933335099999999       0.97266668999999994       0.91537207700000001     
   0.0000000000000000        0.0000000000000000       0.85070538399999995     
   0.0000000000000000        0.0000000000000000        1.1467053759999999     
  -0.0000000000000000       -0.0000000000000000       0.64670538599999994     
   0.0000000000000000        0.0000000000000000        1.3507053469999999     
  0.66666666799999996       0.33333333599999998        1.1840387410000000     
  0.66666666799999996       0.33333333599999998        1.4800387339999999     
  0.66666666799999996       0.33333333599999998       0.98003869799999999     
  0.66666666799999996       0.33333333599999998       0.68403863099999995     
  0.33333333599999998       0.66666667199999996       0.51737212399999999     
  0.33333333599999998       0.66666667199999996       0.81337205499999998     
  0.33333333599999998       0.66666667199999996        1.3133720549999999     
  0.33333333599999998       0.66666667199999996        1.0173720620000000

我知道可以用WRITE命令里的FMT格式，限制输出的数据格式，得到0.748705380。但我不能这么做，因为这只是我中间计算的结果，后期我还要用这个中间计算的结果，做后面的数学处理。也就是说，我想将代码中coor数组里存储的中间数据始终保持0.748705380这类的格式，而不是0.74870537999999998这类格式。

请问要如何做到呢？我的意思是问，要怎样才能保证代码运行过程中的，数组里保存的中间数据永远保持相同的位数格式呢？

能麻烦你再多给我些建议吗？

谢谢啦，盼复。

kyra · 发表于 2024-5-5 17:15:05

0.748705380
0.74870537999999998
这两个数字相比，已经有15位有效数字了。符合你 dp = SELECTED_REAL_KIND(15,14) 的预期。
你还想继续增加有效数字，可以看看你的编译器是否支持双双精度（四精度） dp = SELECTED_REAL_KIND(20,14)

不管怎么样，计算机的资源是有限的，是无法做到无限精度的。
就好像十进制无法准确的表达 1/3 一样，只能用 0.333333333 近似。
你可以了解一下 IEEE 浮点数规则的原理。

Kieran · 发表于 2024-5-8 08:37:06

kyra 发表于 2024-5-5 17:15
0.748705380
0.74870537999999998
这两个数字相比，已经有15位有效数字了。符合你 dp = SELECTED_REAL_KIND ...

谢谢你的回复。

我的目的并不是想增加精确度，也不是想增加输出的数据的位数。

而是希望，让运算的数据，永远保持相同位数的精度。尤其是在程序内部运行过程中的中间数据。比如，coor数组存储的数据，永远保持精确到小数点后第9位。第9位以后的其它位的数据可以忽略不计。这样再用这个数组里的数据，进行其它数学运算得到的中间数据也都是保留精度到小数点后第9位。

当然，最终输出的数据，我可以通过WRITE命令中FMT标签来控制精确保留到小数点后第9位。

但程序运行过程中的变量，数组coor的精确度要如何控制呢？我对这一点就不清楚了。

能麻烦你给些建议吗？

谢谢，盼复。

kyra · 发表于 2024-5-8 13:35:54

没法控制，计算机使用二进制，而不是十进制。
就像三进制计算 1/3 * 3 = 1 是准确的。但十进制就可能是 0.333333 * 3 = 0.999999 不一定准确。

[通用算法] 数据累加误差

新人勋章

美女勋章

元老勋章

热心勋章

规矩勋章

管理勋章