技術メモ

役に立てる技術的な何か、時々自分用の覚書。幅広く色々なことに興味があります。

円周率素数 ネイピア数素数

3.14159265
と言われて何を思い浮かべるだろうか。
そう、円周率である。

円周率と言えば、3.1415...と続く数なのは小学生でも周知の事実だが、
円周率を頭から数えたとき出てくる素数というのはあまり知られていないかもしれない。

まず
3
素数
31
素数
314
素数ではない
31415
素数ではない
314159
素数

こう見てみると、素数って意外と多いんじゃないかと思うかもしれない。
しかし次に素数が現れるのはなんと38桁目の


31415926535897932384626433832795028841
となる。


いやちょっと待てよ、円周率もやったならネイピア数もやってあげないと可哀想じゃないか
と思ったので、ネイピア数(e=2.71828...)でも同じことをやってみた。

まず
2
素数
27
素数ではない
271
素数
2718
素数ではない


...と続けていけば


2718281

2718281828459045235360287471352662497757247093699959574966967627724076630353547594571

と続き、

1781桁目、2780桁目に素数が見つかることがわかる。*1*2


27182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274274663919320030599218174135966290435729003342952605956307381323286279434907632338298807531952510190115738341879307021540891499348841675092447614606680822648001684774118537423454424371075390777449920695517027618386062613313845830007520449338265602976067371132007093287091274437470472306969772093101416928368190255151086574637721112523897844250569536967707854499699679468644549059879316368892300987931277361782154249992295763514822082698951936680331825288693984964651058209392398294887933203625094431173012381970684161403970198376793206832823764648042953118023287825098194558153017567173613320698112509961818815930416903515988885193458072738667385894228792284998920868058257492796104841984443634632449684875602336248270419786232090021609902353043699418491463140934317381436405462531520961836908887070167683964243781405927145635490613031072085103837505101157477041718986106873969655212671546889570350354021234078498193343210681701210056278802351930332247450158539047304199577770935036604169973297250886876966403555707162268447162560798826517871341951246652010305921236677194325278675398558944896970964097545918569563802363701621120477427228364896134225164450781824423529486363721417402388934412479635743702637552944483379980161254922785092577825620926226483262779333865664816277251640191059004916449982893150566047258027786318641551956532442586982946959308019152987211725563475463964479101459040905862984967912874068705048958586717479854667757573205681288459205413340539220001137863009455606881667400169842055804033637953764520304024322566135278369511778838638744396625322498506549958862342818997077332761717839280349465014345588970719425863987727547109629537415211151368350627526023


27182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274274663919320030599218174135966290435729003342952605956307381323286279434907632338298807531952510190115738341879307021540891499348841675092447614606680822648001684774118537423454424371075390777449920695517027618386062613313845830007520449338265602976067371132007093287091274437470472306969772093101416928368190255151086574637721112523897844250569536967707854499699679468644549059879316368892300987931277361782154249992295763514822082698951936680331825288693984964651058209392398294887933203625094431173012381970684161403970198376793206832823764648042953118023287825098194558153017567173613320698112509961818815930416903515988885193458072738667385894228792284998920868058257492796104841984443634632449684875602336248270419786232090021609902353043699418491463140934317381436405462531520961836908887070167683964243781405927145635490613031072085103837505101157477041718986106873969655212671546889570350354021234078498193343210681701210056278802351930332247450158539047304199577770935036604169973297250886876966403555707162268447162560798826517871341951246652010305921236677194325278675398558944896970964097545918569563802363701621120477427228364896134225164450781824423529486363721417402388934412479635743702637552944483379980161254922785092577825620926226483262779333865664816277251640191059004916449982893150566047258027786318641551956532442586982946959308019152987211725563475463964479101459040905862984967912874068705048958586717479854667757573205681288459205413340539220001137863009455606881667400169842055804033637953764520304024322566135278369511778838638744396625322498506549958862342818997077332761717839280349465014345588970719425863987727547109629537415211151368350627526023264847287039207643100595841166120545297030236472549296669381151373227536450988890313602057248176585118063036442812314965507047510254465011727211555194866850800368532281831521960037356252794495158284188294787610852639813955990067376482922443752871846245780361929819713991475644882626039033814418232625150974827987779964373089970388867782271383605772978824125611907176639465070633045279546618550966661856647097113444740160704626215680717481877844371436988218559670959102596862002353718588748569652200050311734392073211390803293634479727355955277349071783793421637012050054513263835440001863239914907054797780566978533580489669062951194324730995876552368128590413832411607226029983305353708761389396391779574540161372236187893652605381558415871869255386061647798340254351284396129460352913325942794904337299085731580290958631382683291477116396337092400316894586360606458459251269946557248391865642097526850823075442545993769170419777800853627309417101634349076964237222943523661255725088147792231519747

