整数論的関数

メモ帳代わりに、Eulerの関数と約数関数を計算するプログラムを載せておきます。

は、n以下でnと疎な整数の総数、は、で定義され、すなわちそれぞれnの約数の総数、総和に当たります。

は、つまり n,m が疎な時、積構造が保存してが成り立ちます。それを計算で確かめてみました。

実行結果

Prime numbers < 100
  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

  n : phi sig0 sig1
   1:   1   1   1
   2:   1   2   3
   3:   2   2   4
   4:   2   3   7
   5:   4   2   6
   6:   2   4  12
   7:   6   2   8
   8:   4   4  15
   9:   6   3  13
  10:   4   4  18
  11:  10   2  12
  12:   4   6  28
  13:  12   2  14
  14:   6   4  24
  15:   8   4  24
  16:   8   5  31
  17:  16   2  18
  18:   6   6  39
  19:  18   2  20
  20:   8   6  42
  21:  12   4  32
  22:  10   4  36
  23:  22   2  24
  24:   8   8  60
  25:  20   3  31
  26:  12   4  42
  27:  18   4  40
  28:  12   6  56
  29:  28   2  30
  30:   8   8  72
  31:  30   2  32
  32:  16   6  63
  33:  20   4  48
  34:  16   4  54
  35:  24   4  48
  36:  12   9  91
  37:  36   2  38
  38:  18   4  60
  39:  24   4  56
  40:  16   8  90
  41:  40   2  42
  42:  12   8  96
  43:  42   2  44
  44:  20   6  84
  45:  24   6  78
  46:  22   4  72
  47:  46   2  48
  48:  16  10 124
  49:  42   3  57
  50:  20   6  93
  51:  32   4  72
  52:  24   6  98
  53:  52   2  54
  54:  18   8 120
  55:  40   4  72
  56:  24   8 120
  57:  36   4  80
  58:  28   4  90
  59:  58   2  60
  60:  16  12 168
  61:  60   2  62
  62:  30   4  96
  63:  36   6 104
  64:  32   7 127
  65:  48   4  84
  66:  20   8 144
  67:  66   2  68
  68:  32   6 126
  69:  44   4  96
  70:  24   8 144
  71:  70   2  72
  72:  24  12 195
  73:  72   2  74
  74:  36   4 114
  75:  40   6 124
  76:  36   6 140
  77:  60   4  96
  78:  24   8 168
  79:  78   2  80
  80:  32  10 186
  81:  54   5 121
  82:  40   4 126
  83:  82   2  84
  84:  24  12 224
  85:  64   4 108
  86:  42   4 132
  87:  56   4 120
  88:  40   8 180
  89:  88   2  90
  90:  24  12 234
  91:  72   4 112
  92:  44   6 168
  93:  60   4 128
  94:  46   4 144
  95:  72   4 120
  96:  32  12 252
  97:  96   2  98
  98:  42   6 171
  99:  60   6 156
 100:  40   9 217

Euler phi function
gcd( 1, 1) = 1: phi( 1) * phi( 1):   1 *   1==   1: phi( 1)
gcd( 1, 2) = 1: phi( 1) * phi( 2):   1 *   1==   1: phi( 2)
gcd( 1, 3) = 1: phi( 1) * phi( 3):   1 *   2==   2: phi( 3)
gcd( 1, 4) = 1: phi( 1) * phi( 4):   1 *   2==   2: phi( 4)
gcd( 1, 5) = 1: phi( 1) * phi( 5):   1 *   4==   4: phi( 5)
gcd( 1, 6) = 1: phi( 1) * phi( 6):   1 *   2==   2: phi( 6)
gcd( 1, 7) = 1: phi( 1) * phi( 7):   1 *   6==   6: phi( 7)
gcd( 1, 8) = 1: phi( 1) * phi( 8):   1 *   4==   4: phi( 8)
gcd( 1, 9) = 1: phi( 1) * phi( 9):   1 *   6==   6: phi( 9)
gcd( 1,10) = 1: phi( 1) * phi(10):   1 *   4==   4: phi(10)
gcd( 2, 3) = 1: phi( 2) * phi( 3):   1 *   2==   2: phi( 6)
gcd( 2, 5) = 1: phi( 2) * phi( 5):   1 *   4==   4: phi(10)
gcd( 2, 7) = 1: phi( 2) * phi( 7):   1 *   6==   6: phi(14)
gcd( 2, 9) = 1: phi( 2) * phi( 9):   1 *   6==   6: phi(18)
gcd( 3, 4) = 1: phi( 3) * phi( 4):   2 *   2==   4: phi(12)
gcd( 3, 5) = 1: phi( 3) * phi( 5):   2 *   4==   8: phi(15)
gcd( 3, 7) = 1: phi( 3) * phi( 7):   2 *   6==  12: phi(21)
gcd( 3, 8) = 1: phi( 3) * phi( 8):   2 *   4==   8: phi(24)
gcd( 3,10) = 1: phi( 3) * phi(10):   2 *   4==   8: phi(30)
gcd( 4, 5) = 1: phi( 4) * phi( 5):   2 *   4==   8: phi(20)
gcd( 4, 7) = 1: phi( 4) * phi( 7):   2 *   6==  12: phi(28)
gcd( 4, 9) = 1: phi( 4) * phi( 9):   2 *   6==  12: phi(36)
gcd( 5, 6) = 1: phi( 5) * phi( 6):   4 *   2==   8: phi(30)
gcd( 5, 7) = 1: phi( 5) * phi( 7):   4 *   6==  24: phi(35)
gcd( 5, 8) = 1: phi( 5) * phi( 8):   4 *   4==  16: phi(40)
gcd( 5, 9) = 1: phi( 5) * phi( 9):   4 *   6==  24: phi(45)
gcd( 6, 7) = 1: phi( 6) * phi( 7):   2 *   6==  12: phi(42)
gcd( 7, 8) = 1: phi( 7) * phi( 8):   6 *   4==  24: phi(56)
gcd( 7, 9) = 1: phi( 7) * phi( 9):   6 *   6==  36: phi(63)
gcd( 7,10) = 1: phi( 7) * phi(10):   6 *   4==  24: phi(70)
gcd( 8, 9) = 1: phi( 8) * phi( 9):   4 *   6==  24: phi(72)
gcd( 9,10) = 1: phi( 9) * phi(10):   6 *   4==  24: phi(90)
続行するには何かキーを押してください . . .

