Fortran Coder

标题: 浮点型精度求解决 [打印本页]

作者: SongAngel 时间: 2022-11-21 11:14
标题: 浮点型精度求解决
各位老师，大家好，本人在Fortran编程中选用一个单精度实型变量进行逐次累加计数，每次都+1，总数目为54000000，但是加至16777216时，该变量的值无法继续增加，一直保持不变。因此，我单独把这一段程序取了出来进行测试，发现计算结果更加迷惑了，发现16777216是个门槛，逐次+1越不过去，但是+2却可以。我也研究了下单精度的计算机存储方式，总数目的大小也没有超过单精度的存储范围，一定是精度的问题，但还是不能想清楚，希望各位老师帮忙解决，测试代码和结果如下：

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

    Program Test

    Real A,B,C,D

    A=16777215

    B=1

    C=2

    D=3

    Write(*,*) A,B,C,D

    Write(*,*) A+B,A+C,A+D

    A=16777216

    B=1

    C=2

    D=3

    Write(*,*) A,B,C,D

    Write(*,*) A+B,A+C,A+D

    End Program Test

测试结果：
  1.6777215E+07 1.000000    2.000000    3.000000
  1.6777216E+07  1.6777216E+07  1.6777218E+07
  1.6777216E+07 1.000000    2.000000    3.000000
  1.6777216E+07  1.6777218E+07  1.6777220E+07
请按任意键继续. . .

作者: SongAngel 时间: 2022-11-21 11:48
其实这个解决方法也很简单，就是将计数变量改为双精度，就不存在上述问题了，但是心里就是想搞清楚导致上述问题出现的原因是什么，搞不清楚太难受了。查看相关资料说是单精度的浮点数的有效位数为6~7位，这个有效位指的具体是什么呢，什么是6位，什么时候是7位，比如说下面代码的输出结果为1.000000 1.000001 1.000001 1.000002 1.000000，这是不是代表具有7为有效数字，而最后一位则是不准确的，是需要按照四舍五入进行估计的呢

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

    Program Test

    Real A,B,C,D,E

    A=1.000000

    B=1.000001

    C=1.0000014

    D=1.0000015

    E=1.0000001

    Write(*,*) A,B,C,D,E

    End Program Test

作者: 楚香饭 时间: 2022-11-21 12:55
本帖最后由楚香饭于 2022-11-21 13:00 编辑

6-7位。
因为计算机使用的时二进制，而不是十进制。
一般计算机用整型来表达整数，用定点数和浮点数来表达小数。（目前大多数系统上用浮点数，极少用定点数）
浮点数，类似二进制的科学计数法。

普遍采用的浮点数标准是 IEEEE 浮点数标准，IEEE 单精度浮点数有 32 位，双精度有 64 位。
　　两者都是分为 s e m 三段，代表符号位，指数位和尾数位。分段规则如下：
   ┏━━━┯━┯━┯━┓
   ┃ 精度 │s │e │m ┃
   ┠───┼─┼─┼─┨
   ┃单精度│1 │8 │23┃
   ┠───┼─┼─┼─┨
   ┃双精度│1 │11│52┃
   ┗━━━┷━┷━┷━┛
