2021-05-24

【メモ帳】unicode 点字で阿部さん

unicode 表示　

Apple の terminal は unicode 文字を表示してくれるので、unicode の点字文字を使ったグラフィックスに挑戦してみることにします。ここで点字文字は 4 行 2 列の点で出来ています。

Windows の DOS 窓や WSL コンソールでは unicode 文字がうまく表示されないのですが、Windows Terminal や Visual Studio Code の Terminal では表示可能のようでしたので、こちらでも挑戦してみました。

unicode 点字

以前 Julia の unicode plots を見て便利そうだなと思っておりました。 docs.juliaplots.org

調べてみると、元々のアイデアは python 用の点字プロットにあるようでした。 github.com

en.wikipedia.org

Fortran と unicode

Fortran は Fortran 90 の時点で、日本語の JIS 漢字を念頭に多バイト文字種を ISO 規格に取り入れるという先進性を持っていたのですが、先進的過ぎて unicode 登場より前に決めた規格なのでその後停滞して未だに unicode は手際よく扱えません。

ここでは、2，3文字ばかり点字フォントを Fortran のソースファイルにコピペしてきて文字コードを読み取り、力ずくで経験的に点字フォント出力を実現しました。文字コードを Little endian 整数としてみると、64文字毎の固まりが 256 だけ間隔を空けて並んでおりました。

これは本来の文字コード表とやや違っているので、何か endian の変換がうまくいっていないようです。もしかしたら 4 バイトを反転させるのではなく 2 バイト２組で反転させるのかもしれません。文字コードについては興味も関心も知識も無いので、今後の課題とします。

unicode に紙テープやパンチカードの穴の表示を入れることを、ここに提案したいと思います。

このフォント出力の問題があり、プログラムの移植はだるかったので、ゼロから作りました。点を打てるようにすれば、以降はこれまで作ってきたルーチンをそのまま使えるので楽です。

とりま阿部さんを描きます。 fortran66.hatenablog.com

Apple M1 Mac での出力

fpm (fortran package manager ) を使ってみたかったのですが x86_64 用しかなく、arm 用の gfortran と組み合わせが悪いようでうまくいきませんでした。今後の課題としたいと思います。

f:id:fortran66:20210524204603p:plain

Windows WSL での出力

fpm を使ってみました。ファイル名・ファイル構成などを縛られるのでメンドイです。

f:id:fortran66:20210524210405p:plain

ウホッ！！いい男たち~ヤマジュン・パーフェクト

作者:山川純一
発売日: 2003/11/01
メディア: コミック

［改訂新版］プログラマのための文字コード技術入門 WEB+DB PRESS plus

作者:矢野啓介
発売日: 2018/12/22
メディア: Kindle版

群書類従第363-365 (国会図書館コレクション)

作者:不明
発売日: 2017/08/19
メディア: Kindle版

ソース・コード

gfortran, intel fortran どちらでも行けます。

描画領域からのはみ出しチェックをしていないので、はみ出すと配列はみだしで死にます。

いちおう、type-bound procedure を使って OOP 風にしてあります。複数 figure を描いてゆけます。

module uniplot
    implicit none
    private
    public :: fig_t
    type :: fig_t
        private
        integer :: nx, ny
        integer, allocatable :: array(:, :)
    contains
        procedure :: init     
        procedure :: point 
        procedure :: show  
        procedure :: line0
        procedure :: line
    end type fig_t
