覚道先生（５）

（覚道先生（４）からの続き）

覚道先生が定年を迎えられて数年後の2001年の３月。当時夜間学部に勤務されていた木立先生から連絡があり（このときはすでに夜間学部主事だったのだろうか，まだ長尾彰夫先生が主事の頃か），明日の卒業式で記念講演される覚道先生の紹介が必要なので，物理方面の業績について一言考えてほしいとのこと。

そのとき，私の送ったメールの記録がある。

覚道先生のご専門は素粒子物理学ですが，初期には「高エネルギーの核子衝突における中間子の多重発生」，また「相対論的なクォークモデルにおけるハドロンの質量の導出」，あるいは「ゲージ相互作用における不変性の問題」など，現象論から基礎理論に至るまでの幅広い研究に携わってこられました。

この卒業式の記念講演には私も出席し，終わってから覚道先生，木立先生と私の３人で四天王寺裏のちゃんこ屋（萩家本場所）で久しぶりの飲み会となった。

私が着任した頃は，前述の館野常昭さんやその当時大阪教育大学の大学院生で，第１研究室の学生さんたちと同室であった原田孝君らとともに，ディラック粒子２個の２体ポテンシャル問題を解いて中間子の質量スペクトルを求めるという仕事をされていた。

覚道先生と私は，スピンが０でない原子核のミューオン捕獲の理論（ミューオン原子の超微細状態からの捕獲と誘導擬スカラー項の導出）について共著の論文を書いたこともあった。

（覚道先生（６）へ続く）

覚道先生（４）

（覚道先生（３）からの続き）

覚道先生は1991年度と1992年度の２年間夜間学部の主事を務められた。そのころはちょうど大学の自己点検評価がうるさくいわれるようになった時代で，覚道先生も自己評価書の作成にたいへんご尽力され，また苦労されていたようであった。1993年には前年の教養学科に続いて教員養成課程も柏原キャンパスに移転したので，最後の１年は研究室も離れ離れになってしまった。

1994年の３月に卒業生やＯＢが集まって大阪府教職員互助組合のたかつガーデンで覚道先生の退官記念パーティを催した。これは先生の最初の学生さんで，大分大学を経て大阪教育大学の理科教育学教室の教授になられた村井護晏先生や，当時覚道先生と共同研究されていた清教学園高等学校の館野常昭先生らが幹事となり，私もお手伝いして開いたものである。３０人くらいの卒業生が集まっただろうか。この日のために幹事の３人で上本町の近鉄百貨店で記念品を探し回ったのものだった。

大阪教育大学物理学教室主催の退職記念お食事会は，天王寺の近鉄百貨店の隣の天王寺都ホテル新館で開いた。こちらの方は私が幹事であった。当日はお花の準備で天王寺都ホテルの日比谷花壇あたりを走り回っていた。いずれも覚道先生と奥様のお二人をご招待している。

覚道先生は定年後しばらく福井工業大学で物理を教えておられた。ここはその後，下村先生や阪大の原子核実験の南園忠則先生なども定年後に勤められている。覚道先生は，通うのがたいへんだけど芦原温泉も近くでなかなかよろしいとおっしゃっていた。

（覚道先生（５）に続く）

覚道先生（３）

（覚道先生（２）からの続き）

天王寺の１研時代にハイキングが企画されたことがある。木立先生が，５月の連休に葛城山のつつじを見に行きましょうとおっしゃって，覚道先生と木立先生と私の三人で御所からケーブルカーで葛城山頂に向かった。学生さんたちにも声をかけていたのだが，当日はどんよりとした曇り日で今にも雨が降りそうな天気。若者は誰もこないのであった。雨が降らないかばかりを気にしていて，山頂からツツジが見えたのかどうか確かな記憶がない。

帰りは山道を北に縦走して，近鉄南大阪線の磐城駅に出て帰った。このときだろうか，木立先生が珍しいものがあるといって，天王寺（またあべの銀座の近くである）の屋台の洋食屋で山羊だか羊だかの脳みそのフライなるものを食べた。うまいでしょう，といわれたがどうしたものか。

大阪教育大学が1993年に柏原キャンパスに統合移転完了した前後に，覚道先生といっしょのハイキングで私が企画したものがある。ご病気になられてからの回復祝いだったのか，1991年から夜間学部主事になられており，その就任のお祝いなのか慰労会なのか，すっかり忘れている。亡くなられた高木義人先生から私に物理教室でハイキングを企画したらどうかとのお話があった。当時すでに奈良県民だったので，大和三山めぐりというコースを考えた。大和八木−耳成山−藤原宮跡−天香久山−畝傍山−神武天皇陵−大和八木とまわるのだ。事前に家族で予習ハイキングをした，子どもたちはまだ小さかったがなんとかクリアできた。さて，教室のハイキングの方は，私が最若年なのにフーフーいっているところ，覚道先生をはじめとして高齢の先生方は皆さん健脚でいたってお元気なのである。終わってから大和八木の居酒屋でちょっと一杯。

（覚道先生（４）に続く）

覚道先生（２）

（覚道先生（１）からの続き）

1982年の４月１日に大阪教育大学に教務員として着任し，その次の年に助手になった。覚道先生の研究室は隣の部屋にあり，鯖田秀樹先生（池田分校所属）の机の奥に座っていらっしゃった。その手前には私より１年先に着任して助手になられたばかりの木立英行先生の机がある。私と併せて４人が物理学第１研究室（通称１研）の教員であった。覚道先生が阪大の素粒子論，私が阪大の原子核理論，木立先生と鯖田先生は京大の非線形動力学研究室（今はＳＮＳでおなじみの佐々真一先生が主宰しているのか）の出身である。次の年には空いた教務員のポストに同じ研究室から古賀眞史先生がこられた。古賀先生はプラズマ物理がご専門の下村昇先生と島田昌敏先生（放電物理・物理教育）の研究室（５研）に入った。

１研の部屋は昔の師範学校時代の小学校の教室の１部屋を２つに変形間仕切りしたものであり，私の方は１研で卒業研究をする学生さんたちと同居なのだった。詳細は，1987年度卒業の齊藤昌久君の書いたターボパスカルマニュアル（ロールプレイングゲーム風）の図をみていただきたい。

天王寺分校２Ｆ南の物理１研のイメージ

図の右が南で天王寺分校と附属中・高の共用グラウンドに面しており，左が廊下で出入り口は別れている。農民たちの村には学生の机があり，当時はどこから持ってきたのやらベッドまで設置されていた。覚道先生や鯖田先生や木立先生の学生もこの部屋にいるので結構な人数であった。南に面していて暖かい部屋だったので覚道先生はよくウツラウツラされていた。自分がこの歳になってようやく，暖かくなくとも毎日の昼食後には必ず眠くなってしまうことが分かった。

私が卒論学生を持ったのは着任してから４年目の1985年からである。覚道先生は1986年に教授になられてまもなく，1988年の教養学科設置の改組のころから二部（夜間学部）の担当になられた。それまでは覚道先生といっしょに卒論生のゼミを行うことがあった。十分予習していない学生が覚道先生の鋭い質問に答えられないと沈黙の時間が10分くらい続く・・・。というわけでゼミではかなり厳しい一面があった。

（覚道先生（３）に続く）

覚道先生（１）

先日，覚道雄次郎先生が昨年の12月に亡くなられたとの報せがあった。私より25歳ほど年上なので90歳くらいだったろうか。阪大理学部の内山研を出られて和歌山大学に数年勤められた後，1966年から大阪教育大学に移られた。覚道先生には本当にお世話になった。

私が大阪教育大学に就職できたのも覚道先生のおかげなのであった。Ｄ３の終わり頃に初めて話があったときに，阪大豊中キャンパスの教養部の山本邦夫先生（覚道先生の同級生）の部屋へ訪ねていってご挨拶をした。日頃着たことのないスーツにネクタイでこちらはかなり緊張していたのだが，優しい声をかけていただきちょっと安心した。ちょっとまだどうなるかわかりませんねぇ，というような雰囲気だったので，冬の中央環状線の横の歩道をとぼとぼと下って蛍池まで帰ったイメージが残っている。

その年は残念ながら就職できなかったのだが，次の年に再び大阪教育大学の物理教室の人事があり（阪田巻蔵先生の学長就任の前後だったような），このときはたいへん幸運なことに採用されることになった。当時の大阪教育大学の物理学教室は天王寺分校と池田分校を併せて15名の先生方がいらっしゃって，なかなかの大所帯であった（いまは３部局あわせてその半分以下だが，今後よりいっそうの削減が予想される）。

採用が決まる最初の面接の後か，あるいは着任前の別の機会だったか忘れたが，暖い春の日のお昼時。おいしい麺を食べに行きませんかと誘われ，天王寺の今はなきあべの銀座の近くの小さなラーメン屋でじゃじゃ麺をごちそうになった。これがうまいんですよとの講釈があったが，こちらは緊張しているのでなんだか味もよくわからなかったかもしれない。

（覚道先生（２）に続く）

草枕：夏目漱石

夏目漱石の草枕は読んだことがない。冒頭の一節は有名なのだけれど。青空文庫の草枕からその部分を引用してみると，

　山路を登りながら、こう考えた。
　智に働けば角が立つ。情に棹させば流される。意地を通せば窮屈だ。とかくに人の世は住みにくい。
　住みにくさが高じると、安い所へ引き越したくなる。どこへ越しても住みにくいと悟った時、詩が生れて、画が出来る。
　人の世を作ったものは神でもなければ鬼でもない。やはり向う三軒両隣にちらちらするただの人である。ただの人が作った人の世が住みにくいからとて、越す国はあるまい。あれば人でなしの国へ行くばかりだ。人でなしの国は人の世よりもなお住みにくかろう。
　越す事のならぬ世が住みにくければ、住みにくい所をどれほどか、寛容て、束の間の命を、束の間でも住みよくせねばならぬ。ここに詩人という天職が出来て、ここに画家という使命が降る。あらゆる芸術の士は人の世を長閑にし、人の心を豊かにするが故に尊い。

草枕は普通の小説ではないようで，物語が進まない気がして，三四郎・それから・門・彼岸過迄・行人・こころ・明暗などの延長線上のテンションでは，読む気力が出なかったのだと思う。この度，youtubeで草枕の朗読をざっくりと聞いてみた。前編と後編をあわせて５時間半ほどある。

蘆雪や若冲の話が出てきたので，ああ，そんな話だったのかと思ってなんだか親しみが湧いてきた。その若冲の部分を青空文庫から引用すると次のようになっている。

