Fortran Coder

查看: 186|回复: 11

[数值问题] 浮点型精度求解决

[复制链接]

19

帖子

5

主题

0

精华

入门

F 币
108 元
贡献
61 点
发表于 2022-11-21 11:14:56 | 显示全部楼层 |阅读模式
各位老师,大家好,本人在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
请按任意键继续. . .

19

帖子

5

主题

0

精华

入门

F 币
108 元
贡献
61 点
 楼主| 发表于 2022-11-21 11:48:05 | 显示全部楼层
其实这个解决方法也很简单,就是将计数变量改为双精度,就不存在上述问题了,但是心里就是想搞清楚导致上述问题出现的原因是什么,搞不清楚太难受了。查看相关资料说是单精度的浮点数的有效位数为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

648

帖子

4

主题

0

精华

大师

农村外出务工人员

F 币
379 元
贡献
186 点

新人勋章爱心勋章水王勋章元老勋章热心勋章

发表于 2022-11-21 12:55:18 | 显示全部楼层
本帖最后由 楚香饭 于 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,则依赖于具体平台。微机平台上翻转,工作站上不翻转。

82

帖子

2

主题

1

精华

大师

Vim

F 币
657 元
贡献
325 点

规矩勋章

发表于 2022-11-21 14:03:42 | 显示全部楼层
本帖最后由 Transpose 于 2022-11-21 14:07 编辑

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

19

帖子

5

主题

0

精华

入门

F 币
108 元
贡献
61 点
 楼主| 发表于 2022-11-21 14:07:48 | 显示全部楼层
楚香饭 发表于 2022-11-21 12:55
6-7位。
因为计算机使用的时二进制,而不是十进制。
一般计算机用整型来表达整数,用定点数和浮点数来表达 ...

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

19

帖子

5

主题

0

精华

入门

F 币
108 元
贡献
61 点
 楼主| 发表于 2022-11-21 14:17:19 | 显示全部楼层
Transpose 发表于 2022-11-21 14:03
再仔细研究一下ieee754你就发现是合理的,这个数是2**24
同理,双精度也会有这个限制,是2**53。
因为单精 ...

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

648

帖子

4

主题

0

精华

大师

农村外出务工人员

F 币
379 元
贡献
186 点

新人勋章爱心勋章水王勋章元老勋章热心勋章

发表于 2022-11-21 14:21:06 | 显示全部楼层
了解一下大数吃小数的问题。

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

1.000e20 + 1.000e16 = 1.0001e20 但只存储三位小数,则还是 1.000e20

1.000e20 + 1.000e17 = 1.001e20 只存储三位小数,结果是 1.001e20

为什么+1e16不行,但+1e17就行了呢?

648

帖子

4

主题

0

精华

大师

农村外出务工人员

F 币
379 元
贡献
186 点

新人勋章爱心勋章水王勋章元老勋章热心勋章

发表于 2022-11-21 14:22:02 | 显示全部楼层
实际上,大量数据的累计,要么用整形(因为整形在有效范围内是完全精确表达的,不存在精度问题)
要么,用乘法。

19

帖子

5

主题

0

精华

入门

F 币
108 元
贡献
61 点
 楼主| 发表于 2022-11-21 16:15:40 | 显示全部楼层
本帖最后由 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

648

帖子

4

主题

0

精华

大师

农村外出务工人员

F 币
379 元
贡献
186 点

新人勋章爱心勋章水王勋章元老勋章热心勋章

发表于 2022-11-21 17:50:27 | 显示全部楼层
你描述的基本对。但是你想复杂了。
都算不上泰勒级数,它就是个等比数列的线性组合。只不过十进制公比是10,而二进制公比是2

十进制可拆成:
(-1)^s10^e\sum_{i=1}^{m}k10^{i}
二进制就拆成
(-1)^s2^e\sum_{i=1}^{m}k2^{i}
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

捐赠本站|Archiver|关于我们 About Us|小黑屋|Fcode ( 京ICP备18005632-2号 )

GMT+8, 2022-12-9 03:11

Powered by Tencent X3.4

© 2013-2022 Tencent

快速回复 返回顶部 返回列表