ARM Fortran 正式版出る！･･･ が、ページがバタ臭くてｗ

ARM が HPC 用コンパイラツール一式を正式出荷したようです。

www.hpcwire.com

ハイライト孫引き引用

• Arm Fortran Compiler,released as a beta compiler in June 2017, is now a fully supported commercial Fortran compiler with Fortran 2003 and prior standards support.
• Arm C/C++ Compileris an LLVM-based commercial compiler with support for C++ 14 standard and tuned for server and HPC workloads on a wide-range of Arm-based platforms.
• The Arm Performance Librariesare processor-optimized math libraries with BLAS, LAPACK and FFT functionality for HPC applications on Arm-based hardware.
• Arm Forge(formerly Allinea Forge) now comprises Arm DDT, the powerful scalable debugger, and Arm MAP, the parallel profiler for debugging, profiling, and optimizing applications.
• Arm Performance Reports(formerly Allinea Performance Reports) is a non-intrusive performance tool that can analyze the applications running on your system to seek out inefficiencies and pinpoint exactly where to focus optimization work.

しか～し、ARM の当該 Fortran ページに行ってみると、写真多用のバタ臭い今風の Web ページで Fortran らしくなく、落ち着かないことこの上ない。眩暈がしてくるｗ

www.arm.com

作動システムに入ってはいないが Raspberry pi で 64bit 版 Linux 載せれば動くのか、動かぬのか。求人柱。

Arm Fortran Compiler supports:

• Fortran 2003 and prior standards
• Partial support for Fortran 2008
• OpenMP 3.1
• 64-bit Arm platforms including Cavium ThunderX2 and Qualcomm Centriq.
• Full support for SVE, an Arm architecture extension suited for HPC
• Leading Linux distributions including Red Hat 7.3+, SLES 12+ and Ubuntu 16.04+

ζ(3) アペリーの定数

中央二項係数でアペリーの定数

中央二項係数の逆数和（その６）
より、
ζ（3）＝Σ１／ｎ^3＝5/2Σ(-1)^(n-1)/n^3(2n,n)
を用いて、ζ（3）アペリーの定数を計算してみます。

実行結果

 Central Binomial Coefficients 2mCm
   1                                       2.
   2                                       6.
   3                                      20.
   4                                      70.
   5                                     252.
   6                                     924.
   7                                    3432.
   8                                   12870.
   9                                   48620.
  10                                  184756.
  11                                  705432.
  12                                 2704156.
  13                                10400600.
  14                                40116600.
  15                               155117520.
  16                               601080390.
  17                              2333606220.
  18                              9075135300.
  19                             35345263800.
  20                            137846528820.
  21                            538257874440.
  22                           2104098963720.
  23                           8233430727600.
  24                          32247603683100.
  25                         126410606437752.
  26                         495918532948104.
  27                        1946939425648112.
  28                        7648690600760440.
  29                       30067266499541040.
  30                      118264581564861424.
  31                      465428353255261088.
  32                     1832624140942590534.
  33                     7219428434016265740.
  34                    28453041475240576740.
  35                   112186277816662845432.
  36                   442512540276836779204.
  37                  1746130564335626209832.
  38                  6892620648693261354600.
  39                 27217014869199032015600.
  40                107507208733336176461620.
  41                424784580848791721628840.
  42               1678910486211891090247320.
  43               6637553085023755473070800.
  44              26248505381684851188961800.
  45             103827421287553411369671120.
  46             410795449442059149332177040.
  47            1625701140345170250548615520.
  48            6435067013866298908421603100.
  49           25477612258980856902730428600.
  50          100891344545564193334812497256.
  51          399608854866744452032002440112.
  52         1583065848125949175357548128136.
  53         6272525058612251449529907677520.
  54        24857784491537440929618523018320.
  55        98527218530093856775578873054432.
  56       390590044887157789360330532465784.
  57      1548655265692941410446222812934512.
  58      6141219157058215937976400809912720.
  59     24356699707654619143838606602026720.
  60     96614908840363322603893139521372656.
 zeta(3) Apery's constant
 s = sum 1/((-1)**(i + 1) * 2iCi * i^2)
   1.20205690315959428539973816151145
続行するには何かキーを押してください . . .

最後の桁まであっています。
ζ(3) = 1.202056903159594285399738161511449990764986292
Apéry's constant - Wikipedia

ソース・プログラム

    program pascal 
      use, intrinsic :: iso_fortran_env
      implicit none
      integer, parameter :: kq = real128
      real(kq), parameter :: pi = 4.0_kq * atan(1.0_kq)
      real(kq), allocatable :: tri(:), xCy(:)
      integer(int64) :: i
      real(kq) :: s 
      xCy = [real(kq)::]
      tri = [1]
      do i = 1, 120
        tri = [tri, 0.0_kq] + [0.0_kq, tri]
        if (mod(i, 2_int64) == 0) xCy = [xCy, tri(i / 2 + 1)] ! 2mCm
      end do  
      !
      print *, 'Central Binomial Coefficients 2mCm'
      print '(i4, 1x, f40.0)', (i, xCy(i), i = 1, size(xCy))
 
      s = 0.0_kq
      do i = 1, size(xCy)
        s = s + (-1)**(i + 1) * 1.0_kq / (xCy(i) * i**3)   
      end do    
      s = s * 5 / 2.0_kq
      print *, 'zeta(3) Apery''s constant'
      print *, 's = sum 1/((-1)**(i + 1) * 2iCi * i^2)', s
    end program pascal 

最近の Fortran 本

第3版 有限要素法による流れのシミュレーション OpenMPに基づくFortranソースコード付

建築構造設計・解析入門 Fortran 解析プログラム付

Introduction to Computational Economics Using Fortran

Computational Economics

Computational Physics: With Worked Out Examples in FORTRAN and MATLAB (De Gruyter Textbook)

Introduction to Molecular Dynamics Simulations: A Practical Guide Using C/C++, FORTRAN, and Python

CRCの本は予告通りの期限に出たためしがないｗ Fortran本が３冊キャンセルされてるし～

Classical Fortran: Programming for Engineering and Scientific Applications, Second Edition

いつの間にか第二版に。

New to the Second Edition

Additional case study on file I/O
More about CPU timing on Pentium processors
More about the g77 compiler and Linux

2017年とは思えない第二版での追加。

Using OpenMP -- The Next Step: Affinity, Accelerators, Tasking, and SIMD (Scientific and Engineering Computation)

【ニュース遅報】Forcheck 社買収されてた

新聞によりますと、Fortran ソースコード静的解析プログラムの Forcheck を開発販売していた 蘭 Forcheck社が今年の１月に買収されたそうです。そしてその買収会社がまた買収された模様です。
www.tomihisa.co.jp

Forcheck は地味に Fortran2008 まで対応しているようです。Fortran90 時代以来忘れていましたが、まだ生きていて驚きｗ　昔、行数限定の試供品が提供されていたような･･･？

www.forcheck.nl

GPL のソース配布がなされていて、販売ページが死んでいるので、オープンソース化したのでしょうか？謎です。
　

中央二項係数 再び

中央二項係数

中央二項係数とは、二項係数の丁度真ん中の値になります。2nＣn とも書き表せます。
以前も計算したのですが、32bit 整数の範囲だったのであまり多くは求められませんでした。
fortran66.hatenablog.com

今回 64bit 整数で計算してもいいのですが、ひとひねりして 128bit 実数の仮数部のビット数が 64bit より長いことを利用して、指数部を使わず仮数部の範囲で求めてみることにします。

また、せっかく求めたので中央二項係数を使って、特殊な数値を計算して検算とします。

ここのページを参考にして、 を求めます。
中央二項係数の逆数和（その２）

実行結果

 Central Binomial Coefficients 2mCm
   1                                       2.
   2                                       6.
   3                                      20.
   4                                      70.
   5                                     252.
   6                                     924.
   7                                    3432.
   8                                   12870.
   9                                   48620.
  10                                  184756.
  11                                  705432.
  12                                 2704156.
  13                                10400600.
  14                                40116600.
  15                               155117520.
  16                               601080390.
  17                              2333606220.
  18                              9075135300.
  19                             35345263800.
  20                            137846528820.
  21                            538257874440.
  22                           2104098963720.
  23                           8233430727600.
  24                          32247603683100.
  25                         126410606437752.
  26                         495918532948104.
  27                        1946939425648112.
  28                        7648690600760440.
  29                       30067266499541040.
  30                      118264581564861424.
  31                      465428353255261088.
  32                     1832624140942590534.
  33                     7219428434016265740.
  34                    28453041475240576740.
  35                   112186277816662845432.
  36                   442512540276836779204.
  37                  1746130564335626209832.
  38                  6892620648693261354600.
  39                 27217014869199032015600.
  40                107507208733336176461620.
  41                424784580848791721628840.
  42               1678910486211891090247320.
  43               6637553085023755473070800.
  44              26248505381684851188961800.
  45             103827421287553411369671120.
  46             410795449442059149332177040.
  47            1625701140345170250548615520.
  48            6435067013866298908421603100.
  49           25477612258980856902730428600.
  50          100891344545564193334812497256.
  51          399608854866744452032002440112.
  52         1583065848125949175357548128136.
  53         6272525058612251449529907677520.
  54        24857784491537440929618523018320.
  55        98527218530093856775578873054432.
  56       390590044887157789360330532465784.
  57      1548655265692941410446222812934512.
  58      6141219157058215937976400809912720.
  59     24356699707654619143838606602026720.
  60     96614908840363322603893139521372656.
 s = sum 1/((-1)**(i + 1) * 2iCi * i^2)
  0.463129641154388784993858142463065
   = 2 * {arcsinh(1/2)}^2
  0.463129641154388784993858142463066
