Fortran Coder

查看: 106|回复: 3

[通用算法] Fortran科学计算中的任意精度问题

[复制链接]

57

帖子

2

主题

0

精华

专家

F 币
391 元
贡献
170 点

规矩勋章

发表于 2020-9-6 05:32:24 | 显示全部楼层 |阅读模式
本帖最后由 风平老涡 于 2020-9-6 05:38 编辑

在许多科学计算中,特别是采用实验数据的计算,IEEE32位浮点运算已经足够。在另一些科学计算中,则需要64位浮点运算。还有些应用如数学物理,实验数学等会需要更高的精度。举几个例子,超新星爆发模拟(32-64位精度),量子化学计算(32-64位精度),多体分子运动模拟(32-120位精度),电磁散射理论(32-100位精度),常微分的泰勒算法(100-600位精度),数学物理中的伊辛积分(100-1000位精度)及实验数学(100-50000位精度及更高)。

Fortran2008标准定义了四个整数精度,即int8, int16, int32及int64,和三个实型精度real32, real64, 及real128。如下表:
Type                          描述                                                             范围                                                             精度
int8                 8位(1字节)整型                                          -128 ~ +127                                          8位(2~3位十进制)
int16             16位(2字节)整型                                      -32768 ~ +32767                                    16位(4~5位十进制)
int32             32位(4字节)整型                         -2,147,483,648 ~ +2,147,483,647                       32位(9~10位十进制)
int64             64位(8字节)整型                              -9,223,372,036,854,775,808 ~                       64位(18~19位十进制)
                                                                                       9,223,372,036,854,775,807
real32           32位(4字节)实型                       +/- (1.1754 x 10^-38 ~ 3.4028 x 10^38)               24位(6~7位十进制)
real64           64位(8字节)实型                    +/- (2.2250 x 10^-308 to 1.7976 x 10^308)             52位(15~16位十进制)
real128       128位(16字节)实型              +/- (3.3621 x 10^-4932 to 1.18973 x 10^4932)         112位(33~36位十进制)


当计算需要大于36位十进制精度时,现有精度就不能满足要求了。

为了解决有限精度问题,一些编程语言及软件提供了任意精度功能。比如Python自带任意精度整型(Bigint),任意精度实型(Bigfloat)可通过第三方的库解决。而Julia自带任意精度的整型和实型及相应的运算功能。软件有诸如Mathematica,SageMath,Maxima,MATLAB 和 Maple等。

David Smith 为 Fortran提供了一个任意精度数学模块FMZM(附件)。FMZM模块包含有任意精度的整型,实型,复型及有理数型(分数),还可以进行四则运算,三角函数,特殊函数,双曲函数,矩阵乘积,矢量内积等运算。FMZM的使用也比较容易,请自查内含的用户说明。这里举两个使用例子:
[Fortran] 纯文本查看 复制代码
program arbitrary_precision
  use fmzm
  implicit none

  type(IM) :: m                              !定义一个任意精度的整型
  type(FM) :: t1, t2, pi                     !定义任意精度的实型
  integer (kind=4) :: i, n
  character(500) :: r                        !定义一个长度为500的字符串

  
  print *, "Input n="
  read (*,*) n
  m = 1                                      !初始赋值
  do i = 1, n
     m = m * i
  end do
  call im_form('i500', m, r)                 !截取任意精度整型变量的500位精度,并转化成字符串
  write(*,'(2/,I3,A,G0)') n,"! = ", r
  call fm_set(500)                           !设定实型精度为500
  t1 = to_fm('1') / to_fm('5')               !转化常数为任意精度的实型,计算并赋值
  t2 = to_fm('1') / to_fm('239')  
  pi = 4 * ( 4 * atan(t1) - atan(t2))        !计算Pi
  call fm_form('f500.499', pi, r)            !格式化任意精度实型变量,并转化成字符串
  write(*,'(2/,A,G0)') "Pi=", r
end program arbitrary_precision


结果如下:

Input n=
253


253! = 51734609926400789218043308997295274695423561272066399607484636163134302903130041238314437828111213744932542876617296316904840977852744354743364096544589631199800576352102197345093407901685444661637384445171444589249826159309289810622514481898739824349965672944938199095203108731528570965561754517676626034976542767771987626709597099937322577683908278497279328468806763572731103332796695726049211496386749680456221513530752014396144012492800000000000000000000000000000000000000000000000000000000000000


Pi=3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588174881520920962829254091715364367892590360011330530548820466521384146951941511609433057270365759591953092186117381932611793105118548074462379962749567351885752724891227938183011949



FM6.zip

37 KB, 下载次数: 0

FM5.zip

145.44 KB, 下载次数: 0

FM4.zip

76.1 KB, 下载次数: 0

FM3.zip

101.49 KB, 下载次数: 0

FM2.zip

58.78 KB, 下载次数: 0

FM1.zip

644.16 KB, 下载次数: 0

回复

使用道具 举报

167

帖子

14

主题

0

精华

大宗师

F 币
4021 元
贡献
823 点
发表于 2020-9-8 10:01:11 | 显示全部楼层
非常感谢推荐,这个应该是很有用的
不过,真正常见的对精度敏感的地方在于矩阵计算上
如果用于矩阵计算的话,没有指令集的支持
这种自定义精度的计算,就算我们自己搭起来类似于blas,lapack的框架
那速度恐怕也完全没办法忍受啊
顺便问一句,这个有没有现成的优化好的blas,lapack框架啊

57

帖子

2

主题

0

精华

专家

F 币
391 元
贡献
170 点

规矩勋章

 楼主| 发表于 2020-9-10 08:46:13 | 显示全部楼层
liudy02 发表于 2020-9-8 10:01
非常感谢推荐,这个应该是很有用的
不过,真正常见的对精度敏感的地方在于矩阵计算上
如果用于矩阵计算的话 ...

有现成的, 不过要写一下接口。http://mplapack.sourceforge.net/

167

帖子

14

主题

0

精华

大宗师

F 币
4021 元
贡献
823 点
发表于 2020-9-10 08:48:53 | 显示全部楼层
风平老涡 发表于 2020-9-10 08:46
有现成的, 不过要写一下接口。http://mplapack.sourceforge.net/

唔,非常感谢啊
有了这个,就可以做很多事情了
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

QQ|捐赠本站|Archiver|关于我们 About Us|群聊|Fcode

GMT+8, 2020-9-27 07:42

Powered by Discuz! X3.2

© 2001-2017 Comsenz Inc.

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