　　最终计算公示为：n = ( -1 ) ** s * m * ( 2 ** e )
因为单精度有23位二进制表示m段，所以转换成十进制的精度为 log10(2**23) = 6.923(也就是6到7位之间）
　　现在详细解释一下各段的含义（以单精度 50.5 为例）：
　　(1) s 段：s 段为 1，表示负数，否则表示正数。50.5 是正数，因此，s 段取为 0。
　　(2) e 段：e 段可以认为是一个整型。比如取值为 10000100 ，即十进制的 132。但参与最后计算的 e 并不是直接的 132。而是 132 于偏移量 127 的差。也就是 5。（所以才能出现负指数）（双精度的偏移量为1023）
　　(3) m 段：m 段的每一位表示一个数值，取为 0 则不计入总数，取为 1 则计入总数。
   第一位表示 1/2，第二位表示 1/4，第三位表示 1/8 ....
   比如取为：10010100000000000000000。则表示 1/2 + 1/16 + 1/64 = 0.578125。在这个基础上 + 1。即 1.578125 最终参与计算。

　　最后的计算 n = ( -1 ) ** 0 * 1.578125 * ( 2 ** 5 ) = 50.5
　　将以上三段写在一起：
   ┏━┯━━━━┯━━━━━━━━━━━━┓
   ┃s │ e │ m          ┃
   ┠─┼────┼────────────┨
   ┃0 │10000100│10010100000000000000000 ┃
   ┗━┷━━━━┷━━━━━━━━━━━━┛
　　最后写在一起：01000010010010100000000000000000
　　每隔 8 位隔开：
   ┏━━━━┯━━━━┯━━━━┯━━━━┯━━━━┓
   ┃ 字节  │01000010│01001010│00000000│00000000┃
   ┠────┼────┼────┼────┼────┨
   ┃十六进制│  42  │  4A  │  00  │  00  ┃
   ┗━━━━┷━━━━┷━━━━┷━━━━┷━━━━┛
　　最终就表示为：42 4A 00 00
　　至于储存是是否翻转为：00 00 4A 42，则依赖于具体平台。微机平台上翻转，工作站上不翻转。

作者: Transpose 时间: 2022-11-21 14:03
本帖最后由 Transpose 于 2022-11-21 14:07 编辑

再仔细研究一下ieee754你就发现是合理的，这个数是2**24
同理，双精度也会有这个限制，是2**53。
因为单精度的最小有效精度是1/2**23，所以超过2**24每个浮点数的间隔就只能是2了，+1自然没有效果，

作者: SongAngel 时间: 2022-11-21 14:07

楚香饭发表于 2022-11-21 12:55
6-7位。
因为计算机使用的时二进制，而不是十进制。
一般计算机用整型来表达整数，用定点数和浮点数来表达 ...

非常感谢您的回复，已经对精度有了一个比较清晰的了解，但是对于第一个累加的问题，应该如何解释呢，还是没有搞清楚，还望再解释下，谢谢了。

作者: SongAngel 时间: 2022-11-21 14:17

Transpose 发表于 2022-11-21 14:03
再仔细研究一下ieee754你就发现是合理的，这个数是2**24
同理，双精度也会有这个限制，是2**53。
因为单精 ...

您好，“超过2**24每个浮点数的间隔就只能是2了，+1自然没有效果”这句话没有理解的太清楚，能否再详细解释下，这个最小精度和间隔有什么具体的联系呢？谢谢了

作者: 楚香饭 时间: 2022-11-21 14:21
了解一下大数吃小数的问题。

用十进制来距离，你会容易理解一些，理解之后自己想到二进制上就行了。

1.000e20 + 1.000e16 = 1.0001e20 但只存储三位小数，则还是 1.000e20
但
1.000e20 + 1.000e17 = 1.001e20 只存储三位小数，结果是 1.001e20

为什么+1e16不行，但+1e17就行了呢？

作者: 楚香饭 时间: 2022-11-21 14:22
实际上，大量数据的累计，要么用整形（因为整形在有效范围内是完全精确表达的，不存在精度问题）
要么，用乘法。

作者: SongAngel 时间: 2022-11-21 16:15
本帖最后由 SongAngel 于 2022-11-21 17:15 编辑

楚香饭发表于 2022-11-21 14:21
了解一下大数吃小数的问题。

用十进制来距离，你会容易理解一些，理解之后自己想到二进制上就行了。

通过您举的这个例子我明白了，根源还是在于精度问题，我的具体理解是这样的：
根据您所给的关于浮点数由二进制转十进制的计算公式n = ( -1 ) ** s * m * ( 2 ** e )
可以看出，对于任意一个浮点数就是首先将其拆分成一种由很多份（1/2**n）的组合，对组合求和，然后再乘以相应的指数。因此，能否对该数据进行准确表示的关键就在于能否找到这样一个组合，具体怎么找有他内部的规则。我觉得这种形式类似于数学上的泰勒级数展开，就是把一个复杂的函数拆分成多个简单的函数进行拟合、拼接和近似，以便于数值计算。
回到这个问题，就是：
16777216（也就是2**24）按照规则换成二进制应为：s段为0；e段十进制为127+24=151，即二进制10010111；m段全部为0。
而对于16777216+1=16777217
16777217（假设能够准确表示：X*2**24）按照规则换成二进制应为：s段为0；e段十进制为127+24=151，即二进制10010111；m段的和X应等于1.000000059604644775390625（共25位），也就是只要能够找到一种（1/2**n）的组合就可以。然而，单精度的最小一份子为1/2**23=0.00000011920928955078125（共24位），可以看出无论怎么组合都无法将最后一位组合出来，只能通过使用更高精度的数据在组合。事实上，1/2**24=0.000000059604644775390625，表明单精度的m段再增加一位，由23变为24就刚刚好。
此时，16777217（1.000000059604644775390625*2**24）按照规则换成二进制应为：s段为0；e段十进制为127+24=151，即二进制10010111；m段全部为000000000000000000000001（共24位）。
但并不存在上述情况，只能使用双精度，其最小精度为1/2**52，远小于16777217所需要的1/2**24，故能够满足。
再次感谢老师的指导和帮助，谢谢。
再补充一句，在单精度下，由于16777216的m段全部为0，故比其稍大一点点点的那个值应为m段的最后一位变为1，e段保持不变，此时m段的和为1+1/2**23，所对应的十进制数为（1+1/2**23）*2**24=16777218，可以看出刚刚好把16777217给跳过去了，所以说+1不行，只能+2才可以。
再补充一句，对于上述跳过的数字有很多，但是对于该值的显示方式，经测试发现应该满足四舍五入法则，具体为判断该跳过值的m段的求和结果中与单精度的最小一份子精度相同的最后一位的大小，大于等于5则按照大值显示，小于5则按照小值显示，具体如下：
0.00000011920928955078125（单精度的最小一份子1/2**23）
1.000000059604644775390625（16777217）小于5，显示为16777216
1.000000178813934326171875（16777219）大于5，显示为16777220
1.000000298023223876953125（16777221）小于5，显示为16777220
1.000000655651092529296875（16777227）大于5，显示为16777228
1.000000774860382080078125（16777229）小于5，显示为16777228

代码及运算结果
16777216.0
16777220.0
16777220.0
16777228.0
16777228.0

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

    Program Test

    Real A,B,C,D,E

    A=16777217

    B=16777219

    C=16777221

    D=16777227

    E=16777229

    Write(*,"(F10.1)") A,B,C,D,E

    End Program Test

作者: 楚香饭 时间: 2022-11-21 17:50
你描述的基本对。但是你想复杂了。
都算不上泰勒级数，它就是个等比数列的线性组合。只不过十进制公比是10，而二进制公比是2

十进制可拆成：
[latex](-1)^s10^e\sum_{i=1}^{m}k10^{i}[/latex]
二进制就拆成
[latex](-1)^s2^e\sum_{i=1}^{m}k2^{i}[/latex]

作者: SongAngel 时间: 2022-11-21 21:50

楚香饭发表于 2022-11-21 17:50
你描述的基本对。但是你想复杂了。
都算不上泰勒级数，它就是个等比数列的线性组合。只不过十进制公比是10 ...

是的，就是等比数列的线性组合

，不过十进制和二进制的公比应该分别是0.1和0.5吧。

作者: 楚香饭 时间: 2022-11-22 08:39

SongAngel 发表于 2022-11-21 21:50
是的，就是等比数列的线性组合，不过十进制和二进制的公比应该分别是0.1和0.5吧。 ...

对对对，不严谨了，不严谨了

欢迎光临 Fortran Coder (http://bbs.fcode.cn/)