Fortran 中完全斐波那契连续的问题-解网

问：

我关注这个网站已经有一段时间了。我只是想问一个关于我的 fortran90 程序输出的问题。以下程序的目标是研究连续是否定义为斐波那契连续的倒数（a（n-1）/a（n））。我应该能够说，如果 abs（a（n-1）/a（n） - g ） < e 与 e 从 1.E-6 到 1.E-20，则 g 是通过将 1 减去黄金比例 phi 或使用解析公式（sqrt（5） -1）/2 获得的完备常数。在这里，它遵循程序。

module precisione
implicit none 
   integer, parameter :: ik = selected_int_kind(33)        
   integer, parameter :: rk = selected_real_kind(33)   
   real (kind=rk),parameter :: g = (sqrt(5.) -1)/real(2)            
   real (kind=rk),parameter :: phi = (1+sqrt(5.))/real(2)  
   real,parameter :: e = 1.E-15
end module precisione

module successione                                         
use precisione
implicit none 
  contains
   recursive subroutine nonfidia(n,f,d,r)                 
    integer(kind=ik), intent(in) :: n                     
    integer(kind=ik), intent(out) :: f                    
    real(kind=rk),intent(out) :: d
    real(kind=rk),intent(out) :: r                           
    integer(kind=ik):: p,s                                
    if(n > 2) then                                      
      call nonfidia(n-1,s,d,r)                            
      call nonfidia(n-2,p,d,r)                            
      f = p+s                                          
      r = real(p)/real(s)                                
      d = abs(r-g)
    else                                                  
      f=1
      r=1
      d=1-g
    end if   
   end subroutine
end module

program calcoli
 use successione                  
 implicit none
   integer(kind=ik) :: n                               
   integer(kind=ik) :: f                               
   real(kind=rk) :: d                                     
   real(kind=rk) :: ratio                                 
   integer(kind=ik) :: i                                 
   read*, n                                                
   do i=1,n                                              
    call nonfidia(i,f,d,ratio)                            
    write(unit=800,fmt=*)i,f,ratio,d,e,g                  
    print*, i,f,ratio,d,e,g                                 
    if (d < e) then                                       
     print*, "inside range , ratio = ", ratio             
     exit                                                
    end if
   end do   
    if (d > e)then                                        
     print*,"not in range, add terms"   
    end if
   
end program

问题我不能超过 21 项，因为在这一点上，差异不知何故变为零，并且它落在第一个 if 条件中，实际上 0.0000 对于我们考虑的每个 epsilon 都是 < e。我不认为问题出在减法取消原因上，通常您可以通过修改 kind 参数使用更高的精度来解决这个问题。这是输出：

程序输出

