Fortran科学计算中的任意精度问题

风平老涡

在许多科学计算中，特别是采用实验数据的计算，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

liudy02

风平老涡 发表于 2020-9-10 08:46
有现成的， 不过要写一下接口。http://mplapack.sourceforge.net/

唔，非常感谢啊
有了这个，就可以做很多事情了

风平老涡

liudy02 发表于 2020-9-8 10:01
非常感谢推荐，这个应该是很有用的
不过，真正常见的对精度敏感的地方在于矩阵计算上
如果用于矩阵计算的话 ...

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

liudy02

非常感谢推荐，这个应该是很有用的
不过，真正常见的对精度敏感的地方在于矩阵计算上
如果用于矩阵计算的话，没有指令集的支持
这种自定义精度的计算，就算我们自己搭起来类似于blas,lapack的框架
那速度恐怕也完全没办法忍受啊
顺便问一句，这个有没有现成的优化好的blas,lapack框架啊