続行するには何かキーを押してください . . .

最後の桁以外は一致しています。

ソース・プログラム

    program pascal 
      use, intrinsic :: iso_fortran_env
      implicit none
      integer, parameter :: kq = real128
      real(kq), parameter :: pi = 4.0_kq * atan(1.0_kq)
      real(kq), allocatable :: tri(:), xCy(:)
      integer(int64) :: i
      real(kq) :: s 
      xCy = [real(kq)::]
      tri = [1]
      do i = 1, 120
        tri = [tri, 0.0_kq] + [0.0_kq, tri]
        if (mod(i, 2) == 0.0) xCy = [xCy, tri(i / 2 + 1)] ! 2mCm
      end do  
      !
      print *, 'Central Binomial Coefficients 2mCm'
      print '(i4, x, f40.0)', (i, xCy(i), i = 1, size(xCy))
 
      s = 0.0_kq
      do i = 1, size(xCy)
        s = s + (-1)**(i + 1) * 1.0_kq / (xCy(i) * i * i)   
      end do    
      print *, 's = sum 1/((-1)**(i + 1) * 2iCi * i^2)', s
      print *, '  = 2 * {arcsinh(1/2)}^2              ', 2 * asinh(0.5_kq)**2
    end program pascal 

Modern Fortran Explained (Numerical Mathematics and Scientific Computation)

Numerical Computing with Modern Fortran (Applied Mathematics)

Introduction to Programming with Fortran: With Coverage of Fortran 90, 95, 2003, 2008 and 77

Fortran 利用状況情報サイト

Companies using Fortran

商利用の情報のサイトのようです。
idatalabs.com

利用者の多くは、科学技術研究機関・重厚長大産業の模様です。

Wilson の定理で素数を求める

Wilson の定理

計算量が多くて役に立たないとされる Wilson の定理をたわむれに計算して見ます。

Wilson の定理とは、素数の時
(p-1)! ≡ -1 (mod p)
が成り立ち、合成数では 0 になります。

実行結果

10^3 以下の素数リストと個数

2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 101 103 107 109 113 127 131 137 139 149 151 157 163 167 173 179 181 191 193 197 199 211 223 227 229 233 239 241 251 257 263 269 271 277 281 283 293 307 311 313 317
331 337 347 349 353 359 367 373 379 383 389 397 401 409 419 421 431 433 439 443 449 457 461 463 467 479 487 491 499 503 509 521 523 541 547 557 563 569 571 577 587 593 599 601 607 613 617 619 631 641 643 647 653 659 661 673 677 683 691
 701 709 719 727 733 739 743 751 757 761 769 773 787 797 809 811 821 823 827 829 839 853 857 859 863 877 881 883 887 907 911 919 929 937 941 947 953 967 971 977 983 991 997

the number of primes less than   1000 =    168
続行するには何かキーを押してください . . .

ソースプログラム

素数の上限が 10**4 を越えたあたりで、4byte 整数のオーバーフローが起こるので、一部 integer(8) :: m で 8byte 整数を使っています。

関数 iwilson は、引数が素数の時はその引数の値を、合成数の時は 0 を返します。
またこの関数は pure elemental なので並列実行出来るはずですが、Thereshold for auto-parallelization の閾値を下げないと、自動並列化してくれません。

    module m_wilson
      implicit none
    contains  
      integer pure elemental function iwilson(k)
        integer, intent(in) :: k
        integer :: i
        integer(8) :: m
        m = 1
        do i = 2, k - 1
          m = mod(i * m, k)     
          ! if (m == 0) exit
        end do  
        if (m /= 0) then 
          iwilson = k 
        else
          iwilson = 0
        end if  
      end function iwilson
    end module m_wilson
    
    program wilson
      use m_wilson
      implicit none
      integer, parameter :: n = 10**3
      integer, allocatable :: m(:), iprim(:)
      integer :: i
      m = iwilson([2, (i, i = 3, n, 2)])
      iprim = pack(m, m /= 0)
      print '(*(g0, x))', iprim
      print '(/, a, i7, a, i7)', 'the number of primes less than', n, ' =', size(iprim)
    end program wilson

Number Theory (Dover Books on Mathematics)

An Introduction to the Theory of Numbers

数論入門―証明を理解しながら学べる (ブルーバックス)