元号ゲームの結果をながめていたら，「弘文」が目に入った。こんな元号あったんじゃないか。調べてみると天皇の名称（弘文天皇）だった。確かに，過去の元号制定の条件（昭和大礼記録）では「元号は本邦はもとより言うを俟たず、支那、朝鮮、南詔、交趾（ベトナム）等の年号、その帝王、后妃、人臣の諡号、名字等及び宮殿、土地の名称等と重複せざるものなるべきこと」となっているので，過去の天皇の名称（漢風諡号）は排除すべきだった。
明日からの旅行の準備に忙しくてもう時間がない。面倒だなと思いつつ，日本の天皇のリストでもつくっておこうとして，天皇名のCSVファイルをつくって，EXCELに読み込ませると，何と文字化け。EXCELはUNICODEのCSVファイルをちゃんと読み込んでくれないのだった。SHIFT-JISに変換するとよいのだが，面倒な話だ。ということで，ここで意欲を失ったので，今日はおしまい。というか明日発表なので，もうおしまい。
とりあえず，データだけあげておく（正しいかどうか確認していないので注意せよ）。
"1","神武","2","綏靖","3","安寧","4","懿徳","5","考昭","6","考安","7","考霊","8","考元","9","開化","10","崇神","11","垂仁","12","景行","13","成務","14","仲哀","15","応神","16","仁徳","17","履中","18","反正","19","允恭","20","安康","21","雄略","22","清寧","23","顕宗","24","仁賢","25","武烈","26","継体","27","安閑","28","宣化","29","欽明","30","敏達","31","用明","32","崇峻","33","推古","34","舒明","35","皇極","36","考徳","37","斉明","38","天智","39","弘文","40","天武","41","持統","42","文武","43","元明","44","元正","45","聖武","46","考謙","47","淳仁","48","称徳","49","光仁","50","桓武","51","平城","52","嵯峨","53","淳和","54","仁明","55","文徳","56","清和","57","陽成","58","光考","59","宇多","60","醍醐","61","朱雀","62","村上","63","冷泉","64","円融","65","花山","66","一条","67","三条","68","後一条","69","後朱雀","70","後冷泉","71","後三条","72","白河","73","堀河","74","鳥羽","75","崇徳","76","近衛","77","後白河","78","二条","79","六条","80","高倉","81","安徳","82","後鳥羽","83","土御門","84","順徳","85","仲恭","86","後堀河","87","四条","88","後嵯峨","89","後深草","90","亀山","91","後宇多","92","伏見","93","後伏見","94","後二条","95","花園","96","後醍醐","97","後村上","98","長慶","99","後亀山","100","後小松","101","称光","102","後花園","103","後土御門","104","後柏原","105","後奈良","106","正親町","107","後陽成","108","後水尾","109","明正","110","後光明","111","後西","112","霊元","113","東山","114","中御門","115","桜町","117","後桜町","118","後桃園","119","光格","120","仁考","121","孝明","122","明治","123","大正","124","昭和","125","今上",
P. S. ミカド文庫（モラロジー研究所由来）というサイトに，「日本の年号候補・未採用文字案」というページを発見した。まじめに取り組むとこういうふうになる。
P. P. S. 明日から１週間ほどBLOG投稿はお休みになると思いますのでよろしく。
2019年3月30日土曜日
元号ゲーム（１）
前回の元号予想プログラムを少しチューニングした。条件は次のとおり。
(1) これまでの日本の元号に登場した漢字の範囲で考える。ただし１文字目と２文字目は区別して考える（４文字元号は２文字元号に読み替える）。
(2) １文字目の漢字と２文字目の漢字のそれぞれの出現確率の積で順序づける。
(3) 既存のアジア（日本，中国，朝鮮，ベトナム，台湾）の元号とそれを入れ替えたものに一致するものは排除する（例えば，「安久」は排除される）。
(4) 慶応，明治，大正，昭和，平成の各漢字は排除する。
(5) 画数の多い漢字，MTSHで始まる漢字などを排除する（「など」のところにプログラマの主観が入る）
昭和とか平成の場合を考えると，このアルゴリズムはよくないことがわかっているが，目をつぶることにしよう（そりゃ，「苦節」とか「余命」とか「結婚」とか「懲役」とか「西暦」とか「卍」とか「タピオカ」のほうがおもしろいに決まっている）。とにかくこの条件で結果抽出されたものが下記のとおり。
(250, 839, 26, 159, 17)
("弘久", 0.020164609053497942)
("永化", 0.019753086419753086)
("元化", 0.018518518518518517)
("万安", 0.016049382716049384)
("万元", 0.014814814814814815)
("文弘", 0.01316872427983539)
("弘文", 0.011522633744855968)
("万久", 0.008641975308641974)
("文万", 0.006584362139917695)
("文同", 0.006584362139917695)
("文字", 0.006584362139917695)
("永同", 0.006584362139917695)
("永字", 0.006584362139917695)
("元万", 0.006172839506172839)
("元同", 0.006172839506172839)
("元字", 0.006172839506172839)
("安化", 0.006172839506172839)
("斉永", 0.005761316872427984)
("朱永", 0.005761316872427984)
("白永", 0.005761316872427984)
("至永", 0.005761316872427984)
("斉安", 0.005349794238683128)
("昌安", 0.005349794238683128)
("朱安", 0.005349794238683128)
("白安", 0.005349794238683128)
("至安", 0.005349794238683128)
("万文", 0.0049382716049382715)
("斉元", 0.0049382716049382715)
("昌元", 0.0049382716049382715)
("朱元", 0.0049382716049382715)
("白元", 0.0049382716049382715)
("万化", 0.003703703703703704)
("弘万", 0.002880658436213992)
("弘同", 0.002880658436213992)
("弘字", 0.002880658436213992)
("斉久", 0.002880658436213992)
("昌久", 0.002880658436213992)
("朱久", 0.002880658436213992)
("白久", 0.002880658436213992)
("至久", 0.002880658436213992)
("万弘", 0.0024691358024691358)
("久化", 0.0024691358024691358)
("安万", 0.00205761316872428)
("安同", 0.00205761316872428)
("安字", 0.00205761316872428)
("久弘", 0.0016460905349794238)
("斉文", 0.0016460905349794238)
("昌文", 0.0016460905349794238)
("朱文", 0.0016460905349794238)
("白文", 0.0016460905349794238)
Juliaが使える方は，これで適当に遊んでください。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
function exch(knj)
#
# knj input kanji string with 2 char
# return exchanged string with 2 char
return knj[4]*knj[1]
end
function gfrq(c,n,r,m)
#
# c input array of kanji char
# n size of c
# r output array of sorted (kanji, count)
# m size of r
#
sort!(c)
k=0
for j=2:n
k=k+1
if c[j]!=c[j-1]
push!(r,(c[j-1],k))
k=0
end
end
sort!(r, by=x->x[2], rev=true)
m=length(r)
# println((m,r))
end
#japan
g1="大化,白雉,朱鳥,大宝,慶雲,和銅,霊亀,養老,神亀,天平,感宝,勝宝,宝字,神護,景雲,宝亀,天応,延暦,大同,弘仁,天長,承和,嘉祥,仁寿,斉衡,天安,貞観,元慶,仁和,寛平,昌泰,延喜,延長,承平,天慶,天暦,天徳,応和,康保,安和,天禄,天延,貞元,天元,永観,寛和,永延,永祚,正暦,長徳,長保,寛弘,長和,寛仁,治安,万寿,長元,長暦,長久,寛徳,永承,天喜,康平,治暦,延久,承保,承暦,永保,応徳,寛治,嘉保,永長,承徳,康和,長治,嘉承,天仁,天永,永久,元永,保安,天治,大治,天承,長承,保延,永治,康治,天養,久安,仁平,久寿,保元,平治,永暦,応保,長寛,永万,仁安,嘉応,承安,安元,治承,養和,寿永,元暦,文治,建久,正治,建仁,元久,建永,承元,建暦,建保,承久,承応,元仁,嘉禄,安貞,寛喜,貞永,天福,文暦,嘉禎,暦仁,延応,仁治,寛元,宝治,建長,康元,正嘉,正元,文応,弘長,文永,建治,弘安,文永,建治,弘安,正応,永仁,正安,乾元,嘉元,徳治,延慶,応長,正和,文保,元応,元亨,正中,嘉暦,元徳,元弘,正慶,建武,延元,興国,正平,建徳,文中,天授,弘和,元中,暦応,康永,貞和,観応,文和,延文,康安,貞治,応安,永和,康暦,永徳,至徳,嘉慶,康応,明徳,応永,正長,永享,嘉吉,文安,宝徳,享徳,康正,長禄,寛正,文正,応仁,文明,長享,延徳,明応,文亀,永正,大永,享禄,天文,弘治,永禄,元亀,天正,文禄,慶長,元和,寛永,正保,慶安,承応,明暦,万治,寛文,延宝,天和,貞享,元禄,宝永,正徳,享保,元文,寛保,延享,寛延,宝暦,明和,安永,天明,寛政,享和,文化,文政,天保,弘化,嘉永,安政,万延,文久,元治,慶応,明治,大正,昭和,平成,"
#china
g2="建元,元光,元朔,元狩,元鼎,元封,太初,天漢,太始,征和,後元,始元,元鳳,元平,本始,地節,元康,神爵,五鳳,甘露,黄龍,初元,永光,建昭,竟寧,建始,河平,陽朔,鴻嘉,永始,元延,綏和,建平,太初,元将,元寿,元始,居摂,初始,始建,天鳳,地皇,更始,建武,建武,中元,永平,建初,元和,章和,永元,元興,延平,永初,元初,永寧,建光,延光,永建,陽嘉,永和,漢安,建康,永憙,本初,建和,和平,元嘉,永興,永寿,延熹,永康,建寧,熹平,光和,中平,光熹,昭寧,永漢,初平,興平,建安,延康,黄初,太和,青龍,景初,正始,嘉平,正元,甘露,景元,咸熙,章武,建興,延熙,景耀,炎興,黄武,黄龍,嘉禾,赤烏,太元,神鳳,建興,五鳳,太平,永安,元興,甘露,宝鼎,建衡,鳳凰,天冊,天璽,天紀,泰始,咸寧,太康,太熙,永熙,永平,元康,永康,永寧,太安,永安,建武,永興,光熙,永嘉,建興,建武,大興,永昌,太寧,咸和,咸康,建元,永和,升平,隆和,興寧,太和,咸安,寧康,太元,隆安,大亨,元興,義熙,元熙,桓玄,天康,永鳳,河瑞,光興,嘉平,建元,麟嘉,漢昌,光初,太和,建平,延熙,建武,太寧,青龍,永寧,永興,元璽,光寿,建熙,建初,建興,晏平,玉衡,玉恒,漢興,太和,嘉寧,建興,和平,升平,皇始,寿光,永興,甘露,建元,太安,太初,延初,元光,燕元,建興,永康,建始,延平,青龍,建平,長楽,光始,建始,燕興,更始,昌平,建明,建平,建武,中興,建平,太上,正始,太平,太興,白雀,建初,皇初,弘始,永和,建義,太初,更始,永康,建弘,永弘,龍昇,鳳翔,昌武,真興,承光,勝光,太安,麟嘉,龍飛,承康,咸寧,神鼎,太初,建和,弘昌,嘉平,庚子,建初,嘉興,永建,神璽,天璽,永安,玄始,承玄,義和,承和,承平,真興,承陽,縁禾,太縁,永初,景平,元嘉,太初,孝建,大明,永光,景和,泰始,泰豫,元徽,昇明,建元,永明,隆昌,延興,建武,永泰,永元,中興,天監,普通,大通,大通,大同,大同,太清,大宝,天正,太始,承聖,天成,紹泰,太平,天啓,大定,天保,広運,永定,天嘉,天康,光大,太建,至徳,禎明,登国,皇始,天興,天賜,永興,神瑞,泰常,始光,神䴥,延和,太延,太平,真君,正平,承平,興安,興光,太安,和平,天安,皇興,延興,承明,太和,景明,正始,永平,延昌,熙平,神亀,正光,孝昌,武泰,建義,永安,建明,普泰,中興,太昌,永興,永熙,天平,元象,興和,武定,大統,天保,乾明,皇建,太寧,河清,天統,武平,隆化,徳昌,承光,武成,保定,天和,建徳,宣政,大成,大象,大定,開皇,仁寿,大業,義寧,皇泰,武徳,貞観,永徽,顕慶,龍朔,麟徳,乾封,総章,咸亨,上元,儀鳳,調露,永隆,開耀,永淳,弘道,嗣聖,文明,光宅,垂拱,永昌,載初,神龍,景龍,唐隆,景雲,太極,延和,先天,開元,天宝,至徳,乾元,上元,宝応,広徳,永泰,大暦,建中,興元,貞元,永貞,元和,長慶,宝暦,大和,開成,会昌,大中,咸通,乾符,広明,中和,光啓,文徳,龍紀,大順,景福,乾寧,光化,天復,天祐,天授,如意,長寿,延載,証聖,天冊,万歳,万歳,登封,万歳,通天,神功,聖暦,久視,大足,長安,仁安,大興,宝暦,中興,正暦,永徳,朱雀,太始,建興,咸和,建初,承平,義熙,甘露,章和,永平,和平,建昌,延昌,延和,義和,重光,延寿,賛普,長寿,見龍,上元,元封,応道,龍興,全義,大豊,保和,天啓,建極,法尭,貞明,大同,嵯耶,中興,安国,始元,天瑞,景星,安和,貞祐,初歴,孝治,天応,尊聖,興聖,大明,鼎新,光聖,文徳,神武,文経,至治,明徳,広徳,順徳,明政,広明,明応,明統,明聖,明治,明啓,乾興,明通,正治,聖明,天明,保安,正安,正徳,保徳,明侯,上徳,広安,上明,保定,建安,天祐,上治,天授,開明,天政,文安,日新,文治,永嘉,保天,広運,永貞,大宝,龍興,盛明,建徳,利貞,盛徳,嘉会,元亨,安定,鳳歴,元寿,天開,天輔,仁寿,道隆,利正,興正,天定,彝泰,開平,乾化,鳳歴,貞明,龍徳,同光,天成,長興,応順,清泰,天福,開運,天福,乾祐,広順,顕徳,天復,天祐,武義,順義,乾貞,大和,天祚,昇元,保大,中興,交泰,天宝,宝大,宝正,龍啓,永和,通文,永隆,天徳,乾亨,白龍,大有,光天,応乾,乾和,大宝,天復,武成,永平,通正,天漢,光天,乾徳,咸康,明徳,広政,乾祐,天会,広運,応天,同慶,天尊,中興,天興,天寿,建隆,乾徳,開宝,太平,興国,雍熙,端拱,淳化,至道,咸平,景徳,大中,祥符,天禧,乾興,天聖,明道,景祐,宝元,康定,慶暦,皇祐,至和,嘉祐,治平,熙寧,元豊,元祐,紹聖,元符,建中,靖国,崇寧,大観,政和,重和,宣和,靖康,建炎,紹興,隆興,乾道,淳熙,紹熙,慶元,嘉泰,開禧,嘉定,宝慶,紹定,端平,嘉熙,淳祐,宝祐,開慶,景定,咸淳,徳祐,景炎,祥興,神冊,天賛,天顕,会同,大同,天禄,応暦,保寧,乾亨,統和,開泰,太平,景福,重熙,清寧,咸雍,太康,太安,寿昌,乾統,天慶,保大,甘露,建福,徳興,神暦,延慶,康国,咸清,紹興,崇福,天禧,収国,天輔,天会,天眷,皇統,天徳,貞元,正隆,大定,明昌,承安,泰和,大安,崇慶,至寧,貞祐,興定,元光,正大,開興,天興,顕道,開運,広運,大慶,天授,礼法,延祚,延嗣,寧国,天祐,垂聖,福聖,承道,奲都,拱化,乾道,天賜,礼盛,国慶,大安,天安,礼定,天儀,治平,天祐,民安,永安,貞観,雍寧,元徳,正徳,大徳,大慶,人慶,天盛,乾祐,天慶,応天,皇建,光定,乾定,宝義,中統,至元,元貞,大徳,至大,皇慶,延祐,至治,泰定,致和,天順,天暦,至順,元統,至元,至正,至正,宣光,天元,洪武,建文,永楽,洪熙,宣徳,正統,景泰,天順,成化,弘治,正徳,嘉靖,隆慶,万暦,泰昌,天啓,崇禎,弘光,隆武,紹武,永暦,天命,天聡,崇徳,順治,康熙,雍正,乾隆,嘉慶,道光,咸豊,祺祥,同治,光緒,保慶,宣統,民国,洪憲,大同,康徳,"
#korea
g3="永楽,延寿,建元,開国,大昌,鴻済,建福,仁平,太和,慶雲,正開,武泰,聖冊,水徳,万歳,政開,天授,光徳,天開,開国,建陽,光武,光武,隆熙,大韓,民国,"
#vietnam
g4="天徳,太平,天福,興統,応天,景瑞,順天,天成,通瑞,乾符,有道,明道,天感,聖武,崇興,大宝,龍瑞,太平,彰聖,嘉慶,龍彰,天嗣,天貺,宝象,神武,太寧,英武,昭勝,広祐,会豊,龍符,会祥,大慶,天符,睿武,天符,慶寿,天順,天彰,宝嗣,紹明,大定,政隆,宝応,天感,至宝,貞符,天資,嘉瑞,天嘉,宝祐,治平,龍応,建嘉,天彰,有道,建中,天応,政平,元豊,紹隆,宝符,紹宝,重興,興隆,大慶,開泰,開祐,紹豊,大治,大定,紹慶,隆慶,昌符,光泰,建新,聖元,紹成,開大,興慶,重光,天慶,順天,紹平,大宝,大和,延寧,天興,光順,洪徳,景統,泰貞,端慶,洪順,光紹,統元,明徳,大正,広和,永定,景暦,光宝,淳福,崇康,延成,端泰,興治,洪寧,武安,宝定,康佑,乾統,隆泰,順徳,元和,順平,天祐,正治,洪福,嘉泰,光興,慎徳,弘定,永祚,徳隆,陽和,福泰,慶徳,盛徳,永寿,万慶,景治,陽徳,徳元,永治,正和,永盛,保泰,永慶,龍徳,永佑,景興,昭統,泰徳,光中,景盛,宝興,嘉隆,明命,紹治,嗣徳,協和,建福,咸宜,同慶,成泰,維新,啓定,保大,"
#taiwan
g5="永暦,康熙,永和,雍正,乾隆,順天,天運,嘉慶,光明,道光,咸豊,祺祥,同治,光緒,永清,大靖,民国,"
#east-asia
g=g1*g2*g3*g4*g5
println((div(length(g1),3),div(length(g2),3),div(length(g3),3),div(length(g4),3),div(length(g5),3)))
m=div(length(g),3)
#
# making array of tupples with numbered era names of east asia
#
c=fill((" ",0),m)
for i=1:7:m*7
c[div(i,7)+1]=(g[i]*g[i+3],div(i,7)+1)
end
#
# the first and the second characters of the Japanese Gengo
#
n=div(length(g1),3)
a=fill(' ',n)
b=fill(' ',n)
for i=1:7:n*7
a[div(i,7)+1]=g1[i]
b[div(i,7)+1]=g1[i+3]
end
p=[ ]
q=[ ]
mp=0
mq=0
gfrq(a,n,p,mp)
gfrq(b,n,q,mq)
mp=length(p)
mq=length(q)
#
r=[ ]
for j in 1:mp
for k in 1:mq
if p[j][1]!=q[k][1]
push!(r,[string(p[j][1],q[k][1]),p[j][2]*q[k][2]/(mp*mq)])
end
end
end
mr=length(r)
#
# exclude previousely used names of era (Japan China Korea Vietnum Taiwan)
# & their exchanges
#
for i in 1:mr, j in 1:m
if r[i][1] == c[j][1] || r[i][1] == exch(c[j][1])
r[i][2]=0
end
#
# exlude recently used characters & complex characters & MTSH characters
#
if occursin(r[i][1][1],"慶応明治大正昭和平成寛嘉亀徳康暦承雉吉祥福禄寿喜建延享亨興景観祚感乾老雲霊国勝武衡護養銅授泰神政天中長禎貞宝保仁")==true ||
occursin(r[i][1][4],"慶応明治大正昭和平成寛嘉亀徳康暦承雉吉祥福禄寿喜建延享亨興景観祚感乾老雲霊国勝武衡護養銅授泰神政天中長禎貞宝保仁")==true
r[i][2]=0
end
end
#
sort!(r, by=x->x[2],rev=true)
for k in 1:50
println((r[k][1],r[k][2]))
end
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
(1) これまでの日本の元号に登場した漢字の範囲で考える。ただし１文字目と２文字目は区別して考える（４文字元号は２文字元号に読み替える）。
(2) １文字目の漢字と２文字目の漢字のそれぞれの出現確率の積で順序づける。
(3) 既存のアジア（日本，中国，朝鮮，ベトナム，台湾）の元号とそれを入れ替えたものに一致するものは排除する（例えば，「安久」は排除される）。
(4) 慶応，明治，大正，昭和，平成の各漢字は排除する。
(5) 画数の多い漢字，MTSHで始まる漢字などを排除する（「など」のところにプログラマの主観が入る）
昭和とか平成の場合を考えると，このアルゴリズムはよくないことがわかっているが，目をつぶることにしよう（そりゃ，「苦節」とか「余命」とか「結婚」とか「懲役」とか「西暦」とか「卍」とか「タピオカ」のほうがおもしろいに決まっている）。とにかくこの条件で結果抽出されたものが下記のとおり。
(250, 839, 26, 159, 17)
("弘久", 0.020164609053497942)
("永化", 0.019753086419753086)
("元化", 0.018518518518518517)
("万安", 0.016049382716049384)
("万元", 0.014814814814814815)
("文弘", 0.01316872427983539)
("弘文", 0.011522633744855968)
("万久", 0.008641975308641974)
("文万", 0.006584362139917695)
("文同", 0.006584362139917695)
("文字", 0.006584362139917695)
("永同", 0.006584362139917695)
("永字", 0.006584362139917695)
("元万", 0.006172839506172839)
("元同", 0.006172839506172839)
("元字", 0.006172839506172839)
("安化", 0.006172839506172839)
("斉永", 0.005761316872427984)
("朱永", 0.005761316872427984)
("白永", 0.005761316872427984)
("至永", 0.005761316872427984)
("斉安", 0.005349794238683128)
("昌安", 0.005349794238683128)
("朱安", 0.005349794238683128)
("白安", 0.005349794238683128)
("至安", 0.005349794238683128)
("万文", 0.0049382716049382715)
("斉元", 0.0049382716049382715)
("昌元", 0.0049382716049382715)
("朱元", 0.0049382716049382715)
("白元", 0.0049382716049382715)
("万化", 0.003703703703703704)
("弘万", 0.002880658436213992)
("弘同", 0.002880658436213992)
("弘字", 0.002880658436213992)
("斉久", 0.002880658436213992)
("昌久", 0.002880658436213992)
("朱久", 0.002880658436213992)
("白久", 0.002880658436213992)
("至久", 0.002880658436213992)
("万弘", 0.0024691358024691358)
("久化", 0.0024691358024691358)
("安万", 0.00205761316872428)
("安同", 0.00205761316872428)
("安字", 0.00205761316872428)
("久弘", 0.0016460905349794238)
("斉文", 0.0016460905349794238)
("昌文", 0.0016460905349794238)
("朱文", 0.0016460905349794238)
("白文", 0.0016460905349794238)
Juliaが使える方は，これで適当に遊んでください。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
function exch(knj)
#
# knj input kanji string with 2 char
# return exchanged string with 2 char
return knj[4]*knj[1]
end
function gfrq(c,n,r,m)
#
# c input array of kanji char
# n size of c
# r output array of sorted (kanji, count)
# m size of r
#
sort!(c)
k=0
for j=2:n
k=k+1
if c[j]!=c[j-1]
push!(r,(c[j-1],k))
k=0
end
end
sort!(r, by=x->x[2], rev=true)
m=length(r)
# println((m,r))
end
#japan
g1="大化,白雉,朱鳥,大宝,慶雲,和銅,霊亀,養老,神亀,天平,感宝,勝宝,宝字,神護,景雲,宝亀,天応,延暦,大同,弘仁,天長,承和,嘉祥,仁寿,斉衡,天安,貞観,元慶,仁和,寛平,昌泰,延喜,延長,承平,天慶,天暦,天徳,応和,康保,安和,天禄,天延,貞元,天元,永観,寛和,永延,永祚,正暦,長徳,長保,寛弘,長和,寛仁,治安,万寿,長元,長暦,長久,寛徳,永承,天喜,康平,治暦,延久,承保,承暦,永保,応徳,寛治,嘉保,永長,承徳,康和,長治,嘉承,天仁,天永,永久,元永,保安,天治,大治,天承,長承,保延,永治,康治,天養,久安,仁平,久寿,保元,平治,永暦,応保,長寛,永万,仁安,嘉応,承安,安元,治承,養和,寿永,元暦,文治,建久,正治,建仁,元久,建永,承元,建暦,建保,承久,承応,元仁,嘉禄,安貞,寛喜,貞永,天福,文暦,嘉禎,暦仁,延応,仁治,寛元,宝治,建長,康元,正嘉,正元,文応,弘長,文永,建治,弘安,文永,建治,弘安,正応,永仁,正安,乾元,嘉元,徳治,延慶,応長,正和,文保,元応,元亨,正中,嘉暦,元徳,元弘,正慶,建武,延元,興国,正平,建徳,文中,天授,弘和,元中,暦応,康永,貞和,観応,文和,延文,康安,貞治,応安,永和,康暦,永徳,至徳,嘉慶,康応,明徳,応永,正長,永享,嘉吉,文安,宝徳,享徳,康正,長禄,寛正,文正,応仁,文明,長享,延徳,明応,文亀,永正,大永,享禄,天文,弘治,永禄,元亀,天正,文禄,慶長,元和,寛永,正保,慶安,承応,明暦,万治,寛文,延宝,天和,貞享,元禄,宝永,正徳,享保,元文,寛保,延享,寛延,宝暦,明和,安永,天明,寛政,享和,文化,文政,天保,弘化,嘉永,安政,万延,文久,元治,慶応,明治,大正,昭和,平成,"
#china
g2="建元,元光,元朔,元狩,元鼎,元封,太初,天漢,太始,征和,後元,始元,元鳳,元平,本始,地節,元康,神爵,五鳳,甘露,黄龍,初元,永光,建昭,竟寧,建始,河平,陽朔,鴻嘉,永始,元延,綏和,建平,太初,元将,元寿,元始,居摂,初始,始建,天鳳,地皇,更始,建武,建武,中元,永平,建初,元和,章和,永元,元興,延平,永初,元初,永寧,建光,延光,永建,陽嘉,永和,漢安,建康,永憙,本初,建和,和平,元嘉,永興,永寿,延熹,永康,建寧,熹平,光和,中平,光熹,昭寧,永漢,初平,興平,建安,延康,黄初,太和,青龍,景初,正始,嘉平,正元,甘露,景元,咸熙,章武,建興,延熙,景耀,炎興,黄武,黄龍,嘉禾,赤烏,太元,神鳳,建興,五鳳,太平,永安,元興,甘露,宝鼎,建衡,鳳凰,天冊,天璽,天紀,泰始,咸寧,太康,太熙,永熙,永平,元康,永康,永寧,太安,永安,建武,永興,光熙,永嘉,建興,建武,大興,永昌,太寧,咸和,咸康,建元,永和,升平,隆和,興寧,太和,咸安,寧康,太元,隆安,大亨,元興,義熙,元熙,桓玄,天康,永鳳,河瑞,光興,嘉平,建元,麟嘉,漢昌,光初,太和,建平,延熙,建武,太寧,青龍,永寧,永興,元璽,光寿,建熙,建初,建興,晏平,玉衡,玉恒,漢興,太和,嘉寧,建興,和平,升平,皇始,寿光,永興,甘露,建元,太安,太初,延初,元光,燕元,建興,永康,建始,延平,青龍,建平,長楽,光始,建始,燕興,更始,昌平,建明,建平,建武,中興,建平,太上,正始,太平,太興,白雀,建初,皇初,弘始,永和,建義,太初,更始,永康,建弘,永弘,龍昇,鳳翔,昌武,真興,承光,勝光,太安,麟嘉,龍飛,承康,咸寧,神鼎,太初,建和,弘昌,嘉平,庚子,建初,嘉興,永建,神璽,天璽,永安,玄始,承玄,義和,承和,承平,真興,承陽,縁禾,太縁,永初,景平,元嘉,太初,孝建,大明,永光,景和,泰始,泰豫,元徽,昇明,建元,永明,隆昌,延興,建武,永泰,永元,中興,天監,普通,大通,大通,大同,大同,太清,大宝,天正,太始,承聖,天成,紹泰,太平,天啓,大定,天保,広運,永定,天嘉,天康,光大,太建,至徳,禎明,登国,皇始,天興,天賜,永興,神瑞,泰常,始光,神䴥,延和,太延,太平,真君,正平,承平,興安,興光,太安,和平,天安,皇興,延興,承明,太和,景明,正始,永平,延昌,熙平,神亀,正光,孝昌,武泰,建義,永安,建明,普泰,中興,太昌,永興,永熙,天平,元象,興和,武定,大統,天保,乾明,皇建,太寧,河清,天統,武平,隆化,徳昌,承光,武成,保定,天和,建徳,宣政,大成,大象,大定,開皇,仁寿,大業,義寧,皇泰,武徳,貞観,永徽,顕慶,龍朔,麟徳,乾封,総章,咸亨,上元,儀鳳,調露,永隆,開耀,永淳,弘道,嗣聖,文明,光宅,垂拱,永昌,載初,神龍,景龍,唐隆,景雲,太極,延和,先天,開元,天宝,至徳,乾元,上元,宝応,広徳,永泰,大暦,建中,興元,貞元,永貞,元和,長慶,宝暦,大和,開成,会昌,大中,咸通,乾符,広明,中和,光啓,文徳,龍紀,大順,景福,乾寧,光化,天復,天祐,天授,如意,長寿,延載,証聖,天冊,万歳,万歳,登封,万歳,通天,神功,聖暦,久視,大足,長安,仁安,大興,宝暦,中興,正暦,永徳,朱雀,太始,建興,咸和,建初,承平,義熙,甘露,章和,永平,和平,建昌,延昌,延和,義和,重光,延寿,賛普,長寿,見龍,上元,元封,応道,龍興,全義,大豊,保和,天啓,建極,法尭,貞明,大同,嵯耶,中興,安国,始元,天瑞,景星,安和,貞祐,初歴,孝治,天応,尊聖,興聖,大明,鼎新,光聖,文徳,神武,文経,至治,明徳,広徳,順徳,明政,広明,明応,明統,明聖,明治,明啓,乾興,明通,正治,聖明,天明,保安,正安,正徳,保徳,明侯,上徳,広安,上明,保定,建安,天祐,上治,天授,開明,天政,文安,日新,文治,永嘉,保天,広運,永貞,大宝,龍興,盛明,建徳,利貞,盛徳,嘉会,元亨,安定,鳳歴,元寿,天開,天輔,仁寿,道隆,利正,興正,天定,彝泰,開平,乾化,鳳歴,貞明,龍徳,同光,天成,長興,応順,清泰,天福,開運,天福,乾祐,広順,顕徳,天復,天祐,武義,順義,乾貞,大和,天祚,昇元,保大,中興,交泰,天宝,宝大,宝正,龍啓,永和,通文,永隆,天徳,乾亨,白龍,大有,光天,応乾,乾和,大宝,天復,武成,永平,通正,天漢,光天,乾徳,咸康,明徳,広政,乾祐,天会,広運,応天,同慶,天尊,中興,天興,天寿,建隆,乾徳,開宝,太平,興国,雍熙,端拱,淳化,至道,咸平,景徳,大中,祥符,天禧,乾興,天聖,明道,景祐,宝元,康定,慶暦,皇祐,至和,嘉祐,治平,熙寧,元豊,元祐,紹聖,元符,建中,靖国,崇寧,大観,政和,重和,宣和,靖康,建炎,紹興,隆興,乾道,淳熙,紹熙,慶元,嘉泰,開禧,嘉定,宝慶,紹定,端平,嘉熙,淳祐,宝祐,開慶,景定,咸淳,徳祐,景炎,祥興,神冊,天賛,天顕,会同,大同,天禄,応暦,保寧,乾亨,統和,開泰,太平,景福,重熙,清寧,咸雍,太康,太安,寿昌,乾統,天慶,保大,甘露,建福,徳興,神暦,延慶,康国,咸清,紹興,崇福,天禧,収国,天輔,天会,天眷,皇統,天徳,貞元,正隆,大定,明昌,承安,泰和,大安,崇慶,至寧,貞祐,興定,元光,正大,開興,天興,顕道,開運,広運,大慶,天授,礼法,延祚,延嗣,寧国,天祐,垂聖,福聖,承道,奲都,拱化,乾道,天賜,礼盛,国慶,大安,天安,礼定,天儀,治平,天祐,民安,永安,貞観,雍寧,元徳,正徳,大徳,大慶,人慶,天盛,乾祐,天慶,応天,皇建,光定,乾定,宝義,中統,至元,元貞,大徳,至大,皇慶,延祐,至治,泰定,致和,天順,天暦,至順,元統,至元,至正,至正,宣光,天元,洪武,建文,永楽,洪熙,宣徳,正統,景泰,天順,成化,弘治,正徳,嘉靖,隆慶,万暦,泰昌,天啓,崇禎,弘光,隆武,紹武,永暦,天命,天聡,崇徳,順治,康熙,雍正,乾隆,嘉慶,道光,咸豊,祺祥,同治,光緒,保慶,宣統,民国,洪憲,大同,康徳,"
#korea
g3="永楽,延寿,建元,開国,大昌,鴻済,建福,仁平,太和,慶雲,正開,武泰,聖冊,水徳,万歳,政開,天授,光徳,天開,開国,建陽,光武,光武,隆熙,大韓,民国,"
#vietnam
g4="天徳,太平,天福,興統,応天,景瑞,順天,天成,通瑞,乾符,有道,明道,天感,聖武,崇興,大宝,龍瑞,太平,彰聖,嘉慶,龍彰,天嗣,天貺,宝象,神武,太寧,英武,昭勝,広祐,会豊,龍符,会祥,大慶,天符,睿武,天符,慶寿,天順,天彰,宝嗣,紹明,大定,政隆,宝応,天感,至宝,貞符,天資,嘉瑞,天嘉,宝祐,治平,龍応,建嘉,天彰,有道,建中,天応,政平,元豊,紹隆,宝符,紹宝,重興,興隆,大慶,開泰,開祐,紹豊,大治,大定,紹慶,隆慶,昌符,光泰,建新,聖元,紹成,開大,興慶,重光,天慶,順天,紹平,大宝,大和,延寧,天興,光順,洪徳,景統,泰貞,端慶,洪順,光紹,統元,明徳,大正,広和,永定,景暦,光宝,淳福,崇康,延成,端泰,興治,洪寧,武安,宝定,康佑,乾統,隆泰,順徳,元和,順平,天祐,正治,洪福,嘉泰,光興,慎徳,弘定,永祚,徳隆,陽和,福泰,慶徳,盛徳,永寿,万慶,景治,陽徳,徳元,永治,正和,永盛,保泰,永慶,龍徳,永佑,景興,昭統,泰徳,光中,景盛,宝興,嘉隆,明命,紹治,嗣徳,協和,建福,咸宜,同慶,成泰,維新,啓定,保大,"
#taiwan
g5="永暦,康熙,永和,雍正,乾隆,順天,天運,嘉慶,光明,道光,咸豊,祺祥,同治,光緒,永清,大靖,民国,"
#east-asia
g=g1*g2*g3*g4*g5
println((div(length(g1),3),div(length(g2),3),div(length(g3),3),div(length(g4),3),div(length(g5),3)))
m=div(length(g),3)
#
# making array of tupples with numbered era names of east asia
#
c=fill((" ",0),m)
for i=1:7:m*7
c[div(i,7)+1]=(g[i]*g[i+3],div(i,7)+1)
end
#
# the first and the second characters of the Japanese Gengo
#
n=div(length(g1),3)
a=fill(' ',n)
b=fill(' ',n)
for i=1:7:n*7
a[div(i,7)+1]=g1[i]
b[div(i,7)+1]=g1[i+3]
end
p=[ ]
q=[ ]
mp=0
mq=0
gfrq(a,n,p,mp)
gfrq(b,n,q,mq)
mp=length(p)
mq=length(q)
#
r=[ ]
for j in 1:mp
for k in 1:mq
if p[j][1]!=q[k][1]
push!(r,[string(p[j][1],q[k][1]),p[j][2]*q[k][2]/(mp*mq)])
end
end
end
mr=length(r)
#
# exclude previousely used names of era (Japan China Korea Vietnum Taiwan)
# & their exchanges
#
for i in 1:mr, j in 1:m
if r[i][1] == c[j][1] || r[i][1] == exch(c[j][1])
r[i][2]=0
end
#
# exlude recently used characters & complex characters & MTSH characters
#
if occursin(r[i][1][1],"慶応明治大正昭和平成寛嘉亀徳康暦承雉吉祥福禄寿喜建延享亨興景観祚感乾老雲霊国勝武衡護養銅授泰神政天中長禎貞宝保仁")==true ||
occursin(r[i][1][4],"慶応明治大正昭和平成寛嘉亀徳康暦承雉吉祥福禄寿喜建延享亨興景観祚感乾老雲霊国勝武衡護養銅授泰神政天中長禎貞宝保仁")==true
r[i][2]=0
end
end
#
sort!(r, by=x->x[2],rev=true)
for k in 1:50
println((r[k][1],r[k][2]))
end
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2019年3月29日金曜日
2019年3月28日木曜日
数学パワーが...orz
AI人材25万人育成の背景にあるのが，これなのか。
経済産業省「理数系人材の産業界での活躍に向けた意見交換会の報告書『数理資本主義の時代～数学パワーが世界を変える～』を取りまとめました」。理数系人材の産業界での活躍に向けた意見交換会の報告書が，「数理資本主義の時代〜数学パワーが世界を変える〜」として2019年3月26日に公表されている。
※1ここで言う「数学」は、純粋数学、応用数学、統計学、確率論、さらには数学的な表現を必要とする量子論、素粒子物理学、宇宙物理学なども含む広範な概念。
という注釈があるのだけれど，物理屋さんはよんでいないのかな。数学者の加藤文元さんは，参考人としてよばれて「NPO法人数理の翼」の話などしているようだけど。
キャリアパスの流れを変えるには，数理資本主義というキーワードも必要だったのだろうが，なんだか気持ち悪いし実効性もないような気がする。この際，グローバル化しきれない経済産業省と右傾化した文部科学省を統合して，経済教育開発省とすると，発展途上国に向かって逆進化の道を進んでいる我が国の未来にふさわしいのではないだろうか。なんで定年退職の最後から２番目の出勤日にこんなくだらない記事を書いているのだろうか。反省することしきりなり。
経済産業省「理数系人材の産業界での活躍に向けた意見交換会の報告書『数理資本主義の時代～数学パワーが世界を変える～』を取りまとめました」。理数系人材の産業界での活躍に向けた意見交換会の報告書が，「数理資本主義の時代〜数学パワーが世界を変える〜」として2019年3月26日に公表されている。
※1ここで言う「数学」は、純粋数学、応用数学、統計学、確率論、さらには数学的な表現を必要とする量子論、素粒子物理学、宇宙物理学なども含む広範な概念。
という注釈があるのだけれど，物理屋さんはよんでいないのかな。数学者の加藤文元さんは，参考人としてよばれて「NPO法人数理の翼」の話などしているようだけど。
キャリアパスの流れを変えるには，数理資本主義というキーワードも必要だったのだろうが，なんだか気持ち悪いし実効性もないような気がする。この際，グローバル化しきれない経済産業省と右傾化した文部科学省を統合して，経済教育開発省とすると，発展途上国に向かって逆進化の道を進んでいる我が国の未来にふさわしいのではないだろうか。なんで定年退職の最後から２番目の出勤日にこんなくだらない記事を書いているのだろうか。反省することしきりなり。
2019年3月27日水曜日
AI人材年25万人育成...orz
うーん，政府の中枢が経産省に乗っ取られると，なかなか面白い政策が次々と実行され，我が国はますますボロボロになっていくような気がガガ。Radio Ga Ga。
日本経済新聞の3月27日の夕刊によりますと，内閣府の統合イノベーション戦略推進会議が29日に公表する「AI戦略」の概要に，人工知能（AI）を使いこなす人材を25万人育てる心目標を掲げるとのこと。全大学生がAIの初等教育を受けるように大学に要請し，社会人向けの専門課程も設置するそうだ。小学校の教育課程へのプログラミング教育導入と見事に対をなしているのだった（経産省起点で文科省に手を突っ込む作戦か）。そのAI戦略の概要は次の通り。
理念：
「数理・データサイエンス・ＡＩ」はすべての国民にとって
「読み・書き・そろばん」と並ぶスキルに
主な取り組み：
・全ての大学生・高専生に初級レベルのＡＩ教育（この対象は６０万人である）
・ＡＩと専門分野のダブルメジャーを促進
・大学に社会人専門コースを設置し，学び直しを支援
かつて，ソフトウェア技術者の不足をソフトウェア危機としてあおり，情報系のバブルをもたらしたまま，後始末におわれた経験を生かしているようには見えないのだけれども。いったいどうなっているのだろうか。必要な人材数の数量的評価も杜撰（ニュースの社会科学的な裏側）なのかもしれない。
日本経済新聞の3月27日の夕刊によりますと，内閣府の統合イノベーション戦略推進会議が29日に公表する「AI戦略」の概要に，人工知能（AI）を使いこなす人材を25万人育てる心目標を掲げるとのこと。全大学生がAIの初等教育を受けるように大学に要請し，社会人向けの専門課程も設置するそうだ。小学校の教育課程へのプログラミング教育導入と見事に対をなしているのだった（経産省起点で文科省に手を突っ込む作戦か）。そのAI戦略の概要は次の通り。
理念：
「数理・データサイエンス・ＡＩ」はすべての国民にとって
「読み・書き・そろばん」と並ぶスキルに
主な取り組み：
・全ての大学生・高専生に初級レベルのＡＩ教育（この対象は６０万人である）
・ＡＩと専門分野のダブルメジャーを促進
・大学に社会人専門コースを設置し，学び直しを支援
かつて，ソフトウェア技術者の不足をソフトウェア危機としてあおり，情報系のバブルをもたらしたまま，後始末におわれた経験を生かしているようには見えないのだけれども。いったいどうなっているのだろうか。必要な人材数の数量的評価も杜撰（ニュースの社会科学的な裏側）なのかもしれない。
2019年3月26日火曜日
Juliaで標準入出力
今朝の日経新聞に，全国統一プログラミング王決定戦本戦の記事が載っていた。新聞に示された例題を見るとそんなに難しそうではない。自分でもできそうなレベルだ。さっそくウェブサイトを訪れると，昨年のマストドンブームでお見かけしたC++の江添亮さんに影響されてC++を勉強するためにバーチャルに通った，会津大学のAOJオンラインプログラミングチャレンジと同じようなフレーバーのサイトだった。
さっそくJuliaで挑戦してみようとしたところ，C++とかPythonはありそうだが，そもそもJuliaなんて受け付けられないのではないか。それはそれとして，Juliaで標準入力をどうするのかは今まで考えてこなかった。この機会に調べてみよう。
(1) Jupyter環境の場合
なんのことはない，readline()とすると stdin> として標準入力が要求された。めでたしめでたし。
(2) REPL環境の場合
標準入力のプロンプトは自分で設定しないと出ないのだが，入力待ちの環境中で入力と改行を繰り返せばよいだけである。入力データの型はStringになっていた。例えばMacOSXの Julia REPL環境でプログラム開発ループを実現するには，まずbashのコマンドラインでJuliaを起動し，REPL環境が立ち上がれば，julia>プロンプトに対し，
edit("/Users/xxxx/yyyy/test.jl")
include("/Users/xxxx/yyyy/test.jl")
を繰り返して，編集とプログラムの読み込みを実行すればよかった。
ここで，ユーザxxxxの作業ディレクトリyyyyにJuliaのテストプログラム test.jl を置いたものと仮定している。editでは標準に設定されているエディタが起動した（ここがまだよくわかっていない）。したがって，これ以降はエディタで編集し，includeを繰り返えせばよい。
［注］Macにインストールしたjuliaは，/Applications/Julia-1.1.app/Contents/Resources/julia/bin/julia にあった。bashでパスを有効にするには，.bashrcにパスを記述して，.bash_profileにでbash起動時に.bashrcを読み込むようにする必要がある。
cat .bash_profile
if [ -f ~/.bashrc ] ; then
. ~/.bashrc
fi
cat .bashrc
export PATH=$PATH:/Applications/Julia-1.1.app/Contents/Resources/julia/bin
(3) コマンドライン環境の場合
この場合はもう少し勉強してコンパイルできるようになってから考えよう。
P. S. At Coder のプログラム言語，Juliaは隅の方にあったけど，0.5.0なのだった...orz
参考：Juliaで競技プログラミング(AtCoder)やるときのTips(Julia 0.5.0用)
さっそくJuliaで挑戦してみようとしたところ，C++とかPythonはありそうだが，そもそもJuliaなんて受け付けられないのではないか。それはそれとして，Juliaで標準入力をどうするのかは今まで考えてこなかった。この機会に調べてみよう。
(1) Jupyter環境の場合
なんのことはない，readline()とすると stdin> として標準入力が要求された。めでたしめでたし。
(2) REPL環境の場合
標準入力のプロンプトは自分で設定しないと出ないのだが，入力待ちの環境中で入力と改行を繰り返せばよいだけである。入力データの型はStringになっていた。例えばMacOSXの Julia REPL環境でプログラム開発ループを実現するには，まずbashのコマンドラインでJuliaを起動し，REPL環境が立ち上がれば，julia>プロンプトに対し，
edit("/Users/xxxx/yyyy/test.jl")
include("/Users/xxxx/yyyy/test.jl")
を繰り返して，編集とプログラムの読み込みを実行すればよかった。
ここで，ユーザxxxxの作業ディレクトリyyyyにJuliaのテストプログラム test.jl を置いたものと仮定している。editでは標準に設定されているエディタが起動した（ここがまだよくわかっていない）。したがって，これ以降はエディタで編集し，includeを繰り返えせばよい。
［注］Macにインストールしたjuliaは，/Applications/Julia-1.1.app/Contents/Resources/julia/bin/julia にあった。bashでパスを有効にするには，.bashrcにパスを記述して，.bash_profileにでbash起動時に.bashrcを読み込むようにする必要がある。
cat .bash_profile
if [ -f ~/.bashrc ] ; then
. ~/.bashrc
fi
cat .bashrc
export PATH=$PATH:/Applications/Julia-1.1.app/Contents/Resources/julia/bin
(3) コマンドライン環境の場合
この場合はもう少し勉強してコンパイルできるようになってから考えよう。
P. S. At Coder のプログラム言語，Juliaは隅の方にあったけど，0.5.0なのだった...orz
参考：Juliaで競技プログラミング(AtCoder)やるときのTips(Julia 0.5.0用)
2019年3月25日月曜日
JupyterLab
「JupyterLabのすすめ」という記事があったので，JupyterLabをさっそくインストールしてみた。
退職後の４月から非常勤講師としてこれまでと同じような授業を担当する。とりあえず，現在使用しているコンピュータ（MacProやMacBookProなど）を明け渡さなければならないので，MacBookAirに新しい環境を構築して移動作業をしている。その結果，現在の自分の環境は，次のようになっている。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
※mbaの設定
(1) homebrewの導入
xcode-select —install （コマンドライン開発ツールの導入）
sudo xcodebuild -license
/usr/bin/ruby -e "＄(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/master/install)"
(2) homebrewの練習と基本環境の設定
brew update
brew upgrade
brew doctor
brew install wget
brew install nkf
brew cask reinstall xquartz
brew install lynx
brew install w3m
brew install xxxx （xxxxをインストール）以下適当に繰り返した。
(3) python3環境（jupyterを含む）の設定
brew install python3
python3 -V
Python 3.7.2
pip3 -V
pip 19.0.2 from /usr/local/lib/python3.7/site-packages/pip (python 3.7)
python -V
Python 2.7.10
・これで，python3がインストールできた。
which python
/usr/bin/python
which python3
/usr/local/bin/python3
・pip3のアップデート方法
pip3 install --upgrade setuptools
Requirement already up-to-date:...
pip3 install --upgrade pip
Collecting pip
...
Successfully installed pip-19.0.3
・環境設定してその中で作業
python3 -m venv env
. env/bin/activate
(env) … python -V
Python 3.7.2
(env) … pip -V
pip 18.1 from /Users/admin/env/lib/python3.7/site-packages/pip (python 3.7)
(env) ... pip install --upgrade pip
(env) ... pip install --upgrade setuptools
(env) ... pip install numpy
(env) ... pip install scipy
(env) ... pip install scikit-learn
(env) ... pip install matplotlib
(env) ... pip install Pillow
(env) ... pip install sympy
(env) ... pip install jupyter
(env) ... deactivate
(4) julia環境の設定
https://julialang.org/downloads/
からダウンロードすればコマンドラインで実行可能に
echo 'export PATH=$PATH:/Applications/
Julia-1.1.app/Contents/Resources/julia/bin' >> ~/.bashrc
source ~/.bashrc
julia
julia> using Pkg
julia> Pkg.add(“IJulia”)
export PATH=＄PATH:/Users/admin/env/bin
jupyter notebook
これで，jupyter が起動して，julia 1.1 が使えるようになった。
versioninfo()
using Pkg
Pkg.update()
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
(5) jupyterlabの導入と起動は簡単であった。
pip3 install jupyterlab
jupyter lab
（付）JupyterLab環境での設定
Settings -> Advanced Settings Editor -> Notebook
codeCellConfigのlineNumbersをTrueにする。
Settings -> Advanced Settings Editor -> Extension Manager -> "enabled": trueにする。
退職後の４月から非常勤講師としてこれまでと同じような授業を担当する。とりあえず，現在使用しているコンピュータ（MacProやMacBookProなど）を明け渡さなければならないので，MacBookAirに新しい環境を構築して移動作業をしている。その結果，現在の自分の環境は，次のようになっている。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
※mbaの設定
(1) homebrewの導入
xcode-select —install （コマンドライン開発ツールの導入）
sudo xcodebuild -license
/usr/bin/ruby -e "＄(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/master/install)"
(2) homebrewの練習と基本環境の設定
brew update
brew upgrade
brew doctor
brew install wget
brew install nkf
brew cask reinstall xquartz
brew install lynx
brew install w3m
brew install xxxx （xxxxをインストール）以下適当に繰り返した。
(3) python3環境（jupyterを含む）の設定
brew install python3
python3 -V
Python 3.7.2
pip3 -V
pip 19.0.2 from /usr/local/lib/python3.7/site-packages/pip (python 3.7)
python -V
Python 2.7.10
・これで，python3がインストールできた。
which python
/usr/bin/python
which python3
/usr/local/bin/python3
・pip3のアップデート方法
pip3 install --upgrade setuptools
Requirement already up-to-date:...
pip3 install --upgrade pip
Collecting pip
...
Successfully installed pip-19.0.3
・環境設定してその中で作業
python3 -m venv env
. env/bin/activate
(env) … python -V
Python 3.7.2
(env) … pip -V
pip 18.1 from /Users/admin/env/lib/python3.7/site-packages/pip (python 3.7)
(env) ... pip install --upgrade pip
(env) ... pip install --upgrade setuptools
(env) ... pip install numpy
(env) ... pip install scipy
(env) ... pip install scikit-learn
(env) ... pip install matplotlib
(env) ... pip install Pillow
(env) ... pip install sympy
(env) ... pip install jupyter
(env) ... deactivate
(4) julia環境の設定
https://julialang.org/downloads/
からダウンロードすればコマンドラインで実行可能に
echo 'export PATH=$PATH:/Applications/
Julia-1.1.app/Contents/Resources/julia/bin' >> ~/.bashrc
source ~/.bashrc
julia
julia> using Pkg
julia> Pkg.add(“IJulia”)
export PATH=＄PATH:/Users/admin/env/bin
jupyter notebook
これで，jupyter が起動して，julia 1.1 が使えるようになった。
versioninfo()
using Pkg
Pkg.update()
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
(5) jupyterlabの導入と起動は簡単であった。
pip3 install jupyterlab
jupyter lab
（付）JupyterLab環境での設定
Settings -> Advanced Settings Editor -> Notebook
codeCellConfigのlineNumbersをTrueにする。
Settings -> Advanced Settings Editor -> Extension Manager -> "enabled": trueにする。
Mathematicaノートブックを長らく使ってきたので，Jupyter Notebook（この記事もWikipedia中文はあるが日本語はない・・・）はすんなりと馴染むことができた。JupyterLabはJupyter Notebookの後継で，今後はこれに移行するのだそうだ。まだ違いがよくわからないが，Jupyter Notebook のことだってそれほどわかっていたわけではない。ぼちぼちやっていくことにしよう。
2019年3月24日日曜日
金沢泉丘高等学校のプール
大乗寺丘陵公園から円光寺の方に下り，泉丘高校に向って歩いた。前身の一中創立70周年から80周年の間で，理数科ができて２年目の1969年に入学した。1984年には新校舎が落成しており，我々の時代の旧校舎はなくなっている。正門のあたりの古い木々や碑などはそのまま残っているが，校舎も制服も生徒も変わってしまった（昔よりみなさんよくできるのである）。グラウンド横の部室裏のプールは当時のままだと思ったが，沿革によると1992年にプールが完成とあるので同じ場所で新しくしたのだろう。校庭の周辺の桜の木とか，相撲場はどうだっただろうか。
小学校５年までは泳げなかったので，水泳の授業（小学校にはまだプールがなかった）では，金沢市営の25mプールに集められた数名の泳げない組で練習していた。というか，先生はあまり教えてくれないのである。同級生の畦地勉君に浮き方をおそわって，ようやく潜った状態でばた足で前進できるようになった。
中学校の体育では25mを泳ぐことが要求された。３年間でようやく25mをなんとか泳げたがほぼ息継ぎのないクロールだったかもしれない。背泳ぎも少しできるようになった。中学の体育の先生は親切だったので，泳ぎの苦手な子に飛び込みを教えてくださった。プールサイドに腰掛けた状態から始まって，徐々に慣らして行き，最後にはスタート台から飛び込むことができるようになった。
高校の体育では200ｍを泳ぐことが要求された。これはちょっと難しかったので，夏休み前に10名程度泳ぎの苦手な連中がプールで特別指導された。というかここでも指導された記憶はあまりないのだが，梅雨も明けきらないやや肌寒いプールでばしゃばしゃしていた。授業本番では背泳ぎで200mをなんとか泳ぐことができた。これ以降は夏休みに市営プールに一人で通って1000mを平泳ぎで泳ぐことができるようになった。
写真 金沢泉丘高等学校のプール（2019.3.22）
中学校の体育では25mを泳ぐことが要求された。３年間でようやく25mをなんとか泳げたがほぼ息継ぎのないクロールだったかもしれない。背泳ぎも少しできるようになった。中学の体育の先生は親切だったので，泳ぎの苦手な子に飛び込みを教えてくださった。プールサイドに腰掛けた状態から始まって，徐々に慣らして行き，最後にはスタート台から飛び込むことができるようになった。
高校の体育では200ｍを泳ぐことが要求された。これはちょっと難しかったので，夏休み前に10名程度泳ぎの苦手な連中がプールで特別指導された。というかここでも指導された記憶はあまりないのだが，梅雨も明けきらないやや肌寒いプールでばしゃばしゃしていた。授業本番では背泳ぎで200mをなんとか泳ぐことができた。これ以降は夏休みに市営プールに一人で通って1000mを平泳ぎで泳ぐことができるようになった。
2019年3月23日土曜日
2019年3月22日金曜日
MathJaxで数式表示（２）
MathJaxの数式表示で行番号がでないのがどうなのかを調べて試しているうちに，調子が悪くなってしまった。古い設定がどこかに残っているためなのか，Safariで数式が小さく表示される現象が生じている。再起動で回復しても，そのページを閲覧した後にはどこでも発症してしまう。どうしたものか。FirefoxやChromeでは問題ない。
今のヘッダーの設定では，標準で連番あり，"¥begin{equation*}"はだめ，"¥[ ¥]"で行番号なし，"¥tag{n}" で番号指定可になっている（¥はバックスラッシュを表す）。
\begin{equation}
\begin{aligned}
\sqrt{x}+\cos{y}\\
\dfrac{\log(x)}{e^x}
\end{aligned}
\end{equation}
\[
\sqrt{\dfrac{z^2+y^2}{x^2}}
\]
\begin{equation}
y= a x^2 + b x + c \tag{4}
\end{equation}
\begin{equation}
m \ddot{\bm{r}} = \bm{F}
\end{equation}
これで"¥bm{}"も使えるようになったのは次のヘッダーである。
参考にしたのは以下のサイトである。
・Mathjaxを使ってBloggerで数式を書く
・おもちゃ箱：Bloggerで数式を表示する
・MathJaxの使い方（黒木玄）
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
<script async='async' src='https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.2/MathJax.js?config=TeX-MML-AM_CHTML' type='text/javascript'/>
<script type='text/x-mathjax-config'>
MathJax.Hub.Config({
HTML: ["input/TeX","output/HTML-CSS"],
TeX: {
Macros: {
bm: ["\\boldsymbol{#1}", 1],
argmax: ["\\mathop{\\rm arg\\,max}\\limits"],
argmin: ["\\mathop{\\rm arg\\,min}\\limits"]},
extensions: ["AMSmath.js","AMSsymbols.js"],
equationNumbers: { autoNumber: "AMS" } },
extensions: ["tex2jax.js"],
jax: ["input/TeX","output/HTML-CSS"],
tex2jax: { inlineMath: [ ['$','$'], ["\\(","\\)"] ],
displayMath: [ ['$$','$$'], ["\\[","\\]"] ],
processEscapes: true },
"HTML-CSS": { availableFonts: ["TeX"],
linebreaks: { automatic: true } }
});
</script>
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
今のヘッダーの設定では，標準で連番あり，"¥begin{equation*}"はだめ，"¥[ ¥]"で行番号なし，"¥tag{n}" で番号指定可になっている（¥はバックスラッシュを表す）。
\begin{equation}
\begin{aligned}
\sqrt{x}+\cos{y}\\
\dfrac{\log(x)}{e^x}
\end{aligned}
\end{equation}
\[
\sqrt{\dfrac{z^2+y^2}{x^2}}
\]
\begin{equation}
y= a x^2 + b x + c \tag{4}
\end{equation}
\begin{equation}
m \ddot{\bm{r}} = \bm{F}
\end{equation}
これで"¥bm{}"も使えるようになったのは次のヘッダーである。
参考にしたのは以下のサイトである。
・Mathjaxを使ってBloggerで数式を書く
・おもちゃ箱：Bloggerで数式を表示する
・MathJaxの使い方（黒木玄）
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
<script async='async' src='https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.2/MathJax.js?config=TeX-MML-AM_CHTML' type='text/javascript'/>
<script type='text/x-mathjax-config'>
MathJax.Hub.Config({
HTML: ["input/TeX","output/HTML-CSS"],
TeX: {
Macros: {
bm: ["\\boldsymbol{#1}", 1],
argmax: ["\\mathop{\\rm arg\\,max}\\limits"],
argmin: ["\\mathop{\\rm arg\\,min}\\limits"]},
extensions: ["AMSmath.js","AMSsymbols.js"],
equationNumbers: { autoNumber: "AMS" } },
extensions: ["tex2jax.js"],
jax: ["input/TeX","output/HTML-CSS"],
tex2jax: { inlineMath: [ ['$','$'], ["\\(","\\)"] ],
displayMath: [ ['$$','$$'], ["\\[","\\]"] ],
processEscapes: true },
"HTML-CSS": { availableFonts: ["TeX"],
linebreaks: { automatic: true } }
});
</script>
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2019年3月21日木曜日
国立大学法人法と学長の任期
以下は国立大学法人法からの学長の任期に係わる部分の抜粋である。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
国立大学法人法
第十条 各国立大学法人に、役員として、その長である学長及び監事二人を置く。
２ 略
第十二条 学長の任命は、国立大学法人の申出に基づいて、文部科学大臣が行う。
２ 前項の申出は、第一号に掲げる委員及び第二号に掲げる委員各同数をもって構成する会議（以下「学長選考会議」という。）の選考により行うものとする。
一 第二十条第二項第三号に掲げる者（注１）の中から同条第一項に規定する経営協議会において選出された者
二 第二十一条第二項第三号又は第四号に掲げる者（注２）の中から同条第一項に規定する教育研究評議会において選出された者
３ 前項各号に掲げる者のほか、学長選考会議の定めるところにより、学長又は理事を学長選考会議の委員に加えることができる。ただし、その数は、学長選考会議の委員の総数の三分の一を超えてはならない。
４ 学長選考会議に議長を置き、委員の互選によってこれを定める。
５ 議長は、学長選考会議を主宰する。
第十五条 学長の任期は、二年以上六年を超えない範囲内において、学長選考会議の議を経て、各国立大学法人の規則で定める。
２ 略
３ 略
４ 役員は、再任されることができる。以下略
（注１）経営協議会構成員から，学長と学長が指名する理事及び職員を除いた「三 当該国立大学法人の役員又は職員以外の者で大学に関し広くかつ高い識見を有するもののうちから，次条第一項に規定する教育研究評議会の意見を聴いて学長が任命するもの」
（注２）教育評議会構成員から，学長と学長が指名する理事を除いた「三 学部、研究科、大学附置の研究所その他の教育研究上の重要な組織の長のうち、教育研究評議会が定める者，四 その他教育研究評議会が定めるところにより学長が指名する職員」
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
国立大学法人法
第十条 各国立大学法人に、役員として、その長である学長及び監事二人を置く。
２ 略
第十二条 学長の任命は、国立大学法人の申出に基づいて、文部科学大臣が行う。
２ 前項の申出は、第一号に掲げる委員及び第二号に掲げる委員各同数をもって構成する会議（以下「学長選考会議」という。）の選考により行うものとする。
一 第二十条第二項第三号に掲げる者（注１）の中から同条第一項に規定する経営協議会において選出された者
二 第二十一条第二項第三号又は第四号に掲げる者（注２）の中から同条第一項に規定する教育研究評議会において選出された者
３ 前項各号に掲げる者のほか、学長選考会議の定めるところにより、学長又は理事を学長選考会議の委員に加えることができる。ただし、その数は、学長選考会議の委員の総数の三分の一を超えてはならない。
４ 学長選考会議に議長を置き、委員の互選によってこれを定める。
５ 議長は、学長選考会議を主宰する。
第十五条 学長の任期は、二年以上六年を超えない範囲内において、学長選考会議の議を経て、各国立大学法人の規則で定める。
２ 略
３ 略
４ 役員は、再任されることができる。以下略
（注１）経営協議会構成員から，学長と学長が指名する理事及び職員を除いた「三 当該国立大学法人の役員又は職員以外の者で大学に関し広くかつ高い識見を有するもののうちから，次条第一項に規定する教育研究評議会の意見を聴いて学長が任命するもの」
（注２）教育評議会構成員から，学長と学長が指名する理事を除いた「三 学部、研究科、大学附置の研究所その他の教育研究上の重要な組織の長のうち、教育研究評議会が定める者，四 その他教育研究評議会が定めるところにより学長が指名する職員」
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
したがって，国立大学法人法は学長が無期限に再任されることを排除しておらず，それは各国立大学法人の学長選考に関る規程に委ねられているようにみえる。
2019年3月20日水曜日
退職挨拶の練習（２）
懺悔の時間がやってまいりました。理科教育講座の越桐です。
1982年に本学に教務員として着任して37年になります。最初の10年ほどは天王寺分校に研
究室がありました。そのころは附属学校のグラウンドの裏にあった池田宿舎に住んでいました。その当時いっしょにすんでいた方は栗林学長と赤松先生だけになりました。
小専理科の実験の授業を池田分校の物理の伊藤太郎先生とともに担当していましたが，授業は３限目からなので，昼頃に池田分校に行き，夕方前のまだ明るいうちに帰宅するという夢のような日々を過ごしていました。ごめんなさい。
先生方，職員の皆様方には長らく大変お世話になりました。大学はこれから大変な時代を迎えますが，どうかお身体にお気をつけ，皆さん力を合わせてこの難関に立ち向かっていただければと思います。本日はどうもありがとうございました。
1982年に本学に教務員として着任して37年になります。最初の10年ほどは天王寺分校に研
究室がありました。そのころは附属学校のグラウンドの裏にあった池田宿舎に住んでいました。その当時いっしょにすんでいた方は栗林学長と赤松先生だけになりました。
小専理科の実験の授業を池田分校の物理の伊藤太郎先生とともに担当していましたが，授業は３限目からなので，昼頃に池田分校に行き，夕方前のまだ明るいうちに帰宅するという夢のような日々を過ごしていました。ごめんなさい。
先生方，職員の皆様方には長らく大変お世話になりました。大学はこれから大変な時代を迎えますが，どうかお身体にお気をつけ，皆さん力を合わせてこの難関に立ち向かっていただければと思います。本日はどうもありがとうございました。
2019年3月19日火曜日
運動物体からの斜方投射
昨日届いた「大学の物理教育（2019 Vol.25 No.1）」の講義室に，茨城大学工学部の坪井一洋さんが，「質点の投射角を再考する 大学の物理教育25(2019)34-37」という記事を書いていた。とてもおもしろかったので，その計算をフォローしてみることにした。
坪井さんの記事の趣旨は次のとおりである。一様重力場での放物運動では，斜方投射の投射角を45度としたときに質点を最も遠くまで投げることができる。しかし現実の現象（ゴルフ，野球，走り幅跳び，砲丸投げ）ではそうなっていない。空気抵抗で説明できる部分もあるがそうでない部分もある。スポーツ科学では初速度の大きさはが初速角に依存するという議論もある。そこで，これをモデル化して物理教育（初等力学）の題材にしようというものである。
原点をOとする慣性系Sを考える。Sのx軸上を一様な速度$u$で正方向に運動している物体がある。この物体がS系の原点Oを通過した瞬間（$t=0$）に，質量$m$の質点を，物体に固定した座標系Mからみて図のように初速度$v$，投射角$\phi$でx-y平面内に投射する。この質点には一様な重力$mg$が y 軸負方向に働いている。放物運動した質点は原点からの距離が$L$であるx軸上の点まで到達したとする。なお，S系からみた質点の速度ベクトル$\boldsymbol{w}$は図のように与えられる。
S系における原点からの距離$L$を最も大きくするにはどうすればよいか。原点における物体の初速度を$(V_x,V_y)$とすれば，高等学校の物理では次の結果を得ている。$L=\dfrac{2V_xV_y}{g}$
ところでこの$V_x$と$V_y$は今のモデルでは次のように与えられる。ただし，$0 < \phi < \pi/2$，$0 < \theta < \pi/2$とする。
\begin{equation}
\begin{aligned}
V_x = w \cos \theta = u + v \cos \phi \\
V_y = w \sin \theta = v \sin \phi
\end{aligned}
\end{equation}
つまり，$u,v$が与えられたときに，$L$を最大にする$\phi$または$\theta$を求めるのがここでの問題となる。なお，$w$は，ベクトルの図に余弦定理をあてはめた次の式に従うため，$\theta$の関数であると考えることにする（$\phi$の関数とすることもできる）。
\begin{equation}
w^2+u^2-2w u \cos\theta = v^2
\end{equation}
$w$についての２次方程式を解いて具体的な$\theta$の関数形を求め，さらに$\dfrac{d w}{d \theta}$を計算しておく。
\begin{equation}
\begin{aligned}
w &= u \cos\theta + \sqrt{v^2-u^2\sin^2\theta}\\
\frac{dw}{d\theta} &= -u\sin\theta + \frac{-u^2\sin\theta \cos\theta}{\sqrt{v^2-u^2\sin^2\theta}}\\
&= -u\sin\theta \ \frac{\sqrt{v^2-u^2\sin^2\theta} + u \cos\theta }{\sqrt{v^2-u^2\sin^2\theta}}\\
&= -\frac{u\sin\theta\ w }{\sqrt{v^2-u^2\sin^2\theta}}
\end{aligned}
\end{equation}
以下では，$\gamma=\dfrac{u}{v}$とおく。
(1) $L(\phi)$を最大にする$\phi$に対する条件を求める。
\begin{equation}
L(\phi)=\frac{2\ (u+v\cos\phi)\ v\sin\phi}{g}\\
\frac{dL}{d\phi}=\frac{2v}{g}\{-v \sin^2\phi +(u + v\cos\phi) \cos\phi \}=0\\
2v \cos^2\phi +u\cos\phi -v=0\\
2 \cos^2\phi + \gamma \cos\phi + 1 = 0\\
\therefore \cos\phi = \frac{\sqrt{\gamma^2+8}-\gamma}{4}
\end{equation}
(2) $L(\theta)$を最大にする$\theta$に対する条件を求める。
\begin{equation}
L(\theta)=\frac{2\ w\cos\theta\ w\sin\theta}{g}=\frac{w^2 \sin 2\theta}{g}\\
\frac{dL}{d\theta}=\frac{2w}{g}\,\frac{d w}{d\theta}\,\sin 2\theta + \frac{2w^2}{g}\,\cos 2\theta =0\\
\frac{-u\sin\theta \ \sin 2\theta}{\sqrt{v^2-u^2\sin^2\theta}}+ \cos 2\theta =0\\
4 u^2 \sin^4\theta\,(1-\sin^2\theta) = (v^2-u^2\sin^2\theta)(1-2\sin^2\theta)^2\\
v(1-2\sin^2\theta) = u\sin\theta\\
2 \sin^2\theta + \gamma \sin\theta + 1 = 0\\
\therefore \sin \theta = \frac{\sqrt{\gamma^2+8}-\gamma}{4}
\end{equation}
従って，最大の$L$を与える$\phi$と$\theta$は，$\phi+\theta=\frac{\pi}{2}$を満足している。
坪井さんの記事の趣旨は次のとおりである。一様重力場での放物運動では，斜方投射の投射角を45度としたときに質点を最も遠くまで投げることができる。しかし現実の現象（ゴルフ，野球，走り幅跳び，砲丸投げ）ではそうなっていない。空気抵抗で説明できる部分もあるがそうでない部分もある。スポーツ科学では初速度の大きさはが初速角に依存するという議論もある。そこで，これをモデル化して物理教育（初等力学）の題材にしようというものである。
図 運動物体からの斜方投射モデル
原点をOとする慣性系Sを考える。Sのx軸上を一様な速度$u$で正方向に運動している物体がある。この物体がS系の原点Oを通過した瞬間（$t=0$）に，質量$m$の質点を，物体に固定した座標系Mからみて図のように初速度$v$，投射角$\phi$でx-y平面内に投射する。この質点には一様な重力$mg$が y 軸負方向に働いている。放物運動した質点は原点からの距離が$L$であるx軸上の点まで到達したとする。なお，S系からみた質点の速度ベクトル$\boldsymbol{w}$は図のように与えられる。
S系における原点からの距離$L$を最も大きくするにはどうすればよいか。原点における物体の初速度を$(V_x,V_y)$とすれば，高等学校の物理では次の結果を得ている。$L=\dfrac{2V_xV_y}{g}$
ところでこの$V_x$と$V_y$は今のモデルでは次のように与えられる。ただし，$0 < \phi < \pi/2$，$0 < \theta < \pi/2$とする。
\begin{equation}
\begin{aligned}
V_x = w \cos \theta = u + v \cos \phi \\
V_y = w \sin \theta = v \sin \phi
\end{aligned}
\end{equation}
つまり，$u,v$が与えられたときに，$L$を最大にする$\phi$または$\theta$を求めるのがここでの問題となる。なお，$w$は，ベクトルの図に余弦定理をあてはめた次の式に従うため，$\theta$の関数であると考えることにする（$\phi$の関数とすることもできる）。
\begin{equation}
w^2+u^2-2w u \cos\theta = v^2
\end{equation}
$w$についての２次方程式を解いて具体的な$\theta$の関数形を求め，さらに$\dfrac{d w}{d \theta}$を計算しておく。
\begin{equation}
\begin{aligned}
w &= u \cos\theta + \sqrt{v^2-u^2\sin^2\theta}\\
\frac{dw}{d\theta} &= -u\sin\theta + \frac{-u^2\sin\theta \cos\theta}{\sqrt{v^2-u^2\sin^2\theta}}\\
&= -u\sin\theta \ \frac{\sqrt{v^2-u^2\sin^2\theta} + u \cos\theta }{\sqrt{v^2-u^2\sin^2\theta}}\\
&= -\frac{u\sin\theta\ w }{\sqrt{v^2-u^2\sin^2\theta}}
\end{aligned}
\end{equation}
以下では，$\gamma=\dfrac{u}{v}$とおく。
(1) $L(\phi)$を最大にする$\phi$に対する条件を求める。
\begin{equation}
L(\phi)=\frac{2\ (u+v\cos\phi)\ v\sin\phi}{g}\\
\frac{dL}{d\phi}=\frac{2v}{g}\{-v \sin^2\phi +(u + v\cos\phi) \cos\phi \}=0\\
2v \cos^2\phi +u\cos\phi -v=0\\
2 \cos^2\phi + \gamma \cos\phi + 1 = 0\\
\therefore \cos\phi = \frac{\sqrt{\gamma^2+8}-\gamma}{4}
\end{equation}
(2) $L(\theta)$を最大にする$\theta$に対する条件を求める。
\begin{equation}
L(\theta)=\frac{2\ w\cos\theta\ w\sin\theta}{g}=\frac{w^2 \sin 2\theta}{g}\\
\frac{dL}{d\theta}=\frac{2w}{g}\,\frac{d w}{d\theta}\,\sin 2\theta + \frac{2w^2}{g}\,\cos 2\theta =0\\
\frac{-u\sin\theta \ \sin 2\theta}{\sqrt{v^2-u^2\sin^2\theta}}+ \cos 2\theta =0\\
4 u^2 \sin^4\theta\,(1-\sin^2\theta) = (v^2-u^2\sin^2\theta)(1-2\sin^2\theta)^2\\
v(1-2\sin^2\theta) = u\sin\theta\\
2 \sin^2\theta + \gamma \sin\theta + 1 = 0\\
\therefore \sin \theta = \frac{\sqrt{\gamma^2+8}-\gamma}{4}
\end{equation}
従って，最大の$L$を与える$\phi$と$\theta$は，$\phi+\theta=\frac{\pi}{2}$を満足している。
2019年3月18日月曜日
バラカ：桐野夏生
まわりからどんどん書店が消えていく。昔は，毎日の仕事の帰りには本屋に立ち寄って，新刊をチェックすることができたが，最近はその習慣が失われてしまった。先日，久しぶりに上本町の近鉄の書店（JUNKUDO）に立ち寄ったときに見つけたのが，桐野夏生の「バラカ（上・下）集英社文庫」だった。
桐野夏生の本は，昔「ＯＵＴ」を図書館で借りて読んだ。その他にいくつか文庫本を買っていたはずだと思って調べると，昨年の部屋の片づけの際に処分していた。確か「柔らかな頬」と「グロテスク」だった。スパーク・ジョイがなかったのである。どんなストーリーだったかを思い出そうとしても印象に残っている部分がない。その点，「ＯＵＴ」にははっきりした読後記憶がある。桐野夏生は1951年に金沢に生まれているけれど，３歳で引っ越しているので，自分の世界線との交わりはほとんどない。
さて，「バラカ」は，福島の原発事故が起こったもう一つの虚構の日本での出来事を描写している。部分的には納得できるのだが，全体として詰めが十分ではないと感じた。自分が期待していたのは，現実の日本を照射するもう一つの世界の精密な描写だった。著者の関心はその物語世界の中の人物造形や関係の方にあって，それはそれでおもしろかったが・・・後半はご都合主義的なまとめ方のような・・・
P. S. 作中でバラカ／おじいちゃんが愛読するのがカズオ・イシグロの「わたしを離さないで」だったので，これは読んでみようと思った。
桐野夏生の本は，昔「ＯＵＴ」を図書館で借りて読んだ。その他にいくつか文庫本を買っていたはずだと思って調べると，昨年の部屋の片づけの際に処分していた。確か「柔らかな頬」と「グロテスク」だった。スパーク・ジョイがなかったのである。どんなストーリーだったかを思い出そうとしても印象に残っている部分がない。その点，「ＯＵＴ」にははっきりした読後記憶がある。桐野夏生は1951年に金沢に生まれているけれど，３歳で引っ越しているので，自分の世界線との交わりはほとんどない。
さて，「バラカ」は，福島の原発事故が起こったもう一つの虚構の日本での出来事を描写している。部分的には納得できるのだが，全体として詰めが十分ではないと感じた。自分が期待していたのは，現実の日本を照射するもう一つの世界の精密な描写だった。著者の関心はその物語世界の中の人物造形や関係の方にあって，それはそれでおもしろかったが・・・後半はご都合主義的なまとめ方のような・・・
P. S. 作中でバラカ／おじいちゃんが愛読するのがカズオ・イシグロの「わたしを離さないで」だったので，これは読んでみようと思った。
2019年3月17日日曜日
素数日
ほぼ毎日（多少前後して編集しているけれど）ブログを書き始めて，今日で100記事になった。いつまで続くのだろう？何のために続けているのだろう？惰性なのだろうか？
西暦（グレゴリオ暦）のyyyy年mm月dd日を8桁の整数 yyyymmdd で表したものが素数であるものを素数日と言い習わしているようだ。あちらこちらで目にします。で，Juliaで素数日を求める関数を書いてみた。あらかじめ，Primesパッケージを読み出すので，自分の仕事はほとんどない。
using Pkg
Pkg.add("Primes")
using Primes
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
function primeday(y)
u=28
if(y%4==0 && y%100!=0 || y%400==0)
u=29
end
# println(u)
mx=(31,u,31,30,31,30,31,31,30,31,30,31)
for m in 1:12, d in 1:mx[m]
if(isprime(y*10000+m*100+d))
println(y*10000+m*100+d)
end
end
end
primeday(2019)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
20190221
20190227
20190301
20190319
20190323
20190421
20190523
20190529
20190601
20190613
20190719
20190811
20190823
20190913
20191009
20191027
20191109
20191117
20191231
西暦（グレゴリオ暦）のyyyy年mm月dd日を8桁の整数 yyyymmdd で表したものが素数であるものを素数日と言い習わしているようだ。あちらこちらで目にします。で，Juliaで素数日を求める関数を書いてみた。あらかじめ，Primesパッケージを読み出すので，自分の仕事はほとんどない。
using Pkg
Pkg.add("Primes")
using Primes
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
function primeday(y)
u=28
if(y%4==0 && y%100!=0 || y%400==0)
u=29
end
# println(u)
mx=(31,u,31,30,31,30,31,31,30,31,30,31)
for m in 1:12, d in 1:mx[m]
if(isprime(y*10000+m*100+d))
println(y*10000+m*100+d)
end
end
end
primeday(2019)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
20190221
20190227
20190301
20190319
20190323
20190421
20190523
20190529
20190601
20190613
20190719
20190811
20190823
20190913
20191009
20191027
20191109
20191117
20191231
2019年3月16日土曜日
機械学習と物理
長らく物理学会には足を向けていない。今年の第74回年会は九州大学の伊都キャンパスで開催されている。ちょうど今の時間に，シンポジウム「機械学習と物理」が行われている。こんな感じ。
１（一般シンポジウム講演）はじめに
阪大理・物理，橋本幸士
２（一般シンポジウム講演）物性物理のグランドチャレンジに対する重回帰分析と機械学習
東大院工，今田正俊
３（一般シンポジウム講演）強い相互作用の最難問 — 中性子星の状態方程式
東大理・物理・原子核理論，福嶋健二
４（一般シンポジウム講演）データ駆動手法による相関物質の予測と理解
産総研 CD-FMat，三宅隆
５（一般シンポジウム講演）機械学習によるマルコフ連鎖モンテカルロ法の高速化へ向けて
理研(AIP/iTHEMS), 慶應大・数理，田中章詞
６（一般シンポジウム講演）機械学習による特徴抽出と，繰り込み群や熱力学との関係
OIST，船井正太郎
７（一般シンポジウム講演）広域撮像宇宙サーベイによるビッグデータ宇宙論
東大理・物理・宇宙理論，吉田直紀
８（一般シンポジウム講演）量子力学と機械学習の数理
東北大院情報科学，大関真之
ちなみに，昨年の2018年度日本物理学会科学セミナーのテーマも「ＡＩ（人工知能）と物理学（東京大学駒場キャンパス 数理科学研究棟 大講義室）」であった。そのプログラムは以下の通り（ちょっと被っている）。
8月11日（土・祝）10:00-16:30
１ はじめに
日本物理学会会長 川村 光
２ 情報処理技術としてのAI
中島 秀之（札幌市立大学 学長）
３ 思考力を競うゲームの人工知能技術発展の歴史と現状
保木 邦仁（電気通信大学大学院情報理工学研究科 准教授）
４ 広域宇宙撮像データを用いたビッグデータ宇宙論
吉田 直紀（東京大学大学院理学系研究科 教授）
５ 量子コンピュータが人工知能を加速する
大関 真之（東北大学大学院情報科学研究科 准教授）
６ 深層学習と時空
橋本 幸士 (大阪大学大学院理学研究科 教授)
8月12日（日）10:00-16:40
７ 深層学習とは何か、そしてどんなことが出来るようになっているのか
瀧 雅人（理化学研究所数理創造プログラム(iTHEMS) 上級研究員）
８ 多層畳み込みニューラルネットワークで求めた量子相転移の相図
大槻 東巳（上智大学理工学部機能創造理工学科 教授）
９ 人工知能と脳科学
甘利 俊一（理化学研究所脳神経科学研究センター 特別顧問）
10 量子力学の問題をニューラルネットワークで解く
斎藤 弘樹（電気通信大学大学院情報理工学研究科 教授）
11 AI は物理において何の役に立つか？
寺倉 清之（物質・材料研究機構 名誉フェロー・エグゼクティブアドバイザー）
12 おわりに
日本物理学会科学セミナー担当理事 迫田 和彰
日本物理学会誌の方でも「シリーズ 人工知能と物理学」という特集が開始され，産総研の神嶌敏弘さんの「変わりゆく機械学習と変わらない機械学習」が読める。そのうち日本物理教育学会誌にも，「機械学習と物理教育」とか「人工知能と物理教育」などの特集が組まれる時代がくるのだろうか。５年後？１０年後？
参考：Quantum Machine Learning（Wikipedia）
List of Acronyms
ANN: Artificial neural network
BM: Boltzmann machine
BN: Bayesian network
CDL: Classical deep learning
CML: Classical machine learning
HMM: Hidden Markov model
HQMM: Hidden quantum Markov model
k-NN: k-nearest neighbours
NMR: Nuclear magnetic resonance
PCA: Principal component analysis
QBN: Quantum Bayesian network
QDL: Quantum deep learning
QML: Quantum machine learning
QPCA: Quantum principal component analysis
QRAM: Quantum random access memory
RAM: Random access memory
SQW: Stochastic quantum walk
SVM: Support vector machine
WNN: Weightless neural network
１（一般シンポジウム講演）はじめに
阪大理・物理，橋本幸士
２（一般シンポジウム講演）物性物理のグランドチャレンジに対する重回帰分析と機械学習
東大院工，今田正俊
３（一般シンポジウム講演）強い相互作用の最難問 — 中性子星の状態方程式
東大理・物理・原子核理論，福嶋健二
４（一般シンポジウム講演）データ駆動手法による相関物質の予測と理解
産総研 CD-FMat，三宅隆
５（一般シンポジウム講演）機械学習によるマルコフ連鎖モンテカルロ法の高速化へ向けて
理研(AIP/iTHEMS), 慶應大・数理，田中章詞
６（一般シンポジウム講演）機械学習による特徴抽出と，繰り込み群や熱力学との関係
OIST，船井正太郎
７（一般シンポジウム講演）広域撮像宇宙サーベイによるビッグデータ宇宙論
東大理・物理・宇宙理論，吉田直紀
８（一般シンポジウム講演）量子力学と機械学習の数理
東北大院情報科学，大関真之
ちなみに，昨年の2018年度日本物理学会科学セミナーのテーマも「ＡＩ（人工知能）と物理学（東京大学駒場キャンパス 数理科学研究棟 大講義室）」であった。そのプログラムは以下の通り（ちょっと被っている）。
8月11日（土・祝）10:00-16:30
１ はじめに
日本物理学会会長 川村 光
２ 情報処理技術としてのAI
中島 秀之（札幌市立大学 学長）
３ 思考力を競うゲームの人工知能技術発展の歴史と現状
保木 邦仁（電気通信大学大学院情報理工学研究科 准教授）
４ 広域宇宙撮像データを用いたビッグデータ宇宙論
吉田 直紀（東京大学大学院理学系研究科 教授）
５ 量子コンピュータが人工知能を加速する
大関 真之（東北大学大学院情報科学研究科 准教授）
６ 深層学習と時空
橋本 幸士 (大阪大学大学院理学研究科 教授)
8月12日（日）10:00-16:40
７ 深層学習とは何か、そしてどんなことが出来るようになっているのか
瀧 雅人（理化学研究所数理創造プログラム(iTHEMS) 上級研究員）
８ 多層畳み込みニューラルネットワークで求めた量子相転移の相図
大槻 東巳（上智大学理工学部機能創造理工学科 教授）
９ 人工知能と脳科学
甘利 俊一（理化学研究所脳神経科学研究センター 特別顧問）
10 量子力学の問題をニューラルネットワークで解く
斎藤 弘樹（電気通信大学大学院情報理工学研究科 教授）
11 AI は物理において何の役に立つか？
寺倉 清之（物質・材料研究機構 名誉フェロー・エグゼクティブアドバイザー）
12 おわりに
日本物理学会科学セミナー担当理事 迫田 和彰
参考：Quantum Machine Learning（Wikipedia）
List of Acronyms
ANN: Artificial neural network
BM: Boltzmann machine
BN: Bayesian network
CDL: Classical deep learning
CML: Classical machine learning
HMM: Hidden Markov model
HQMM: Hidden quantum Markov model
k-NN: k-nearest neighbours
NMR: Nuclear magnetic resonance
PCA: Principal component analysis
QBN: Quantum Bayesian network
QDL: Quantum deep learning
QML: Quantum machine learning
QPCA: Quantum principal component analysis
QRAM: Quantum random access memory
RAM: Random access memory
SQW: Stochastic quantum walk
SVM: Support vector machine
WNN: Weightless neural network
2019年3月15日金曜日
円周率の分数近似
円周率の日の余震が続いている。3月14日，Googleは円周率の計算の世界記録を更新したことを発表したとの報道があった。31.4兆桁（31,415,926,535,897桁）までの計算を111.8日（start to end で121.1日）かけて計算したそうだ。このプロジェクトを主導したのは，筑波大学出身の技術者 Emma Haruka Iwao （岩尾・エマ・はるか）さん。今年の1月には終わっていたが，このよき日に発表された。つぎに記録に挑戦するグループは，314兆桁を越えないといけないのか。
さて，だいぶ昔（2011.6.30）のWolfram Blogで，Jon McLooneが，"All Rational Approximations of Pi Are Useless"という記事を書いている。
どういうことかというと・・・円周率のよくある近似値22/7を例にとってみる。22/7 = 3.1426... は小数点以下2桁まで正しいので，小数点も含め，3 . 1 4 の4文字が表現されたことになる。このために必要な数は 2 2 / 7 の4文字である。これじゃ文字数が節約できていない。どちらを覚えてもいっしょじゃないか，ということらしい。
もっとやってみると，355/113 = 3.141592... これは小数点以下6桁すなわち8文字まで正しい円周率に対応している。一方，355/113は7文字，さきほどよりちょっとだけましだ。
McLooneの記事にあるMathematicaの関数を用いて，小数点以下50桁のπとその分数表示を求めると，
\begin{equation}
\frac{26151465932107044561886949}{8324270144388272579650158} = \\
\\
3.14159265358979323846264338327950288419716939937510...
\end{equation}
ということで共に52文字必要としている。でこぼこはあるが，おおむねこの比が1で推移することから，分数表示意味ないじゃんということになったようだ。
必要なMathematica関数は以下のようになっている。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
In[1]:= characterCount[expr_] := StringLength[
StringReplace[ToString[Unevaluated[expr], InputForm],
" " -> ""]]; SetAttribute[characterCount, HoldAll];
In[2]:= matchingCharacters[expr_, target_] :=
Ceiling[-Log10[Abs[expr - target]]] + 1
In[3]:= $MaxExtraPrecision = Infinity;
In[4]:= candidates = Table[Rationalize[Pi, 10^(-n)], {n, 1001}];
In[5]:= a = Table[candidates[[k]], {k, 1, 10}]
Out[5]= {22/7, 22/7, 201/64, 333/106, 355/113, 355/113,
75948/24175, 100798/32085, 103993/33102, 312689/99532}
In[6]:= b =
Table[characterCount[candidates[[k]]], {k, 1, 1001, 50}]
Out[6]= {4, 54, 103, 154, 203, 253, 303, 353, 404, 455, 504,
553, 603, 652, 704, 753, 803, 853, 904, 953, 1003}
In[7]:= c =
Table[matchingCharacters[candidates[[k]], Pi], {k, 1, 1001, 50}]
Out[7]= {4, 53, 103, 153, 204, 253, 303, 353, 404, 453, 503,
553, 604, 653, 704, 753, 803, 853, 903, 953, 1003}
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
さて，だいぶ昔（2011.6.30）のWolfram Blogで，Jon McLooneが，"All Rational Approximations of Pi Are Useless"という記事を書いている。
どういうことかというと・・・円周率のよくある近似値22/7を例にとってみる。22/7 = 3.1426... は小数点以下2桁まで正しいので，小数点も含め，3 . 1 4 の4文字が表現されたことになる。このために必要な数は 2 2 / 7 の4文字である。これじゃ文字数が節約できていない。どちらを覚えてもいっしょじゃないか，ということらしい。
もっとやってみると，355/113 = 3.141592... これは小数点以下6桁すなわち8文字まで正しい円周率に対応している。一方，355/113は7文字，さきほどよりちょっとだけましだ。
McLooneの記事にあるMathematicaの関数を用いて，小数点以下50桁のπとその分数表示を求めると，
\begin{equation}
\frac{26151465932107044561886949}{8324270144388272579650158} = \\
\\
3.14159265358979323846264338327950288419716939937510...
\end{equation}
ということで共に52文字必要としている。でこぼこはあるが，おおむねこの比が1で推移することから，分数表示意味ないじゃんということになったようだ。
必要なMathematica関数は以下のようになっている。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
In[1]:= characterCount[expr_] := StringLength[
StringReplace[ToString[Unevaluated[expr], InputForm],
" " -> ""]]; SetAttribute[characterCount, HoldAll];
In[2]:= matchingCharacters[expr_, target_] :=
Ceiling[-Log10[Abs[expr - target]]] + 1
In[3]:= $MaxExtraPrecision = Infinity;
In[4]:= candidates = Table[Rationalize[Pi, 10^(-n)], {n, 1001}];
In[5]:= a = Table[candidates[[k]], {k, 1, 10}]
Out[5]= {22/7, 22/7, 201/64, 333/106, 355/113, 355/113,
75948/24175, 100798/32085, 103993/33102, 312689/99532}
In[6]:= b =
Table[characterCount[candidates[[k]]], {k, 1, 1001, 50}]
Out[6]= {4, 54, 103, 154, 203, 253, 303, 353, 404, 455, 504,
553, 603, 652, 704, 753, 803, 853, 904, 953, 1003}
In[7]:= c =
Table[matchingCharacters[candidates[[k]], Pi], {k, 1, 1001, 50}]
Out[7]= {4, 53, 103, 153, 204, 253, 303, 353, 404, 453, 503,
553, 604, 653, 704, 753, 803, 853, 903, 953, 1003}
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2019年3月14日木曜日
ディオファントス方程式 $x^3+y^3+z^3=n$（２）
（ディオファントス方程式 $x^3+y^3+z^3=n$（１）からの続き）
というわけで，今日，３月14日は円周率の日なのかホワイトデイなのかアインシュタインの誕生日なのかホーキングの命日なのか巷では協議が続いている。
x^3+y^3+z^3=n が n=9k+4, 9k+5 では解を持たないことは簡単に示すことができる。xを9k+mとして９で割った余りをmとする。m=0,1,... 8である。このときx^3を９で割った余りは，m^3を９で割った余りになるので，(0,1,8,27,64,125,216,343,512) ≡ (0,1,8,0,1,8,0,1,8)。
そこで，この剰余(0,1,8)の３つの和のすべての組み合わせを９で割った余りの可能性を調べれば良い。(0,0,0)≡0, (0,0,1)≡1, (0,0,8)≡8, (0,1,1)≡2, (0,1,8)≡0, (0,8,8)≡7, (1,1,1)≡3, (1,1,8)≡1, (1,8,8)≡8, (8,8,8)≡6であり，ここから4と5は出てこない。
少しだけ前回のプログラムを高速化してみる。k^3の配列の次元をmとして，前回のものはO(m^3)，今回のものはO(m^2)である。
（注）Juliaは0の配列添字を許さないので，この度もx,y,zから0は除いて考えている。0を含めれば，91=0^3+3^3+4^3 など簡単なものがあると思うけれど。
（注）Juliaは0の配列添字を許さないので，この度もx,y,zから0は除いて考えている。0を含めれば，91=0^3+3^3+4^3 など簡単なものがあると思うけれど。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
function pre(m,t)
for l in 1:m
t[l]=l^3
end
end
function cinv(x)
return Int(ceil(cbrt(abs(x))))
end
function dpe(n,m,t)
l=0
if(n%9 ==4 || n%9 ==5)
return (n, n%9, 0)
end
for i in 1:m
d=cinv(n-t[i])
if(n>=t[i])
for j in 1:d
k=cinv(n-t[i]-t[j])
l=l+1
if k !=0 && (t[i]+t[j]+t[k]==n || t[i]+t[j]+t[k]==-n)
return (n,l)
end
end
for j in d+1:i
k=cinv(n-t[i]-t[j])
l=l+1
if k !=0 && (t[i]+t[j]-t[k]==n || t[i]+t[j]-t[k]==-n)
return (n,l)
end
end
else
for j in d+1:i
k=cinv(t[i]-n-t[j])
l=l+1
if k !=0 && (t[i]-t[j]+t[k]==n || t[i]-t[j]+t[k]==-n)
return (n,l)
end
end
for j in 1:d
k=cinv(t[i]-n-t[j])
l=l+1
if k !=0 && (t[i]-t[j]-t[k]==n || t[i]-t[j]-t[k]==-n)
return (n,l)
end
end
end
end
return (m,l)
end
function pst(l,m,t)
for n in 1:l
println(dpe(n,m,t))
end
end
m=10000
t=ones(Int64,m)
pre(m,t)
@time pst(99,m,t)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
(1, 1)
(2, 26)
(3, 1)
(4, 4, 0)
(5, 5, 0)
(6, 1)
(7, 5491)
(8, 1)
(9, 23487)
(10, 1)
(11, 4)
(12, 61)
(13, 4, 0)
(14, 5, 0)
(15, 1057)
(16, 100)
(17, 1)
(18, 4)
(19, 185)
(20, 1561)
(21, 131)
(22, 4, 0)
(23, 5, 0)
(24, 2)
(25, 2)
(26, 48674)
(27, 1)
(28, 146)
(29, 1)
(10000, 50004997)
(31, 4, 0)
(32, 5, 0)
(10000, 50004997)
(34, 10)
(35, 96)
(36, 1)
(37, 1574)
(38, 364)
(10000, 50004997)
(40, 4, 0)
(41, 5, 0)
(10000, 50004997)
(43, 2)
(44, 30)
(45, 152329)
(46, 422)
(47, 31)
(48, 6)
(49, 4, 0)
(50, 5, 0)
(51, 317009)
(10000, 50004997)
(53, 9)
(54, 70)
(55, 1)
(56, 240)
(57, 726)
(58, 4, 0)
(59, 5, 0)
(60, 9)
(61, 466761)
(62, 2)
(63, 23)
(64, 5)
(65, 6188)
(66, 5)
(67, 4, 0)
(68, 5, 0)
(69, 342)
(70, 219)
(71, 12)
(72, 50)
(73, 5)
(10000, 50004998)
(10000, 50004998)
(76, 4, 0)
(77, 5, 0)
(78, 1509)
(79, 612)
(80, 6)
(81, 3)
(82, 96)
(83, 11)
(10000, 50004994)
(85, 4, 0)
(86, 5, 0)
(87, 9120551)
(88, 13)
(89, 884141)
(90, 7)
(91, 72576)
(92, 1)
(93, 20)
(94, 4, 0)
(95, 5, 0)
(96, 239)
(97, 6)
(98, 108)
(99, 2)
14.962059 seconds (103.76 k allocations: 4.814 MiB)
function pre(m,t)
for l in 1:m
t[l]=l^3
end
end
function cinv(x)
return Int(ceil(cbrt(abs(x))))
end
function dpe(n,m,t)
l=0
if(n%9 ==4 || n%9 ==5)
return (n, n%9, 0)
end
for i in 1:m
d=cinv(n-t[i])
if(n>=t[i])
for j in 1:d
k=cinv(n-t[i]-t[j])
l=l+1
if k !=0 && (t[i]+t[j]+t[k]==n || t[i]+t[j]+t[k]==-n)
return (n,l)
end
end
for j in d+1:i
k=cinv(n-t[i]-t[j])
l=l+1
if k !=0 && (t[i]+t[j]-t[k]==n || t[i]+t[j]-t[k]==-n)
return (n,l)
end
end
else
for j in d+1:i
k=cinv(t[i]-n-t[j])
l=l+1
if k !=0 && (t[i]-t[j]+t[k]==n || t[i]-t[j]+t[k]==-n)
return (n,l)
end
end
for j in 1:d
k=cinv(t[i]-n-t[j])
l=l+1
if k !=0 && (t[i]-t[j]-t[k]==n || t[i]-t[j]-t[k]==-n)
return (n,l)
end
end
end
end
return (m,l)
end
function pst(l,m,t)
for n in 1:l
println(dpe(n,m,t))
end
end
m=10000
t=ones(Int64,m)
pre(m,t)
@time pst(99,m,t)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
(1, 1)
(2, 26)
(3, 1)
(4, 4, 0)
(5, 5, 0)
(6, 1)
(7, 5491)
(8, 1)
(9, 23487)
(10, 1)
(11, 4)
(12, 61)
(13, 4, 0)
(14, 5, 0)
(15, 1057)
(16, 100)
(17, 1)
(18, 4)
(19, 185)
(20, 1561)
(21, 131)
(22, 4, 0)
(23, 5, 0)
(24, 2)
(25, 2)
(26, 48674)
(27, 1)
(28, 146)
(29, 1)
(10000, 50004997)
(31, 4, 0)
(32, 5, 0)
(10000, 50004997)
(34, 10)
(35, 96)
(36, 1)
(37, 1574)
(38, 364)
(10000, 50004997)
(40, 4, 0)
(41, 5, 0)
(10000, 50004997)
(43, 2)
(44, 30)
(45, 152329)
(46, 422)
(47, 31)
(48, 6)
(49, 4, 0)
(50, 5, 0)
(51, 317009)
(10000, 50004997)
(53, 9)
(54, 70)
(55, 1)
(56, 240)
(57, 726)
(58, 4, 0)
(59, 5, 0)
(60, 9)
(61, 466761)
(62, 2)
(63, 23)
(64, 5)
(65, 6188)
(66, 5)
(67, 4, 0)
(68, 5, 0)
(69, 342)
(70, 219)
(71, 12)
(72, 50)
(73, 5)
(10000, 50004998)
(10000, 50004998)
(76, 4, 0)
(77, 5, 0)
(78, 1509)
(79, 612)
(80, 6)
(81, 3)
(82, 96)
(83, 11)
(10000, 50004994)
(85, 4, 0)
(86, 5, 0)
(87, 9120551)
(88, 13)
(89, 884141)
(90, 7)
(91, 72576)
(92, 1)
(93, 20)
(94, 4, 0)
(95, 5, 0)
(96, 239)
(97, 6)
(98, 108)
(99, 2)
14.962059 seconds (103.76 k allocations: 4.814 MiB)
2019年3月13日水曜日
ディオファントス方程式 $x^3+y^3+z^3=n$（１）
昨日，ディオファントス方程式，$x^3+y^3+z^3=n$ の $n=33$ の解が見つかったというニュースが飛び込んできた。ここで，$n$は自然数，$x, y, z$ は整数である。例えば，$(x,y,z)=(1,1,-1)$ ならば $n=1$であり，$(x,y,z)=(7,-6,-5)$ならば $n=2$であり，$(x,y,z)=(1, 1, 1)$ならば $n=3$である。
\begin{equation}
33= (8,866,128,975,287,528)^3 \\
+(-8,778,405,442,862,239)^3 \\
+(-2,736,111,468,807,040)^3
\end{equation}
20年ほどまえの，KOYAMA, TSURUOKA & SEKIGAWA の論文では，nが1000以下では，
n= 30, 33, 42, 52, 74, 75, 110, 114, 156, 165, 195, 290, 318, 366, 390, 420, 435, 444, 452, 501, 530, 534, 564, 579, 588, 600, 606, 609, 618, 627, 633, 732, 735, 758, 767, 786, 789, 795, 830, 834, 861, 894, 903, 906, 912, 921, 933, 948, 964, 969, 975.
が未発見となっている。
簡単なJuliaのプログラムをつくってみた。$1^3$から$1000^3$までの組み合わせでつくれない２桁以下の$n$は，30,33,39,42,52,74,75,84,87の９つである。『数学者の密室（三島久典さん）』のページにはさらに詳しい情報があり，$0 \le n \le 99$, not equal 4 or 5 (mod 9) , $0 \le |x| \le |y| \le |z| \le 10^{10}$の範囲で，未解決なものはその記事の執筆時点で 33, 42, 74となっていた。
つまり，
30=(-283059965)^3+(-2218888517)^3+(2220422932)^3
33= ? → 上述のとおり解決
39=(117367)^3+(134476)^3+(-159380)^3
42= ?
52=(23961292454)^3+(60702901317)^3+(-61922712865)^3
74= ?
75=(4381159)^3+(435203083)^3+(-435203231)^3
84=(-8241191)^3+(-41531726)^3+(41639611)^3
87=(-1972)^3+(-4126)^3+(4271)^3
である。
（注）オンライン整数列大事典では，このあたり（A173515）。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c=zeros(Int64,100)
function dpt(n,c)
a=ones(Int128,n)
for m in 1:n
a[m]=m^3
end
for i in 1:n, j in -n:n, k in j:n
if(i*j*k != 0 && j+k !=0 && i+j !=0 && i+k !=0)
b=a[i]+sign(j)*a[abs(j)]+sign(k)*a[abs(k)]
if(b>0 && b<100)
c[b]=b
end
end
end
end
for m in 10:10
println(m)
@time dft(c,100*m)
println(c,count(x->x==0,c))
end
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
10
（ディオファントス方程式 $x^3+y^3+z^3=n$（２）に続く）
\begin{equation}
33= (8,866,128,975,287,528)^3 \\
+(-8,778,405,442,862,239)^3 \\
+(-2,736,111,468,807,040)^3
\end{equation}
20年ほどまえの，KOYAMA, TSURUOKA & SEKIGAWA の論文では，nが1000以下では，
n= 30, 33, 42, 52, 74, 75, 110, 114, 156, 165, 195, 290, 318, 366, 390, 420, 435, 444, 452, 501, 530, 534, 564, 579, 588, 600, 606, 609, 618, 627, 633, 732, 735, 758, 767, 786, 789, 795, 830, 834, 861, 894, 903, 906, 912, 921, 933, 948, 964, 969, 975.
が未発見となっている。
簡単なJuliaのプログラムをつくってみた。$1^3$から$1000^3$までの組み合わせでつくれない２桁以下の$n$は，30,33,39,42,52,74,75,84,87の９つである。『数学者の密室（三島久典さん）』のページにはさらに詳しい情報があり，$0 \le n \le 99$, not equal 4 or 5 (mod 9) , $0 \le |x| \le |y| \le |z| \le 10^{10}$の範囲で，未解決なものはその記事の執筆時点で 33, 42, 74となっていた。
つまり，
30=(-283059965)^3+(-2218888517)^3+(2220422932)^3
33= ? → 上述のとおり解決
39=(117367)^3+(134476)^3+(-159380)^3
42= ?
52=(23961292454)^3+(60702901317)^3+(-61922712865)^3
74= ?
75=(4381159)^3+(435203083)^3+(-435203231)^3
84=(-8241191)^3+(-41531726)^3+(41639611)^3
87=(-1972)^3+(-4126)^3+(4271)^3
である。
（注）オンライン整数列大事典では，このあたり（A173515）。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c=zeros(Int64,100)
function dpt(n,c)
a=ones(Int128,n)
for m in 1:n
a[m]=m^3
end
for i in 1:n, j in -n:n, k in j:n
if(i*j*k != 0 && j+k !=0 && i+j !=0 && i+k !=0)
b=a[i]+sign(j)*a[abs(j)]+sign(k)*a[abs(k)]
if(b>0 && b<100)
c[b]=b
end
end
end
end
for m in 10:10
println(m)
@time dft(c,100*m)
println(c,count(x->x==0,c))
end
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
10
（ディオファントス方程式 $x^3+y^3+z^3=n$（２）に続く）
2019年3月12日火曜日
World Wide Web
本日2019年3月12日のグーグルは，ワールド・ワイド・ウェブ(WWW)誕生30周年のデザインになっていた。ワールド・ワイド・ウェブの歴史が平成の歴史とほぼ重なっているのか。1989年の3月12日に欧州原子核研究機構（CERN=Conseil Europeen pour la Recherche Nucleaire）のティム・バーナーズ＝リーが，情報管理システムの提案をしたのが出発点となっている。1990年の12月には，NeXTコンピュータ上にシステムを構築し，1991年には，世界最初のウェブサイト http://info.cern.ch/ が公開されている。
その2年後の1993年に，画像とテキストが統合されたブラウザ，NCSAモザイクが登場する。そのシステムを開発したマーク・アンドリーセンはNetscape社を立ち上げ，しばらくは，ネットスケープ・ナビゲータが標準ブラウザの位置を占めていた。まもなく，それはMicrosoftのインターネット・エクスプローラに取って代わられてしまった。
ウェブサーバのソフトとしては，ティム・バーナーズ＝リーによって最初に作られたCERN-httpd，NCSAモザイクに対応した NCSA HTTPdなどがあった。20年ほど前に，最初にサーバにインストールしたのはNCSA HTTPdだったが，これが，Apache HTTP Serverに引き継がれたので切り替えた。10年前ごろには忙しくなったのか歳を取ったせいか，サーバーのお守りもできなくなってしまった。
このあたりの歴史は，TheWebKANZAKIのウェブブラウザ小史2001に詳しい。
その2年後の1993年に，画像とテキストが統合されたブラウザ，NCSAモザイクが登場する。そのシステムを開発したマーク・アンドリーセンはNetscape社を立ち上げ，しばらくは，ネットスケープ・ナビゲータが標準ブラウザの位置を占めていた。まもなく，それはMicrosoftのインターネット・エクスプローラに取って代わられてしまった。
ウェブサーバのソフトとしては，ティム・バーナーズ＝リーによって最初に作られたCERN-httpd，NCSAモザイクに対応した NCSA HTTPdなどがあった。20年ほど前に，最初にサーバにインストールしたのはNCSA HTTPdだったが，これが，Apache HTTP Serverに引き継がれたので切り替えた。10年前ごろには忙しくなったのか歳を取ったせいか，サーバーのお守りもできなくなってしまった。
このあたりの歴史は，TheWebKANZAKIのウェブブラウザ小史2001に詳しい。
写真：NCSA MOSAIC beta（ウェブアーカイブから）
2019年3月11日月曜日
第二種ソフトウェア危機とマクロユーザサービス
標題の文書が発掘された。たぶん，1993年前後のものではないかと思われる。大阪大学大型計算機センターの何かによせて寄稿することを想定した文書だが（利用相談員の自己紹介かな），どこかの段階でボツになったのかもしれない。とりあえず，参考のためにテキストを起こして再現してみる。第二種というのは，研究テーマだった原子核の弱い相互作用に存在するがどうかが問題となっていたＧパリティ非保存の「第二種カレント」にかけているのだが，そんなもの誰にも分からない。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
第二種ソフトウェア危機とマクロユーザサービス
−ユーザのオタク化による情報環境破壊の進行と対策−
大阪教育大学 物理学教室 越桐國雄 (メールアドレス略)
自己紹介： 私の専門分野は原子核理論（電弱相互作用と核構造）です。センターのACOSをTSSで使い始めたのは1980年ごろで，それまでは阪大核物理研究センターのTOSBACを使っていました。最近は主にワークステーションを使っています。スーパーコンピュータでUNIXが使えれば，また変わるかもしれません。相談分野は次のとおりです。
機種：ACOS，FACOM，NeXTstation，SPARCstation，Macintosh，PC-9801，FM-TOWNS．言語：FORTRAN，Mathematica．
（#ただ，どれもほとんどあてになりませんのでよろしく。）
大阪教育大の紹介： 大教大では柏原市への移転統合を機会に，情報処理センターが新たに設置されました。ホストはFACOM-M770で，ワークステーション30台，パーソナルコンピュータ50台等の構成です。各研究室からはTCP/IPで学内LANにアクセスします。N1ネットワークへは64kbpsの専用回線で接続し，SINETによる学外へのIP接続も予定しています。しかしあまりの急激な変化（これまではホスト無し，TSS端末数台とMacII40台等があっただけ）に戸惑っているような次第です。（#ほんとにもうたいへんなのです。）
センターへの希望： ダウンサイジングとネットワークング及びフリーソフトウェアの浸透により，スタッフのそろわない中小規模大学や情報系以外の研究室では，ユーザのシステム管理やソフトウェア管理の負担が年々重くなっています。また，ネットワークの発展に伴うコンピュータのメディア化が我々の情報処理能力の限界を越えて進もうとしています（ニュースとメールを読み，ファイル転送してmakeすると１日が終わってしまう!?）。こうしてあふれた情報ゴミがユーザのオタク化を促進しています。いわゆる「ソフトウェア危機」＝ソフトウェア生産場面での需給ギャップに対して，これらはソフトウェア消費場面での需給ギャップといえるかもしれません。これを，「第二種ソフトウェア危機」と呼ぶことにします。第二種ソフトウェア危機は人間の基本的な欲求に密接に関係しているため，より本質的で深刻です。これを回避するための一助として，研究室，教室あるいは小規模大学のユーザグループ（＝マクロユーザ）に対する適切な情報提供，事例紹介や講習等のサポートをセンターにお願いしたいと考えています。（#でも難しいですよね。中途半端な結論で終わってしまった・・・）
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
(注) FM-TOWNSは，日本で最初に発売されたCD-ROMドライブ付きパソコンだった。情報処理センターの副センター長だった定金晃三先生らと，これいいよね（デジタル教材メディア活用のはしりだぁ），とかいいながらネットワーク整備予算を使って一教室分導入したものである。そのうちビジネス用のPCにもCD-ROMがつくのが当たり前になってしまった。そして，これからはCD-ROMドライブが付かないのものが普通になるようだ。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
第二種ソフトウェア危機とマクロユーザサービス
−ユーザのオタク化による情報環境破壊の進行と対策−
大阪教育大学 物理学教室 越桐國雄 (メールアドレス略)
自己紹介： 私の専門分野は原子核理論（電弱相互作用と核構造）です。センターのACOSをTSSで使い始めたのは1980年ごろで，それまでは阪大核物理研究センターのTOSBACを使っていました。最近は主にワークステーションを使っています。スーパーコンピュータでUNIXが使えれば，また変わるかもしれません。相談分野は次のとおりです。
機種：ACOS，FACOM，NeXTstation，SPARCstation，Macintosh，PC-9801，FM-TOWNS．言語：FORTRAN，Mathematica．
（#ただ，どれもほとんどあてになりませんのでよろしく。）
大阪教育大の紹介： 大教大では柏原市への移転統合を機会に，情報処理センターが新たに設置されました。ホストはFACOM-M770で，ワークステーション30台，パーソナルコンピュータ50台等の構成です。各研究室からはTCP/IPで学内LANにアクセスします。N1ネットワークへは64kbpsの専用回線で接続し，SINETによる学外へのIP接続も予定しています。しかしあまりの急激な変化（これまではホスト無し，TSS端末数台とMacII40台等があっただけ）に戸惑っているような次第です。（#ほんとにもうたいへんなのです。）
センターへの希望： ダウンサイジングとネットワークング及びフリーソフトウェアの浸透により，スタッフのそろわない中小規模大学や情報系以外の研究室では，ユーザのシステム管理やソフトウェア管理の負担が年々重くなっています。また，ネットワークの発展に伴うコンピュータのメディア化が我々の情報処理能力の限界を越えて進もうとしています（ニュースとメールを読み，ファイル転送してmakeすると１日が終わってしまう!?）。こうしてあふれた情報ゴミがユーザのオタク化を促進しています。いわゆる「ソフトウェア危機」＝ソフトウェア生産場面での需給ギャップに対して，これらはソフトウェア消費場面での需給ギャップといえるかもしれません。これを，「第二種ソフトウェア危機」と呼ぶことにします。第二種ソフトウェア危機は人間の基本的な欲求に密接に関係しているため，より本質的で深刻です。これを回避するための一助として，研究室，教室あるいは小規模大学のユーザグループ（＝マクロユーザ）に対する適切な情報提供，事例紹介や講習等のサポートをセンターにお願いしたいと考えています。（#でも難しいですよね。中途半端な結論で終わってしまった・・・）
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
(注) FM-TOWNSは，日本で最初に発売されたCD-ROMドライブ付きパソコンだった。情報処理センターの副センター長だった定金晃三先生らと，これいいよね（デジタル教材メディア活用のはしりだぁ），とかいいながらネットワーク整備予算を使って一教室分導入したものである。そのうちビジネス用のPCにもCD-ROMがつくのが当たり前になってしまった。そして，これからはCD-ROMドライブが付かないのものが普通になるようだ。
2019年3月10日日曜日
JuliaのDoubleFloats
Juliaの任意精度実数としてはBigFloatがあるが，拡張精度として，DoubleFloatsのモジュールを見つけた。Double64はBigFloatより高速であるという説明があったので，ためしてみると，sqrtにバグがあって動かないような気がする。どうしたものか。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
using DoubleFloats
z=Double64(0.5)
#println((sqrt(z),cbrt(z)))
println((sin(z),cos(z),tan(z)))
#println((asin(z),acos(z),atan(z)))
println((csc(z),sec(z),cot(z)))
#println((acsc(1/z),asec(1/z),acot(1/z)))
println((sinh(z),cosh(z),tanh(z)))
#println((asinh(z),acosh(z),atanh(z)))
println((log(z),exp(z),log2(z),log10(z)))
sqrt(z)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
(0.479425538604203, 0.8775825618903728, 0.5463024898437905)
(2.085829642933488, 1.139493927324549, 1.830487721712452)
(0.5210953054937474, 1.1276259652063807, 0.46211715726000974)
(-0.6931471805599453, 1.6487212707001282, -1.0, -0.3010299956639812)
UndefVarError: v not defined
Stacktrace:
[1] mul_ddfp_dd at /Users/xxxxx/.julia/packages/DoubleFloats/obU7N/src/math/ops/op_ddfp_dd.jl:16 [inlined]
[2] sqrt_dd_dd(::Tuple{Float64,Float64}) at /Users/xxxxx/.julia/packages/DoubleFloats/obU7N/src/math/ops/op_dd_dd.jl:74
[3] sqrt_db_db at /Users/koshi/.julia/packages/DoubleFloats/obU7N/src/math/ops/op_db_db.jl:26 [inlined]
[4] sqrt(::DoubleFloat{Float64}) at /Users/xxxxx/.julia/packages/DoubleFloats/obU7N/src/math/ops/arith.jl:9
[5] top-level scope at In[1]:11
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
using DoubleFloats
z=Double64(0.5)
#println((sqrt(z),cbrt(z)))
println((sin(z),cos(z),tan(z)))
#println((asin(z),acos(z),atan(z)))
println((csc(z),sec(z),cot(z)))
#println((acsc(1/z),asec(1/z),acot(1/z)))
println((sinh(z),cosh(z),tanh(z)))
#println((asinh(z),acosh(z),atanh(z)))
println((log(z),exp(z),log2(z),log10(z)))
sqrt(z)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
(0.479425538604203, 0.8775825618903728, 0.5463024898437905)
(2.085829642933488, 1.139493927324549, 1.830487721712452)
(0.5210953054937474, 1.1276259652063807, 0.46211715726000974)
(-0.6931471805599453, 1.6487212707001282, -1.0, -0.3010299956639812)
UndefVarError: v not defined
Stacktrace:
[1] mul_ddfp_dd at /Users/xxxxx/.julia/packages/DoubleFloats/obU7N/src/math/ops/op_ddfp_dd.jl:16 [inlined]
[2] sqrt_dd_dd(::Tuple{Float64,Float64}) at /Users/xxxxx/.julia/packages/DoubleFloats/obU7N/src/math/ops/op_dd_dd.jl:74
[3] sqrt_db_db at /Users/koshi/.julia/packages/DoubleFloats/obU7N/src/math/ops/op_db_db.jl:26 [inlined]
[4] sqrt(::DoubleFloat{Float64}) at /Users/xxxxx/.julia/packages/DoubleFloats/obU7N/src/math/ops/arith.jl:9
[5] top-level scope at In[1]:11
2019年3月9日土曜日
Juliaの配列添字範囲
Juliaでは配列の添字が１から始まっている。Cなどでは配列の添字が０から始まっている。Fortranでは，配列の添字範囲を指定できた。そこで，Juliaでも配列の添字範囲を指定したい。
調べてみると，OffsetArrays.jl というモジュールがあった。
using Pkg
Pkg.add("OffsetArrays")
using OffsetArrays
とすればよい。
使用例は以下の通り。昔作ったFortranプログラムをJuliaにしたもの。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
setprecision(BigFloat,2^6)
o=BigFloat(0)
e=BigFloat(1)
d=BigFloat(2)
s=e/d
N = 300
fc = OffsetArray{BigFloat}(undef,-N:N);
fd = OffsetArray{BigFloat}(undef,-N:N);
fj = OffsetArray{BigFloat}(undef,-N:N);
bn = OffsetArray{BigFloat}(undef,-N:N,-N:N);
function factorials(fc,fd,fj,bn,n)
# k! = k*(k-1)*...*2*1
# fc[k] = sqrt(k!); fc[-k] = k!
# k!! = k*(k-2)*...*3*1 or k*(k-2)*...*4*2
# fd[k] = sqrt((2*k-1)!!/2^k); fd[-k] = fd[k]^2
# fj[k] = sqrt(k+1); fj[-k] = 1/sqrt(k+1)
#
fc[0]=e
fd[0]=e
fj[0]=e
for k in 1:n
kk=BigFloat(k)
fc[k]=sqrt(kk)*fc[k-1]
fc[-k]=kk*fc[-1+k]
fd[k]=sqrt(kk-s)*fd[k-1]
fd[-k]=(kk-s)*fd[-1+k]
fj[k]=sqrt(kk+e)
fj[-k]=e/sqrt(kk+e)
end
bn[0,0]=e
for k in 1:n
bn[k,0]=e
bn[0,k]=e
bn[k,k]=e
pd=e
for l in 1:div(k,2)
j=k-l
pd=pd*BigFloat(j+1)/BigFloat(l)
sd=sqrt(pd)
bn[k,l]=sd
bn[k,j]=sd
bn[l,k]=pd
bn[j,k]=pd
end
end
println(fc[n])
println(fc[-10*div(n,10)])
println((bn[45,90]))
end
@time factorials(fc,fd,fj,bn,300)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
調べてみると，OffsetArrays.jl というモジュールがあった。
using Pkg
Pkg.add("OffsetArrays")
using OffsetArrays
とすればよい。
使用例は以下の通り。昔作ったFortranプログラムをJuliaにしたもの。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
setprecision(BigFloat,2^6)
o=BigFloat(0)
e=BigFloat(1)
d=BigFloat(2)
s=e/d
N = 300
fc = OffsetArray{BigFloat}(undef,-N:N);
fd = OffsetArray{BigFloat}(undef,-N:N);
fj = OffsetArray{BigFloat}(undef,-N:N);
bn = OffsetArray{BigFloat}(undef,-N:N,-N:N);
function factorials(fc,fd,fj,bn,n)
# k! = k*(k-1)*...*2*1
# fc[k] = sqrt(k!); fc[-k] = k!
# k!! = k*(k-2)*...*3*1 or k*(k-2)*...*4*2
# fd[k] = sqrt((2*k-1)!!/2^k); fd[-k] = fd[k]^2
# fj[k] = sqrt(k+1); fj[-k] = 1/sqrt(k+1)
#
fc[0]=e
fd[0]=e
fj[0]=e
for k in 1:n
kk=BigFloat(k)
fc[k]=sqrt(kk)*fc[k-1]
fc[-k]=kk*fc[-1+k]
fd[k]=sqrt(kk-s)*fd[k-1]
fd[-k]=(kk-s)*fd[-1+k]
fj[k]=sqrt(kk+e)
fj[-k]=e/sqrt(kk+e)
end
bn[0,0]=e
for k in 1:n
bn[k,0]=e
bn[0,k]=e
bn[k,k]=e
pd=e
for l in 1:div(k,2)
j=k-l
pd=pd*BigFloat(j+1)/BigFloat(l)
sd=sqrt(pd)
bn[k,l]=sd
bn[k,j]=sd
bn[l,k]=pd
bn[j,k]=pd
end
end
println(fc[n])
println(fc[-10*div(n,10)])
println((bn[45,90]))
end
@time factorials(fc,fd,fj,bn,300)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1.74944994845934542383e+307
3.03013619603034000008e+308
1.0382742128755341137e+26
0.732669 seconds (1.29 M allocations: 58.421 MiB, 10.08% gc time)
2019年3月8日金曜日
退職挨拶の練習（１）
理科教育講座の越桐です。1982年に本学に着任してこの３月末で37年になります。
本学に在席していた時間を２つに分けると，20世紀が19年，21世紀が18年ということで，20世紀のほうが長い旧世紀型人間のひとりなので，最近の社会や大学の変化の速さにはとてもついていけない今日この頃です。
さて，退職の挨拶とは，これまで言えなかったことを告白する懺悔の時間のようなものなので，１つだけ紹介させてください。昔の天王寺分校にはデータステーションという，情報処理センターの前身がありました。雨の日には雨漏りするのでコンピュータにシートをかけていました。２階の私の研究室から徒歩30秒ところです。
インターネットの普及の直前の1991年に，本学ははじめてインターネットにつながりました。当時，キャンパスネットワークはないので，利用するにはデータステーションまで行って，コンピュータを使う必要がありました。
そこで，データステーションの垣本先生と結託して，キャンパスネットワークの実験をしようということで，日本橋（にっぽんばし）までいって，私費で10BASE2のケーブルを30メートル買い，研究室のNeXTステーションにヤミで接続したのは私です。ごめんなさい。これが，本学におけるキャンパスＬＡＮのはしりでした。
今では無線LANがどこでも使える夢の世界が実現したのですが，いいのか悪いのか・・・
そんなわけで，先生方，職員の皆様方には長らく大変お世話になりました。大学はこれから大変な時代を迎えますが，どうかお身体にお気をつけ，皆さん仲たがいせずにこの難関に立ち向かっていただければと思います。本日はどうもありがとうございました。
本学に在席していた時間を２つに分けると，20世紀が19年，21世紀が18年ということで，20世紀のほうが長い旧世紀型人間のひとりなので，最近の社会や大学の変化の速さにはとてもついていけない今日この頃です。
さて，退職の挨拶とは，これまで言えなかったことを告白する懺悔の時間のようなものなので，１つだけ紹介させてください。昔の天王寺分校にはデータステーションという，情報処理センターの前身がありました。雨の日には雨漏りするのでコンピュータにシートをかけていました。２階の私の研究室から徒歩30秒ところです。
インターネットの普及の直前の1991年に，本学ははじめてインターネットにつながりました。当時，キャンパスネットワークはないので，利用するにはデータステーションまで行って，コンピュータを使う必要がありました。
そこで，データステーションの垣本先生と結託して，キャンパスネットワークの実験をしようということで，日本橋（にっぽんばし）までいって，私費で10BASE2のケーブルを30メートル買い，研究室のNeXTステーションにヤミで接続したのは私です。ごめんなさい。これが，本学におけるキャンパスＬＡＮのはしりでした。
今では無線LANがどこでも使える夢の世界が実現したのですが，いいのか悪いのか・・・
そんなわけで，先生方，職員の皆様方には長らく大変お世話になりました。大学はこれから大変な時代を迎えますが，どうかお身体にお気をつけ，皆さん仲たがいせずにこの難関に立ち向かっていただければと思います。本日はどうもありがとうございました。
2019年3月7日木曜日
Juliaの数学的定数と精度
（ネイピア数の計算からの続き）
円周率やネイピア数の計算で，JuliaのBigFloatの精度が$10^{-77}$であったが，これは，
$2^{-255}$に対応していた。すなわち，BigFloatの精度は次のようになる。
\begin{equation}
{\rm eps(BigFloat)}=2^{-255}=1.7272337110188889250772703725600\\
79914223200072887256277004740694033718360632485\ {\rm e}^{-77}
\end{equation}
円周率πやネイピア数eの有効桁数がこれによって決まるが，BigFloatの精度は制御することができる。このためには，setprecision(BigFloat,n) とすればよい。デフォルト値はn=256に対応しており，m = log10(BigInt(2)^(-n+1)) から，10進法における有効桁数mが求まることになる。
その様子は次のプログラムでわかる。Juliaの数学的定数として，円周率πやネイピア数ℯに加えて，カタラン数，オイラーのガンマγ，黄金比φなどが Base.MathConstants に準備されている。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
function cnst(n)
setprecision(BigFloat,2^n)
println(" ")
println("- - - - - - - - - - - - - - - - - - - - - - - - - - - -")
println(2^n," eps=",eps(BigFloat))
println("- - - - - - - - - - - - - - - - - - - - - - - - - - - -")
println(" π=",BigFloat(Base.MathConstants.π))
println(" =",BigFloat(Base.MathConstants.ℯ))
println(" c=",BigFloat(Base.MathConstants.catalan))
println(" γ=",BigFloat(Base.MathConstants.γ))
println(" φ=",BigFloat(Base.MathConstants.φ))
end
for k in 4:2:10
@time cnst(k)
end
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - - - - - - - - -
16 eps=3.05176e-05
- - - - - - - - - -
π=3.1416
ℯ=2.71826
c=0.91597
γ=0.577209
φ=1.61804
0.000738 seconds (475 allocations: 36.148 KiB)
- - - - - - - - - -
64 eps=1.08420217248550443401e-19
- - - - - - - - - -
π=3.14159265358979323851
ℯ=2.71828182845904523543
c=0.915965594177219015048
γ=0.577215664901532860616
φ=1.61803398874989484821
0.000528 seconds (471 allocations: 37.184 KiB)
- - - - - - - - - -
256 eps=1.727233711018888925077270372560079914223200072887256277004740694033718360632485e-77
- - - - - - - - - -
π=3.141592653589793238462643383279502884197169399375105820974944592307816406286198
ℯ=2.718281828459045235360287471352662497757247093699959574966967627724076630353555
c=0.9159655941772190150546035149323841107741493742816721342664981196217630197762584
γ=0.5772156649015328606065120900824024310421593359399235988057672348848677267776685
φ=1.61803398874989484820458683436563811772030917980576286213544862270526046281891
0.000546 seconds (471 allocations: 40.613 KiB)
- - - - - - - - - -
1024 eps=1.112536929253600691545116358666202032109607990231165915276663708443602217406959097927141579506255510282033669865517905502576217080776730054428006192688859410565388996766001165239805073721291818035960782523471251867104187625403325308329079474360245589984295819824250317954385059152437399890443876874974725790226e-308
- - - - - - - - - -
π=3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724586997
ℯ=2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427427466391932003059921817413596629043572900334295260595630738132328627943490763233829880753195251019011573834187930702154089149934884167509244761460668082264800168477411853742345442437107539077744992069551702766
c=0.9159655941772190150546035149323841107741493742816721342664981196217630197762547694793565129261151062485744226191961995790358988033258590594315947374811584069953320287733194605190387274781640878659090247064841521630002287276409423882599577415088163974702524820115607076448838078733704899008647751132259971343386
γ=0.5772156649015328606065120900824024310421593359399235988057672348848677267776646709369470632917467495146314472498070824809605040144865428362241739976449235362535003337429373377376739427925952582470949160087352039481656708532331517766115286211995015079847937450857057400299213547861466940296043254215190587755369
φ=1.618033988749894848204586834365638117720309179805762862135448622705260462818902449707207204189391137484754088075386891752126633862223536931793180060766726354433389086595939582905638322661319928290267880675208766892501711696207032221043216269548626296313614438149758701220340805887954454749246185695364864449242
0.000710 seconds (481 allocations: 56.924 KiB)
円周率やネイピア数の計算で，JuliaのBigFloatの精度が$10^{-77}$であったが，これは，
$2^{-255}$に対応していた。すなわち，BigFloatの精度は次のようになる。
\begin{equation}
{\rm eps(BigFloat)}=2^{-255}=1.7272337110188889250772703725600\\
79914223200072887256277004740694033718360632485\ {\rm e}^{-77}
\end{equation}
円周率πやネイピア数eの有効桁数がこれによって決まるが，BigFloatの精度は制御することができる。このためには，setprecision(BigFloat,n) とすればよい。デフォルト値はn=256に対応しており，m = log10(BigInt(2)^(-n+1)) から，10進法における有効桁数mが求まることになる。
その様子は次のプログラムでわかる。Juliaの数学的定数として，円周率πやネイピア数ℯに加えて，カタラン数，オイラーのガンマγ，黄金比φなどが Base.MathConstants に準備されている。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
function cnst(n)
setprecision(BigFloat,2^n)
println(" ")
println("- - - - - - - - - - - - - - - - - - - - - - - - - - - -")
println(2^n," eps=",eps(BigFloat))
println("- - - - - - - - - - - - - - - - - - - - - - - - - - - -")
println(" π=",BigFloat(Base.MathConstants.π))
println(" =",BigFloat(Base.MathConstants.ℯ))
println(" c=",BigFloat(Base.MathConstants.catalan))
println(" γ=",BigFloat(Base.MathConstants.γ))
println(" φ=",BigFloat(Base.MathConstants.φ))
end
for k in 4:2:10
@time cnst(k)
end
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - - - - - - - - -
16 eps=3.05176e-05
- - - - - - - - - -
π=3.1416
ℯ=2.71826
c=0.91597
γ=0.577209
φ=1.61804
0.000738 seconds (475 allocations: 36.148 KiB)
- - - - - - - - - -
64 eps=1.08420217248550443401e-19
- - - - - - - - - -
π=3.14159265358979323851
ℯ=2.71828182845904523543
c=0.915965594177219015048
γ=0.577215664901532860616
φ=1.61803398874989484821
0.000528 seconds (471 allocations: 37.184 KiB)
- - - - - - - - - -
256 eps=1.727233711018888925077270372560079914223200072887256277004740694033718360632485e-77
- - - - - - - - - -
π=3.141592653589793238462643383279502884197169399375105820974944592307816406286198
ℯ=2.718281828459045235360287471352662497757247093699959574966967627724076630353555
c=0.9159655941772190150546035149323841107741493742816721342664981196217630197762584
γ=0.5772156649015328606065120900824024310421593359399235988057672348848677267776685
φ=1.61803398874989484820458683436563811772030917980576286213544862270526046281891
0.000546 seconds (471 allocations: 40.613 KiB)
- - - - - - - - - -
1024 eps=1.112536929253600691545116358666202032109607990231165915276663708443602217406959097927141579506255510282033669865517905502576217080776730054428006192688859410565388996766001165239805073721291818035960782523471251867104187625403325308329079474360245589984295819824250317954385059152437399890443876874974725790226e-308
- - - - - - - - - -
π=3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724586997
ℯ=2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427427466391932003059921817413596629043572900334295260595630738132328627943490763233829880753195251019011573834187930702154089149934884167509244761460668082264800168477411853742345442437107539077744992069551702766
c=0.9159655941772190150546035149323841107741493742816721342664981196217630197762547694793565129261151062485744226191961995790358988033258590594315947374811584069953320287733194605190387274781640878659090247064841521630002287276409423882599577415088163974702524820115607076448838078733704899008647751132259971343386
γ=0.5772156649015328606065120900824024310421593359399235988057672348848677267776646709369470632917467495146314472498070824809605040144865428362241739976449235362535003337429373377376739427925952582470949160087352039481656708532331517766115286211995015079847937450857057400299213547861466940296043254215190587755369
φ=1.618033988749894848204586834365638117720309179805762862135448622705260462818902449707207204189391137484754088075386891752126633862223536931793180060766726354433389086595939582905638322661319928290267880675208766892501711696207032221043216269548626296313614438149758701220340805887954454749246185695364864449242
0.000710 seconds (481 allocations: 56.924 KiB)
2019年3月6日水曜日
ネイピア数の計算
（円周率の計算からの続き）
ネイピア数は自然対数の底 $e$ のこと。欧米では Euler's Number とよばれることが多いとあるが，Eulerが関わる数はたくさんあってややこしいので，ここではネイピア数とする。
円周率の計算方法というか表現方法はたくさんあって，その中には与えられた級数の少数の項で良い精度の近似を与えるものがあった。ところで，ネイピア数の表式をざっと調べても，なんだかぱっとしないのである。
\begin{equation}
\begin{aligned}
e^x &= 1 + x + \frac{x^2}{2!}+\frac{x^3}{3!}+\frac{x^4}{4!}+\cdots\\
e &= 1 + 1 + \frac{1}{2!}+\frac{1}{3!}+\frac{1}{4!}+\cdots
\end{aligned}
\end{equation}
を多少変形したものか，連分数で表現したものはあるけれど，ラマヌジャンの公式のようにはっとするものが見当たらない。単に調査不足だけかもしれないし，そもそも指数関数を微分や積分しても元に戻るだけなのが原因かもしれないが，よくわからないのだった。
５項くらいで精度が70桁程度になる公式はないものなのだろうか？Improving the Convergence of Newton's Series Approximation for $e$ とかあったけど，目的のものではない。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
f=ones(BigFloat,1000)
function fact(f,n)
# calculation of factorial 1!..n!
for i in 2:n
f[i]=BigFloat(i)*f[i-1]
end
# println(f[n])
end
function cf(x,a,b,n)
# continued fraction
for k = n:-1:1
x = a[k] + b[k]/x
end
return x
end
function napier1(n)
# by series expansion
o=BigFloat(1)
sum=o
for k in 1:n
sum=sum+o/f[k]
end
println(sum-exp(o))
end
function napier2(n)
# by series expansion
o=BigFloat(1)
sum=o
for k in 1:n
sum=sum+o/((BigFloat(4)^k)*f[k])
end
println(sum*sum*sum*sum-exp(o))
end
function napier3(n)
# by continued fraction
o=BigFloat(1)
z=BigFloat(2)
a=ones(BigFloat,1000)
for k in 1:1000
a[k]=z*(z*k-z-o)
end
a[1]=o
a[2]=o
b=ones(BigFloat,1000)
b[1]=z
println(cf(z,a,b,n)-exp(o))
end
function napier4(n)
# sqrt(e) by series expansion
o=BigFloat(1)
z=BigFloat(2)
sum=BigFloat(0)
for k in 0:n
sum=sum+BigFloat(4*k+3)/(z^(2*k+1)*f[2*k+1])
end
println(sum*sum-exp(o))
end
fact(f,160)
@time napier1(57)
@time napier2(41)
@time napier3(26)
@time napier4(24)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-1.381786968815111140061816298048063931378560058309805021603792555226974688505988e-76
0.018460 seconds (30.73 k allocations: 1.488 MiB)
-2.763573937630222280123632596096127862757120116619610043207585110453949377011977e-76
0.029598 seconds (49.70 k allocations: 2.358 MiB)
-1.381786968815111140061816298048063931378560058309805021603792555226974688505988e-76
0.054272 seconds (100.99 k allocations: 5.090 MiB)
1.727233711018888925077270372560079914223200072887256277004740694033718360632485e-76
0.032685 seconds (72.55 k allocations: 3.431 MiB)
（Juliaの数学的定数と精度へ続く）
ネイピア数は自然対数の底 $e$ のこと。欧米では Euler's Number とよばれることが多いとあるが，Eulerが関わる数はたくさんあってややこしいので，ここではネイピア数とする。
円周率の計算方法というか表現方法はたくさんあって，その中には与えられた級数の少数の項で良い精度の近似を与えるものがあった。ところで，ネイピア数の表式をざっと調べても，なんだかぱっとしないのである。
\begin{equation}
\begin{aligned}
e^x &= 1 + x + \frac{x^2}{2!}+\frac{x^3}{3!}+\frac{x^4}{4!}+\cdots\\
e &= 1 + 1 + \frac{1}{2!}+\frac{1}{3!}+\frac{1}{4!}+\cdots
\end{aligned}
\end{equation}
を多少変形したものか，連分数で表現したものはあるけれど，ラマヌジャンの公式のようにはっとするものが見当たらない。単に調査不足だけかもしれないし，そもそも指数関数を微分や積分しても元に戻るだけなのが原因かもしれないが，よくわからないのだった。
５項くらいで精度が70桁程度になる公式はないものなのだろうか？Improving the Convergence of Newton's Series Approximation for $e$ とかあったけど，目的のものではない。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
f=ones(BigFloat,1000)
function fact(f,n)
# calculation of factorial 1!..n!
for i in 2:n
f[i]=BigFloat(i)*f[i-1]
end
# println(f[n])
end
function cf(x,a,b,n)
# continued fraction
for k = n:-1:1
x = a[k] + b[k]/x
end
return x
end
function napier1(n)
# by series expansion
o=BigFloat(1)
sum=o
for k in 1:n
sum=sum+o/f[k]
end
println(sum-exp(o))
end
function napier2(n)
# by series expansion
o=BigFloat(1)
sum=o
for k in 1:n
sum=sum+o/((BigFloat(4)^k)*f[k])
end
println(sum*sum*sum*sum-exp(o))
end
function napier3(n)
# by continued fraction
o=BigFloat(1)
z=BigFloat(2)
a=ones(BigFloat,1000)
for k in 1:1000
a[k]=z*(z*k-z-o)
end
a[1]=o
a[2]=o
b=ones(BigFloat,1000)
b[1]=z
println(cf(z,a,b,n)-exp(o))
end
function napier4(n)
# sqrt(e) by series expansion
o=BigFloat(1)
z=BigFloat(2)
sum=BigFloat(0)
for k in 0:n
sum=sum+BigFloat(4*k+3)/(z^(2*k+1)*f[2*k+1])
end
println(sum*sum-exp(o))
end
fact(f,160)
@time napier1(57)
@time napier2(41)
@time napier3(26)
@time napier4(24)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-1.381786968815111140061816298048063931378560058309805021603792555226974688505988e-76
0.018460 seconds (30.73 k allocations: 1.488 MiB)
-2.763573937630222280123632596096127862757120116619610043207585110453949377011977e-76
0.029598 seconds (49.70 k allocations: 2.358 MiB)
-1.381786968815111140061816298048063931378560058309805021603792555226974688505988e-76
0.054272 seconds (100.99 k allocations: 5.090 MiB)
1.727233711018888925077270372560079914223200072887256277004740694033718360632485e-76
0.032685 seconds (72.55 k allocations: 3.431 MiB)
（Juliaの数学的定数と精度へ続く）
登録:
投稿 (Atom)