素数って楽しいですね。

*1:確率的素数判定法を使っているので厳密には「素数である可能性が非常に高い」ということです

*2:合成数なのに素数と判定されてしまう確率はおよそ4^{-50}

Nが現れる素数(N=1,2,3,4)

2が現れる素数という面白い素数が紹介されていた。
2が現れる素数 - INTEGERS

昔せっかく高速素数判定器を作ったので、どうせならNが現れる素数を見つけてやろう!と思い立った。

プログラム

(※プログラムはpython(2.7.12)で動作します)

ルールとしては
①四隅のみの数字を変える(もちろん先頭は1以上の数字)
②四隅の数字はN以外の数字にする
としています。

なので、それぞれ5832(8*9*9*9)個の数字の中から素数を探すことになります。

高速素数判定のプログラム(再掲)

primechecker.pyという名前で保存

import random
import numpy as np

class PrimeChecker:
    def __init__(self, list_limit = pow(10,3)):
        if list_limit < 5: list_limit = 5
        self.prime_list = self.set_prime_list(list_limit)
        self.list_limit = list_limit

    @staticmethod
    def set_prime_list(n):
        def mark(primes,x):
            for i in xrange(2*x, n+1, x):
                primes[i-2] = False
        primes = [(i%2!=0 and i%3!=0 and i%5!=0) for i in range(2,n+1)]
        primes[0] = primes[1] = primes[3] = True
        for x in xrange(7, int(np.sqrt(n+1))+1):
           if primes[x-2]: mark(primes, x)
        prime_list = []
        for i in range(n-1):
            if primes[i]: prime_list.append(i+2)
        return prime_list

    def is_prime(self, n):
        k = 50 # number of iteration times
        if n < self.list_limit: return True if n in self.prime_list else False
        for prime in self.prime_list:
            if n%prime == 0: return False
        else:
            s, d = 0, n-1
            while d%2==0: s,d = s+1, d/2
            while k > 0:
                k = k-1
                a = random.randint(1,n-1)
                t, y = d, pow(a,d,n)
                while t != n-1 and y != 1 and y != n-1:
                    y = pow(y,2,n)
                    t <<= 1
                if y != n-1 and t%2 == 0:
                    return False
            return True

if __name__=="__main__":
    pc = PrimeChecker()
    pc.is_prime(58427) # True

Nが現れる素数を見つけるためのプログラム

from primechecker import PrimeChecker

def display_letter(p):
    """
    display letter on commandline for checking
    """
    print("-"*18)
    for i in range(0,216,18)[::-1]:
        print(str(p/pow(10,i) + pow(10,18))[1:])
        p = p%pow(10,i)
    print("-"*18)

def tex_command(p,num,f):
    """
    output tex command to file
    """
    f.write("\\begin{align} \n")
    for i in range(0,216,18)[::-1]:
        f.write("&")
        string = str(p/pow(10,i) + pow(10,18))[1:]
        string = string.replace(str(num), "\color{red}{%d}"%num)
        f.write(string)
        p = p%pow(10,i)
        f.write("\\\\ \n")
    f.write("\\end{align} \n")

