ウラムの螺旋

ウラムの螺旋とは、らせん状に整数を書いて、素数を選んでゆくと、奇妙に水平・垂直、斜め４５度方向に並んでいることが多いというものです。

31 * 31 = 961 個の整数のらせん状の配置

201*201 = 40401 までの素数

アスペクト比が狂っているので注意が必要です。

素数の生成式がよく、二次式で表わされることと関係があるようです。この図の意味は、素数が全くランダムに並んでいるのでなく、何らかの規則性をもって配列されていることを示唆していることのようです。

ソース・プログラム

素朴に整数を配列しました。
素数はエラトステネスのふるいで求めます。

    module m_prim
      implicit none
    contains
      function eratos(n) result(ires)
        integer, intent(in) :: n
        integer, allocatable :: ires(:)
        integer :: i, k
        allocate(ires(n))
        forall(i = 1:n) ires(i) = i
        ires(1) = 0
        do i = 2, int(sqrt(real(n)))
          if (ires(i) /= 0) ires(ires(i)**2::ires(i)) = 0
        end do  
      end function eratos  
      
      subroutine mk_map(imap, n)
        integer, intent(out) :: imap(-n:, -n:)
        integer, intent(in ) :: n
        integer :: i, k, m, ix, iy
        do i = 0, n
          k = 2 * i + 1
          m = k * k
          ix = i
          iy = -i
          do ix = i, -i, -1
            imap(ix, iy) = m
            m = m - 1
          end do 
          ix = -i
          do iy = -i + 1, i
            imap(ix, iy) = m
            m = m - 1
          end do 
          iy = i
          do ix = -i + 1, i 
            imap(ix, iy) = m
            m = m - 1
          end do
          ix = i
          do iy = i - 1, -i + 1, -1
            imap(ix, iy) = m
            m = m - 1
          end do  
        end do  
      end subroutine mk_map
    end module m_prim
  
    program ulam
      use m_prim
      implicit none
      integer, parameter :: n = 15
      integer :: imap(-n:n, -n:n)
      character :: cmap(-n:n, -n:n)
      integer :: ix, iy
      integer, allocatable :: ip(:)
      ip = eratos( (2 * n + 1)**2 )
      ip(1) = 1
      !
      call mk_map(imap, n)
      !
      cmap = ' '
      forall (ix = -n:n, iy = -n:n, ip(imap(ix, iy)) /= 0) cmap(ix, iy) = '*'
      do iy = n, -n, -1
      !  print '(201a1)', cmap(:, iy)
        print '(201i4)', imap(:, iy)
      end do  
    end program ulam