ソース・プログラム

効率を考えず、Fortran の命令を使うことを考えました。(20150101 補足修正 sigmak)

    program main
      implicit none
      integer :: i, j
      
      print '(a)', 'Prime numbers < 100'
      print '(25i3)', eratos(100) 
      
      print '(/a)', '  n : phi sig0 sig1'
      do i = 1, 100
        print '(i4, a, 3i4)', i, ':', phi(i), sigma0(i), sigma1(i)    
      end do    
      
      print '(/a)', 'Euler phi function'
      do i = 1, 10
        do j = i, 10  
          if (phi(i) * phi(j) == phi(i * j)) & 
              print '(5(a, i2), 3(a, i4), a, i2, a)', 'gcd(', i, ',', j, ') =', gcd(i, j), &
                         ': phi(', i, ') * phi(', j, '):', phi(i), ' *', phi(j), '==', phi(i * j), ': phi(', i * j, ')' 
        end do   
      end do
       
    contains
      integer function gcd(n, m)
        integer, value :: n, m
        do  
          n = mod(n, m) 
          gcd = m
          if (n == 0) exit
          m = n
          n = gcd
        end do
      end function gcd
    
      function eratos(n) result(ires)
        integer, intent(in) :: n
        integer, allocatable :: ires(:)
        integer :: i, iwk(n)
        iwk(1) = 0
        forall(i = 2:n) iwk(i) = i
        forall(i = 2:int(sqrt(real(n)))) iwk(i*i::i) = 0    
        ires = pack(iwk, iwk /= 0)
      end function eratos  
        
      function ifactor(n) result(ires)
        integer, intent(in) ::n
        integer :: i
        integer, allocatable :: ires(:), iwk(:), prim(:)
        prim = eratos(n)
        allocate(ires(size(prim)), source = 0)
        allocate(iwk (size(prim)), source = n)
        ires = 0
        do while( any(iwk /= 1) )
          where(mod(iwk, prim) == 0)
            ires = ires + 1
            iwk = iwk / prim
          elsewhere
            iwk = 1
          end where  
        end do
      end function ifactor

      integer function phi(n)
        integer, intent(in) :: n
        integer, allocatable :: ifact(:), iprim(:), ifc(:), ipr(:)       
        iprim = eratos(n)
        ifact = ifactor(n)
        ipr = pack(iprim, ifact /= 0)
        ifc = pack(ifact, ifact /= 0) - 1
        phi = product(ipr**ifc * (ipr - 1))
      end function phi
     
      integer function sigma0(n)
        integer, intent(in) ::n
        integer :: i, iwk(n)
        forall(i = 1:n) iwk(i) = i ! iwk = [1:n]
        sigma0 = count(mod(n, iwk) == 0)
      end function sigma0  
        
      integer function sigma1(n)
        integer, intent(in) ::n
        integer :: i, iwk(n)
        forall(i = 1:n) iwk(i) = i ! iwk = [1:n]
        sigma1 = sum(iwk, mod(n, iwk) == 0)
      end function sigma1  

      integer function sigmak(n, k)
        integer, intent(in) ::n, k
        integer :: i, iwk(n)
        forall(i = 1:n) iwk(i) = i ! iwk = [1:n]
        sigmak = sum(iwk**k, mod(n, iwk) == 0)
      end function sigmak  
    end program main