整数論的関数２ 位数 原始根 指数表

メモ帳代わりに。ソースだけ。

    program order
      implicit none
      integer, allocatable :: iprm(:), irtp(:), iwk(:)
      integer :: i, k, m
      iprm = eratos(100)
      do k = 1, size(iprm)
        print *, '   ', iprm(k)  
        do i = 2, iprm(k) - 1
          print *, i, ':', ord(iprm(k), i)   
        end do    
        print *
      end do    
      
      allocate(irtp(size(iprm)))
      do k = 1, size(iprm)
        do i = 2, iprm(k) - 1
          if ( ord(iprm(k), i) == iprm(k) - 1 .and. any(i == iprm) ) exit    
        end do    
        irtp(k) = i
      end do    

      do k = 2, size(iprm)
        print *, '   ', iprm(k)  
        print *, '   ', irtp(k)
        allocate(iwk(iprm(k) - 1))
        m = 1
        iwk(1) = 0
        do i = 1, iprm(k) - 2
          m = mod(m * irtp(k), iprm(k))  
          iwk(m) = i 
        end do    
       ! print '(2i3, a, *(i3))', iprm(k), irtp(k), ':', iwk
        do i = 1, iprm(k) - 1
         print *, i, ':', iwk(i)   
        end do    
        print *
        deallocate(iwk)
      end do    
      
    contains
       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  
        
      integer function ord(ip, m)
        integer, intent(in) :: ip, m
        integer :: k
        k = m
        do ord = 2, ip - 1
          k = mod(k * m, ip)
          if (k == 1) exit
        end do
      end function ord
    end program order