contains
    subroutine init(fig, nx, ny)
        class(fig_t), intent(out) :: fig 
        integer, intent(in) :: nx, ny
        fig%nx = nx
        fig%ny = ny 
        allocate(fig%array(0:(nx+1)/2, 0:(ny+3)/4) )
    end subroutine init  

    subroutine point(fig, ix, iy)
        class(fig_t), intent(in out) :: fig 
        integer, intent(in) :: ix, iy
        integer :: iax, iay
        iax = ix / 2
        iay = iy / 4
        fig%array(iax, iay) = ior(fig%array(iax, iay), icode(mod(ix, 2), mod(iy, 4))) 
    end subroutine point

    pure elemental integer function icode(kx, ky)
        integer, intent(in) :: kx, ky
        if (ky == 3) then
            icode = 64 + 64 * kx
        else ! 0, 1, 2
            icode = 2**(ky + 3*kx)  
        end if
    end function icode

    subroutine line0(fig, ix0, iy0, ix1, iy1)
        class(fig_t), intent(in out) :: fig 
        integer, intent(in) :: ix0, iy0, ix1, iy1
        integer :: i, ix, iy, nx, ny
        real :: d
        nx = ix1 - ix0
        ny = iy1 - iy0  
        if (abs(nx) < abs(ny)) then
            d = nx / real(ny)
            do i = 0, abs(ny)
                ix = nint(ix0 + d * sign(i, ny))
                iy = iy0 + sign(i, ny)
                call fig%point(ix, iy)
              end do  
        else
            d = ny / real(nx)
            do i = 0, abs(nx)
                iy = nint(iy0 + d * sign(i, nx))
                ix = ix0 + sign(i, nx)
                call fig%point(ix, iy)
            end do  
        end if  
    end subroutine line0

    subroutine show(fig)
        class(fig_t), intent(in) :: fig 
        integer :: iy
        do iy = 0, ubound(fig%array, 2)
            print '(*(a))', reverse_endian(shift_code(fig%array(:, iy)))
        end do
    end subroutine show    

    pure elemental integer function shift_code(k)
        integer, intent(in) :: k
        integer, parameter :: n0 = Z'E2A080' !14852224
        shift_code = n0 + 256 * (k /64) + mod(k, 64)  !E2A180, E2A280, E2A380      
    end function shift_code    
    
    pure elemental character(len = 4) function reverse_endian(i)
        integer, intent(in) :: i
        character:: tmp(4)
        tmp = transfer(i, ' ', size = 4)
        reverse_endian = transfer(tmp(4:1:-1), '    ')  !array 4 to len 4
    end function reverse_endian
    
    
    subroutine line(fig, x, y, ipen)
        class(fig_t), intent(in out) :: fig
        real, intent(in) :: x, y
        integer, intent(in) :: ipen
        integer, save :: ix0 = 0, iy0 = 0 
        integer :: ix, iy
        real, parameter :: xn = 80.0, yn = 100.0, fx = 1.0, fy = 0.85
        ix = nint( fx * x + xn)      
        iy = nint(-fy * y + yn)
        if (ipen == 1) call fig%line0(ix0, iy0, ix, iy)
        ix0 = ix
        iy0 = iy
    end subroutine line

end module uniplot