我不知道如何解决这个问题。从 1.E-6 开始，它确实有效，输出在第 17 个周期，1.E-7 在第 19 个周期，它是连贯的，但从 1.E-8 开始，我不断收到相同的错误消息。请帮忙

 100
 1 1   1.00000000000000000000000000000000000        0.381965994834899902343750000000000000         1.00000000E-15  0.618034005165100097656250000000000000      
 2 1   1.00000000000000000000000000000000000        0.381965994834899902343750000000000000         1.00000000E-15  0.618034005165100097656250000000000000      
 3 2   1.00000000000000000000000000000000000        0.381965994834899902343750000000000000         1.00000000E-15  0.618034005165100097656250000000000000      
 4 3  0.500000000000000000000000000000000000        0.118034005165100097656250000000000000         1.00000000E-15  0.618034005165100097656250000000000000      
 5 5  0.666666686534881591796875000000000000         4.86326813697814941406250000000000000E-0002   1.00000000E-15  0.618034005165100097656250000000000000      
 6 8  0.600000023841857910156250000000000000         1.80339813232421875000000000000000000E-0002   1.00000000E-15  0.618034005165100097656250000000000000      
 7 13  0.625000000000000000000000000000000000         6.96599483489990234375000000000000000E-0003   1.00000000E-15  0.618034005165100097656250000000000000      
 8 21  0.615384638309478759765625000000000000         2.64936685562133789062500000000000000E-0003   1.00000000E-15  0.618034005165100097656250000000000000      
 9 34  0.619047641754150390625000000000000000         1.01363658905029296875000000000000000E-0003   1.00000000E-15  0.618034005165100097656250000000000000      
 10 55  0.617647051811218261718750000000000000         3.86953353881835937500000000000000000E-0004   1.00000000E-15  0.618034005165100097656250000000000000      
 11 89  0.618181824684143066406250000000000000         1.47819519042968750000000000000000000E-0004   1.00000000E-15  0.618034005165100097656250000000000000      
 12 144  0.617977499961853027343750000000000000         5.65052032470703125000000000000000000E-0005   1.00000000E-15  0.618034005165100097656250000000000000      
 13 233  0.618055582046508789062500000000000000         2.15768814086914062500000000000000000E-0005   1.00000000E-15  0.618034005165100097656250000000000000      
 14 377  0.618025779724121093750000000000000000         8.22544097900390625000000000000000000E-0006   1.00000000E-15  0.618034005165100097656250000000000000      
 15 610  0.618037164211273193359375000000000000         3.15904617309570312500000000000000000E-0006   1.00000000E-15  0.618034005165100097656250000000000000      
 16 987  0.618032813072204589843750000000000000         1.19209289550781250000000000000000000E-0006   1.00000000E-15  0.618034005165100097656250000000000000      
 17 1597  0.618034422397613525390625000000000000         4.17232513427734375000000000000000000E-0007   1.00000000E-15  0.618034005165100097656250000000000000      
 18 2584  0.618033826351165771484375000000000000         1.78813934326171875000000000000000000E-0007   1.00000000E-15  0.618034005165100097656250000000000000      
 19 4181  0.618034064769744873046875000000000000         5.96046447753906250000000000000000000E-0008   1.00000000E-15  0.618034005165100097656250000000000000      
 20 6765  0.618033945560455322265625000000000000         5.96046447753906250000000000000000000E-0008   1.00000000E-15  0.618034005165100097656250000000000000      
 21 10946  0.618034005165100097656250000000000000         0.00000000000000000000000000000000000         1.00000000E-15  0.618034005165100097656250000000000000      
 inside range , ratio =   0.618034005165100097656250000000000000

Fortran 精度

答：

2赞 Ian Bush 7/31/2020 #1

你对实际种类的使用不一致，特别是很多算术运算都是用默认种类完成的，这可能比你通过使用参数请求的精度要低。下面是一个始终如一地使用更高精度的程序。注意rk

a）常量也有一个种类，如果不指定任何内容，它们将是默认的种类。如果你想要一个不同的类型，你需要明确地说出来

b）除非您指定额外的可选参数，否则 real 函数也将转换为默认类型

c）我无法说服代码和列表在 Stackoverflow 上正确格式化在一起......

ian@eris:~/work/stack$ cat fab.f90
Module precisione
  Implicit None 
  Integer, Parameter :: ik = selected_int_Kind(33)        
  Integer, Parameter :: rk = selected_real_Kind(33)   
  Real (kind=rk),Parameter :: g = (Sqrt(5.0_rk) -1.0_rk)/2.0_rk   
  Real (kind=rk),Parameter :: phi = (1.0_rk+Sqrt(5.0_rk))/2.0_rk  
  Real (kind=rk),Parameter :: e = 1.E-15_rk
End Module precisione

Module successione                                         
  Use precisione
  Implicit None 
Contains
  Recursive Subroutine nonfidia(n,f,d,r)                 
    Integer(kind=ik), Intent(in) :: n                     
    Integer(kind=ik), Intent(out) :: f                    
    Real(kind=rk),Intent(out) :: d
    Real(kind=rk),Intent(out) :: r                           
    Integer(kind=ik):: p,s                                
    If(n > 2) Then                                      
       Call nonfidia(n-1,s,d,r)                            
       Call nonfidia(n-2,p,d,r)                            
       f = p+s                                          
       r = Real(p, Kind = Kind( r ) )/Real(s, Kind = Kind( r ) )                                
       d = Abs(r-g)
    Else                                                  
       f=1
       r=1.0_rk
       d=1.0_rk-g
    End If
  End Subroutine nonfidia
End Module successione

Program calcoli
  Use successione                  
  Implicit None
  Integer(kind=ik) :: n                               
  Integer(kind=ik) :: f                               
  Real(kind=rk) :: d                                     
  Real(kind=rk) :: ratio                                 
  Integer(kind=ik) :: i                                 
  Read*, n                                                
  Do i=1,n                                              
     Call nonfidia(i,f,d,ratio)                            
     Write(unit=800,fmt=*)i,f,ratio,d,e,g                  
     Print*, i,f,ratio,d,e,g                                 
     If (d < e) Then                                       
        Print*, "inside range , ratio = ", ratio             
        Exit                                                
     End If
  End Do
  If (d > e)Then                                        
     Print*,"not in range, add terms"   
  End If