letters = ["""
000000000000000000
000000001111000000
000000011111000000
000000111111000000
000001101111000000
000000001111000000
000000001111000000
000000001111000000
000000001111000000
000000001111000000
000001111111111000
000000000000000000
""","""
000000000000000000
000000222222000000
000022222222220000
000222000002220000
000000000022220000
000000000222200000
000000002222000000
000000022220000000
000002222000000000
000222222222222000
000222222222222000
000000000000000000
""","""
000000000000000000
000000333333000000
000033333333330000
000333000003333000
000000000033333000
000000033333300000
000000033333300000
000000000033333000
000333000003333000
000033333333330000
000000333333000000
000000000000000000
""","""
000000000000000000
000000004444400000
000000044444400000
000000444444400000
000004440444400000
000044400444400000
000444000444400000
004444444444444400
004444444444444400
000000000444400000
000000000444400000
000000000000000000
"""]

pc = PrimeChecker()

f = open("result.txt","w")
for n, letter in enumerate(letters):
    num = n+1
    prime_letter = int(letter.replace("\n",""))
    for i in range(1000,10000)[::-1]:
        if (i/1000)!=num and ((i%1000)/100)!=num and ((i%100)/10)!=num and i%10 != num:
            pl = prime_letter
            pl += pow(10,215) * (i/1000)
            pl += pow(10,198) * ((i%1000)/100)
            pl += pow(10,17) * ((i%100)/10)
            pl += i%10
            if pc.is_prime(pl):
                display_letter(pl)
                tex_command(pl,num,f)
                break
f.close()

実際の素数

1

900000000000000007000000001111000000000000011111000000000000111111000000000001101111000000000000001111000000000000001111000000000000001111000000000000001111000000000000001111000000000001111111111000400000000000000003
\begin{align}
&900000000000000007\\
&00000000\color{red}{1}\color{red}{1}\color{red}{1}\color{red}{1}000000\\
&0000000\color{red}{1}\color{red}{1}\color{red}{1}\color{red}{1}\color{red}{1}000000\\
&000000\color{red}{1}\color{red}{1}\color{red}{1}\color{red}{1}\color{red}{1}\color{red}{1}000000\\
&00000\color{red}{1}\color{red}{1}0\color{red}{1}\color{red}{1}\color{red}{1}\color{red}{1}000000\\
&00000000\color{red}{1}\color{red}{1}\color{red}{1}\color{red}{1}000000\\
&00000000\color{red}{1}\color{red}{1}\color{red}{1}\color{red}{1}000000\\
&00000000\color{red}{1}\color{red}{1}\color{red}{1}\color{red}{1}000000\\
&00000000\color{red}{1}\color{red}{1}\color{red}{1}\color{red}{1}000000\\
&00000000\color{red}{1}\color{red}{1}\color{red}{1}\color{red}{1}000000\\
&00000\color{red}{1}\color{red}{1}\color{red}{1}\color{red}{1}\color{red}{1}\color{red}{1}\color{red}{1}\color{red}{1}\color{red}{1}\color{red}{1}000\\
&400000000000000003\\
\end{align}

2

900000000000000004000000222222000000000022222222220000000222000002220000000000000022220000000000000222200000000000002222000000000000022220000000000002222000000000000222222222222000000222222222222000400000000000000009
\begin{align}
&900000000000000004\\
&000000\color{red}{2}\color{red}{2}\color{red}{2}\color{red}{2}\color{red}{2}\color{red}{2}000000\\
&0000\color{red}{2}\color{red}{2}\color{red}{2}\color{red}{2}\color{red}{2}\color{red}{2}\color{red}{2}\color{red}{2}\color{red}{2}\color{red}{2}0000\\
&000\color{red}{2}\color{red}{2}\color{red}{2}00000\color{red}{2}\color{red}{2}\color{red}{2}0000\\
&0000000000\color{red}{2}\color{red}{2}\color{red}{2}\color{red}{2}0000\\
&000000000\color{red}{2}\color{red}{2}\color{red}{2}\color{red}{2}00000\\
&00000000\color{red}{2}\color{red}{2}\color{red}{2}\color{red}{2}000000\\
&0000000\color{red}{2}\color{red}{2}\color{red}{2}\color{red}{2}0000000\\
&00000\color{red}{2}\color{red}{2}\color{red}{2}\color{red}{2}000000000\\
&000\color{red}{2}\color{red}{2}\color{red}{2}\color{red}{2}\color{red}{2}\color{red}{2}\color{red}{2}\color{red}{2}\color{red}{2}\color{red}{2}\color{red}{2}\color{red}{2}000\\
&000\color{red}{2}\color{red}{2}\color{red}{2}\color{red}{2}\color{red}{2}\color{red}{2}\color{red}{2}\color{red}{2}\color{red}{2}\color{red}{2}\color{red}{2}\color{red}{2}000\\
&400000000000000009\\
\end{align}

3

800000000000000002000000333333000000000033333333330000000333000003333000000000000033333000000000033333300000000000033333300000000000000033333000000333000003333000000033333333330000000000333333000000400000000000000009
\begin{align}
&800000000000000002\\
&000000\color{red}{3}\color{red}{3}\color{red}{3}\color{red}{3}\color{red}{3}\color{red}{3}000000\\
&0000\color{red}{3}\color{red}{3}\color{red}{3}\color{red}{3}\color{red}{3}\color{red}{3}\color{red}{3}\color{red}{3}\color{red}{3}\color{red}{3}0000\\
&000\color{red}{3}\color{red}{3}\color{red}{3}00000\color{red}{3}\color{red}{3}\color{red}{3}\color{red}{3}000\\
&0000000000\color{red}{3}\color{red}{3}\color{red}{3}\color{red}{3}\color{red}{3}000\\
&0000000\color{red}{3}\color{red}{3}\color{red}{3}\color{red}{3}\color{red}{3}\color{red}{3}00000\\
&0000000\color{red}{3}\color{red}{3}\color{red}{3}\color{red}{3}\color{red}{3}\color{red}{3}00000\\
&0000000000\color{red}{3}\color{red}{3}\color{red}{3}\color{red}{3}\color{red}{3}000\\
&000\color{red}{3}\color{red}{3}\color{red}{3}00000\color{red}{3}\color{red}{3}\color{red}{3}\color{red}{3}000\\
&0000\color{red}{3}\color{red}{3}\color{red}{3}\color{red}{3}\color{red}{3}\color{red}{3}\color{red}{3}\color{red}{3}\color{red}{3}\color{red}{3}0000\\
&000000\color{red}{3}\color{red}{3}\color{red}{3}\color{red}{3}\color{red}{3}\color{red}{3}000000\\
&400000000000000009\\
\end{align}

4

900000000000000009000000004444400000000000044444400000000000444444400000000004440444400000000044400444400000000444000444400000004444444444444400004444444444444400000000000444400000000000000444400000500000000000000003
\begin{align}
&900000000000000009\\
&00000000\color{red}{4}\color{red}{4}\color{red}{4}\color{red}{4}\color{red}{4}00000\\
&0000000\color{red}{4}\color{red}{4}\color{red}{4}\color{red}{4}\color{red}{4}\color{red}{4}00000\\
&000000\color{red}{4}\color{red}{4}\color{red}{4}\color{red}{4}\color{red}{4}\color{red}{4}\color{red}{4}00000\\
&00000\color{red}{4}\color{red}{4}\color{red}{4}0\color{red}{4}\color{red}{4}\color{red}{4}\color{red}{4}00000\\
&0000\color{red}{4}\color{red}{4}\color{red}{4}00\color{red}{4}\color{red}{4}\color{red}{4}\color{red}{4}00000\\
&000\color{red}{4}\color{red}{4}\color{red}{4}000\color{red}{4}\color{red}{4}\color{red}{4}\color{red}{4}00000\\
&00\color{red}{4}\color{red}{4}\color{red}{4}\color{red}{4}\color{red}{4}\color{red}{4}\color{red}{4}\color{red}{4}\color{red}{4}\color{red}{4}\color{red}{4}\color{red}{4}\color{red}{4}\color{red}{4}00\\
&00\color{red}{4}\color{red}{4}\color{red}{4}\color{red}{4}\color{red}{4}\color{red}{4}\color{red}{4}\color{red}{4}\color{red}{4}\color{red}{4}\color{red}{4}\color{red}{4}\color{red}{4}\color{red}{4}00\\
&000000000\color{red}{4}\color{red}{4}\color{red}{4}\color{red}{4}00000\\
&000000000\color{red}{4}\color{red}{4}\color{red}{4}\color{red}{4}00000\\
&500000000000000003\\
\end{align}

本当は5,6,7,8,9と求めたかったけれど、文字を作るのが面倒で断念してしまいました。
誰か頑張って作ってください。。。

(追記)

そこまで珍しいことだとは思わず、一例しか挙げていませんでした。
存在自体が不思議という話があったので、数を数えてみたのですが
1: 14個
2: 15個
3: 4個
4: 17個
と、そこまで珍しくないことがわかります*1

素数定理に従えば、216桁以下の素数の存在確率は1/log(10^217)≒0.2%らしいです。
リーマン予想の意味,素数分布との関係 | 高校数学の美しい物語
5832*0.002 = 11.6個
だいたいそれくらいですね。

5832個の数字が独立に0.2%の確率で素数かどうか決まるとすると、
各Nに対してそのような素数が存在しない確率は実に0.000850%と非常に低い確率であることがわかります。

ASCIIコード素数を作りたい

それでは、画素数を増やしたいとすればどうなるのでしょうか。
ここでASCIIコードの文字(93文字分)を出来るだけ大きい画素数素数で表現したいとします。
「各文字に対して99.9%の確率でそのような素数が作れます」と言われれば、まあ使い物になりそうだなと言えそうですね。

素数定理に従えば、桁数をnとした時、
\displaystyle
\left( 1 - \frac{1}{nlog10} \right)^{5832} > 0.001
を満たす最大のnを求めればよいことになります。
これを計算するとn=366桁で、
「ASCIIコード素数を作るには19*19くらいの画素数まで大きくできそうだ」
ということが分かりますね。

是非誰か19×19の素数でASCIIコード素数を作ってみてください

*1:4~17/5832を珍しいと考えるかどうかは人によるけれど

2次元配列の特定の条件を満たす要素の中からランダムに選択【numpy】

「0~1の値を取る2次元配列があって、0.7以上の要素を持つインデックスの中からランダムにインデックスを取ってくる。」
こんなことをしたい時にnumpyを駆使すれば綺麗に書けたのでメモしておく。


例えばこんな配列(numpyのarray)があったとする

import numpy as np
a = np.random.randint(0,10,size=[8,6])
array([[6, 9, 4, 1, 2, 4],
       [4, 9, 3, 8, 9, 7],
       [9, 6, 6, 8, 9, 7],
       [6, 6, 2, 7, 2, 3],
       [9, 9, 3, 2, 1, 0],
       [9, 9, 0, 3, 8, 2],
       [0, 5, 4, 4, 0, 7],
       [1, 8, 6, 0, 9, 4]])

この中から8以上の値を持つ要素は

a > 7
array([[False,  True, False, False, False, False],
       [False,  True, False,  True,  True, False],
       [ True, False, False,  True,  True, False],
       [False, False, False, False, False, False],
       [ True,  True, False, False, False, False],
       [ True,  True, False, False,  True, False],
       [False, False, False, False, False, False],
       [False,  True, False, False,  True, False]], dtype=bool)

これだけある。
このTrueの中からランダムに一つインデックスを選びたい。
例えば、[0,1]とかだ。

これをnumpyを使わずに書こうとするとこんな風になるだろう。

import random
index_list = []
for i in range(len(a)):
    for j in range(len(a[0])):
        if a[i][j] > 7:
            index_list.append([i,j])
coord = random.choice(index_list)

randomの便利な関数のおかげでこれでも割とすっきりしているが、
numpyを使えばさらにすっきり書ける。

index = np.random.choice(np.where(a.reshape(-1,) > 7)[0])
coord = [index/a.shape[1], index%a.shape[1]]

たったこれだけだ。
numpy凄い。



※解説

  • reshape[-1,] ... 2次元配列を1次元配列に変換できる
  • np.where ... Trueの箇所のインデックスを返してくれる
  • np.random.choice ... 配列の中から一つ選び取ってくれる