program uniplot_main
    use :: uniplot
    implicit none
    type(fig_t) :: fig1

    call fig1%init(140, 140)
    ! chin chin
    call fig1%line(  0.0,  11.0, 0)
    call fig1%line(  0.0,   8.0, 1)    

    call fig1%line( -8.0, -26.5, 0)
    call fig1%line(-11.0, -24.0, 1)
    call fig1%line(  8.0, -26.5, 0)
    call fig1%line( 11.0, -24.0, 1)    

    call fig1%line(  5.5, -26.0, 0)
    call fig1%line(  2.0, -36.0, 1)
    call fig1%line( -5.5, -26.0, 0)
    call fig1%line( -2.0, -36.0, 1)    

    call fig1%line( 20.0, -18.0, 0)
    call fig1%line(  8.5, -19.0, 1)
    call fig1%line(  4.0, -21.5, 1)
    call fig1%line(  0.0, -23.0, 1)
    call fig1%line( -4.0, -21.5, 1)
    call fig1%line( -8.5, -19.0, 1)
    call fig1%line(-20.0, -18.0, 1)    

    call fig1%line( 12.0, -16.0, 0)
    call fig1%line( 22.0,  14.0, 1)
    call fig1%line(-12.0, -16.0, 0)
    call fig1%line(-22.0,  14.0, 1)
     
    call fig1%line( 53.0,  -9.0, 0)
    call fig1%line( 28.5,   1.0, 1)
    call fig1%line( 28.5, -14.0, 0)
    call fig1%line( 28.5,  25.0, 1)
    call fig1%line( 28.5,  33.0, 0)
    call fig1%line( 28.5,  76.0, 1)
    call fig1%line( -2.5,  76.0, 1)
    call fig1%line( -2.5,  72.0, 1)
    call fig1%line( -0.5,  68.0, 1)
    call fig1%line(  1.0,  66.0, 1)
    call fig1%line( -1.5,  67.0, 1)
    call fig1%line( -4.0,  72.0, 1)
    call fig1%line( -4.0,  76.0, 1)
    call fig1%line(-28.5,  76.0, 1)
    call fig1%line(-28.5,  33.0, 1)    

    call fig1%line(-28.5,  25.0, 0)
    call fig1%line(-28.5, -14.0, 1)
    call fig1%line(-53.0,  -9.0, 0)
    call fig1%line(-28.5,   1.0, 1)     

    call fig1%line(  0.0,   0.0, 0)
    call fig1%line(  6.5,   0.0, 1)
    call fig1%line( 10.0,   3.0, 1)
    call fig1%line( 15.0,   7.0, 1)
    call fig1%line( 22.0,  14.0, 1)
    call fig1%line( 28.5,  25.0, 1)
    call fig1%line( 31.0,  26.0, 1)
    call fig1%line( 35.0,  34.0, 1)
    call fig1%line( 38.0,  44.0, 1)
    call fig1%line( 38.0,  53.0, 1)
    call fig1%line( 36.0,  55.0, 1)
    call fig1%line( 32.0,  55.0, 1)
    call fig1%line( 28.5,  51.0, 1)    

    call fig1%line(  0.0,   0.0, 0)
    call fig1%line( -6.5,   0.0, 1)
    call fig1%line(-10.0,   3.0, 1)
    call fig1%line(-15.0,   7.0, 1)
    call fig1%line(-22.0,  14.0, 1)
    call fig1%line(-28.5,  25.0, 1)
    call fig1%line(-31.0,  26.0, 1)
    call fig1%line(-35.0,  34.0, 1)
    call fig1%line(-38.0,  44.0, 1)
    call fig1%line(-38.0,  53.0, 1)
    call fig1%line(-36.0,  55.0, 1)
    call fig1%line(-32.0,  55.0, 1)
    call fig1%line(-28.5,  51.0, 1)    

    call fig1%line( 34.0,  55.0, 0)
    call fig1%line( 34.0,  76.0, 1)
    call fig1%line( 31.0,  82.0, 1)
    call fig1%line( 26.0,  87.0, 1)
    call fig1%line( 21.0,  91.0, 1)
    call fig1%line( 15.0,  95.0, 1)
    call fig1%line(  0.0,  95.0, 1)    

    call fig1%line(-34.0,  55.0, 0)
    call fig1%line(-34.0,  76.0, 1)
    call fig1%line(-31.0,  82.0, 1)
    call fig1%line(-26.0,  87.0, 1)
    call fig1%line(-21.0,  91.0, 1)
    call fig1%line(-15.0,  95.0, 1)
    call fig1%line(  0.0,  95.0, 1)    

    ! nose
    call fig1%line( -5.0, 41.0, 0)
    call fig1%line( -4.0, 40.0, 1)
    call fig1%line(  0.0, 40.0, 1)
    call fig1%line(  2.0, 42.0, 1)
    call fig1%line(  3.0, 42.0, 1)
    call fig1%line(  5.5, 37.0, 1)
    call fig1%line(  5.0, 37.0, 1)
    call fig1%line(  4.0, 38.5, 1)
    call fig1%line(  0.5, 38.5, 1)
    call fig1%line(  0.0, 37.0, 1)
    call fig1%line( -2.0, 37.0, 1)
    call fig1%line( -3.0, 38.5, 1)
    call fig1%line( -7.0, 38.5, 1)
    call fig1%line( -8.0, 37.0, 1)
    call fig1%line( -8.5, 40.0, 1)
    call fig1%line( -4.0, 45.0, 1)
    call fig1%line( -4.0, 54.0, 1)
    call fig1%line( -5.0, 55.0, 1)    

    call fig1%line( -6.0, 53.0, 0)
    call fig1%line( -5.0, 53.0, 1)
    call fig1%line( -5.0, 47.0, 1)
    call fig1%line( -7.5, 46.0, 1)
    call fig1%line( -6.0, 53.0, 1)    

    ! left eye
    call fig1%line(-24.0, 55.0, 0)
    call fig1%line(-22.0, 53.0, 1)
    call fig1%line(-17.0, 54.5, 1)
    call fig1%line( -8.0, 55.0, 1)
    call fig1%line( -7.0, 55.5, 1)
    call fig1%line( -8.5, 56.5, 1)
    call fig1%line(-24.0, 55.0, 1)    

    call fig1%line( -8.0, 54.5, 0)
    call fig1%line(-12.0, 52.5, 1)    

    call fig1%line(-23.0, 56.0, 0)
    call fig1%line(-21.5, 57.0, 1)
    call fig1%line(-10.0, 58.0, 1)
    call fig1%line( -9.0, 57.0, 1)    

    ! left eyebrow
    call fig1%line(-27.5, 56.5, 0)
    call fig1%line(-24.0, 59.0, 1)
    call fig1%line(-11.0, 59.5, 1)
    call fig1%line( -7.0, 61.0, 1)
    call fig1%line( -4.0, 65.0, 1)
    call fig1%line( -9.0, 63.0, 1)
    call fig1%line(-25.0, 62.0, 1)
    call fig1%line(-27.5, 56.5, 1)    

    ! right eyebrow
    call fig1%line( 27.5, 56.5, 0)
    call fig1%line( 24.0, 59.0, 1)
    call fig1%line( 11.0, 59.5, 1)
    call fig1%line(  7.0, 61.0, 1)
    call fig1%line(  4.0, 65.0, 1)
    call fig1%line(  9.0, 63.0, 1)
    call fig1%line( 25.0, 62.0, 1)
    call fig1%line( 27.5, 56.5, 1)    

    ! right eye
    call fig1%line( 19.0, 53.0, 0)
    call fig1%line( 23.0, 55.0, 1)
    call fig1%line( 16.0, 55.0, 1)
    call fig1%line(  9.0, 56.0, 1)
    call fig1%line(  9.5, 55.0, 1)
    call fig1%line( 19.0, 53.0, 1)    

    call fig1%line(  9.0, 58.0, 0)
    call fig1%line( 12.0, 58.0, 1)
    call fig1%line( 19.0, 56.0, 1)
    call fig1%line( 21.5, 56.0, 1)    

    call fig1%line(  0.0, 29.0, 0)
    call fig1%line(  5.0, 29.0, 1)
    call fig1%line( 11.0, 27.0, 1)
    call fig1%line(  6.0, 32.0, 1)
    call fig1%line(  0.0, 30.0, 1)
    call fig1%line( -6.0, 32.0, 1)
    call fig1%line(-11.0, 27.0, 1)
    call fig1%line( -5.0, 29.0, 1)
    call fig1%line(  0.0, 29.0, 1)    

    call fig1%line(-6.5, 21.5, 0)
    call fig1%line(-3.5, 20.0, 1)
    call fig1%line( 3.5, 20.0, 1)
    call fig1%line( 6.5, 21.5, 1)
    call fig1%line(-5.0, 21.0, 0)
    call fig1%line( 5.0, 21.0, 1)
    call fig1%show()
  end program uniplot_main