横を向く。床にかかっている若冲の鶴の図が目につく。これは商売柄だけに、部屋に這入った時、すでに逸品と認めた。若冲の図は大抵精緻な彩色ものが多いが、この鶴は世間に気兼なしの一筆がきで、一本足ですらりと立った上に、卵形の胴がふわっと乗かっている様子は、はなはだ吾意を得て、飄逸の趣は、長い嘴のさきまで籠っている。

ところで，松岡正剛の千夜千冊の583夜が「草枕」である。そこには，「そんなことを思いながら横に目を凝らすと若冲の鶏がいる。複製の商品だ。」とある。若冲といえば鶏というのは普通はそうだが，ここは精緻な彩色の鶏ではなくて，一筆書きの卵形の鶴でなければならないのである。複製の商品というのはどこからでてきたのだ。まあ，そんな瑣末なことはどうでもよいが，千夜千冊の記事の中で松岡正剛は漱石を巡る重要な指摘をしているように思った。

P. S. 草枕の那美のモデルが，前田卓（つな）なのだそうで，ここからも歴史がつながる。

「春立つや子規より手紙漱石へ」　（榎本好宏　1937-）

統計の日の標語

総務省が2019年度の統計の日の標語を募集している。募集期間は２月１日から３月３１日まで。総務省統計局のページはいたって冷静なのでこの話題はでてこない。担当が違うのであろう。

このタイミングでこの募集なので，twitterでは大喜利状態が発生している。そこでさっそくその中から川柳形式のものを選んでみたところ，３日でこのありさまである。募集要領では誰でも応募できるが１人５編以内であり，自作で未発表のものに限れられるため，これらの作品が選ばれる可能性は低い（そうでなくてもね…orz）。

一方，２月１日には柴山文部科学大臣が，大学改革案としてデータサイエンス必修化を打ち出している。もう浮き足立って支離滅裂な日本なのだった。初等中等教育から特別な教科の道徳もプログラミング教育も外し，大学改革案からデータサイエンスを外し，これらすべてを官僚と政治家の必須科目として設定したほうがよいのではないか。

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
統計の日の標語募集によせて（twitterから36選）

・統計に　隠れて見えぬ　世の姿
・国家への　信頼無視した　その統計
・マイナスが　プラスに見える　政府かな
・偽統計　国の破滅へ　まっしぐら
・為政者の　望み通りに　数いじる
・偽統計　みんなで作れば　怖くない
・偽データ　集めて早し　偽統計
・統計は　政府に任せちゃ　いけないよ
・統計は　いまや出世の　一里塚
・統計の　真の目的　出世かな
・統計と　言えば信じる　愚か者
・統計の　捏造改竄　誰のため
・統計が　民を誘う　夢の国
・忖度で　統計自在に　改竄し
・官邸の　顔色うかがい　手を入れる
・官邸の　意のままになす　数合わせ
・官邸の　為なら資料　ほっ統計〜
・大東亜　築く力だ　この改ざん
・統計の　不正で作れ　好景気
・統計の　操作が作る　独裁者
・統計を　疑う目から　民主主義
・統計を　とってくれろと　泣く子かな
・妄想に　あわせてつくれ　偽統計
・合わぬなら　作ってしまえ　偽統計
・ない数字　作ってしまえ　ホトトギス
・秀吉は　検地を始めた　非国民
・絶対に　組織ぐるみと　認めない
・計（かず）数え　頭鶏（あたまにわとり）　計（かず）忘れ
・騙せれば　調査方法　ほっ統計
・統計は　どうせ不正だ　ほっとうけい
・意図的に　抽出省略　記入ミス
・結論へ　鉛筆舐めて　数合わせ
・統計は　答えを先に　決めてから
・統計は　伸び縮みする　物差しだ
・統計の　せいでベクトル　数三Ｃ
・あの方の　笑顔がみたくて　この数字
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

能村登四郎

能村登四郎という名前やイメージが頭のどこかに引っかかっていたのは，NHKの俳句入門か何かで見かけたからかもしれない。その名前から，石川県にゆかりのある人かなと思ったら，やはり「我が父の生地」という前書きで「汗ばみて加賀強情の血ありけり」という句が「増殖する俳句歳時記（清水哲男）」にあった。

さらに調べると，Wikipediaの社会性俳句の項に，富山出身の沢木欣一（細見綾子の夫）と並んでその名前を見つけた。

能村登四郎の句にそれほど引かれるというわけではない。むしろ老いを突き付けられてしまうことに対する軽い拒否感があるのかもしれない。それでも清水哲男の文章がいいなと思ってしまうのが次の句の解説なので以下に引用する。

「老境や空ほたる籠朱房垂れ」　（能村登四郎　1911-2001）

老境にはいまだしの私が言うのも生意気だが、年齢を重ねるに連れて、たしかに事物は大きく鮮明に見えてくるのだと思う。事物の大小は相対的に感じるわけで、視線を活発に動かす若年時には、大小の区別は社会常識の範囲内に納まっている。ところが、身体的にも精神的にも目配りが不活発になってくる（活発に動かす必要もなくなる）と、相対化が徹底しないので、突然のように心はある一つのものを拡大してとらえるようになる。単なる細くて赤い紐が、大きな朱房に見えたりする。そのように目に見えるというよりも、そのように見たくなるというべきか。生理的な衰えとともに、人は事物を相対化せず、個としていつくしむ「レンズ」を育てていくようである。（清水哲男　June 18 2000）

初めての旅

数日前に，日本映画製作者連盟の新年記者発表があり，クイーンのボヘミアン・ラプソディが2018年度の興行収入トップで100億円を越えた，とテレビのニュースでやっていた。東映グループ会長ですっかりおっさんになってしまった岡田裕介が映っていた。

岡田裕介といえば，映画の「赤頭巾ちゃん気をつけて」の薫クンであり，同じ森谷司郎監督の「初めての旅」の主人公である。両作品ともに森和代が主人公の相手役として出演している。高校時代に米島誠二君が森和代を絶賛しており，今はなくなってしまった金沢劇場（東宝の配給館，ゴジラシリーズも日本のいちばん長い日もここで観たのだ）に一人で行ったのが48年前の高校２年の冬か。さらば青春を含む映画音楽は小椋佳だった（青春−砂漠の少年− 岡田裕介と森和代のナレーション入り）。そして，その小椋佳が東大出身の銀行員だとわかったときに，二人でかなりがっかりしていたのもなつかしい。

庄司薫の薫君シリーズでは「白鳥の歌なんか聞こえない」が，NHKの銀河テレビ小説でドラマ化され，主演が荒谷公之と仁科明子だったのだが，やはり岡田裕介−森和代ペアには勝てなかった。岡田裕介は70年代半ばには映画プロデューサに転じ，その後父親の岡田茂の後をついで東映の社長になる。

P. S. 何年か後に，「初めての旅」は曽野綾子が原作だということに気がついた。あの曽野綾子だ…orz。まあ，ストーリーから考えて “むべなるかな” なのであった（岡田裕介が岡田茂の息子であるという状況から考えると２乗で効いているのか，現実の権力関係は青春を凌駕する…orz）。

Flickrの短縮URL

Twitterで奥村晴彦さんがつぶやいている。これはたぶんPythonかなと思ったらそのようである。以下引用。

全然関係ないけどFlickrの短縮URLの理屈がわかった。Base58だ
n = 8348087410
a = "123456789abcdefghijkmnopqrstuvwxyzABCDEFGHJKLMNPQRSTUVWXYZ"
e = ""
while n != 0:
    e = a[n % 58] + e
    n //= 58
print("https://flic.kr/p/ " + e)

後学のために，Juliaに翻訳してみた。

n = 8348087410

a = "123456789abcdefghijkmnopqrstuvwxyzABCDEFGHJKLMNPQRSTUVWXYZ"

e = ""

while n != 0

    m = n % 58 +1

    e = a[m:m] * e

    n = n ÷ 58

end

print("flic.kr/p/" * e)

その過程で，Pythonの文法を調べようとしたところ，KEKの山本昇さんのPython入門（初級編）リリース0.1がひっかかった。山本昇さんは阪大理学部の森田研の隣の内山研に在席しており，私の２学年下かな。素粒子の理論の研究室では当時かなり珍しいことであったのだが，修論でコンピュータを使ったシミュレーションを行っていたのが印象深い。Lattice QCDなどがはやる前です。

さて，Juliaに翻訳する際の主な注意

(1) Juliaでは文字列定数は文字定数の配列になるので，データを１文字の文字列として取り出したい場合は，m:m のような範囲指定をする必要がある。

(2) Juliaでは文字列定数の配列は添字が１から始まる（Pythonでは0から）。そのために ｎを５８で割った剰余に１を加えている。

(3) Juliaでは文字列の結合には * を用いる（Pythonでは + ）。

(4) Juliaでは整数の除算には÷を用いる（Pythonでは // ）。

(5) JuliaではWhileなどの制御文は end で閉じる。字下げも不要である。

ついでに，逆変換のプログラムも考える。

(6) Juliaでは文字列eの頭を削るには strip(e, e[1])とする。

(7) Juliaでは文字列eの頭の文字を切り出すには first(e)とする。

(8) Juliaでは文字列aと文字cのマッチする添字を得るには findfirst(isequal(c),a)とする。

これらから，

function frickr(key)
  e = " "*key
  a = "123456789abcdefghijkmnopqrstuvwxyzABCDEFGHJKLMNPQRSTUVWXYZ"
  n = 0
  while e != ""
    n = n * 58
    e = strip(e,e[1]) 
    m = findfirst(isequal(first(e)),a)
    n = n + m - 1
#    println((m,n,e))
    if length(e)==1
      return n
    end
  end
end
    
frickr("dHG8xy")

8348087410

天牛書店

ニュースによると，天牛堺書店が経営破綻したとのこと。しばらく前に大阪市営地下鉄なんば駅の改札を出てすぐの ekimo なんば店が閉店していて残念だったのだが（と，調べてみると2016年9月・・・もう２年半も過ぎてしまったのか），このたび破産手続を開始したようだ。

所用で和泉中央駅に降りたときにも見かけたような気がする。あれはいったい何の用事だったのだろう。凄い勢いで記憶が失われている・・・orz（P. S. …桃山学院大学でした…）。というわけで，自分の活動範囲とはあまり重なっていない天牛堺書店を利用したことは少なかった。

天牛という名前がついているので，てっきり道頓堀にあった天牛書店の親戚か何かだと思っていた。天牛堺書店は，1963年に第１号店がオープンし，1968年にこの商号になった別会社で天牛書店とはまったく関係がないらしい。

天牛書店には，1972年に阪大に入学し大阪に住むようになってからたいへんお世話になった。当時は２〜３ヶ月に１度くらいは，梅田や桜橋を経由しながら，道頓堀から日本橋にかけての古本屋をハシゴするのが習慣になっていた。角座前店の方だろうか，奥に眼鏡をかけたおじいさんが鎮座していたが，あれが天牛新一郎さんだったのかもしれない。そのころの道頓堀は今に比べるとずっと落ち着いた（半分寂れかけたというのはいいすぎだろうが）繁華街であった。最近よくお世話になっている松竹座の他に，まだ，中座も角座も文楽の朝日座もあった時代だが，そちらの方にはあまり縁がなかった。朝日座に文楽ではなく映画を見に入ったように思う。

1970年代の後半に，出張で大阪に出てきた父に連れられて漫才を観たことがある。あれは梅田コマ劇場だったろうか，夢路いとし・喜美こいしの「もしもし鈴木さん（これでいろは四十七文字を大学生になってから覚えた）」とか「売り声（さおだけー）」がおもしろかった。梅田コマ劇場では，田村正和の眠狂四郎の円月殺法も。まあ，子どもの人生教育の一貫だったのかもしれない。

天牛書店は，その後，アメリカ村から緑地公園に移ってしまう。それぞれ１，２回は訪れたが，ちょっと遠いのだ。梅田の古本屋はカッパ横丁から紀伊國屋書店東横のうめ茶小路に移転し，理工古書の太田書店も入ったのは有り難いことだ。阪大豊中キャンパス下交差点近くの太田書店には石橋に住んでいる頃はたいへんお世話になった（学生時代は蛍池に住んでいたのでちょっと行きにくかったのだけれど）。

Juliaでパズル（５）

（Juliaでパズル（４）から続く）

結論：失敗でした。７時間以上かかるプログラムを高速化するために，重複計算部分を配列に落とし込もうとしたが，プログラムの見通しが非常に悪くなって，しかも時間も短縮できないという残念な結果に終わってしまった。

30年くらい前，スーパーコンピュータを使って原子核の殻模型計算を高速化するための試みについて，大阪大学の大型計算機センターニュースに投稿したことがある（スーパーコンピュータと原子核殻模型計算 1987）。導入されて間もなかったNECのベクトル型計算機の最初のモデルのSX-1を使い，$sd$殻の11粒子系である${}^{29}{\rm Si}$の基底状態（$J=\frac{1}{2}, T=\frac{1}{2}$，2360次元）を求める計算を高速化しようというものである（西村道明さんとの共同研究のための予備計算，大坪久夫先生にRCNPの計算費を確保していただいた）。

結果は，ベクトル化よりも重複計算を間接アドレッシングやメモリに蓄えて呼び出すことで避ける効果が高いというなさけないものであったが，その時から重複計算を避けようというのがスローガンとしてたたき込まれていたのであった（しかし，ベクトルマシンそのものが陳腐化して消えてしまうわけで，世の中の栄枯盛衰，諸行無常なり）。

今回の結果は，そもそも重複計算がそんなにあったのかというところからして問題であり，コードが複雑になって保守性にも欠けるという結末で終わってしまった。もともと重複だと思っている関数呼び出しも大したコストがかかっておらず，多次元の配列呼び出しとあまりかわらなかったのかもしれない。とりあえず，悪戦苦闘した（残念な）プログラムは以下の通り

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
function sqr(x)
  return sqrt(abs(x))
end

function cbr(x)
  return cbrt(abs(x))
end

function rog(x)
  return log(abs(x))
end

function cc(a,b)
  if (a==1 && b==1) || (a==11 && b==1)
    return 10*a+b
  else
    return 0
  end
end

function ne(op, op1)
  expr = Expr(:call, op, op1)
  return expr
end

function me(op, op1, op2)
  expr = Expr(:call, op, op1, op2)
  return expr
end

op1=[:+,:sqr,:cbr,:exp,:rog]
op2=[:+,:-,:*,:/,:cc]
ne1=Array{Any}(undef, 5)
ne2=Array{Any}(undef, 5, 5, 5)
ne11=Array{Any}(undef, 5, 5, 5, 5)
nex1=Array{Any}(undef, 5, 5, 5, 5, 5, 5, 5)
ne1x=Array{Any}(undef, 5, 5, 5, 5, 5, 5, 5)
mexy11=Array{Any}(undef, 5, 5, 5, 5, 5, 5, 5, 5, 5)
mey1x1=Array{Any}(undef, 5, 5, 5, 5, 5, 5, 5, 5, 5)
mey11x=Array{Any}(undef, 5, 5, 5, 5, 5, 5, 5, 5, 5)
me1yx1=Array{Any}(undef, 5, 5, 5, 5, 5, 5, 5, 5, 5)
me1y1x=Array{Any}(undef, 5, 5, 5, 5, 5, 5, 5, 5, 5)

function stor(ne1,ne11,ne1x,nex1,mexy11,mey1x1,mey11x,me1yx1,me1y1x,op1,op2)
  for i in 1:5
    ne1[i]=Expr(:call,op1[i],1)
  end
  for i in 1:5, j in 1:5, k in 1:5
    ne2[i,j,k]=Expr(:call,op2[i],ne1[j],ne1[k])
  end
  for i in 1:5, j in 1:5, k in 1:5, l in 1:5
    ne11[i,j,k,l]=Expr(:call,op1[i],ne2[j,k,l])
  end
  for i in 1:5, j in 1:5, k in 1:5, l in 1:5, m in 1:5, n in 1:5, o in 1:5
    nex1[i,j,k,l,m,n,o]=Expr(:call,op1[i],Expr(:call,op2[j],ne11[k,l,m,n],ne1[o]))
  end
  for i in 1:5, j in 1:5, k in 1:5, l in 1:5, m in 1:5, n in 1:5, o in 1:5
    ne1x[i,j,k,l,m,n,o]=Expr(:call,op1[i],Expr(:call,op2[j],ne1[k],ne11[l,m,n,o]))
  end
  for i in 1:5, j in 1:5, k in 1:5, l in 1:5, m in 1:5, n in 1:5, o in 1:5, p in 1:5, q in 1:5
    mexy11[i,j,k,l,m,n,o,p,q]=Expr(:call,op2[i],ne11[j,k,l,m],ne11[n,o,p,q])
  end
  for i in 1:5, j in 1:5, k in 1:5, l in 1:5, m in 1:5, n in 1:5, o in 1:5, p in 1:5, q in 1:5
    mey1x1[i,j,k,l,m,n,o,p,q]=Expr(:call,op2[i],nex1[j,k,l,m,n,o,p],ne1[q])
  end
  for i in 1:5, j in 1:5, k in 1:5, l in 1:5, m in 1:5, n in 1:5, o in 1:5, p in 1:5, q in 1:5
    mey11x[i,j,k,l,m,n,o,p,q]=Expr(:call,op2[i],ne1x[j,k,l,m,n,o,p],ne1[q])
  end
  for i in 1:5, j in 1:5, k in 1:5, l in 1:5, m in 1:5, n in 1:5, o in 1:5, p in 1:5, q in 1:5
    me1yx1[i,j,k,l,m,n,o,p,q]=Expr(:call,op2[i],ne1[j],nex1[k,l,m,n,o,p,q])
  end
  for i in 1:5, j in 1:5, k in 1:5, l in 1:5, m in 1:5, n in 1:5, o in 1:5, p in 1:5, q in 1:5
    me1y1x[i,j,k,l,m,n,o,p,q]=Expr(:call,op2[i],ne1[j],ne1x[k,l,m,n,o,p,q])
  end
end

@time stor(ne1,ne11,ne1x,nex1,mexy11,mey1x1,mey11x,me1yx1,me1y1x,op1,op2)
 
function doe1(a,b,c,d)
# ((1,1)_a,(1,1)_c)_b
#   x=ne11[d[5],a,d[1],d[2]]
#   y=ne11[d[6],c,d[3],d[4]]
    z=ne(op1[d[7]],mexy11[b,d[5],a,d[1],d[2],d[6],c,d[3],d[4]])
    if b==5
      z=:0
    end
    return (z,eval(z))
end

function doe2(a,b,c,d)
# (((1,1)_a,1)_b,1)_c
#   x=ne11[d[3],a,d[1],d[2]]
#   y=nex1[d[5],b,d[3],a,d[1],d[2],d[4]]
    z=ne(op1[d[7]],mey1x1[c,d[5],b,d[3],a,d[1],d[2],d[4],d[6]])
    if  (a!=5 && b==5) || ((a!=5 || b!=5) && c==5)
      z=:0
    end
    return (z,eval(z))
end

function doe3(a,b,c,d)
# ((a,(b,c)),d)
#   x=me(b,ne(d[1],1),ne(d[2],1))
#   y=ne1x[d[5],a,d[3],d[4],b,d[1],d[2]]
    z=ne(op1[d[7]],mey11x[c,d[5],a,d[3],d[4],b,d[1],d[2],d[6]])
    if a==5 || c==5
      z=:0
    end
    return (z,eval(z))
end

function doe4(a,b,c,d)
# (a,((b,c),d))
#   x=me(b,ne(d[1],1),ne(d[2],1))
#   y=nex1[d[5],c,d[3],b,d[1],d[2],d[4]]
    z=ne(op1[d[7]],me1yx1[a,d[6],d[5],c,d[3],b,d[1],d[2],d[4]])
    if a==5 || (b!=5 && c==5)
      z=:0
    end
    return (z,eval(z))
end

function doe5(a,b,c,d)
# (a,(b,(c,d)))
#   x=me(c,ne(d[1],1),ne(d[2],1))
#   y=ne1x[d[5],b,d[3],d[4],c,d[1],d[2]]
    z=ne(op1[d[7]],me1y1x[a,d[6],d[5],b,d[3],d[4],c,d[1],d[2]])
    if a==5 || b==5
      z=:0
    end
    return (z,eval(z))
end

function main(pz)
  d=[1,1,1,1,1,1,1]

for a in 1:5, b in 1:5, c in 1:5
  for d[1] in 1:5, d[2] in 1:5, d[3] in 1:5, d[4] in 1:5, d[5] in 1:1, d[6] in 1:1, d[7] in 1:5
    e,f = doe1(a,b,c,d)
    if f>0.5 && f<200
      push!(pz,(e,f))
    end
    e,f = doe2(a,b,c,d)
    if f>0.5 && f<200
      push!(pz,(e,f))
    end
    e,f = doe3(a,b,c,d)
    if f>0.5 && f<200
      push!(pz,(e,f))
    end
    e,f = doe4(a,b,c,d)
    if f>0.5 && f<200
      push!(pz,(e,f))
    end
    e,f = doe5(a,b,c,d)
    if f>0.5 && f<200
      push!(pz,(e,f))
    end
  end
  end
end

function fine(pz)
  println(length(pz))
  qz=[]
  sort!(pz, by=x->x[2])
  pl=length(pz)
  for i in 1:pl
    push!(qz,Int(ceil(pz[i][2])))
  end

  for i in 1:121
    if in(i,qz)==true
      j=findfirst(isequal(i),qz)
      rz=repr(pz[j][1])
      rz=replace(rz,":"=>"")
      rz=replace(rz,"cc(+1, +1)"=>"11")
      rz=replace(rz,"cc(+(11), +1)"=>"111")
      rz=replace(rz,"+1"=>"1")
      println(i," -> ",rz," = ",pz[j][2])
    end
  end
end

pz=[]
@time main(pz)
@time fine(pz)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 18.116161 seconds (20.53 M allocations: 1.368 GiB, 85.01% gc time)
420.117500 seconds (115.60 M allocations: 6.676 GiB, 4.81% gc time)
740539
1 -> (exp(1 + +((rog(1) - rog(exp(1) + exp(1)))))) = 0.5000000000000001
2 -> (cbr(1 + +(1 / exp(11)))) = 1.000005567202603
3 -> (cbr(1 - +(exp(1) * sqr(11)))) = 2.0012925727345556
4 -> (+((1 + +((exp(1) - cbr(1 / exp(1))))))) = 3.001750517885256
5 -> (exp(1 + +((exp(1) - sqr(exp(1) + exp(1)))))) = 4.0013741789600505
6 -> (+(1 * +(exp(1) / rog(1 - exp(1))))) = 5.021535230269061
7 -> (exp(1 + +(1 / cbr(1 + 1)))) = 6.011657650957131
8 -> (exp(+((1 + cbr(exp(1) * exp(1)))) - 1)) = 7.012778894114402
9 -> (+((+(exp(1) * sqr(11)) - 1))) = 8.015520899439874
10 -> (exp(1 + +((rog(1) + cbr(1 - exp(1)))))) = 9.004695764350682
11 -> (rog(1 + +(exp(11) / exp(1)))) = 10.000045398899218
12 -> (exp(+((1 + rog(11))) - 1)) = 11.000000000000002
13 -> (rog(1 + +(exp(1) * exp(11)))) = 12.000006144193478
14 -> (exp(+(rog(1 + exp(1)) * exp(1)) - 1)) = 13.063412466658708
15 -> (exp(+((sqr(exp(1) + rog(1)) + rog(1))) + 1)) = 14.135951028495938
16 -> (sqr(1 + +((exp(1) - exp(exp(1) + exp(1)))))) = 15.031080541832813
17 -> (+((1 + +(exp(1) / exp(1 - exp(1)))))) = 16.15426224147926
18 -> (+((+((1 + exp(exp(1) + rog(1)))) + 1))) = 17.154262241479262
19 -> (exp(+((exp(1) + exp(1 - exp(1)))) * 1)) = 18.131593378402822
20 -> (exp(1 + +((rog(1) + cbr(exp(1) * exp(1)))))) = 19.062709434872296
21 -> (+(+((exp(1) + rog(1))) * +(exp(1) * exp(1)))) = 20.085536923187664
22 -> (exp(+(+(11) / exp(1)) - 1)) = 21.045228353448717
23 -> (exp(1 + +((exp(1) - sqr(1 / exp(1)))))) = 22.46034183113636
24 -> (exp(1 + +(exp(1) / cbr(1 + 1)))) = 23.511783406204675
25 -> (exp(1 / +((rog(1 + exp(1)) - 1)))) = 24.34239023161335
26 -> (sqr(1 + +(exp(exp(1) + exp(1)) * exp(1)))) = 25.005158374895853
27 -> (exp(1 / +((1 - rog(1 + 1))))) = 26.020673411168517
28 -> (exp(1 + +((1 + sqr(1 - exp(1)))))) = 27.40793292860369
29 -> (cbr(1 - +(exp(11) / exp(1)))) = 28.031200676839138
30 -> (+(+((exp(1) + exp(1))) * +((exp(1) + exp(1))))) = 29.556224395722598
31 -> (exp(+((exp(1) + rog(exp(1) + exp(1)))) - 1)) = 30.308524482958514
32 -> (exp(+((exp(1) + cbr(1 / exp(1)))) * 1)) = 31.025614547492058
33 -> (exp(+((cbr(exp(1) + exp(1)) + exp(1))) - 1)) = 32.350947201581
35 -> (exp(1 / +((1 - cbr(1 / exp(1)))))) = 34.046473983702036
36 -> (exp(+(sqr(1 - exp(1)) * exp(1)) * 1)) = 35.27632820034533
37 -> (exp(+(rog(exp(1) + exp(1)) * exp(1)) - 1)) = 36.68805458104296
38 -> (exp(+(11) - +(exp(1) * exp(1)))) = 37.00096158422464
39 -> (exp(+((sqr(1 + exp(1)) + exp(1))) - 1)) = 38.34279034771416
40 -> (exp(+((exp(1) + cbr(exp(1) * exp(1)))) - 1)) = 39.09583226076508
41 -> (+(+((exp(1) + exp(1))) * +(exp(1) * exp(1)))) = 40.17107384637533
42 -> (exp(+((+((1 + 1)) + exp(1))) - 1)) = 41.1935556747161
43 -> (+((+((exp(1 + exp(1)) + rog(1))) + 1))) = 42.193555674716116
44 -> (+((+((1 + exp(1 + exp(1)))) + 1))) = 43.193555674716116
45 -> (exp(1 + +((exp(1) + exp(rog(1) - exp(1)))))) = 44.00353027860415
47 -> (exp(+(sqr(1 + 1) * exp(1)) * 1)) = 46.72274206040535
48 -> (exp(1 + +((exp(exp(1) - 1) - exp(1))))) = 47.30706718593521
50 -> (exp(1 + +((exp(1) + exp(1 - exp(1)))))) = 49.286780801520734
51 -> (exp(+((exp(1) + cbr(1 - exp(1)))) * 1)) = 50.20065233451701
52 -> (exp(+((exp(1) + cbr(11))) - 1)) = 51.5350376483723
54 -> (exp(+((cbr(1 + 1) + exp(1))) * 1)) = 53.42094397564407
55 -> (cbr(1 - +(exp(1) * exp(11)))) = 54.598038212039256
56 -> (exp(+(exp(1) / rog(1 - exp(1))) - 1)) = 55.78668552594455
57 -> (exp(+((exp(1) + sqr(1 - exp(1)))) * 1)) = 56.21110430560197
58 -> (exp(+(1 / exp(1)) * +(11))) = 57.20686180895067
60 -> (exp(1 + +((exp(1) + exp(rog(1) - 1))))) = 59.511005963978846
62 -> (exp(+((cbr(exp(1) + rog(1)) + exp(1))) * 1)) = 61.18452227596258
63 -> (exp(+((sqr(1 + 1) + exp(1))) * 1)) = 62.33327490494041
65 -> (exp(+((exp(1) + exp(1 / exp(1)))) * 1)) = 64.2607927021827
67 -> (sqr(1 - +(exp(exp(1) * exp(1)) * exp(1)))) = 66.31488392984271
68 -> (exp(+(cbr(1 + exp(1)) * exp(1)) * 1)) = 67.43919204479707
69 -> (exp(1 + +((1 + cbr(11))))) = 68.30480329758417
70 -> (exp(+(sqr(1 + exp(1)) * exp(1)) - 1)) = 69.52046855397886
71 -> (exp(1 + +(cbr(1 - exp(1)) * exp(1)))) = 70.51403099154743
72 -> (exp(+((cbr(1 + exp(1)) + exp(1))) * 1)) = 71.34344142876819
74 -> (exp(+(cbr(exp(1) * exp(1)) * exp(1)) - 1)) = 73.2948282500825
75 -> (exp(1 + +((rog(1) + sqr(11))))) = 74.93527870314952
76 -> (exp(1 + +((exp(1) + sqr(1 / exp(1)))))) = 75.55134476063684
77 -> (exp(+((1 + sqr(exp(1) + exp(1)))) + 1)) = 76.06924026135655
79 -> (exp(+((sqr(exp(1) + rog(1)) + exp(1))) * 1)) = 78.80710038074831
82 -> (+(+(exp(1) * exp(1)) * +(11))) = 81.27961708823715
83 -> (exp(+((exp(1) + rog(1 + 1))) + 1)) = 82.3871113494322
84 -> (+((+(exp(exp(1) + exp(1)) / exp(1)) - 1))) = 83.48412584713863
85 -> (exp(1 + +((exp(1) + cbr(1 / exp(1)))))) = 84.33636424122227
86 -> (+((1 + +(exp(exp(1) + exp(1)) / exp(1))))) = 85.48412584713863
88 -> (exp(+((cbr(exp(1) + exp(1)) + exp(1))) * 1)) = 87.93899191149563
89 -> (exp(+(sqr(exp(1) + rog(1)) * exp(1)) * 1)) = 88.38383317988601
96 -> (exp(1 + +(sqr(1 - exp(1)) * exp(1)))) = 95.89100192175613
97 -> (exp(1 + +(rog(1 + exp(1)) * exp(1)))) = 96.52628755961122
98 -> (exp(+(1 / exp(1 - exp(1))) - 1)) = 97.02236556502673
100 -> (exp(+(rog(exp(1) + exp(1)) * exp(1)) * 1)) = 99.72847208916271
105 -> (exp(+((sqr(1 + exp(1)) + exp(1))) * 1)) = 104.22651025460627
107 -> (exp(+((exp(1) + cbr(exp(1) * exp(1)))) * 1)) = 106.27349040292063
110 -> (+((+(111) - 1))) = 110
111 -> (+((+(exp(1 + exp(1)) * exp(1)) - 1))) = 110.9756938401968
112 -> (exp(+((1 + 1)) + +((exp(1) + rog(1))))) = 111.97569384019675
113 -> (+((1 + +(exp(1 + exp(1)) * exp(1))))) = 112.9756938401968
120 -> (exp(+(cbr(exp(1) + exp(1)) * exp(1)) * 1)) = 119.07124802191112
121 -> (exp(1 + +(cbr(exp(1) + rog(1)) * exp(1)))) = 120.74343163832631
 10.825410 seconds (63.96 M allocations: 993.083 MiB, 36.17% gc time)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

地獄からの復活劇

新聞によりますと，岡山の小学校５年生が，地獄−極楽の移動時間を求めたとのこと。移動時間は12.7年とのことで，移動に用いられる蜘蛛の糸の強度も割り出したそうだ。

これは，是非とも調べてみなければいけません。一般財団法人理数教育研究所（代表理事岡本和夫）が主催している「算数・数学の自由研究」作品コンクールの審査員特別賞の受賞作の１つであった。「地獄からの復活劇 〜御釈迦様からの試練〜岡山県 岡山大学教育学部附属小学校５年末光由季」

地獄と極楽の距離は，何万里という芥川の本文の表現の下限から１万里≒４万ｋｍとしている。静止軌道半径に近いので，極楽は静止軌道上にあるのかもしれない。ということは蜘蛛の糸は宇宙エレベータの暗示だったのか・・・

すばらしいのは，蜘蛛の糸を登る時間を自分で実験しているところだ。近所の公園の登り棒を３ｍ昇る時間の平均から，５秒で１ｍ（0.2ｍ／秒）を得ている。最近は，危険だということで公園から次々と遊具が撤去されている。うちの近くの公園でも子どもたちが遊んでいたブランコや回転遊具は何年か前に消えてしまい，もうすぐ鉄棒までなくなるそうだ。組み体操をやめるのは賛成だが，公園の遊具は残してほしかった。

話がそれた。結局，４万ｋｍ÷ 0.2 ｍ／秒＝２億秒＝ 6.4 年・・・あれ？ファクター２違うな。２億秒まではあっている。で，１年≒ π × 10^７秒ですよね。２億秒から時間に換算する計算が２倍になっているようだ。１日の半分は休憩とすればつじつまがあうということにしておこう。

さらに蜘蛛の糸の強度まで計算してピアノ線と比較しているところがすごい。最後のコメントも含蓄がある。宇宙エレベータ協会の皆様は是非，彼女を顧問に迎えられたい。

軌道エレベータの基部の想像図（Wikipediaより）

Juliaでパズル（４）

（Juliaでパズル（３）からの続き）

目的はJuliaプログラミングの習得ということなのでこの話題がだらだたと続いている。高速が話題のJuliaで初めて実行時間が１時間以上のプログラムを作ったが，アルゴリズムがポンコツだといくらでも遅いプログラムができるのであった・・・orz

昨日の２項演算だけの組み合わせで整数を作るものから，単項演算を加えた組み合わせに拡張する。$5 \times 5^3 = 625$回から$5 \times 5^3 \times 5^7 ≒ 5000$万回に増えて７時間で終わらなかったので，次は時間短縮の方法を考えなければならない。

単項演算は１つのエントリーポイントで１回だけ作用するものとしている。
以下は，７重ループを５重ループに落として７分程度で実行したもの。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
function sqr(x)
  return sqrt(abs(x))
end

function cbr(x)
  return cbrt(abs(x))
end

function rog(x)
  return log(abs(x))
end

function cc(a,b)
  if (a==1 && b==1) || (a==11 && b==1)
    return 10*a+b
  else
    return 0
  end
end

function ne(op, op1)
  expr = Expr(:call,Symbol(op), op1)
  return expr
end

function me(op, op1, op2)
  expr = Expr(:call,Symbol(op), op1, op2)
  return expr
end

function doe1(a,b,c,d)
# ((1,1)_a,(1,1)_c)_b
    x=me(a,ne(d[1],1),ne(d[2],1))
    y=me(c,ne(d[3],1),ne(d[4],1))
    z=ne(d[7],me(b,ne(d[5],x),ne(d[6],y)))
    if b=="cc"
      z=:0
    end
    return (z,eval(z))
end

function doe2(a,b,c,d)
# (((1,1)_a,1)_b,1)_c
    x=me(a,ne(d[1],1),ne(d[2],1))
    y=me(b,ne(d[3],x),ne(d[4],1))
    z=ne(d[7],me(c,ne(d[5],y),ne(d[6],1)))
    if  (a!="cc" && b=="cc") || ((a!="cc" || b!="cc") && c=="cc")
      z=:0
    end
    return (z,eval(z))
end

function doe3(a,b,c,d)
# ((1,(1,1)_b)_a,1)_c
    x=me(b,ne(d[1],1),ne(d[2],1))
    y=me(a,ne(d[3],1),ne(d[4],x))
    z=ne(d[7],me(c,ne(d[5],y),ne(d[6],1)))
    if a=="cc" || c=="cc"
      z=:0
    end
    return (z,eval(z))
end

function doe4(a,b,c,d)
# (1,((1,1)_b,1)_c)_a
    x=me(b,ne(d[1],1),ne(d[2],1))
    y=me(c,ne(d[3],x),ne(d[4],1))
    z=ne(d[7],me(a,ne(d[6],1),ne(d[5],y)))
    if a=="cc" || (b!="cc" && c=="cc")
      z=:0
    end
    return (z,eval(z))
end

function doe5(a,b,c,d)
# (1,(1,(1,1)_c)_b)_a
    x=me(c,ne(d[1],1),ne(d[2],1))
    y=me(b,ne(d[3],1),ne(d[4],x))
    z=ne(d[7],me(a,ne(d[6],1),ne(d[5],y)))
    if a=="cc" || b=="cc"
      z=:0
    end
    return (z,eval(z))
end

function main(pz)
  uno=["+","-","*","/","cc"]
  dos=["+","sqr","cbr","exp","rog"]
  d=["+","+","+","+","+","+","+"]
for a in uno, b in uno, c in uno
#a="*"; b="-"; c="/";
  for d[1] in dos, d[2] in dos, d[3] in dos, d[4] in dos, d[7] in dos
#, d[5] in dos, d[6] in dos, d[7] in dos
    e,f = doe1(a,b,c,d)
    if f>0.5 && f<200
      push!(pz,(e,f))
    end
    e,f = doe2(a,b,c,d)
    if f>0.5 && f<200
      push!(pz,(e,f))
    end
    e,f = doe3(a,b,c,d)
    if f>0.5 && f<200
      push!(pz,(e,f))
    end
    e,f = doe4(a,b,c,d)
    if f>0.5 && f<200
      push!(pz,(e,f))
    end
    e,f = doe5(a,b,c,d)
    if f>0.5 && f<200
      push!(pz,(e,f))
    end
  end
  end
end

function fine(pz)
  println(length(pz))
  qz=[]
  sort!(pz, by=x->x[2])
  pl=length(pz)
  for i in 1:pl
    push!(qz,Int(ceil(pz[i][2])))
  end

  for i in 1:121
    if in(i,qz)==true
      j=findfirst(isequal(i),qz)
      rz=repr(pz[j][1])
      rz=replace(rz,":"=>"")
      rz=replace(rz,"cc(+1, +1)"=>"11")
      rz=replace(rz,"cc(+(11), +1)"=>"111")
      rz=replace(rz,"+1"=>"1")
      println(i," -> ",rz," = ",pz[j][2])
    end
  end
end

pz=[]
@time main(pz)
@time fine(pz)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

408.936226 seconds (144.22 M allocations: 8.591 GiB, 5.96% gc time)

751191

1 -> (exp(+((rog(1) - rog(exp(1) + exp(1)))) + 1)) = 0.5000000000000001

2 -> (cbr(+(1 / exp(11)) + 1)) = 1.000005567202603

3 -> (cbr(+(exp(1) * sqr(11)) - 1)) = 2.0012925727345556

4 -> (+((+((exp(1) - cbr(1 / exp(1)))) + 1))) = 3.001750517885256

5 -> (exp(+((exp(1) - sqr(exp(1) + exp(1)))) + 1)) = 4.0013741789600505

6 -> (+(+(exp(1) / rog(1 - exp(1))) * 1)) = 5.021535230269061

7 -> (exp(+(1 / cbr(1 + 1)) + 1)) = 6.011657650957131

8 -> (exp(+((1 + cbr(exp(1) * exp(1)))) - 1)) = 7.012778894114402

9 -> (+((+(exp(1) * sqr(11)) - 1))) = 8.015520899439874

10 -> (exp(+((rog(1) + cbr(1 - exp(1)))) + 1)) = 9.004695764350682

11 -> (rog(+(exp(11) / exp(1)) + 1)) = 10.000045398899218

12 -> (exp(+((1 + rog(11))) - 1)) = 11.000000000000002

13 -> (rog(+(exp(1) * exp(11)) + 1)) = 12.000006144193478

14 -> (exp(+(rog(1 + exp(1)) * exp(1)) - 1)) = 13.063412466658708

15 -> (exp(+((sqr(exp(1) + rog(1)) + rog(1))) + 1)) = 14.135951028495938

16 -> (sqr(+((exp(1) - exp(exp(1) + exp(1)))) + 1)) = 15.031080541832813

17 -> (+((+(exp(1) / exp(1 - exp(1))) + 1))) = 16.15426224147926

18 -> (+((+((1 + exp(exp(1) + rog(1)))) + 1))) = 17.154262241479262

19 -> (exp(+((exp(1) + exp(1 - exp(1)))) * 1)) = 18.131593378402822

20 -> (exp(+((rog(1) + cbr(exp(1) * exp(1)))) + 1)) = 19.062709434872296

21 -> (+(+((exp(1) + rog(1))) * +(exp(1) * exp(1)))) = 20.085536923187664

22 -> (exp(+(+(11) / exp(1)) - 1)) = 21.045228353448717

23 -> (exp(+((exp(1) - sqr(1 / exp(1)))) + 1)) = 22.46034183113636

24 -> (exp(+(exp(1) / cbr(1 + 1)) + 1)) = 23.511783406204675

25 -> (sqr(+(exp(exp(1) * exp(1)) / exp(1)) - 1)) = 24.37815447036041

26 -> (sqr(+(exp(exp(1) + exp(1)) * exp(1)) + 1)) = 25.005158374895853

27 -> (exp(+((exp(1) + rog(1 - exp(1)))) * 1)) = 26.039293433236857

28 -> (exp(+((1 + sqr(1 - exp(1)))) + 1)) = 27.40793292860369

29 -> (cbr(+(exp(11) / exp(1)) - 1)) = 28.031200676839138

30 -> (+(+((exp(1) + exp(1))) * +((exp(1) + exp(1))))) = 29.556224395722598

31 -> (exp(+((exp(1) + rog(exp(1) + exp(1)))) - 1)) = 30.308524482958514

32 -> (exp(+((exp(1) + cbr(1 / exp(1)))) * 1)) = 31.025614547492058

33 -> (exp(+((cbr(exp(1) + exp(1)) + exp(1))) - 1)) = 32.350947201581

35 -> (exp(+((exp(1) - exp(1 - exp(1)))) + 1)) = 34.42929324111139

36 -> (exp(+(sqr(1 - exp(1)) * exp(1)) * 1)) = 35.27632820034533

37 -> (exp(+(rog(exp(1) + exp(1)) * exp(1)) - 1)) = 36.68805458104296

38 -> (exp(+(11) - +(exp(1) * exp(1)))) = 37.00096158422464

39 -> (exp(+((sqr(1 + exp(1)) + exp(1))) - 1)) = 38.34279034771416

40 -> (exp(+((exp(1) + cbr(exp(1) * exp(1)))) - 1)) = 39.09583226076508

41 -> (+(+((exp(1) + exp(1))) * +(exp(1) * exp(1)))) = 40.17107384637533

42 -> (exp(+((+((1 + 1)) + exp(1))) - 1)) = 41.1935556747161

43 -> (+((+((exp(1 + exp(1)) + rog(1))) + 1))) = 42.193555674716116

44 -> (+((+((1 + exp(1 + exp(1)))) + 1))) = 43.193555674716116

45 -> (exp(+((exp(1) + exp(rog(1) - exp(1)))) + 1)) = 44.00353027860415

47 -> (exp(+(sqr(1 + 1) * exp(1)) * 1)) = 46.72274206040535

48 -> (exp(+((exp(exp(1) - 1) - exp(1))) + 1)) = 47.30706718593521

50 -> (exp(+((exp(1) + exp(1 - exp(1)))) + 1)) = 49.286780801520734

51 -> (exp(+((exp(1) + cbr(1 - exp(1)))) * 1)) = 50.20065233451701

52 -> (exp(+((exp(1) + cbr(11))) - 1)) = 51.5350376483723

54 -> (exp(+((cbr(1 + 1) + exp(1))) * 1)) = 53.42094397564407

55 -> (cbr(+(exp(1) * exp(11)) - 1)) = 54.598038212039256

56 -> (exp(+(exp(1) / rog(1 - exp(1))) - 1)) = 55.78668552594455

57 -> (exp(+((exp(1) + sqr(1 - exp(1)))) * 1)) = 56.21110430560197

58 -> (exp(+(1 / exp(1)) * +(11))) = 57.20686180895067

60 -> (exp(+((exp(1) + exp(rog(1) - 1))) + 1)) = 59.511005963978846

62 -> (exp(+((cbr(exp(1) + rog(1)) + exp(1))) * 1)) = 61.18452227596258

63 -> (exp(+((sqr(1 + 1) + exp(1))) * 1)) = 62.33327490494041

65 -> (exp(+((exp(1) + exp(1 / exp(1)))) * 1)) = 64.2607927021827

67 -> (sqr(+(exp(exp(1) * exp(1)) * exp(1)) - 1)) = 66.31488392984271

68 -> (exp(+(cbr(1 + exp(1)) * exp(1)) * 1)) = 67.43919204479707

69 -> (exp(+((1 + cbr(11))) + 1)) = 68.30480329758417

70 -> (exp(+(sqr(1 + exp(1)) * exp(1)) - 1)) = 69.52046855397886

71 -> (exp(+(cbr(1 - exp(1)) * exp(1)) + 1)) = 70.51403099154743

72 -> (exp(+((cbr(1 + exp(1)) + exp(1))) * 1)) = 71.34344142876819

74 -> (exp(+(cbr(exp(1) * exp(1)) * exp(1)) - 1)) = 73.2948282500825

75 -> (exp(+((rog(1) + sqr(11))) + 1)) = 74.93527870314952

76 -> (exp(+((exp(1) + sqr(1 / exp(1)))) + 1)) = 75.55134476063684

77 -> (exp(+((1 + sqr(exp(1) + exp(1)))) + 1)) = 76.06924026135655

79 -> (exp(+((sqr(exp(1) + rog(1)) + exp(1))) * 1)) = 78.80710038074831

82 -> (+(+(exp(1) * exp(1)) * +(11))) = 81.27961708823715

83 -> (exp(+((exp(1) + rog(1 + 1))) + 1)) = 82.3871113494322

84 -> (+((+(exp(exp(1) + exp(1)) / exp(1)) - 1))) = 83.48412584713863

85 -> (exp(+((exp(1) + cbr(1 / exp(1)))) + 1)) = 84.33636424122227

86 -> (+((+(exp(exp(1) + exp(1)) / exp(1)) + 1))) = 85.48412584713863

88 -> (exp(+((cbr(exp(1) + exp(1)) + exp(1))) * 1)) = 87.93899191149563

89 -> (exp(+(sqr(exp(1) + rog(1)) * exp(1)) * 1)) = 88.38383317988601

96 -> (exp(+(sqr(1 - exp(1)) * exp(1)) + 1)) = 95.89100192175613

97 -> (exp(+(rog(1 + exp(1)) * exp(1)) + 1)) = 96.52628755961122

98 -> (exp(+(1 / exp(1 - exp(1))) - 1)) = 97.02236556502673

100 -> (exp(+(rog(exp(1) + exp(1)) * exp(1)) * 1)) = 99.72847208916271

105 -> (exp(+((sqr(1 + exp(1)) + exp(1))) * 1)) = 104.22651025460627

107 -> (exp(+((exp(1) + cbr(exp(1) * exp(1)))) * 1)) = 106.27349040292063

110 -> (+((+(cc(+(11), 1)) - 1))) = 110

111 -> (+((+(exp(1 + exp(1)) * exp(1)) - 1))) = 110.9756938401968

112 -> (exp(+((1 + 1)) + +((exp(1) + rog(1))))) = 111.97569384019675

113 -> (+((+(exp(1 + exp(1)) * exp(1)) + 1))) = 112.9756938401968

120 -> (exp(+(cbr(exp(1) + exp(1)) * exp(1)) * 1)) = 119.07124802191112

121 -> (exp(+(cbr(exp(1) + rog(1)) * exp(1)) + 1)) = 120.74343163832631

  6.644247 seconds (64.72 M allocations: 1004.683 MiB, 1.93% gc time)

Juliaでパズル（３）

（Juliaでパズル（２）からの続き）

世の中で流布している数学パズルだと，４つの数字と加減乗除を組み合わせて10をつくる Make 10 が有名らしい。もちろん，この手の問題では計算結果はすべて整数となるのが普通であるが，今回は，実数の範囲で計算して最後に整数化できるものを探している。

とりあえず，２項演算だけを組み合わせた結果をすべて探すプログラムを作ってみた。

(1) findfirst()のように，ある目的の数に一致する数をコレクション（配列またはタプル）の中から探すとともに，その位置の添字を返すことができれば有り難い。いまのところ「タプルの配列に対して１コマンドである数字を見つける方法」がわからない。

(2) ２次元配列ならば融通がききそうなのであるが，「タプルの配列を１コマンドで多次元配列に展開する方法」がわからない。もちろん地道にやればできることなのだが。

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
function cc(a,b)
  return 10*a+b
end

function me(op, op1, op2)
  expr = Expr(:call,Symbol(op), op1, op2)
  return expr
end

function dop1(a,b,c)
# ((a,b),(cd))
    x=me(a,1,1)
    y=me(c,1,1)
    z=me(b,x,y)
    if b=="cc"
      z=:0
    end
    return (z,eval(z))
end

function dop2(a,b,c)
# (((a,b),c),d)
    x=me(a,1,1)
    y=me(b,x,1)
    z=me(c,y,1)
    if  (a!="cc" && b=="cc") || ((a!="cc" || b!="cc") && c=="cc")
      z=:0
    end
    return (z,eval(z))
end

function dop3(a,b,c)
# ((a,(b,c)),d)
    x=me(b,1,1)
    y=me(a,1,x)
    z=me(c,y,1)
    if a=="cc" || c=="cc"
      z=:0
    end
    return (z,eval(z))
end

function dop4(a,b,c)
# (a,((b,c),d))
    x=me(b,1,1)
    y=me(c,x,1)
    z=me(a,1,y)
    if a=="cc" || (b!="cc" && c=="cc")
      z=:0
    end
    return (z,eval(z))
end

function dop5(a,b,c)
# (a,(b,(c,d)))
    x=me(c,1,1)
    y=me(b,1,x)
    z=me(a,1,y)
    if a=="cc" || b=="cc"
      z=:0
    end
    return (z,eval(z))
end

pz=[]
uno=["+","-","*","/","cc"]
for a in uno, b in uno, c in uno
    d,e = dop1(a,b,c)
    if e>0.5 && e<Inf
      push!(pz,(d,e))
    end
    d,e = dop2(a,b,c)
    if e>0.5 && e<Inf
      push!(pz,(d,e))
    end
    d,e = dop3(a,b,c)
    if e>0.5 && e<Inf
      push!(pz,(d,e))
    end
    d,e = dop4(a,b,c)
    if e>0.5 && e<Inf
      push!(pz,(d,e))
    end
    d,e = dop5(a,b,c)
    if e>0.5 && e<Inf
      push!(pz,(d,e))
    end
end

qz=[]
sort!(pz, by=x->x[2])
pl=length(pz)
for i in 1:pl
    push!(qz,Int(ceil(pz[i][2])))
end

for i in 1:121
  if in(i,qz)==true
    j=findfirst(isequal(i),qz)
    rz=repr(pz[j][1])
    rz=replace(rz,":"=>"")
    rz=replace(rz,"cc(1, 1)"=>"11")
    rz=replace(rz,"cc(11, 1)"=>"111")
    println(i," -> ",rz," = ",pz[j][2])
  end
end
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

1 -> (1 - 1 / 11) = 0.9090909090909091

2 -> (1 + 1 / 11) = 1.0909090909090908

3 -> ((1 + 1) + 1 * 1) = 3

4 -> ((1 + 1) + (1 + 1)) = 4

6 -> (11 / (1 + 1)) = 5.5

9 -> (11 - (1 + 1)) = 9

10 -> (1 * 11 - 1) = 10

11 -> ((1 + 11) - 1) = 11

12 -> (1 + 1 * 11) = 12

13 -> ((1 + 1) + 11) = 13

22 -> ((1 + 1) * 11) = 22

110 -> (111 - 1) = 110

111 -> (1 * 111) = 111

112 -> (1 + 111) = 112

121 -> (11 * 11) = 121

「冬麗の不思議をにぎる赤ン坊」（野澤節子　1920-1995）

（Juliaでパズル（４）に続く）

Juliaでパズル（２）

（Juliaでパズル（１）からの続き）

４つの数字から整数をつくるパズルをコンピュータで処理させるためには，式に対応する文字列を，プログラム内で式として評価する必要がある。これは，JuliaドキュメントのMetaprogrammingのところに詳しいが，ここでのポイントは，Juliaでは，Juliaプログラムの中で，Juliaコードを生成したり操作したりすることができるということである。具体例としては次のような例文が示されている。

julia> function math_expr(op, op1, op2)
           expr = Expr(:call, op, op1, op2)
           return expr
       end
math_expr (generic function with 1 method)

julia>  ex = math_expr(:+, 1, Expr(:call, :*, 4, 5))
:(1 + 4 * 5)

julia> eval(ex)
21

そこで（1, 1, 1, 1）パズルで上記を用い，２項演算子を３つ与えた場合の５つのパターンの構文とその値を出力する関数を書いてみる。cc(a, b) = 10*a + bは２つの数字を並べて２桁の数（例えば a=1, b=1 なら cc(1,1) = 11である）を作る関数である。funcion me() においてSymbol(op)は文字列として与えた２項演算子を Expr型に変換する関数である（注：Ver 0.6以前の資料では symbol と小文字で始まっているので，Ver. 1.0 以降では大文字にすること）。

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
function cc(a,b)
  return 10*a+b
end

function me(op, op1, op2)
  expr = Expr(:call,Symbol(op), op1, op2)
  return expr
end

function dop1(a,b,c)
# ((a,b),(cd))
    x=me(a,1,1)
    y=me(c,1,1)
    z=me(b,x,y)
    println((z,eval(z)))
end

function dop2(a,b,c)
# (((a,b),c),d)
    x=me(a,1,1)
    y=me(b,x,1)
    z=me(c,y,1)
    println((z,eval(z)))
end

function dop3(a,b,c)
# ((a,(b,c)),d)
    x=me(b,1,1)
    y=me(a,1,x)
    z=me(c,y,1)
    println((z,eval(z)))
end

function dop4(a,b,c)
# (a,((b,c),d))
    x=me(b,1,1)
    y=me(c,x,1)
    z=me(a,1,y)
    println((z,eval(z)))
end

function dop5(a,b,c)
# (a,(b,(c,d)))
    x=me(c,1,1)
    y=me(b,1,x)
    z=me(a,1,y)
    println((z,eval(z)))
end
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

dop1("cc","/","+")
dop2("cc","+","*")
dop3("-","*","/")
dop4("+","\\","^")
dop5("+","-","*")

(:(cc(1, 1) / (1 + 1)), 5.5)
(:((cc(1, 1) + 1) * 1), 12)
(:((1 - 1 * 1) / 1), 0.0)
(:(1 + (1 \ 1) ^ 1), 2.0)
(:(1 + (1 - 1 * 1)), 1)

（Juliaでパズル（３）に続く）

Juliaでパズル（１）

４つの同じ数字に演算をほどこして，１から順に整数を作るというパズルがある。
例えば，6÷6+6-6 = 1,　 6÷6+6÷6=2,　 (6+6+6)/6 = 3,　 などなど。
これを，Juliaの基本演算子や基本関数だけを使って実行してみると，

1　1+1-1*1
2　1+1-1+1
3　1+1*1+1
4　1+1+1+1
5　Int(floor(11/(1+1)))
6　Int(ceil(11/(1+1)))
7　Int(floor(exp(1*1+1*1)))
8　Int(ceil(exp(1*1+1*1)))
9　11-1-1
10　11-1*1
11　11+1-1
12　11+1*1
13　11+1+1
14　Int(ceil(sqrt(sqrt(exp(11)))-1-1))
15　Int(floor(exp(1+sqrt(1+1+1))))
16　Int(ceil(exp(1+sqrt(1+1+1))))
17　Int(floor(exp((1+1)*sqrt(1+1))))
18　Int(ceil(exp((1+1)*sqrt(1+1))))
19　Int(floor(exp(1+1+1)-1))
20　Int(floor(exp(1+1*1+1)))
21　Int(floor(exp(1+1+1)+1))
22　Int(ceil(exp(1+1+1)+1))
23　Int(floor(exp(sqrt(1+1)*cbrt(11))))
24　Int(ceil(exp(sqrt(1+1)*cbrt(11))))
25　Int(floor(exp(sqrt(11))-1-1))
26　Int(floor(exp(sqrt(11))-1*1))
27　Int(floor(exp(sqrt(11))+1-1))
28　Int(ceil(exp(sqrt(11))+1-1))
29　Int(ceil(exp(sqrt(11))+1*1))
30　Int(ceil(exp(sqrt(11))+1+1))
31　Int(floor(sqrt(sqrt(exp(11)))*(1+1)))
32　Int(floor(cbrt(exp(sqrt(111)))-1))
33　Int(floor(cbrt(exp(sqrt(111)))*1))
34　Int(floor(cbrt(exp(sqrt(111)))+1))
35　Int(ceil(cbrt(exp(sqrt(111)))+1))
36　Int(floor(exp(sqrt(11+1+1))))
37　Int(floor(cbrt(exp(11))-1-1))
38　Int(ceil(cbrt(exp(11))-1-1))
39　Int(ceil(cbrt(exp(11))-1*1))
40　Int(ceil(cbrt(exp(11))+1-1))
41　Int(ceil(cbrt(exp(11))+1*1))
42　Int(ceil(cbrt(exp(11))+1+1))
43　Int(floor(cbrt(sqrt(1+1)*exp(11))))
44　Int(ceil(cbrt(sqrt(1+1)*exp(11))))
45　Int(ceil(exp(cbrt(111)-1)))
46　Int(ceil(exp(sqrt(sqrt(11))+1+1)))
47　Int(ceil(sqrt(exp(11-sqrt(11)))))
などなど

（１）もし，ユーザ定義関数が使えるならば，任意の n に対して，すべての整数を作り出すことが簡単にできるのであった・・・

f(n1,n2,n3,n4)=Int(n1/n2+n3-n4)
f(n)=Int(ceil(n+1/n))

f(1,1,1,1)
1
f(f(1,1,1,1))
2
f(f(f(1,1,1,1)))
3

（２）juliaドキュメントの算術的演算と基本関数から主な演算と関数を選んで使うことにする。例えば，２項演算ならば，+, -, *, / の４種とする（ここでは，整数ゼロ除算をさけるために÷と%は含めない）。数字を２〜４桁並べる操作が必要なのでこれに対応する２項演算を１つ付け加えると５種になる。この２項演算は cc(a,b)=10*a+b で実現できる（例えば３桁並べるには，cc(cc(a, b), c)とすればよい），ただし適応できる位置に制限があることに注意する。次に単項演算を選ぶ。ここでは，+, sqrt(abs()), qbrt(abs()), exp(), log(abs()) の５種とする（+は恒等演算子）。

（３）４つの数字（a, b, c, d,  a = b = c = d）に対応する２項演算のパターンは５つある。
２項演算を（a, b）などで表し，ふつうのカッコの優先順位で計算するものとして，
（（a, b）,（c, d）），（（（a, b）, c）, d），（a,（b,（c,  d））），
（（a,（b, c））, d），（a,（（b, c）, d））の５パターンである。１パターンあたりの２項演算の数は３であるから，パターンあたりの総数は$5^3$通りである。

（４）上記の各数字や２項演算の結果に対して５種の単項演算を任意の階数だけ作用することができる。単項演算のエントリーポイントは１パターンにつき$n$個あるとする。単項演算子の種類を$m$個とし，エントリーポイントに$k$個まで入るならば，パターンあたりの配位の数は，$\sum_{j=0}^k  {}_{n+j-1}C_j m^k$となる。

したがって，この他に条件をつけなければ，総数で $5 \times 5^3 \times (5^k \dfrac{n^k}{k!}+\cdots) $個の式が得られる。これをしらみつぶしに調べるのだろうか・・・

（Juliaでパズル（２）に続く）

円周率が22/7より小さいことの証明

というタイトルのWikipedia記事がある（円周率が22/7より小さいことの証明）。

証明は次の通りである。
\begin{equation}
\begin{aligned}
0 &< \int_0^1 \dfrac{x^4(1-x)^4}{1+x^2} dx \\
& = \int_0^1\dfrac{x^4-4x^5+6x^6-4x^7+x^8}{1+x^2} dx\\
&= \int_0^1x^6 -4x^5 + 5x^4 -4x^2 +4 -\dfrac{4}{1+x^2} dx\\
&=\frac{x^7}{7}-\frac{2 x^6}{3}+x^5-\frac{4 x^3}{3}+4 x-4 \tan^{-1} x\  \Big|_0^1\\
&= \frac{1}{7}-\frac{2}{3}+1-\frac{4}{3}+4-\pi\\
&=\frac{22}{7}-\pi
\end{aligned}
\end{equation}

同様にして円周率が355/113より小さいことも証明される。
\begin{equation}
\begin{aligned}
0 &< \int_0^1 \dfrac{x^8(1-x)^8(25+816x^2)}{3164(1+x^2)} dx \\
&= \cdots\\
&= \frac{12 x^{17}}{791}
-\frac{102 x^{16}}{791}
+\frac{3151 x^{15}}{6780}
-\frac{703 x^{14}}{791}
+\frac{393 x^{13}}{452}
-\frac{23 x^{12}}{113}\\
&-\frac{145 x^{11}}{452}
-\frac{5 x^{10}}{791}
+\frac{1409 x^9}{3164}
-\frac{4 x^7}{7}
+\frac{4 x^5}{5}
-\frac{4 x^3}{3}
+4 x−4 \tan^{-1}x \Bigr|_0^1 \\
&= \cdots\\
&=\frac{355}{113}-\pi
\end{aligned}
\end{equation}

 えーっ，3164と816と25はどこから連れてきたのだ・・・という方は上記記事から参考文献をたどってみよう。この世界にはまだまだ自分が知らないことが沢山埋もれている。

角運動量の合成への道（４）

（角運動量の合成への道（３）からの続き）

次に，スピン１の粒子の軌道角運動量$\boldsymbol{L}$とスピン角運動量$\boldsymbol{T}$の合成を考える。$\boldsymbol{J}=\boldsymbol{L}+\boldsymbol{T}$であり，前回までの手順と同様に，$(\boldsymbol{J}^2, J_z) | J M \rangle = (J(J+1)\hbar^2, M\hbar ) | J M \rangle$ となる状態は，$| \ell m \rangle | t \rangle$ の線形結合で表されるとする。

（１）$J_z$を$| \ell m \rangle | t \rangle$に作用させる。
$J_z | \ell m \rangle | t \rangle = (m + t) \hbar  | \ell m \rangle | t \rangle = M  | \ell m \rangle | t \rangle$。したがって，$m=M-t$となる。

（２）$\boldsymbol{J}^2$を$| \ell M-t \rangle | t \rangle$に作用させる。
\begin{equation}
\begin{aligned}
(\boldsymbol{L}^2+\boldsymbol{T}^2+2L_zT_z &+L_{+}T_{-}+L_{-}T_{+})| \ell M-t \rangle | t \rangle \\
= & \sqrt{\ell(\ell+1)-m(m-1)} \sqrt{2-t(t+1)}\ \hbar^2 | \ell m-1 \rangle | t+1 \rangle\\
+ & (\ell(\ell+1)+2+2 m t )\ \hbar^2 | \ell m \rangle | t \rangle \\
+ & \sqrt{\ell(\ell+1)-m(m+1)}\sqrt{2-t(t-1)}\ \hbar^2   | \ell m+1 \rangle | t-1 \rangle
\end{aligned}
\end{equation}
ここで，$\boldsymbol{J}^2$の固有値を決める固有値方程式は次のようになる。
\begin{equation}
\boldsymbol{J}^2 \sum_{t=-1}^1 \alpha_t  | \ell m-t \rangle | t \rangle
= \lambda \sum_{t=-1}^1 \alpha_t  | \ell m-t \rangle | t \rangle
\end{equation}
ここで$\ell(\ell+1)=j$, $m(m\pm 1)=\mu_{\pm}$と略記する。
\begin{equation}
\begin{pmatrix}j+2m & \sqrt{2(j-\mu_{-})} & 0 \\
\sqrt{2(2(j-\mu_{-})} & j+2 & \sqrt{2(j-\mu_{+})}\\
0& \sqrt{2(j-\mu_{+})}&j-2m \end{pmatrix}
\begin{pmatrix} \alpha_{+1} \\ \alpha_0 \\ \alpha_{-1}\end{pmatrix}
= \lambda \begin{pmatrix} \alpha_{+1} \\ \alpha_0 \\ \alpha_{-1} \end{pmatrix}
\end{equation}
この固有値方程式を解くと，$\lambda=(\ell+1)(\ell+2), \ell(\ell+1), \ell(\ell-1)$となる。
また，固有ベクトルの係数は，
\begin{array}{c|ccc}
  \lambda  & \alpha_{+1} & \alpha_{0} & \alpha_{-1} \\
  \hline
 \ell(\ell-1)
& \sqrt{\frac{(\ell-m)(\ell-m+1)}{2\ell(2\ell+1)}}
& -\sqrt{\frac{(\ell+m)(\ell-m)}{\ell(2\ell+1)}}
& \sqrt{\frac{(\ell+m)(\ell+m+1)}{2\ell(2\ell+1)}}  \\
  \hline
 \ell(\ell+1)
& -\sqrt{\frac{(\ell+m)(\ell-m+1)}{2\ell(\ell+1)}}
& \sqrt{\frac{m}{\ell(\ell+1)}}
&\sqrt{\frac{(\ell-m)(\ell+m+1)}{2\ell(\ell+1)}} \\
 \hline
 (\ell+1)(\ell+2)
& \sqrt{\frac{(\ell+m)(\ell+m+1)}{2(\ell+1)(2\ell+1)}}
& \sqrt{\frac{(\ell+m+1)(\ell-m+1)}{(\ell+1)(2\ell+1)}}
&\sqrt{\frac{(\ell-m)(\ell-m+1)}{2(\ell+1)(2\ell+1)}} \\
  \hline
\end{array}

角運動量の合成への道（３）

（角運動量の合成への道（２）からの続き）

２つの独立な自由度からなる系の別の例として，１つの粒子が軌道角運動量$\boldsymbol{L}$とスピン角運動量$\boldsymbol{S}$を持つ場合を考える。これらを合成した全角運動量を $\boldsymbol{J}=\boldsymbol{L}+\boldsymbol{S}$とする（２粒子系の場合と同様，より正確には通常の３次元空間とスピン空間に作用する演算子のテンソル積 $\boldsymbol{J}=\boldsymbol{L}\otimes \boldsymbol{1}_S+\boldsymbol{1}_L\otimes\boldsymbol{S}$）。

軌道角運動量$\boldsymbol{L}$の固有状態を$|\ell m \rangle$，スピン角運動量$\boldsymbol{S}$の固有状態を$| s \rangle$とし，その積$|\ell m \rangle | s \rangle$を考える。このとき，次の関係が成り立っている。
\begin{equation}
\begin{aligned}
\boldsymbol{L}^2 |\ell m\rangle = \ell(\ell+1)\hbar^2 |\ell m\rangle, \quad &
\boldsymbol{S}^2 | s\rangle = \frac{3}{4} \hbar^2 |s \rangle, \\
L_z |\ell m\rangle = m \hbar |\ell m\rangle, \quad &
S_z | s\rangle = s \hbar |s \rangle, \\
L_\pm |\ell m\rangle = \sqrt{\ell(\ell+1)-m(m\pm 1)} \ \hbar |\ell m \pm 1 \rangle, \quad &
S_\pm | s\rangle = \sqrt{\frac{3}{4}-s(s\pm 1)} \ \hbar |s \pm 1 \rangle
\end{aligned}
\end{equation}
また，合成角運動量$\boldsymbol{J}$も一般の角運動量の交換関係，$[J_i, J_j]=i\hbar \epsilon_{ijk} J_k$を満足するので，$\boldsymbol{J}^2$と$J_z$の同時固有状態を$|j\mu \rangle$とすると，次の関係が成り立つ。
\begin{equation}
\boldsymbol{J}^2 |j \mu \rangle = j(j+1)\hbar^2 |j \mu \rangle, \quad
J_z |j \mu\rangle = \mu \hbar |j \mu\rangle
\end{equation}
$[\boldsymbol{J}^2,\boldsymbol{L}^2]=0, [\boldsymbol{J}^2,\boldsymbol{S}^2]=0, [J_z,L_z]=0,[J_z,S_z]=0$などが成立するので，$|j\mu\rangle$を$|\ell m \rangle | s \rangle$から構成することができる（$[\boldsymbol{J}^2,L_z]\neq 0$や$[\boldsymbol{J}^2,S_z]\neq 0$なので，一般には単独の$|\ell m \rangle | s \rangle$は$\boldsymbol{J}^2$の固有状態ではない）。

（１）$J_z$を$|\ell m \rangle | s \rangle$に作用させる。
\begin{equation}
J_z |\ell m \rangle | s \rangle = (L_z+S_z) |\ell m \rangle | s \rangle = (m+s)\hbar |\ell m \rangle | s \rangle
\end{equation}
したがって，$|\ell m \rangle | s \rangle$は，$J_z$の固有値$\mu\hbar = (m+s)\hbar$の固有状態である。以下，$m=\mu-s$と表記することもある。

（２）$\boldsymbol{J}^2$を$|\ell \mu-s \rangle | s \rangle$に作用させる。
\begin{equation}
\begin{aligned}
\boldsymbol{J}^2 &= \boldsymbol{L}^2+ \boldsymbol{S}^2 + 2 \boldsymbol{L}\cdot\boldsymbol{S} = \boldsymbol{L}^2+ \boldsymbol{S}^2 + 2 L_z S_z + L_{+} S_{-} + L_{-} S_{+} \\
\boldsymbol{J}^2 |\ell \mu-s \rangle |s \rangle &= \{ \ell(\ell+1)+\frac{3}{4} + 2 (\mu-s) s \} \hbar^2 |\ell \mu-s \rangle |s \rangle \\
&+ \sqrt{\ell(\ell+1)-(\mu-s)(\mu+s)} \ \hbar^2 |\ell \mu+s \rangle | -s \rangle\\
&= \{ j^2+2s\mu \}\ \hbar^2 |\ell \mu-s \rangle |s \rangle
+ \sqrt{j^2-\mu^2}\ \hbar^2 |\ell \mu+s \rangle | -s \rangle\\
\end{aligned}
\end{equation}
ここで，$j=\ell+\frac{1}{2}$とし，$s^2=\frac{1}{4}$を用いた。これから，$\boldsymbol{J}^2$の固有状態を構成するために，$|\mathscr{j} \mu\rangle = \alpha |\ell \mu-\frac{1}{2} \rangle |\frac{1}{2} \rangle + \beta |\ell \mu+\frac{1}{2} \rangle | -\frac{1}{2} \rangle$ とする。ただし規格化条件より$|\alpha|^2+|\beta|^2 = 1$である。この２つの独立な状態に対する固有値方程式は $\boldsymbol{J}^2 | \mathscr{j} \mu\rangle = \lambda \hbar^2  | \mathscr{j} \mu\rangle$ であり，$\alpha, \beta$を用いると次のように行列形式で表される。
\begin{equation}
\begin{pmatrix} j^2+\mu & \sqrt{j^2-\mu^2} \\ \sqrt{j^2-\mu^2} & j^2-\mu \end{pmatrix}
\begin{pmatrix} \alpha \\ \beta \end{pmatrix}
= \lambda \begin{pmatrix} \alpha \\ \beta \end{pmatrix}
\end{equation}
この$\alpha, \beta$に対する連立方程式が自明でない解を持つための条件は，右辺を移行して0にしたときの左辺の行列式が0になることであり，$(j^2+\mu -\lambda)(j^2-\mu -\lambda)-(j^2-\mu^2)=0$という$\lambda$の２次方程式になる。$\therefore (\lambda -j(j+1))(\lambda -j(j-1)) = 0$から$\boldsymbol{J}^2$の固有値$\lambda \hbar^2$は，$j(j+1)\hbar^2$と$j(j-1)\hbar^2$である（したがって，$\mathscr{j}=j$ と $j-1$）。

（３）固有状態の構成
上記の$\alpha$と$\beta$の連立方程式に，固有値$\lambda$を代入して，規格化条件とともに$\alpha$と$\beta$を求めると次のようになる。
\begin{array}{c|cc||cc}
  \lambda  & \alpha & \beta & \alpha & \beta\\
  \hline
 j(j+1) & \sqrt{\dfrac{j+\mu}{2j}} & \sqrt{\dfrac{j-\mu}{2j}}
 &\sqrt{\dfrac{\ell+m}{2\ell+1}}  & \sqrt{\dfrac{\ell-m+1}{2\ell+1}} \\
  \hline
 j(j-1)  & \sqrt{\dfrac{j-\mu}{2j}} & -\sqrt{\dfrac{j+\mu}{2j}}
&\sqrt{\dfrac{\ell-m+1}{2\ell+1}}  & -\sqrt{\dfrac{\ell+m}{2\ell+1}}
\end{array}
まとめると，
\begin{equation}
\begin{array}{l}
  \begin{array}{|l|}
  \hline

  \quad |\ j\quad\ \mu\rangle = \sqrt{\dfrac{j+\mu}{2j}} |\ell \mu-\frac{1}{2} \rangle |\frac{1}{2} \rangle+\sqrt{\dfrac{j-\mu}{2j}} |\ell \mu+\frac{1}{2} \rangle | -\frac{1}{2} \rangle \quad (\lambda = j(j+1)) \quad \\
  \quad |j-1\mu\rangle = \sqrt{\dfrac{j-\mu}{2j}} |\ell \mu-\frac{1}{2} \rangle |\frac{1}{2} \rangle-\sqrt{\dfrac{j+\mu}{2j}} |\ell \mu+\frac{1}{2} \rangle | -\frac{1}{2} \rangle \quad (\lambda = j(j-1) ) \quad \\

  \hline
  \end{array} \\
\end{array}
\end{equation}

（角運動量の合成へ道（４）に続く）