Fortran科学计算中的任意精度问题
本帖最后由 风平老涡 于 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的使用也比较容易，请自查内含的用户说明。这里举两个使用例子: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
非常感谢推荐，这个应该是很有用的
不过，真正常见的对精度敏感的地方在于矩阵计算上
如果用于矩阵计算的话，没有指令集的支持
这种自定义精度的计算，就算我们自己搭起来类似于blas,lapack的框架
那速度恐怕也完全没办法忍受啊
顺便问一句，这个有没有现成的优化好的blas,lapack框架啊
liudy02 发表于 2020-9-8 10:01
非常感谢推荐，这个应该是很有用的
不过，真正常见的对精度敏感的地方在于矩阵计算上
如果用于矩阵计算的话 ...
有现成的， 不过要写一下接口。http://mplapack.sourceforge.net/ 风平老涡 发表于 2020-9-10 08:46
有现成的， 不过要写一下接口。http://mplapack.sourceforge.net/
唔，非常感谢啊
有了这个，就可以做很多事情了
页:
[1]