Modern Fortran Explained: Incorporating Fortran 2018 (Numerical Mathematics and Scientific Computation)

作者:Metcalf, Michael,Reid, John,Cohen, Malcolm
発売日: 2018/11/06
メディア: ハードカバー

Modern Fortran: Building efficient parallel applications

作者:Curcic, Milan
発売日: 2020/11/24
メディア: ペーパーバック

2021-05-23

【メモ帳】intel fortran 次世代構想 roadmap

intel compiler 講演

intel の Fortran/C++ compiler は backend を llvm ベースに変えつつありますが、その移行期についての解説がなされています。

旧来の intel backend 版は intel fortran classic (ifort)となり、次世代 llvm backend 版が intel fortran (ifx) の名を襲うようです。

techdecoded.intel.io

Fortran compiler は、C++ に比べて若干遅い進行で 2022 年初頭に β 版から脱却する見込みのようです。llvm 版が intel backend 版に匹敵するようになるのは、さらに遅れて 2023 年を見込んでいるようです。

f:id:fortran66:20210523020145p:plain

文法の対応状況としては、現在 F2003 におおむね対応した状況のようです。使い道が限定される割に文法的に直交性の悪い parameterized derived type と、使い道があるけど色々めんどくさい derived type I/O が特に重大な未対応要素のようです。詳しくは下記のページにあります。

software.intel.com

F2008/18 の対応状況はまだこれからのようです。

visual studio からの ifx の利用

上記の講演では、windows visual studio からの ifx の利用方法も紹介されていました。こんなところにメニューが隠れていようとは･･･

追記：言い忘れましたが、target が x64 版でないと出ません。

f:id:fortran66:20210523021747p:plain

intel fortran classic の f2018 対応は済んだことになっていますが、文法エラーは出なくても中身は未実装とか、結構ボロボロです。

Modern Fortran Explained: Incorporating Fortran 2018 (Numerical Mathematics and Scientific Computation)

作者:Metcalf, Michael,Reid, John,Cohen, Malcolm
Oxford University Press, USA

Amazon

Modern Fortran: Building efficient parallel applications

作者:Curcic, Milan
Manning Publications

Amazon

How to learn Applied Mathematics through modern FORTRAN

作者:Hernandez Ramos, Juan Antonio,Escoto Lopez, Javier
Independently published

Amazon

訃報　富永一朗先生死去

ご冥福をお祈り申し上げます。

www.sankei.com

hochi.news

2021-05-22

【メモ帳】分割数計算　改訂

分割数

Hardy-Ramanujan-Rademacher の公式の exp の複素数和の所を、どうせ虚部は x 軸に対して対称の位置に出て打ち消しあうので、cos に置き換えて実数のみで計算します。

Fortran と Julia の以前のプログラムを改訂します。

Hardy-Ramanujan-Rademacher の公式

$\displaystyle{ p(n)=\sum_{k=1}^\infty{\sqrt{k}\over\pi\sqrt{2}}\sum_{0\lt h\leq k\\ (h,k)=1}e^{{-2\pi i n h\over k}+\pi i\sum_{j=1}^{k-1}{j\over k}({h j\over k}-\lfloor{h j\over k}\rfloor-{1\over2})} \left. {d\over dz}{\sinh\left({\pi\over k}\sqrt{{2\over3}(z-{1\over24})}\right)\over \sqrt{z-{1\over24}} }\right|_{z=n} }$

