[R]Gaussian Primes in R

wikipedia のユークリッドの互除法のページでは  Gauss 素数の可視化(Mathmatica)がされていた。
ユークリッドの互除法のステップ数につつき、Gauss 素数の可視化も R&ggplot2 でも実装してみた。

まずはガウス整数の定義から。

# Gaussian Integer

Gaussian Integer は
Z[i] = {a + bi | a, b in Z}
の形の整数

ノルムは次のように定義されている。
N(a + bi) = a^2 + b^2

# Gaussian Prime
Gaussian Integer a + bi が Gaussian Prime となる必要十分条件は

  • one of a, b is zero and the other is a prime number of the form 4n + 3 (with n a nonnegative integer) or its negative − (4n + 3), or
  • both are nonzero and a2 + b2 is a prime number (which will not be of the form 4n + 3).

※以上の定義は wikipedia から

—–

# Implementation in R

R で定義通りに実装すると、次のようになる

library(ggplot2)
library(gmp)

NORM <- 500
N <- ceiling(sqrt(NORM))
re_val <- numeric(0)
im_val <- numeric(0)
for (i in 0:N) {
   for (j in 0:N) {
     if (i**2 + j**2 > NORM) {
        # do nothing
    } else if (i == 0 & isprime(j) & (abs(j) %% 4 == 3)) {
        re_val<- append(re_val, i)
        im_val <- append(im_val, j)
    } else if (j == 0 & isprime(i) & (abs(i) %% 4 == 3)) {
        re_val<- append(re_val, i)
        im_val <- append(im_val, j)
    } else if (isprime(i ** 2 + j ** 2)){
        re_val<- append(re_val, i)
        im_val <- append(im_val, j)
    }
  }
}

df <- data.frame(re_val, im_val)
df <- rbind(df
           ,transform(df, re_val = -re_val)
           ,transform(df, im_val = -im_val)
           ,transform(df, re_val = -re_val, im_val = -im_val))

p <- ggplot(df, aes(re_val, im_val))
p + geom_tile()

# MEMO
素数の判定には isprime{gmp} を利用。
整数を与えると、次のような値を返す

0 Number is not prime
1 Number is probably prime
2 Number is prime

—–

ここから不要なX/Y軸名や背景色を取っ払うと次のようになる。

p <- ggplot(df, aes(re_val, im_val))
p <- p + geom_tile()
p + opts(title            = 'Gaussian Primes with norm less than 500',
         axis.title.x     = theme_blank(),
         axis.title.y     = theme_blank(),
         axis.text.x      = theme_blank(),
         axis.text.y      = theme_blank(),
         axis.ticks       = theme_blank(),
         panel.background = theme_blank()

geom_tile を見つけるまでは geom_point を利用していた。元画像との違いを見比べると、自分は点でプロットしており、wikipediaでは正方形でプロットしていた。
可視化対象が整数という離散的なデータのため、正方形の方がきれいに可視化できる。
geom_tile ではなく geom_point を使うと、次のようになる。

p <- ggplot(df, aes(re_val, im_val))
p + geom_point(aes(shape=15))
Tagged with: ,
Posted in algorithm, R

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Archives
  • RT @__apf__: How to write a research paper: a guide for software engineers & practitioners. docs.google.com/presentation/d… /cc @inwyrd 6 days ago
  • RT @HayatoChiba: 昔、自然と対話しながら数学に打ち込んだら何かを悟れるのではと思いたち、専門書1つだけ持ってパワースポットで名高い奈良の山奥に1週間籠ったことがある。しかし泊まった民宿にドカベンが全巻揃っていたため、水島新司と対話しただけで1週間過ぎた。 それ… 3 weeks ago
  • RT @googlecloud: Ever wonder what underwater fiber optic internet cables look like? Look no further than this deep dive w/ @NatAndLo: https… 3 weeks ago
  • @ijin UTC+01:00 な時間帯で生活しています、、、 6 months ago
  • RT @mattcutts: Google's world-class Site Reliability Engineering team wrote a new book: amazon.com/Site-Reliabili… It's about managing produc… 9 months ago
%d bloggers like this: