ToDo:
Function[]に包まれないSlot[]は、
In[1]:= {#}
???General::slot: Undefined Slot #1:
{#}
^
???-FFS-Error-?Undefined command or element: {#}
のように、きちんとsemantics errorになるが、 Function[]に包まれないSlotSequence[]は、
In[1]:= {##}
セグメントエラー(coreを出力しました)
のように、segmentation faultになる
この結果は、あなたの予想通りですか？
FFS;
l = {
11.11111111111,
11.111111111111,
11.1111111111111,
11.11111111111111,
11.111111111111111,
11.1111111111111111,
11.11111111111111111,
11.111111111111111111,
11.1111111111111111111,
11.11111111111111111111,
11.111111111111111111111,
11.1111111111111111111111,
11.11111111111111111111111,
11.111111111111111111111111,
11.1111111111111111111111111,
11.11111111111111111111111111,
11.111111111111111111111111111,
11.1111111111111111111111111111,
11.11111111111111111111111111111,
11.111111111111111111111111111111,
11.1111111111111111111111111111111,
11.11111111111111111111111111111111,
Null[]};
Scan[With[{a = #}, Print[Map[SameQ[a, #]&, l]]]&, l]; Exit[];
{1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}
{0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}
{0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}
{0,0,0,1,0,0,0,1,1,1,1,1,0,1,0,0,1,1,1,1,1,1}
{0,0,0,0,1,1,1,0,0,0,0,0,1,0,1,1,0,0,0,0,0,0}
{0,0,0,0,1,1,1,0,0,0,0,0,1,0,1,1,0,0,0,0,0,0}
{0,0,0,0,1,1,1,0,0,0,0,0,1,0,1,1,0,0,0,0,0,0}
{0,0,0,1,0,0,0,1,1,1,1,1,0,1,0,0,1,1,1,1,1,1}
{0,0,0,1,0,0,0,1,1,1,1,1,0,1,0,0,1,1,1,1,1,1}
{0,0,0,1,0,0,0,1,1,1,1,1,0,1,0,0,1,1,1,1,1,1}
{0,0,0,1,0,0,0,1,1,1,1,1,0,1,0,0,1,1,1,1,1,1}
{0,0,0,1,0,0,0,1,1,1,1,1,0,1,0,0,1,1,1,1,1,1}
{0,0,0,0,1,1,1,0,0,0,0,0,1,0,1,1,0,0,0,0,0,0}
{0,0,0,1,0,0,0,1,1,1,1,1,0,1,0,0,1,1,1,1,1,1}
{0,0,0,0,1,1,1,0,0,0,0,0,1,0,1,1,0,0,0,0,0,0}
{0,0,0,0,1,1,1,0,0,0,0,0,1,0,1,1,0,0,0,0,0,0}
{0,0,0,1,0,0,0,1,1,1,1,1,0,1,0,0,1,1,1,1,1,1}
{0,0,0,1,0,0,0,1,1,1,1,1,0,1,0,0,1,1,1,1,1,1}
{0,0,0,1,0,0,0,1,1,1,1,1,0,1,0,0,1,1,1,1,1,1}
{0,0,0,1,0,0,0,1,1,1,1,1,0,1,0,0,1,1,1,1,1,1}
{0,0,0,1,0,0,0,1,1,1,1,1,0,1,0,0,1,1,1,1,1,1}
{0,0,0,1,0,0,0,1,1,1,1,1,0,1,0,0,1,1,1,1,1,1}
何が問題か分からないかった人は次のサンプルをご覧有れ！
String2RealConvTest[max_Real:40] := Module[{reference, min, loop, s, i},
reference = 1 / 11;
s = "0.";
loop = "09";
min = 1;
Table[(s = s//loop[1 + Mod[i - 1, StringLength[loop]]];
Print[StringFill[s, " ", max + 2], " - 1 / 9 -> ",
ToExpression[s] - reference]), {i, min, max}]];
String2RealConvTest[]; Exit[];
0.0 - 1 / 9 -> -.090909090909091
0.09 - 1 / 9 -> -9.090909090909E-4
0.090 - 1 / 9 -> -9.090909090909E-4
0.0909 - 1 / 9 -> -9.090909090917E-6
0.09090 - 1 / 9 -> -9.090909090917E-6
0.090909 - 1 / 9 -> -9.090909090792E-8
0.0909090 - 1 / 9 -> -9.090909090792E-8
0.09090909 - 1 / 9 -> -9.09090913659E-10
0.090909090 - 1 / 9 -> -9.09090913659E-10
0.0909090909 - 1 / 9 -> -9.09090858148E-12
0.09090909090 - 1 / 9 -> -9.09090858148E-12
0.090909090909 - 1 / 9 -> -9.0913387929E-14
0.0909090909090 - 1 / 9 -> -9.0913387929E-14
0.09090909090909 - 1 / 9 -> -9.15933995316E-16
0.090909090909090 - 1 / 9 -> -9.15933995316E-16
0.0909090909090909 - 1 / 9 -> -1.38777878078E-17
0.09090909090909090 - 1 / 9 -> -1.38777878078E-17
0.090909090909090909 - 1 / 9 -> 0
0.0909090909090909090 - 1 / 9 -> 0
0.09090909090909090909 - 1 / 9 -> 0
0.090909090909090909090 - 1 / 9 -> 1.387778780781E-17
0.0909090909090909090909 - 1 / 9 -> 1.387778780781E-17
0.09090909090909090909090 - 1 / 9 -> 1.387778780781E-17
0.090909090909090909090909 - 1 / 9 -> 1.387778780781E-17
0.0909090909090909090909090 - 1 / 9 -> 2.775557561563E-17
0.09090909090909090909090909 - 1 / 9 -> 2.775557561563E-17
0.090909090909090909090909090 - 1 / 9 -> 2.775557561563E-17
0.0909090909090909090909090909 - 1 / 9 -> 1.387778780781E-17
0.09090909090909090909090909090 - 1 / 9 -> 2.775557561563E-17
0.090909090909090909090909090909 - 1 / 9 -> 2.775557561563E-17
0.0909090909090909090909090909090 - 1 / 9 -> 2.775557561563E-17
0.09090909090909090909090909090909 - 1 / 9 -> 2.775557561563E-17
0.090909090909090909090909090909090 - 1 / 9 -> 1.387778780781E-17
0.0909090909090909090909090909090909 - 1 / 9 -> 1.387778780781E-17
0.09090909090909090909090909090909090 - 1 / 9 -> 1.387778780781E-17
0.090909090909090909090909090909090909 - 1 / 9 -> 2.775557561563E-17
0.0909090909090909090909090909090909090 - 1 / 9 -> 1.387778780781E-17
0.09090909090909090909090909090909090909 - 1 / 9 -> 2.775557561563E-17
0.090909090909090909090909090909090909090 - 1 / 9 -> 2.775557561563E-17
0.0909090909090909090909090909090909090909 - 1 / 9 -> 1.387778780781E-17
つまり、十進数で記述の精度を真の値に向けて増大させたときに、 SAD内部表現であるReal型の値が真の値に収束するわけではない ということ
例えば、実数のデータ列をWrite[]で書き出した結果をRead[]してくると Real型としては元の値との同一性が保証されないということ。 実例は以下のサンプル参照のこと。
例えば、数値計算の精度をテストするユニットテストの記述時に 標準値をユニットテスト内に記述した段階で、その値が 元の計算値から変化する可能性が有ります。
String2RealConvTest[max_Real:10] := Module[{x, s, i}, StandardForm[
x = 1 - 1e-3;
Table[(x = (1 + x) / 2;
s = ToString[x];
Print[(" "//ToString[i])[-3,-1],
": x=", x, " s=", s, " diff=", ToExpression[s] - x];
), {i, 1, max}];
]];
String2RealConvTest[45]; Exit[0];
1: x=.9995 s=".9995" diff=0
2: x=.99975 s=".99975" diff=0
3: x=.999875 s=".999875" diff=-1.11022302463E-16
4: x=.9999375 s=".9999375" diff=0
5: x=.99996875 s=".99996875" diff=-1.11022302463E-16
6: x=.999984375 s=".999984375" diff=0
7: x=.9999921875 s=".9999921875" diff=0
8: x=.99999609375 s=".99999609375" diff=0
9: x=.999998046875 s=".999998046875" diff=0
10: x=.9999990234375 s=".9999990234375" diff=0
11: x=.99999951171875 s=".99999951171875" diff=0
12: x=.999999755859375 s=".999999755859375" diff=0
13: x=.999999877929688 s=".999999877929688" diff=5.551115123126E-16
14: x=.999999938964844 s=".999999938964844" diff=2.22044604925E-16
15: x=.999999969482422 s=".999999969482422" diff=1.110223024625E-16
16: x=.999999984741211 s=".999999984741211" diff=1.110223024625E-16
17: x=.999999992370605 s=".999999992370605" diff=-4.4408920985E-16
18: x=.999999996185303 s=".999999996185303" diff=2.22044604925E-16
19: x=.999999998092651 s=".999999998092651" diff=-3.33066907388E-16
20: x=.999999999046326 s=".999999999046326" diff=2.22044604925E-16
21: x=.999999999523163 s=".999999999523163" diff=1.110223024625E-16
22: x=.999999999761581 s=".999999999761581" diff=-3.33066907388E-16
23: x=.999999999880791 s=".999999999880791" diff=3.330669073875E-16
24: x=.999999999940395 s=".999999999940395" diff=-3.33066907388E-16
25: x=.999999999970198 s=".999999999970198" diff=3.330669073875E-16
26: x=.999999999985099 s=".999999999985099" diff=2.22044604925E-16
27: x=.999999999992549 s=".999999999992549" diff=-4.4408920985E-16
28: x=.999999999996275 s=".999999999996275" diff=2.22044604925E-16
29: x=.999999999998137 s=".999999999998137" diff=-3.33066907388E-16
30: x=.999999999999069 s=".999999999999069" diff=2.22044604925E-16
31: x=.999999999999534 s=".999999999999534" diff=-3.33066907388E-16
32: x=.999999999999767 s=".999999999999767" diff=-2.22044604925E-16
33: x=.999999999999884 s=".999999999999884" diff=3.330669073875E-16
34: x=.999999999999942 s=".999999999999942" diff=2.22044604925E-16
35: x=.999999999999971 s=".999999999999971" diff=1.110223024625E-16
36: x=.999999999999985 s=".999999999999985" diff=-4.4408920985E-16
37: x=.999999999999993 s=".999999999999993" diff=3.330669073875E-16
38: x=.999999999999996 s=".999999999999996" diff=-3.33066907388E-16
39: x=.999999999999998 s=".999999999999998" diff=-2.22044604925E-16
40: x=.999999999999999 s=".999999999999999" diff=-1.11022302463E-16
41: x=1 s="1" diff=4.440892098501E-16
42: x=1 s="1" diff=2.22044604925E-16
43: x=1 s="1" diff=1.110223024625E-16
44: x=1 s="1" diff=0
45: x=1 s="1" diff=0
LHC Collision OpticsでHorizontal Dispersionを条件に入れると、 CELLモードでなんかうまくマッチングが取れない。
Corssing Angle無しのモデルだと問題なくできるので、 IPの Vertical Crossing Angleが諸悪の根源のようです。
IPの Vertical Bump付近に Couplingも残っているようなので それも含めると自由度足りないが、Couplingパラメータは拘束していない はずなのだが...
Couplingの変化がIRの外に伝搬して境界での Twissパラメータの拘束条件に 悪影響を与えてるのかな？
コミット数の調査ツール(contrib/tools/CommitActivityChecker)を svkによるミラーに対応させたので、手元に有る amorita branchの Subversion repositoryのミラーを調べてみる。
要するに、CVS repositoryでのoldsadモジュール相当の部分
| Committer | 2005 | 2006 | 2007 | 2008 | Total |
| amorita | 191 | 219 | 383 | 326 | 1119 |
| Total | 191 | 219 | 383 | 326 | 1119 |
(2008/04/02 12:30 +0200現在)
拡張モジュール開発用の保管庫
| Committer | 2005 | 2006 | 2007 | 2008 | Total |
| amorita | 17 | 6 | 140 | 125 | 288 |
| Total | 17 | 6 | 140 | 125 | 288 |
(2008/04/02 12:30 +0200現在)
だそうです。
通販や法人営業は残るそうなので、ぷらっとホームと同様に 利益率の下がった店頭販売から撤退ということのようです。
これで、秋葉原からTyanとかSupermicroの現物を見られる店が死滅か...
割高だけど、ECC付きメモリーが確実に在庫しているとか、SMPや RAID5/6等の システムに対する技術情報の蓄積も多く、新しいシステム組む際に相談するには、 良い所だったのですが、これも時代の流れですかね... (コンピュータがコモディティ化しすぎた)
GCC 4.3.0がリリースされ、GCC4系列の開発は 4.4へ移行したことですし... amorita branchでは Revision 1700をもって、 GCC 3.4.6のEnd of Supportを宣言します。
別にconfig/GCC.specから GCC 3.4.6向けのコードを 削除するわけでは有りませんが、何も指定がない場合は ビルドフレームワークはGCC 4.3.0を想定します。 また、amorita branchでのBuild Testの対象から GCC 3.4.6が外れます。
1.0.10.2.4b以降のバックポートがだいぶ溜まったので、 切りの良いところでCVS MAIN trunkのバージョンを 1.0.10.2.5aへバンプしてSAD_1_0_10_2_5aタグを 打ち込みました。
基本的には、一件一件コードを精読して修正を繰り返す作業。 新規のコードマージ時に査読をきちんと行えれば、 既存コードの修正は地道にやれば確実に終わる。
参照例が実在し、実装の方向性は固まっているので、 APIを定義してコードを書くだけ。 ただし、需要動向は不明(SAD Workshop 2006以降、問い合わせは無い)
関連作業として、src/inc/TMACRO1.incのC側からの参照手段の整備も必要。 (Cで、組み込みルーチンを記述する場合)
内容的には、Beam Line Element定義の動的組み込みに近いが ドライバーフレームワーク側に新規にホックを仕込むので 設計に際して注意が必要
長期的な保守を考えれば、これ以上 SPAC系ルーチンが増える前に実装すべき案件
最近の PC事情では、4 Core/8 Coreの Workstationが導入されることも 多いく、動作クロックの上昇が停滞している状況から並列実行性能は ターンアラウンドタイム短縮に重要。
可動範囲とか、各種の副作用検証用にパネル形式のテストコードを組み立てる。
Dispersion matchingまわりは、妙な副作用が多いので最初の実装では 取り込まずにすませることに
要するに誰も使用しない関数のみで構成されるファイルの一斉処分
一応、処分前に保守状況を調査(このまえ作った 調査用 Subversionレポジトリが活躍してくれました)
保守が必要でないであろう行列操作の基本関数等を除いて、 ある程度以上の複雑さを備えたコードが一度も保守されていないということは 1995年当時から実質死んでいたと思われる
1回コミットされているコード。 修正内容は、implicit文や変数の型宣言の機械的な修正のように見える 前後数ヶ月のコミットログから察するに HP-UXへのポーティング作業に 付随する修正と思われる。(コンパイルエラー等の解消目的だとすると コードが実際には利用されていなかった可能性がある)
CVS repositoryが導入された 1995〜96年代にかけて実際に TWISSの計算に 使われていたいたが、途中から使用されなくなったものと思われる。 また、2000/08/19の段階で、呼び出し元は存在していない。
C.K.Allenの実装した MatrixLog[]なのだが、精度が悪い＆収束性がよろしく ないので、直せないかゴソゴソ調査していたのだが...
結論、アルゴリズム的に数値不安定性が有って無理！
と言うわけで、まったく別のアルゴリズムで再実装しました。
実数の場合と同様に Logarithmを Exponentialの逆関数で定義する。
行列$\mathbb{F}$の Logarithmである行列 $\mathbb{A}$は、 $MatrixExp[\mathbb{A}] = \mathbb{F}$を満足する行列。
Exponentialは実数の場合の定義をそのまま拡張して定義する。 例えば、積展開の極限
$MatrixExp[\mathbb{A}] = \lim_{n\rightarrow\infty}\left(1 + \frac{\mathbb{A}}{n}\right)^n$
もしくは、Taylor級数展開
$MatrixExp[\mathbb{A}] = \sum_{n=0}^{\infty}\frac{\mathbb{A}^n}{n!}$
実際の計算では、任意の精度で計算を打ち切るために Taylor級数展開が 使われます。重要なのは、このTaylor級数展開においては 行列のノルム $||\mathbb{A}||$が有界ならば級数が絶対収束し、 級数展開の収束半径は無限大であることです。
複素数まで拡大して、一般化すると有界な複素正方行列$\mathbb{A}$の Exponentialは上記の Taylor展開で一意に定義できます。
一般に逆関数を解くのは骨なのでExponentialと同様に Logarithmを級数展開してみると、 Taylor級数展開では
$MatrixLog[\mathbb{I} - \mathbb{G}] = -\sum_{n=1}^{\infty}\frac{\mathbb{G}^n}{n}, \mathbb{G} = \mathbb{I} - \mathbb{F}$
となります。この例では収束半径は、$||\mathbb{G}|| < 1$です。
つまり、$\mathbb{F}$の全ての固有値が複素平面上で1を中心とする単位円に 入っている必要があります。 つまり、どう頑張っても任意行列$\mathbb{A}$のExponential $MatrixExp[\mathbb{A}]$に対してLogarithmの値を定めるには Taylor展開では役不足と言うことです。
ただし、$\mathbb{F}$のすべての固有値が右半面にあれば、 適当なスケール変換を施した行列$\lambda\mathbb{F}$は 収束半径に入るので次の式で級数展開による定義が可能になります。
$MatrixLog[\mathbb{F}] = MatrixLog[\lambda\mathbb{F}] - Log[\lambda]\mathbb{I}$
したがって、一般の行列$\mathbb{F}$に対するLogarithmの計算では、
という方針になります。
(1)の方針に相当し、コード中ではquadratic modeなどと 書かれているあたりが該当します。
行列$\mathbb{F}_s$のLogarithmの近似値$\mathbb{A}_n$を 次の処方で改良しています。
$\mathbb{X}_n = MatrixExp[-\mathbb{A}_n].\mathbb{F}_s$
$MatrixExp[\mathbb{A}_{n+1}] = MatrixExp[\mathbb{A}_n].X_n$
もちろん、このままでは$MatrixExp[\mathbb{A}_{n+1}]$の逆関数が 解けないので、箱の鍵は箱の中状態です。 ここで、近似が十分に良い($\mathbb{x}_n$が単位行列$\mathbb{I}$に近い)と 仮定すると、$\mathbb{X}_n$のLogarithmを Taylor級数展開で与えること が可能になり、
$\mathbb{W}_n = MatrixLogTaylor[\mathbb{X}_n]$
$MatrixExp[\mathbb{A}_{n+1}] = MatrixExp[\mathbb{A}_n].MatrixExp[\mathbb{W}_n]$
を得ます。 ここで、行列$\mathbb{X}_n$が単位行列に十分近ければ、 行列$\mathbb{W}_n$は零行列に近いので、 Baker-Campbell-Hausdorff級数展開(Zassenhaus式)を用いて
$\mathbb{A}_{n+1} = \mathbb{A}_n + \mathbb{W}_n + \frac{1}{2}\{\mathbb{A}_n, \mathbb{W}_n\} + \frac{1}{12}(\{\mathbb{A}_n, \{\mathbb{A}_n, \mathbb{W}_n\}\} + \{\{\mathbb{A}_n, \mathbb{W}_n\}, \mathbb{W}_n\}) - \frac{1}{24}\{\mathbb{A}_n, \{\mathbb{W}_n, \{\mathbb{A}_n, \mathbb{W}_n\}\}\} + ...$
が得られます。
で、実際のデータ(乱数生成した行列のExponential)を流して $MatrixLogTaylor[\mathbb{X}_n]$の入力を見ていると 2次収束していた$\mathbb{X}_n$が突然発散に転じてあっという間に 無限大になってしまうとい現象に遭遇することになります。 (まあ、それ以前にオリジナルのコードにはバグが有るのだが、それとは独立な現象)
結論から言うと、この手法では解は2次収束なのだが、 数値的に不安定なので演算誤差や打ち切り誤差を種に 指数関数的に発散が発生している。 つまり、数列 $\mathbb{A}_1,\mathbb{A}_2,\mathbb{A}_3,...$に 指数関数的増大列が含まれることが諸悪の根源なので、直せません。Orz (発散の直前で止めることで、それなりに有意な解を得られるが、 精度がそこで制限されてしまう)
(2)の方針に相当し、$\mathbb{F}$を変形して単位行列に近づけます。
アイデア自身は、MatrixExp[]の実装でノルムの大きな行列を計算する際に 使っている次の恒等式そのものです。(ExponentialのTaylor級数展開は 絶対収束だが、ノルムの小さな行列ほど収束が早い)
$MatrixExp[\mathbb{A}] = MatrixExp[\frac{\mathbb{A}}{n}]^n$
実際のコードでは、$n$に$2^k$を使っています。
で、これに$\mathbb{F}=MatrixExp[\mathbb{A}]$を代入してLogarithmをとると
$\mathbb{F} = n MatrixLog[\mathbb{F}^{\frac{1}{n}}]$
が得られます。 ここで、例えば$n = 2$を選ぶと
$\mathbb{F} = 2 MatrixLog[MatrixSqrt[\mathbb{F}]]$
となります。 ここで注目すべきは、平方根の主値をうまく選ぶと $MatrixSqrt[\mathbb{F}]$の固有値分布では$\mathbb{F}$の固有値分布と 比べて各固有値の偏角$\theta$が半分になる点です。 つまり、複素平面全域にあった固有値が右半面に集めることが出来ます。 ようは、平方根を複数回作用させて固有値分布を右半面実軸近傍に 集めれば、Taylor展開を使ってLogarithmを計算できるわけです。
ここで問題となるのは平方根の実装で、ExponentialやLogarithmを 使わずに求めなければ、意味が有りません。 これを愚直に次のような Newton法で実装してみました。
$\mathbb{A} = MatrixSqrt[\mathbb{B}]$
$\mathbb{A}_{n+1}=\frac{1}{2}(\mathbb{A_n}+\mathbb{A_n}^{-1}\mathbb{B})$
$\mathbb{A} = \lim_{n\rightarrow\infty}\mathbb{A}_n$
これによる前処理で、だいぶマシな状態になるのですが、この漸化式自体が 数値的不安定性を持っているという罠が有りました。
つまり、精度を上げようと反復を増やすと発散してしまう問題は 本質的には改善していないと言うことです。
で、何か手段はないかと Wikipedia(英語版)などをあさっていると... Denman-Beavers square root iterationsという 数値的に安定な行列の開平法を見つけました。
$\mathbb{M}_{n+1} = \frac{1}{2}\left(\mathbb{I} + \frac{\mathbb{M}_n + \mathbb{M}_n^{-1}}{2}\right), \mathbb{M}_0 = \mathbb{A}$
$\mathbb{Y}_{n+1} = \frac{1}{2}\mathbb{Y}_n(\mathbb{I}+\mathbb{M}_k^{-1}), \mathbb{Y}_0 = \mathbb{A}$
$\mathbb{Z}_{n+1} = \frac{1}{2}\mathbb{Z}_n(\mathbb{I}+\mathbb{M}_n^{-1}), \mathbb{Z}_0 = \mathbb{I}$
この漸化式は数値的に安定に2次収束し、その極限は
$\lim_{n\rightarrow\infty} \mathbb{M}_n = \mathbb{I}$
$\lim_{n\rightarrow\infty} \mathbb{Y}_n = \mathbb{A}^{\frac{1}{2}}$
$\lim_{n\rightarrow\infty} \mathbb{Z}_n = \mathbb{A}^{-\frac{1}{2}}$
に収束します。
これの周辺論文を調べていたら、応用数学分野の2001年ごろの論文で そのものズバリな下記の論文を見つけました。
APPROXIMATING THE LOGARITHM OF A MATRIX TO SPECIFIED ACCURACY,
S.H.Cheng, N.J.Higham, C.S.Kenney, and A.J.Laub,
SIAM J. Matrix Anal. Appl. Vol.22 No.4 pp.1112-1125
で、現在開発版でのMatrixLog[]はこの論文に基づいて実装しています。
最後に行う単位行列に近い行列のLogarithmも論文と同様に Taylor級数展開ではなく、Pade近似による有理級数展開で実装しています。 (2/2次のPade近似よりも旧実装の MatrixLog[]を突っ込んだ方が精度が悪かった罠)
現在、内部実装のチューニングと精度評価部の調整中 (論文のアルゴリズム自身は、任意精度で計算できるが実装は 倍精度浮動小数点のマシンイプシロンで制約されている)
関数$f(x)$を次の形の有理和で近似する。
$r_{mm}(x) = \sum_{j=1}^{m}\frac{\alpha_j^{(m)}x}{1+\beta_j^{(m)}x}$
$\log{(1-x)}$の近似式で低次のものは
$r_{11}(x) = \frac{-x}{1-\frac{1}{2}x}$
$r_{22}(x) = \frac{-\frac{1}{2}x}{1-\frac{1}{2}\left(1+\frac{1}{\sqrt{3}}\right)x} + \frac{-\frac{1}{2}x}{1-\frac{1}{2}\left(1-\frac{1}{\sqrt{3}}\right)x}$
$r_{33}(x) = \frac{-\frac{4}{9}x}{1-\frac{1}{2}x} + \frac{-\frac{5}{18}x}{1-\frac{1}{2}\left(1 - \sqrt{\frac{3}{5}}\right)x} + \frac{-\frac{5}{18}x}{1-\frac{1}{2}\left(1 + \sqrt{\frac{3}{5}}\right)x}$
Y氏のコメントに有るように、LLVM-GFORTRANの equivalence文の実装が直ったようなので、ベンチマークをとってみました。
使用した llvmは、Revision 49485です。
| Function | Optics | Tracking | Matching | Overall |
| 4.979576766490936 | 7.273632526397705 | 10.365879058837891 | 6.395431518554688 | .3689253422476 |
| Function | Optics | Tracking | Matching | Overall |
| 4.9791539311409 | 7.322125911712646 | 10.631559371948242 | 6.909873962402344 | .377546692887942 |
-msse2/3をつけるとSegmentation faultする。llvmの問題？
| Function | Optics | Tracking | Matching | Overall |
| 4.251880943775177 | 7.825789451599121 | 9.968093872070313 | 7.023014068603516 | .359744608402252 |
なんか、あまり差が出ない...Orz
残り警告数 116個まで減少
| ソース | 警告数 |
| src/phsrot.f | 1 |
| src/qdbend.f | 14 |
| src/qdcell.f | 4 |
| src/qdthin.f | 2 |
| src/qdtwis.f | 3 |
| src/qgettr.f | 1 |
| src/qtwiss.f | 9 |
| src/tbende.f | 5 |
| src/tcsvdm.f | 2 |
| src/temit.f | 3 |
| src/tffsmatch.f | 10 |
| src/tmulta.f | 4 |
| src/tmulte.f | 9 |
| src/tmulti.f | 19 |
| src/tquade.f | 1 |
| src/tsgeo.f | 1 |
| src/tsolque.f | 24 |
| src/tsteee.f | 2 |
| src/tthine.f | 1 |
| src/tturn.f | 1 |
残り警告数 93個まで減少
| ソース | 警告数 |
| src/phsrot.f | 1 |
| src/qdbend.f | 14 |
| src/qdcell.f | 4 |
| src/tcsvdm.f | 2 |
| src/temit.f | 3 |
| src/tffsmatch.f | 10 |
| src/tmulta.f | 4 |
| src/tmulte.f | 9 |
| src/tmulti.f | 19 |
| src/tsolque.f | 24 |
| src/tsteee.f | 2 |
| src/tturn.f | 1 |
あと、tturn1 @ src/tturn1.fも SOLエレメント絡みの動作が 怪しそうなので、残りの調査が済んだらTrackParticles[]で 実際の動作を検証してみる。
残り警告数 63個まで減少、残存する未調査領域は 4ファイル(56警告)
| ソース | 警告数 |
| src/phsrot.f | 1 |
| src/tcsvdm.f | 2 |
| src/temit.f | 3 |
| src/tmulta.f | 4 |
| src/tmulte.f | 9 |
| src/tmulti.f | 19 |
| src/tsolque.f | 24 |
| src/tturn.f | 1 |
Fail-Safe Cという メモリ安全なANSI Cの実装が有るらしい。
warning: 'foo' may be used uninitialized in this function警告の調査完了
| ソース | 警告数 |
| src/phsrot.f | 1 |
| src/temit.f | 3 |
| src/tturn.f | 1 |
まずは、バックポートを処理するかな
-finit-real=[zero|nan|inf|-inf]なるオプションが有る！
gfrotranビルド時は、-finit-real=nanでもつけておくか？
コメントにsrc/bb.fのコピーっぽいので i_p0は533ではないかという意見が出ているが、 BEAMBEAMエレメントのblist構造は長さ600だが、 PHSROTエレメントのblist構造は長さ240しか無いので、 明らかに違う。
正しいi_p0は、1から240までの整数のはず。 また、phsrot関数内の定数nblist(おそらく、blistの 長さを与える定数として導入されたと思われる)の値が200で、 src/tpara.fのtpara関数のラベル3700で確保している 配列長が240であることから見ても、以前の考察通り、 1995年に CVS repositoryが導入された時点で保守されていなかった＆ 死んでいたという仮説が正しいと思われる。
次の作業はタグ打って 数学/物理定数の整理と物理定数更新(1627...1691) かな？
extensions/Standard/SpaceCharge/Scheffに関しては需要なさそうなんで、 バックポートせずに塩漬けの予定
昨今、いろいろなところでネットワーク上に流れるデジタル化された著作物の 著作権がらみの問題が報じられているが、昔から私が持っている疑問が有ります。
「表現をデジタル化(量子化)した著作物のあり得る数は 加算無限個に制限されるのでは？」と言うものです。
議論を簡単にするために、 全ての「US-ASCIIで符号化された有限長の文書」を含む 集合Xを考えます。
まず、集合の個々の要素は有限長の符号列で符号体系はで離散符号なので、 適切な定義を行えば集合は全順序集合になります。
集合Xの元xの長さをLength[x]'、 長さN以下の符号のみで出来たXの部分集合をX(N)として、 関数Order(x)を部分集合X(Length[x])を短い符号列が前に くるようにstrcmpで整列させた時にxの前から数えた順位を 返す関数として定義します。
したがって、集合Xの任意の元x,yに対して
が成り立ちます。ここで、Nullは長さ0の符号列です。
集合Xの任意の元は自然数と一対一に対応付けられます。 したがって、集合Xの濃度は加算無限(整数や有理数と同じ)です。
この著作物符号化空間内での著作権の侵害の可能性を考えてみましょう。 議論の単純化のために、著作権の消滅などは考えずに、純粋に 互いの著作物が侵害する恐れが有るか無いかを判定する 関数Conflict(x,y)を 「xがyを侵害する可能性がある場合にtrueを返し、 そうでない場合にfalseを返す」と定義します。
これを使って何をやるかというと、早い話が「エネラウスの篩」です。 まず、既知の著作物集合A_0を用意します。 ここで、「∀a∈A_n, Conflict(x_n, a) = false」満足する Xの元x_nを用いて、 集合A_nの漸化式「A_{n+1} = A_n ∪ {x_n} 」を定義します。 さて、Xの元は加算無限個しか存在しないので、元x_nは 順位の小さい源から順に試していけば良いわけです。 しかも、A_nの漸化式の定義から「Order(x_n) < Order(x_{n+1})」が 満たされますから、一度検査した元を再検査する必要は有りません。 漸化式が、途中で止まるならもう新しい著作権侵害の恐れのない著作物が ないということです。
つまり、「任意の符号化著作物に対して著作権侵害を判定可能」 であるなら、衝突判定関数Conflictを構成可能で、 これを使って既知の著作物空間と衝突しない元からなる 部分集合X ∩ !A_0から発見済みの元(既知の著作物を含む)に 衝突しない元x_nを次々と選び出すことができると言うことです。
現実のリソースで表現可能な範囲は、集合Xの部分集合でしかなく、 ある十分に大きな自然数Mに対して現実のリソースで表現可能な Xの部分集合X*は、「X* ⊂ X(M)」を満足します。 X(M)は、たかがか有限集合なので たかがか有限の計算機リソースを たかだか有限時間稼働させることで、未発見の元x_nを すべて発見出来ます。
つまり、極論すれば 「任意の符号化著作物に対して著作権侵害を判定可能」であれば、 「有限リソースで表現可能な符号化著作物空間を自動的に埋め尽くすことが可能」 であるということ。
"\\r"の扱いが変わっているのは、腐れ構文も受け入れていた MAIN Levelの式パーサーを直したからだと思うのですが...
この変更(src/calc.[chy]と src/yylex_.c)を「バグだ！修正しろ！」 と言う人には論拠となる仕様書とか実装設計書とかを見せてもらいたい。
バグ報告が上がってるんだけど、開発環境とか KEKB標準のコンパイル環境(PowerMac G5/Mac Pro)では再現しない。
多分、報告者の環境固有の問題臭いんですけど、 どんな環境で実行したのか書いてないから誰も追試出来ない罠。
つーか、ソース配布なんだからビルド環境の問題とソース側の問題の 切り分けが出来ていないなら、ビルド環境の状態を報告してくれないと 追試が出来ない。 もっとも、予想通り MacOSXだとしたら手元に環境を自由にいじれる 試験機が無いので調べようがない。
原因が分からなければ、DarwinにBROKEN付けて終了かな？
報告者のコンパイル環境の不良が原因だったらしぃ
なんか、最近のSADのバグ報告は GCCに含まれる Fortranコンパイラの 特定環境での不良絡みが上がってくるのですが、Cコンパイラの場合は 自分で自分をビルドするので、自己生成の結果を世代間で比較することで、 コンパイラ自身のバグが発覚しやすいのですが... 自己生成でないコンパイラは、テストセットを走らせない限り ほとんど発覚する機会が無い罠
みんなコンパイラツールチェーンをきちんと検証しようぜ(特に、Fortran)！
クエンチを避けるため7TeVではなく、5TeV衝突で立ち上げるようですが、 最近LCU(LHC Commissioning & Upgrade)ミーティングでも 5TeV時の アパーチャや光学系パラメータを検討した話が話題になっています。
暴走するというバグレポート＆調べろ指令が来てたんで、 KEK内のホストに Xvnc Server環境を構築してそこ経由で 調査を3時間ぐらい実施。
報告内容は、「特定ホストalsad3(系列機が存在しないので 特定アーキテクチャとも言える)で特定のSADパネル (X Clientなので、そのまま CERNへ連れてくると画面開くだけで 10分以上かかる)で、KBFOpenDialog[]するとスタックする」 と言うものだったのですが、聞いたこともない＆発生率100%という かなり変な症状です。(確率的な障害なら、レースコンディションを 含んだ実装等に思い当たる節が無いわけでは無いのですが...)
あからさまに変です。なんで大量に fork(2)してるんだよ...Orz
まず、fork(2)システムコールがおかしかったら OSがまともに 動くはずも無いので、この可能性は除外というか、SADの責任じゃないよ！
KBFOpenDialog[]内で連続してFork[]が呼び出されるのは、 コードフロー的にあり得ないというか、そういう事態が起きるとすれば、 インタープリターは既に正常に動作していない状況なわけで... 明らかに実行モジュールsad1.exeが、 内部で問題(多分メモリ系)を起こしているか、腐っているかという線が濃厚。
で、他人の作った実行モジュールにプローブ差し込んで調べるのは 面倒なので、当該環境でクリーンビルドな実行モジュールを生成して 検証してみたのですが...障害が再現しないよ Orz
いつものように amorita branchで検証してたんで、CVS MAIN trunk固有の 障害を疑って CVS HEADでも検証してみるも、再現せず。なにゆぇ...
問題の実行モジュールのVersion Stringから その当時のソースツリー構成を 推定(前後数日の誤差有り)して構築してみるも、再現せず。
この実行モジュールって、どうやって作ったんだ！
Sticky Tagで意図的にソースを固定せずに、cvs checkoutしたままの ソースツリーで作られたとすれば、再現したソースツリーの構成は かなり近いはず(関与するコンポーネントという意味で)なので、 コンパイル環境かコンパイル設定の不良とかストレージの不良とかが 怪しそうなんですけど... Orz
なんか、3時間ぐらい検証作業やってたけど虚しすぎ...
バグ報告の前に、環境をクリーンビルドして障害が再現するか 試験して欲しかった
LHCB1の IP1/5/8(for ATLAS/CMS/LHCb)での基本動作を確認。 遅いという致命的な問題が有るけど、ノブの可動範囲とかηyへの副作用 (ノブ領域にV-Dispersion Sourceがあるが、補正用の自由変数が無いので 制御できない)の評価には十分なので、しばらくはこの実装で進める。 次は、LHCB2の定義(Crossing Bumpの定義)を調べて対応させる作業かな？
KEKBだと、Closed Bumpで使う Steeringのリストと拘束条件だけ 記述してその場でマッチングするので、見ればどの Magnetを 自由変数にしているか分かるんですけど、LHCのラティス定義 (シーケンスファイル)にはマグネットの K値が名前付きの変数を入れてあって、 strファイル(おそらく strengthファイルの意)側に ノブのオンオフを記述する変数(0と 1が入る)と それを使ったパラメータの関係式(一次式)で記述されているので、 いちいち、シーケンスファイルと strファイルを行ったり来たり しないとClosed Bumpに使った自由変数が分からないのが難。 おまけに、リング全体をマッチングする台本は repositoryに 無いようなので、各種拘束条件も自分で決めないといけない罠。 (Frankにも聞いたけど、どうやらリング全体をマッチングする 公式台本は無いらしい。設計グループは内部にはそれに相当するものが 有るはずだが...)
file system mount回りの問題で動かないようです。
諦めて、AFSからコントロール室の計算機は必要なファイルを すべてコピーして試験しようという話になったのだが、 AFS Serverが落ちていてファイルが取ってこれないようです。 というとで、本日の試験終了〜
Fortranの規格外部分を標準規格で置き換える or 外部ライブラリ化する 作業に着手。FNUM/FGETC/FSEEKに関しては、Fortran I/O絡みなので コンパイラ非依存に外部ライブラリ化するのは無理なので後回し。
gfortran 4.4.0で Fortran2008モードでコンパイルしてエラーや 警告メッセージを調査中
EQUIVALENCE絡みはもろに違反しているような...
GCC上のビルドだと見かけ上は計測さえれた CPU時間の桁数が 減って見えますが、有効精度自体は向上しています。
Revision 1886以前のコードでは、積算CPU時間の取得に 非標準のSECOND関数を使用していましたが、これの返り値は 単精度実数であり差分を取ってマイクロ秒まで見ようとすると 有効数字が不足していました。
これに対して、Revision 1887からは Fortran95規格で 標準化されたCPU_TIMEサブルーチンを使用して 倍精度実数で積算CPU時間を取得するようになったので、 差分を取った際にマイクロ秒以下へ桁落ちに伴う誤差が 現れにくくなったので、結果的に表示される桁数が 減少しています。 (積算CPU時間を得るためのgetrusage(2)はマイクロ秒精度です)
カテゴリー: Admin | Emacs | EPICS | Fortran | FreeBSD | GCC | hgsubversion | IPv6 | KEKB | LHC | Lisp | LLVM | MADX | Ryzen | SAD | samba | tDiary | unix | WWW | YaSAI | お仕事 | イベント | 出張 | 宴会 | 数学 | 艦これ | 買いもの | 追記 | 雑記
Before...
_ wallaroo [覚えてるかわからんけど久々に。 ゆーざーずは実物見に行ける（＆相談できる）唯一のショップだっただけに閉店は残念。 ..]
_ A. Morita [うわ、めちゃめちゃお久しぶりです>wallaroo 事前に歳末セールもやっていたらしいですからねぇ]
_ wallaroo＠Y.Yokose [今日行って来たよ〜。 CPU関連は10%OFFで微妙。40or60%OFFは元々安いマウスとRAID板、FANと細..]