← PC版は別ページ

≪いっしょにPython≫　プログラマーの実験ノート

　このページは，プログラム言語Pythonをこれから学ぼうと考えている筆者の備忘録です．
　そこそこ調べて書いていますが，仕様を精密に調べたものではありません．どちらかと言えば，実装されたPython 3を使ってみてどうなったかということを中心に書いたものです．いわば，実験ノートの記録です．（2019.5.10，Python3.7.3に基づいて作成）

Pythonの練習には，
の３通りが使えるが，以下は(2)のNew Fileから，テキストエディタを使って，プログラムを書く場合を考える．

• 初等整数論には，問題の意味は小中学生でも分かるが，専門家でも証明できていない未解決問題が多くあります．3n+1問題（コラッツ予想）と呼ばれる問題もその１つです．
• 3n+1問題（コラッツ予想）とは
　1930年代にコラッツは，どんな正の整数から始めても上記の操作を繰り返せば１になると予想したが，この予想は証明も反証もできていない未解決問題となっている．（3n+1問題とも呼ばれる）
【例】
≪3から始めると≫
3→10→5→16→8→4→2→1 （7回の操作で1になる）
≪11から始めると≫
11→34→17→52→26→13→40→20→10→5→16
→8→4→2→1 （14回の操作で1になる）
≪35から始めると≫
35→106→53→160→80→40→20→10→5→16
→8→4→2→1 （13回の操作で1になる）
　中には非常に長い経過をたどるものがあり，27から始めると111回で１になる．
≪27から始めると≫
27 →82 →41 →124 →62 →31 →94 →47 →142 →71 →214 →107 →322 →161 →484 →242 →121 →364 →182 →91 →274 →137 →412 →206 →103 →310 →155 →466 →233 →700 →350 →175 →526 →263 →790 →395 →1186 →593 →1780 →890 →445 →1336 →668 →334 →167 →502 →251 →754 →377 →1132 →566 →283 →850 →425 →1276 →638 →319 →958 →479 →1438 →719 →2158 →1079 →3238 →1619 →4858 →2429 →7288 →3644 →1822 →911 →2734 →1367 →4102 →2051 →6154 →3077 →9232 →4616 →2308 →1154 →577 →1732 →866 →433 →1300 →650 →325 →976 →488 →244 →122 →61 →184 →92 →46 →23 →70 →35 →106 →53 →160 →80 →40 →20 →10 →5 →16 →8 →4 →2 →1
本当に証明したいことは，次の(1)(2)(3)にまとめられる．
(1) どんな正の整数から始めても，上記の(�T)(�U)の操作を有限回繰り返せば1になる．
(2) どんな正の整数から始めても，同じ数字に戻って来ることはない．
(3) どんな正の整数から始めても，限りなく大きくなることはない．

== 10までの数で確認実験をする ==
【例 1.1】

import datetime as dt #1
dt1 = dt.datetime
tn1 = dt1.now() #2

data1 = []
def t3_p1(num1):
    count2 = 0 #3
    count3 = 0 #4
    tmp = num1
    while tmp != 1:
        if tmp % 2 == 1:
            tmp = tmp * 3 +1
            count3 += 1            
        else:
            tmp //= 2
            count2 += 1
    return num1, count2, count3, count2+count3 #5
for n in range(2, 10): #6
    data1.append(t3_p1(n))
tn2 = dt1.now() #7
print(data1)
print(tn2-tn1)

→
[(2, 1, 0, 1), (3, 5, 2, 7), (4, 2, 0, 2), (5, 4, 1, 5), (6, 6, 2, 8), (7, 11, 5, 16), (8, 3, 0, 3), (9, 13, 6, 19)]
0:00:00
#1 所要時間を測定するために，datetimeモジュールをインポートする．
#2 datetimeオブジェクトのnow()メソッドで，それを呼び出した日時をマイクロ秒まで測定できる．#7との差が所要時間になる．
#3 2で割った回数を数える．
#4 3を掛けて1加えた回数を調べる．
#5 t3_p1()関数は，引数num1, 2で割った回数，3を掛けて1加えた回数, 合計操作回数のタプルを返す．
例えば，3からスタートすれば，(3, 5, 2, 7)のように返される．
• 所要時間は0秒と表示されるが，これは測定できないほどの一瞬だという意味

• あまり多くのデータが出力されても困るので，1000以下，10000以下，100000以下の整数について，合計操作回数が多いものベストテンを見つけることにする．

== 回数の多いものベストテンを表示する ==
【例 1.2】

import datetime as dt
dt1 = dt.datetime
tn1 = dt1.now()

data1 = []
def t3_p1(num1):
    count2 = 0
    count3 = 0
    tmp = num1
    while tmp != 1:
        if tmp % 2 == 1:
            tmp = tmp * 3 +1
            count3 += 1            
        else:
            tmp //= 2
            count2 += 1
    return num1, count2, count3, count2+count3
for n in range(2, 1000):###
    total1 = t3_p1(n)
    if total1[3] > 100:#1
        data1.append(t3_p1(n))#2
data2 = sorted(data1, key= lambda tp:tp[3],\
 reverse=True)#3
tn2 = dt1.now()
disp_count = min(10, len(data2))#4
for k in range(disp_count):
    print(data2[k])
print(tn2-tn1)

→
(871, 113, 65, 178)
(937, 110, 63, 173)
(703, 108, 62, 170)
(763, 97, 55, 152)
(775, 97, 55, 152)
(859, 94, 53, 147)
(865, 94, 53, 147)
(873, 94, 53, 147)
(879, 94, 53, 147)
(889, 94, 53, 147)
0:00:00.040999
• 途中の操作関数は，前記と全く同じ．データが多くなり過ぎないように，結果の処理を工夫する．
#1,#2 合計回数が100回よりも多いものだけ記録する．
#3 記録されたdeta1を合計回数の多いもの順にソートする．
#4 ベストテンだけを選んで表示する．
• ###の第2引数を変えれば上限を変更できる．1000までの所要時間は0.04秒．
• 10000までで0.63秒，最長は6171の261回．100000までで8.98秒，最長は77031の350回．これ以上になると所要時間が長くなる．（所要時間が長くなるのは，###においてその数以下のすべての数を調べるためで，例えば5×2⁶⁰（19桁の整数）を単発で調べるのは0.02秒程度でできる）

• ペル方程式とは，
• Dが1けたの整数である場合は，x, yに様々な値を試してみると，短時間で１つの解が見つかる．
【例】
• 正の整数解のうちでxの値が最小のものx=a, y=bの組を最小解と呼ぶとき，一般に他の解は次の連立漸化式を満たす数列の一般項として表せることが知られている．

　したがって，与えられたDの値に対して，最小解x=a, y=bを求めることが関心事となるが，Dが2桁以上の整数になると，最小解を見つけることが困難になって来る．
　例えば，D=61の場合には，最小解のxは10桁の整数，yは9桁の整数となり，これよりも小さな解はない．
• このように，最小解を見つけること自体は，試し打ちになるので，コンピュータを使って網羅的に調査するのに適した仕事になる．
　D=61, 97, 109のように途方もない計算になるものが含まれるので，初めはその数を避けるとよい．

== １つの方程式で3秒以上要するものは，TimeOutとして放棄する ==
#1,#2 所要時間が3秒を超えるものはTimeOverとする
#3 だから，yが与えられたとき，xはの整数部分と見てよい
#4 解が見つかれば (D, x, y)の形のタプルでリストに追加する
#5 Dの値の範囲を変更するには，#5の第2引数を変更する
#6 Dが平方数であればペル方程式にならないので除外する．とその整数部分の差が0.0ならば平方数と見なす．
【例 2.1】

import math as ma
import datetime as dt
dt1 = dt.datetime
sol_list = []
def pell1(D):
    tn1 = dt1.now()
    y = 1
    while True:
        tn2 = dt1.now()
        if str(tn2-tn1) > '0:00:03': #1
           sol_list.append((D,'TimeOver')) #2
           return
        x = int(ma.sqrt(D)*y+1) #3
        if x**2 - D*(y**2) == 1:
            sol_list.append((D,x,y)) #4
            return
        else:
            y += 1
            continue
for D in range(2, 200): #5
    if ma.sqrt(D)-int(ma.sqrt(D)) == 0.0: #6
        continue;
    else: 
       pell1(D)
print(sol_list)

→
[(2, 3, 2), (3, 2, 1), (5, 9, 4), (6, 5, 2), (7, 8, 3), …
(59, 530, 69), (60, 31, 4), (61, 'TimeOver'), (62, 63, 8), …
(96, 49, 5), (97, 'TimeOver'), (98, 99, 10), (99, 10, 1),…
(106, 'TimeOver'), (107, 962, 93), (108, 1351, 130), (109, 'TimeOver'), …]
• D<200でTimeOverとなる値は
D=61, 97, 106, 109, 137, 139, 149, 151, 157, 163, 166, 172, 181, 191, 193, 199
• D<200でxが5桁以上になるものは，(D, x, y)の順に

(46, 24335, 3588)
(53, 66249, 9100)
(58, 19603, 2574)
(67, 48842, 5967)
(73, 2281249, 267000)
(76, 57799, 6630)
(85, 285769, 30996)
(86, 10405, 1122)
(89, 500001, 53000)
(93, 12151, 1260)
(94, 2143295, 221064)
(103, 227528, 22419)
(113, 1204353, 113296)
(118, 306917, 28254)
(124, 4620799, 414960)
(125, 930249, 83204)
(127, 4730624, 419775)
(129, 16855, 1484)
(131, 10610, 927)
(133, 2588599, 224460)
(134, 145925, 12606)
(154, 21295, 1716)
(161, 11775, 928)
(162, 19601, 1540)
(173, 2499849, 190060)
(177, 62423, 4692)
(179, 4190210, 313191)
(184, 24335, 1794)
(190, 52021, 3774)

• 高校数学Bでは，自然数の累乗の和について，次のような公式を習う．



以上から
…�@
が成り立つ．
• 高校では，これ以上の話を習ったことも，教えたこともないが，自然数のｋ乗の和を

で表すとき，との間の関係式を�@と類似の形で求めてみようという作戦です．
• たとえば，は6次式では4次式なので，
• 連立方程式は，Pythonのsympyをインポートして，solve([])とすれば解ける．

との関係式を求める（必要条件で係数a,b,cを求める）
と仮定する
【例 3.1】

import sympy as sp
a = sp.Symbol('a')
b = sp.Symbol('b')
c = sp.Symbol('c')
sol1 = sp.solve([a+b+c-1, 12*a+6*b+3*c-11,\
 27*a+9*b+3*c-23])
print(sol1)

→
{a: 2/3, b: 2/3, c: -1/3}
したがって
　
が必要条件（十分性の証明は数学的帰納法による．ここでは略）

との関係式を求める（必要条件で係数a,b,cを求める）
と仮定する
【例 3.2】

import sympy as sp
a = sp.Symbol('a')
b = sp.Symbol('b')
c = sp.Symbol('c')
sol1 = sp.solve([a+b+c-1, 20*a+10*b+5*c-17,\
 9*a+3*b+c-7])
print(sol1)

→
{a: 3/5, b: 3/5, c: -1/5}
したがって
　
が必要条件（十分性の証明は数学的帰納法による．ここでは略）

との関係式を求める（必要条件で係数a,b,c,dを求める）
と仮定する
【例 3.3】

import sympy as sp
a = sp.Symbol('a')
b = sp.Symbol('b')
c = sp.Symbol('c')
d = sp.Symbol('d')
sol1 = sp.solve([a+b+c+d-1, 24*a+12*b+6*c+3*d-17,\
 81*a+27*b+9*c+3*d-49,320*a+80*b+20*c+5*d-177])
print(sol1)

→
{a: 2/5, b: 3/5, c: 1/15, d: -1/15}
したがって
　

が必要条件（十分性の証明は数学的帰納法による．ここでは略）

• ここで述べる内容は，他の項目とは異なり，高校の教科書に書かれている証明済みのものです．
• プログラミングの入門として，再帰的に繰り返す方法とループを使って繰り返す方法があり，参考になりそうなので掲載する．解読のしやすさや動作の速さなどを大きな値を使って比較できる．

【ユークリッドの互除法とは】
• ２つの正整数a, bの最大公約数Gと最小公倍数Lを求めたいとき，小中学校では次の(A)の方法，中高校では(B)の方法を使う．
(A) 2つの正整数を共通に割れるだけ割っていく．共通に割れるものがなくなったら止める．
共通な数を掛けたものが最大公約数G，共通でない数も掛けたものが最小公倍数L
(B) 2つの正整数を素因数分解し，指数の低い方を集めたものが最大公約数G，指数の高い方を集めたものが最小公倍数L．（最小の指数を集めたら最大公約数，最大の指数を集めたら最小公倍数，言葉と表すものが逆なので要注意）
• これらの方法(A)(B)はいずれも，共通因数（何で割り切れるか）が分かることを前提としており，そもそも何で割り切れるか分からない次の例のような場合には，どちらも使えない．

• このように，そもそも何で割り切れるかが分からない場合でも，最大公約数を求められる方法がユークリッドの互除法で，高校数学Ａで整数問題の分野を選択したら習う．その内容は次の通り． • ≪例≫ (18, 12)→(12, 6)→(6, 0)→G=6
• ≪証明≫
　a, bの最大公約数をGとおくと，a=Ga’, b=b’（a’, b’は互いに素）とおける．
　このとき，a=bq+r (0≦r<b)ならばGa’=Gb’q+rとなるから
r=G(a’−b’q)となって，rもGで割り切れる．したがって，b, rはGで割り切れる．
　もしも，b, rの最大公約数がH(>G)となるならば，a=Hb”+H”=H(b”+”)となって仮定に反するから，Gよりも大きな公約数はない．
したがって，Gがb, rの最大公約数．
　以上により，(a, b)の最大公約数は(b, r)の最大公約数に等しいから，(a, b)の最大公約数を求める代わりに(b, r)の最大公約数を求めたらよい．
• a>bのとき，r<bだから，有限の整数a, bからスタートする限り，毎回の割り算でrは必ずbよりも小さくなるから，この操作は有限回で終わる．

• 最大公約数Gが求まれば，最小公倍数は次の公式によって求められる．

• 基数（radix）または基底（base）とは，各位の数に使える文字の個数で，これよりも大きな数は上の桁を使って表す．
• 10進法は0,1,2,3,4,5,6,7,8,9の10種類の数字を使って数を表し，
例えば，をと書く．
　2進法は0,1の2種類の数字を使って数を表し，
例えば，をと書く．
　3進法は0,1,3の2種類の数字を使って数を表し，
例えば，をと書く．
• 11進法以上は，1文字で10以上の数を表せるように英字も使って表すのが普通である．
10個の数字（0123456789）と26文字の英小文字（abcdefghijklmnopqrstuvwxyz）を使えば，10+26=36進法まで表示できるが，さらに26文字の英大文字（ABCDEFGHIJKLMNOPQRSTUVWXYZ）も区別して使えば，62進法まで表せる．区切りの良い所で64進法まで表示するには，コンピュータシステムごとに差し支えない2文字，例えば（_~）を追加して64文字にするとよい．
• 16進法は0〜9,a〜fの16種類の文字を使って表し，
例えば，をと書く．
　38進法は0〜9,a〜z, A〜Bの38種類の文字を使って表し，
例えば，をと書く．
（64進法というのは，「アスキー文字だけで表示しよう」と考えた場合の上限になっているだけで，コンピュータの内部処理としては，32768進法でも可能なことに注意．例えば各要素が0≦x<32767となる整数からなる要素数nのリストを作ると，32768進数でn桁の整数を表すことができる．）
• 10進法の整数を，2進法，8進法，16進法に直す場合は，Pythonのビルドイン関数，bin(), oct(), hex()が利用でき, 任意のm進法の整数を10進法に変換するにはint()関数が利用できるので，自作関数の点検に使える．以下においては，任意のm進法から任意のn進法に変換する場合を想定する．
　※余談であるが，16進法，32進法という場合の基数は，16，32と書き，g進法やv進法とは書かない．

m進法で書かれた正の整数をn進法に直す
【例 4.2】

N_set ='0123456789abcdefghijklmnopqrstuvwxyz\
ABCDEFGHIJKLMNOPQRSTUVWXYZ_~'#1
str_n = ''#2
def is_m(str1, m):
    for x in str1:
        if N_set.find(x, 0, m) == -1:#3
            return False
    return True
def m2n(str1, radix_m, radix_n):
    global str_n
    if is_m(str1, radix_m) == False:
        return print('元の基数が範囲外')
    if radix_n > len(N_set):
        return print('変換先の基数が範囲外')
    L1 = list(str1)
    L2 = list(reversed(L1))#4
    N1 = 0
    for k in range(len(L2)):
        N1 += N_set.find(L2[k], 0) * (radix_m**k)
    L3 = []
    while N1 > 0:
        r1 = N1 % radix_n#5
        L3.append(r1)
        N1 = (N1 - r1)//radix_n #6
    L4 = list(reversed(L3))#7
    for x in L4:
        str_n += N_set[x]#8
    print(str_n)
    
m2n('ab90e', 16, 10)#9

→
702734
#1 10よりも基数が大きい場合に対応するため，アルファベットの大文字と小文字，特殊記号2文字を使って64進法まで表す．（なお，ビルドイン関数のint()では，大文字と小文字を区別していないので，37進法以上を指定するとエラーになる．）
#2 変換先を文字列として出力するための入れ物
#3 文字列.find(sub, start, end)により，与えらえた文字列において，文字列subをstart（省略すれば0）からend（省略すれば末尾）まで探し，最初に見つかった位置（先頭0からの番号）を返す．ただし，endに書かれた数字の前まで調べ，end自身は入らないので注意
　関数is_m(str1, m)は，str1に書かれた文字がm進法で使用可能かどうかチェックする．矛盾があれば，以後の計算は行われない．
#4 文字列で与えらえた元の数をリストにして，逆順のリストに直す．
#5, #6がこの変換プログラムの核心．変換先の基数で割った余りが第1位の数．#6により商が0になるまで繰り返す．
#7 できたリストを逆順にする．#8 基数の約束に従って文字を当てはめる．
#9 具体的には m2n('変換元の文字列', 変換元の基数, 変換先の基数)のように書くと，結果が表示される．