無限和の内必要な項の見積もりの式があったのですが、失念しました。仕方ないので、適当に大きめに入れています。

Introduction to Analytic Number Theory (Undergraduate Texts in Mathematics)

作者:Apostol, Tom M.
発売日: 2010/12/04
メディア: ペーパーバック

Modular Functions and Dirichlet Series in Number Theory (Graduate Texts in Mathematics) (Graduate Texts in Mathematics, 41)

作者:Apostol, Tom M.
発売日: 2012/12/21
メディア: ペーパーバック

Integer Partitions

作者:Andrews, George
発売日: 2010/07/12
メディア: ペーパーバック

The Theory of Partitions (Encyclopedia of Mathematics and its Applications, Series Number 2)

作者:Andrews, George E.
発売日: 2008/01/12
メディア: ペーパーバック

整数の分割

作者:ジョージ・アンドリュース,キムモ・エリクソン
発売日: 2006/05/01
メディア: 単行本

Fortran 版

P(721) を計算し、無限和の各項の大きさの変化を見ます。

fortran66.hatenablog.com

    module partition_m
        use, intrinsic :: iso_fortran_env
        implicit none
        integer , parameter :: kd = real128
        real(kd), parameter :: pi = 4 * atan(1.0_kd)
    contains  

        impure elemental real(kd) function partition(n)
            integer, intent(in) :: n
            integer, parameter :: kmax = 21
            integer :: k
            partition = 0.0_kd
            do k = 1, kmax
                partition = partition + sqrt(k / 2.0_kd) / pi * fa(n, k) * fs(k, n) 
                print '(i6, 2f40.5)', k, sqrt(k / 2.0_kd) / pi * fa(n, k) * fs(k, n) 
            end do         
        end function partition
    
        real(kd) pure function fs(k, n) ! d/dx sinh( f * xs ) / xs ; xs = sqrt(x - 1/24)
            integer, intent(in) :: k, n 
            real(kd) :: xs, f
            xs = sqrt( n - 1 / 24.0_kd )
            f  = pi / k * sqrt(2 / 3.0_kd)
            fs = ( f * cosh(f * xs) - sinh(f * xs) / xs ) / (2 * xs * xs) 
        end function fs
        
        pure real(kd) function fa(n, k)
            integer, intent(in) :: n, k
            real(kd) :: s
            integer :: m
            fa = 0.0_kd
            do m = 1, k
                if (igcd(m, k) == 1) then
                    s = -real(2 * n * m, kd) / k + fb(m, k)
                    fa = fa + cos(pi * s) 
                end if
            end do
        end function fa
        
        pure real(kd) function fb(h, k)
            integer, intent(in) :: h, k
            integer :: j
            fb = 0.0_kd
            do j = 1, k - 1
                fb = fb + j * (h * j / real(k, kd) - floor(h * j / real(k, kd)) - 0.5_kd) 
            end do
            fb = fb / k
        end function fb
  
        pure integer function igcd(ia, ib)
            integer, value :: ia, ib
            do 
                igcd = ib
                ib = mod(ia, ib)
                if (ib == 0) exit
                ia = igcd
            end do
        end function igcd  
    end module partition_m
  
    program Pn
        use partition_m
        implicit none
        integer :: i
        real(kd) :: r       
!        print *, nint(partition([1:10]))
        
        do i = 721, 721 
            r = partition(i)
            print *
            print '(a, i3, a, f40.5)', 'p(', i, ')', r
        end do  
    end program Pn

     1       161061755750279601828302117.84821
     2                     -124192062781.96844
     3                           -706763.61926
     4                              2169.16830
     5                                 0.00000
     6                                14.20704
     7                                 6.07837
     8                                 0.18926
     9                                 0.04914
    10                                 0.00000
    11                                 0.08814
    12                                -0.03525
    13                                 0.03248
    14                                -0.00687
    15                                 0.00000
    16                                -0.01133
    17                                 0.00000
    18                                -0.00554
    19                                -0.00860
    20                                -0.00000
    21                                -0.00526
 
p(721)       161061755750279477635534762.00040

Julia 版

百までの分割数を計算します。

fortran66.hatenablog.com

function fs(k, z)
    xs = sqrt(z-1/24)
    fs = √(2/3) * π / k  
    (fs * cosh(fs * xs) - sinh(fs * xs) / xs) / 2xs^2
end

function fb(k, h)
    s = 0
    for j = 1:k-1
        s += j * (h * j / k - floor(h * j / k) - 1/2)
    end 
    s / k
end