End Program calcoli
ian@eris:~/work/stack$ gfortran -O -fcheck=all -Wall -Wextra -std=f2008 fab.f90
ian@eris:~/work/stack$ ./a.out
100
 1 1   1.00000000000000000000000000000000000        0.381966011250105151795413165634361840         9.99999999999999999999999999999999987E-0016  0.618033988749894848204586834365638160      
 2 1   1.00000000000000000000000000000000000        0.381966011250105151795413165634361840         9.99999999999999999999999999999999987E-0016  0.618033988749894848204586834365638160      
 3 2   1.00000000000000000000000000000000000        0.381966011250105151795413165634361840         9.99999999999999999999999999999999987E-0016  0.618033988749894848204586834365638160      
 4 3  0.500000000000000000000000000000000000        0.118033988749894848204586834365638160         9.99999999999999999999999999999999987E-0016  0.618033988749894848204586834365638160      
 5 5  0.666666666666666666666666666666666635         4.86326779167718184620798323010284749E-0002   9.99999999999999999999999999999999987E-0016  0.618033988749894848204586834365638160      
 6 8  0.599999999999999999999999999999999981         1.80339887498948482045868343656381790E-0002   9.99999999999999999999999999999999987E-0016  0.618033988749894848204586834365638160      
 7 13  0.625000000000000000000000000000000000         6.96601125010515179541316563436184030E-0003   9.99999999999999999999999999999999987E-0016  0.618033988749894848204586834365638160      
 8 21  0.615384615384615384615384615384615414         2.64937336527946358920221898102274545E-0003   9.99999999999999999999999999999999987E-0016  0.618033988749894848204586834365638160      
 9 34  0.619047619047619047619047619047619066         1.01363029772419941446078468198090626E-0003   9.99999999999999999999999999999999987E-0016  0.618033988749894848204586834365638160      
 10 55  0.617647058823529411764705882352941154         3.86929926365436439880952012697005884E-0004   9.99999999999999999999999999999999987E-0016  0.618033988749894848204586834365638160      
 11 89  0.618181818181818181818181818181818168         1.47829431923333613594983816180008114E-0004   9.99999999999999999999999999999999987E-0016  0.618033988749894848204586834365638160      
 12 144  0.617977528089887640449438202247191017         5.64606600072077551486321184471430512E-0005   9.99999999999999999999999999999999987E-0016  0.618033988749894848204586834365638160      
 13 233  0.618055555555555555555555555555555577         2.15668056607073509687211899174172578E-0005   9.99999999999999999999999999999999987E-0016  0.618033988749894848204586834365638160      
 14 377  0.618025751072961373390557939914163083         8.23767693347481402889445147507700748E-0006   9.99999999999999999999999999999999987E-0016  0.618033988749894848204586834365638160      
 15 610  0.618037135278514588859416445623342130         3.14652861974065482961125770397041390E-0006   9.99999999999999999999999999999999987E-0016  0.618033988749894848204586834365638160      
 16 987  0.618032786885245901639344262295081946         1.20186464894656524257207055621363756E-0006   9.99999999999999999999999999999999987E-0016  0.618033988749894848204586834365638160      
 17 1597  0.618034447821681864235055724417426583         4.59071787016030468890051788423244255E-0007   9.99999999999999999999999999999999987E-0016  0.618033988749894848204586834365638160      
 18 2584  0.618033813400125234815278647463994977         1.75349769613389308186901643182656767E-0007   9.99999999999999999999999999999999987E-0016  0.618033988749894848204586834365638160      
 19 4181  0.618034055727554179566563467492260061         6.69776593313619766331266219010299554E-0008   9.99999999999999999999999999999999987E-0016  0.618033988749894848204586834365638160      
 20 6765  0.618033963166706529538387945467591504         2.55831883186661988888980466554972085E-0008   9.99999999999999999999999999999999987E-0016  0.618033988749894848204586834365638160      
 21 10946  0.618033998521803399852180339985218040         9.77190855164759350561957988060738922E-0009   9.99999999999999999999999999999999987E-0016  0.618033988749894848204586834365638160      
 22 17711  0.618033985017357938973140873378403072         3.73253690923144596098723508783041383E-0009   9.99999999999999999999999999999999987E-0016  0.618033988749894848204586834365638160      
 23 28657  0.618033990175597086556377392580881950         1.42570223835179055821524379070670731E-0009   9.99999999999999999999999999999999987E-0016  0.618033988749894848204586834365638160      
 24 46368  0.618033988205325051470844819764804385         5.44569796733742014600833774707823291E-0010   9.99999999999999999999999999999999987E-0016  0.618033988749894848204586834365638160      
 25 75025  0.618033988957902001380262249827467253         2.08007153175675415461829093100635152E-0010   9.99999999999999999999999999999999987E-0016  0.618033988749894848204586834365638160      
 26 121393  0.618033988670443185604798400533155613         7.94516625997884338324825462573575097E-0011   9.99999999999999999999999999999999987E-0016  0.618033988749894848204586834365638160      
 27 196418  0.618033988780242682856507376866870418         3.03478346519205425012322582668039137E-0011   9.99999999999999999999999999999999987E-0016  0.618033988749894848204586834365638160      
 28 317811  0.618033988738303006852732437963934105         1.15918413518543964017040549735541358E-0011   9.99999999999999999999999999999999987E-0016  0.618033988749894848204586834365638160      
 29 514229  0.618033988754322537608830405492572618         4.42768940424357112693445786493866695E-0012   9.99999999999999999999999999999999987E-0016  0.618033988749894848204586834365638160      
 30 832040  0.618033988748203621343798191078293911         1.69122686078864328734424877133339579E-0012   9.99999999999999999999999999999999987E-0016  0.618033988749894848204586834365638160      
 31 1346269  0.618033988750540839382721984519974993         6.45991178135150154336833610760961439E-0013   9.99999999999999999999999999999999987E-0016  0.618033988749894848204586834365638160      
 32 2178309  0.618033988749648101530971893432887521         2.46746673614940932750639110506896385E-0013   9.99999999999999999999999999999999987E-0016  0.618033988749894848204586834365638160      
 33 3524578  0.618033988749989097047296779290725048         9.42488427099449250868884736278446728E-0014   9.99999999999999999999999999999999987E-0016  0.618033988749894848204586834365638160      
 34 5702887  0.618033988749858848350071980248415539         3.59998545148541172226208047535336351E-0014   9.99999999999999999999999999999999987E-0016  0.618033988749894848204586834365638160      
 35 9227465  0.618033988749908598925421457588060227         1.37507208346232224220669885957802422E-0014   9.99999999999999999999999999999999987E-0016  0.618033988749894848204586834365638160      
 36 14930352  0.618033988749889595896597819661196270         5.25230798901470444188988345745259482E-0015   9.99999999999999999999999999999999987E-0016  0.618033988749894848204586834365638160      
 37 24157817  0.618033988749896854407719255379913368         2.00620313242101427520822532300926455E-0015   9.99999999999999999999999999999999987E-0016  0.618033988749894848204586834365638160

注意：我在程序完成之前杀死了它

0赞 zizouoziz 7/31/2020

哇，太快了！但是你在b点是什么意思？对于其他一切，我都明白了。完全忘记了_rk附加到参数等。谢谢很多

0赞 Ian Bush 7/31/2020

Real（ a ）返回默认类型的实数值。Real（ a， rk ）返回一个 rk 类型的实数，假设 rk 是实数的有效类型。它们不一定相同。

0赞 JAlex 7/31/2020

好奇的是，当每个递归调用将一帧存储到堆栈中时，Fortran 是否会导致递归函数（如基于 C 的语言）的堆栈溢出，在 30-40 次递归后耗尽空间。

0赞 Ian Bush 7/31/2020

取决于实现。该标准没有说明实现如何处理这个问题，或者据我所知，Fortran 中的任何其他内容都没有说明。没有什么能阻止它变成迭代。我的信念是 C 是一样的，但我的知识在那里不太好。

上一个：x86_64 和 ARMv8.2-A 之间的浮点计算结果不同

下一个：Fortran 中数值运算中自动类型/类型转换的可移植性

Fortran 中完全斐波那契连续的问题

Problem with the complemetary Fibonacci succession in Fortran

评论

评论