function fa(n, k)
    s = 0
    for h = 1:k 
        if (gcd(h,k)==1) 
           z = -2n*h/k + fb(k, h)
           s += cos(π*z)
        end    
    end 
    s
end

function p(n)
    s = 0
    for k=1:10
        s += sqrt(k) * fa(n, k) * fs(k, n)
    end 
    s / √2pi  
end

using Printf
for i = 1:100
    pn = p(i)
    @printf( "%10d %20d %20.7f \n", i, round(Int64, real(pn)), real(pn)) 
end

         1                    1            1.0063342 
         2                    2            2.0044611 
         3                    3            3.0023039 
         4                    5            4.9960379 
         5                    7            6.9972722 
         6                   11           10.9886684 
         7                   15           14.9997219 
         8                   22           21.9944150 
         9                   30           30.0038956 
        10                   42           42.0078618 
        11                   56           56.0028080 
        12                   77           76.9987022 
        13                  101          101.0035727 
        14                  135          135.0035003 
        15                  176          175.9936875 
        16                  231          231.0001458 
        17                  297          296.9912596 
        18                  385          385.0031500 
        19                  490          489.9937920 
        20                  627          627.0072493 
        21                  792          792.0080395 
        22                 1002         1001.9931672 
        23                 1255         1254.9993226 
        24                 1575         1575.0069153 
        25                 1958         1958.0040403 
        26                 2436         2435.9892205 
        27                 3010         3010.0028149 
        28                 3718         3717.9938783 
        29                 4565         4564.9961219 
        30                 5604         5604.0040496 
        31                 6842         6842.0039095 
        32                 8349         8349.0100898 
        33                10143        10142.9938646 
        34                12310        12309.9901594 
        35                14883        14883.0072610 
        36                17977        17977.0032812 
        37                21637        21637.0036938 
        38                26015        26014.9961219 
        39                31185        31184.9842031 
        40                37338        37338.0070623 
        41                44583        44582.9966134 
        42                53174        53174.0087179 
        43                63261        63261.0087796 
        44                75175        75174.9936763 
        45                89134        89134.0016629 
        46               105558       105557.9980458 
        47               124754       124753.9977149 
        48               147273       147272.9904516 
        49               173525       173525.0070202 
        50               204226       204226.0036898 
        51               239943       239943.0021134 
        52               281589       281588.9910221 
        53               329931       329931.0021026 
        54               386155       386155.0101807 
        55               451276       451275.9887966 
        56               526823       526823.0079577 
        57               614154       614153.9990722 
        58               715220       715220.0009011 
        59               831820       831820.0001180 
        60               966467       966466.9958295 
        61              1121505      1121504.9930073 
        62              1300156      1300155.9942735 
        63              1505499      1505499.0088466 
        64              1741630      1741630.0099512 
        65              2012558      2012558.0033675 
        66              2323520      2323519.9945812 
        67              2679689      2679689.0022615 
        68              3087735      3087734.9963524 
        69              3554345      3554344.9974198 
        70              4087968      4087967.9953353 
        71              4697205      4697204.9976341 
        72              5392783      5392782.9993796 
        73              6185689      6185689.0052464 
        74              7089500      7089500.0034465 
        75              8118264      8118264.0027078 
        76              9289091      9289091.0111846 
        77             10619863     10619862.9921330 
        78             12132164     12132163.9936887 
        79             13848650     13848649.9941123 
        80             15796476     15796476.0006777 
        81             18004327     18004326.9991289 
        82             20506255     20506255.0020917 
        83             23338469     23338469.0030397 
        84             26543660     26543659.9989915 
        85             30167357     30167357.0013871 
        86             34262962     34262962.0040397 
        87             38887673     38887673.0087649 
        88             44108109     44108108.9878901 
        89             49995925     49995924.9983584 
        90             56634173     56634173.0039072 
        91             64112359     64112358.9928125 
        92             72533807     72533806.9986574 
        93             82010177     82010177.0021016 
        94             92669720     92669719.9965452 
        95            104651419    104651419.0057203 
        96            118114304    118114304.0000427 
        97            133230930    133230930.0105914 
        98            150198136    150198136.0060634 
        99            169229875    169229874.9836169 
       100            190569292    190569292.0010892

2021-05-21

【メモ帳】補遺

ZDNet 記事の翻訳

先日紹介した記事の日本語訳が出て、少し話題になっておりました。 fortran66.hatenablog.com

japan.zdnet.com

Pyodide

以前、SciPy など Fortran 部分が WebAssembly 化できなくて困っていたのが一応解決したようです。

CBLAS/CLAPACK なども実は f2c で変換しているだけで、文字列の扱いや alternate return の扱いなど、色々細々した所で問題があるようです。

なお、Mozilla は金欠で Iodide などのプロジェクトは放棄されたようです。

fortran66.hatenablog.com

以下のページから実際に試して見られます。

pyodide.org

Guide to Fortran 2008 Programming

作者:Brainerd, Walter S. S.
発売日: 2016/10/22
メディア: ペーパーバック

Modern Fortran Explained: Incorporating Fortran 2018 (Numerical Mathematics and Scientific Computation)

作者:Metcalf, Michael,Reid, John,Cohen, Malcolm
発売日: 2018/11/06
メディア: ハードカバー

Modern Fortran: Style and Usage

作者:Clerman, Norman S
発売日: 2012/02/23
メディア: ペーパーバック

Fortran ハンドブック

作者:田口俊弘
発売日: 2015/07/22
メディア: Kindle版

2021-05-14

【寝言】夕空に月水金！

三日月水星金星

夕空に三日月が出ているのを眺めていると、地平線近くに金星も見えていました。水星の方が高度は高かったですが、金星の方がはるかに見やすかったです。光度は水星が 0.1 等、金星が -3.9 等程のようです。

f:id:fortran66:20210514203753p:plain

May, West

f:id:fortran66:20210514203819j:plain

earthsky.org

Mae West

メイ・ウエスト

www.youtube.com

Go West Young Man by Mae West

メディア: DVD

2021-05-13

【メモ帳】MATLAB の Moler の面白講演

Star Trek 映画に登場！

LINPACK, EISPACK や MATLAB で有名な Moler の講演が結構面白いです。

www.youtube.com

（以前 EISPACK の FORTRAN への移植についての Moler の講演を紹介しました。） fortran66.hatenablog.com

MATLAB は元々は Moler によって FORTRAN66 で書かれていました。その時 N. Wirth の Algorithms + Data Structures = Programs 中の PL/0 コンパイラ実装の章を読んでインタープリンタを作ったそうです。なお MATLAB のインデックスが 1 始まりになっているのは、Wirth のこの本の PASCAL の影響であり、数学では 1 始まりだからだと言っています。

Algorithms + Data Structures = Programs (Prentice-Hall Series in Automatic Computation)

作者:Wirth, Niklaus
発売日: 1976/02/01
メディア: ハードカバー

Moler がスタンフォード大学で MATLAB を使った集中講義をやった時、数学やコンピュータ学科の学生は興味を見せなかったけれど、行列計算を実用していた工学部の学生が食いついて来たそうです。そのとき自分で FORTRAN のプログラムを組むのを止めて MATLAB を使いだした学生の中に、後に MATLAB を商用にして売ることを考えたジャック・リトルがいて、当時キラキラ言語だった C に移植して商売をはじめ、今に続くようです。

Moler が当時最先端のコンピュータ・グラフィックスを用いて SVD 分解の啓蒙映画を Los Alamos 研究所で作ったところ、宇宙大作戦 (Star Trek) の初代の劇場版映画に引用されたそうです。ミスター・スポックの背後に行列対角化中の映像がたしかに映っています。

www.youtube.com

f:id:fortran66:20210513204425p:plain

f:id:fortran66:20210513204435p:plain

【追記】１０年前の記事にもあったｗ eetimes.jp

なお初代『スター・トレック』劇場版は、本邦においては尼さんマニアにとても高評価だった映画です。

www.agonybooth.com

スター・トレックⅠ(字幕版)

発売日: 2013/11/26
メディア: Prime Video

スター・トレックⅠ (吹替版)

発売日: 2013/11/26
メディア: Prime Video

FORTRAN 版 MATLAB

そんな FORTRAN 版の MATLAB が The Cyber Vanguard: Lazy Reading で紹介されていました。

http://cyber.dabamos.de/blog/#article-2021-05-07 から

www.urbanjost.altervista.org

元々 FORTRAN66 で書かれていたようですが、FORTRAN77 で動くようになっているようで、現今の Fortran でもコンパイルできます。

パーサー・インタープリタの所は、プログラムで再帰を実現しているようです。

無利用の捨て引数として配列にスカラーを渡したりしているので (FORTRAN は call by reference なので番地が渡りさえすればいい)、argument check を切らないと動きません。

HELP ファイルを実行ファイルと同じディレクトリに入れておけば HELP も使えます。

FORTRAN 版 MATLAB のマニュアルの例題が面白くて、過去の米国国勢調査の結果から 1980 年の米人口を多項式近似（Lagrange近似）で見積もっているのですが、Vandermonde 行列を使って外挿して解いて、指数関数などと比較してよく合うとみなしています。

ひとしきり計算した後に、多項式を方程式として解いてみせ 1984 年に解があるのでこの年にアメリカの人口が 0 になると結論してみせていますｗ多項式近似はしばしば隅っこの方で飛び跳ねるのを使った面白ネタです。内挿にしておけとｗ

ifort

ifort matlab88.f90 -o matlab -nocheck

gfortran-10

gfortran matlab88.f90 -o matlab -std=legacy -fallow-invalid-boz

適宜 -O2 など最適化オプションを与えてください。

f:id:fortran66:20210513214637p:plain

Moler の別動画に見易いスクリーンショットがあります。 f:id:fortran66:20210513215012p:plain

www.youtube.com

ほとんど同じ持ちネタを話していますが、Star Trek は出てきません。

質疑の所で Python/SciPy などとの競合について問われています。Moler 答えて曰く MATLAB の値段が高すぎると言われるが、ではいくらの値段か？と問い返すと皆答えられない、結局ソフトウェアには金を払えないという思想の問題だと言い切っています。

またサポートの無いコロコロ変わるつぎはぎソフトウェアより安定が求められているとも言っています。トヨタが MathWorks の最大の顧客だとも言っていて、興味深いところです。

Numerical Computing with MATLAB, Revised Reprint by Cleve B. Moler(2008-07-25)

作者:Cleve B. Moler
メディア: 文庫

Computer Solution of Linear Algebraic Systems (Automatic Computation S.)

作者:Forsythe, George E.,Moler, C.B.
発売日: 1967/09/01
メディア: テキスト

LINPACK Users' Guide

作者:Dongarra, J. J.,Bunch, J. R.,Moler, G. B.,Stewart, G. W.
発売日: 1987/01/01
メディア: ペーパーバック

Matrix Eigensystem Routines - EISPACK Guide Extension (Lecture Notes in Computer Science (51))

作者:Dongarra, J. J.,Moler, C. B.,Garbow, B. S.,Boyle, J. M.
発売日: 2008/06/13
メディア: ペーパーバック

Matrix Eigensystem Routines - EISPACK Guide (Lecture Notes in Computer Science, 6)

作者:Klema, V. C.,Dongarra, J. J.,Smith, B. T.,Moler, C. B.,Ikebe, Y.
発売日: 2008/06/13
メディア: ペーパーバック

LAPACK利用の手引―行列計算パッケージ

発売日: 1995/07/01
メディア: 単行本

LAPACK Users' Guide (Software, Environments and Tools) by E. Anderson Z. Bai C. Bischof S. Blackford J. Demmel J. Dongarra J. Du Croz A. Greenbaum S. Hammarling A. McKenney D. Sorensen(1987-01-01)

作者:E. Anderson Z. Bai C. Bischof S. Blackford J. Demmel J. Dongarra J. Du Croz A. Greenbaum S. Hammarling A. McKenney D. Sorensen
メディア: 文庫

2021-05-12

【寝言】水星が東方最大離角近し

夕空に水星

19時も大分過ぎると空も暗くなってきて肉眼でもよく見えました。

www.nao.ac.jp

スマホのデジタルズームなのでガバガバです。 f:id:fortran66:20210512004515j:plain

美少女戦士セーラームーンコスプレ水野亜美風男性用 Lサイズ

メディア:

男性用 L サイズ・・・

MR. バッド・ガイ(スペシャル・エディション)(SHM-CD)

アーティスト:フレディ・マーキュリー
発売日: 2019/10/11
メディア: CD

fortran66のブログ

fortran について書きます。

【メモ帳】unicode 点字で阿部さん

unicode 表示

unicode 点字

Fortran と unicode

Apple M1 Mac での出力

Windows WSL での出力

ソース・コード

【メモ帳】intel fortran 次世代構想 roadmap

intel compiler 講演

visual studio からの ifx の利用

訃報　富永一朗先生死去

【メモ帳】分割数計算　改訂

分割数

Hardy-Ramanujan-Rademacher の公式

Fortran 版

Julia 版

【メモ帳】補遺

ZDNet 記事の翻訳

Pyodide

【寝言】夕空に月水金！

三日月水星金星

May, West

Mae West

【メモ帳】MATLAB の Moler の面白講演

Star Trek 映画に登場！

FORTRAN 版 MATLAB

ifort

gfortran-10

【寝言】水星が東方最大離角近し

夕空に水星

unicode 表示

Fortran と unicode

Apple M1 Mac での出力

Windows WSL での出力

ソース・コード

intel compiler 講演

visual studio からの ifx の利用

訃報 富永一朗先生死去

分割数

Hardy-Ramanujan-Rademacher の公式

Fortran 版

Julia 版

ZDNet 記事の翻訳

Pyodide

三日月 水星 金星

May, West

Mae West

Star Trek 映画に登場！

FORTRAN 版 MATLAB

ifort

gfortran-10

夕空に水星

unicode 表示　

訃報　富永一朗先生死去

三日月水星金星