ポケモンGO 個体値計算

| プログラムメモ | 2017/03/28 Tue 00:18
ポケモンGOでは、多くのトレーナーが最強のポケモンを探しています。
しかし、100%個体値と呼ばれる最強ポケモンが手に入ることは、めったにありません。
なので、90%や80%で我慢するわけですが、この％で表せられる個体値とは何だろうと思い、調べてみました。

各ポケモンには、基本となる種族値（攻撃力・防御力・体力）が設定されていて、
それぞれの個体には、＋αとなる、攻撃力・防御力・体力が与えられる。
ゲーム画面上で表示される値は、強さを表すCPと体力を表すHPだけが表示されている。





個体特有値(Individual)は、攻撃力・防御力・体力それぞれに、0〜15の範囲で決定される。

ポケモンには、1〜40のレベルがあり、0.5キザミで指定される。
このレベルは、強化することによって、0.5上げることができる。
ポケモンレベルにより、CP補正値が一意に定まる。

ポケモン毎の種族値や、CP補正値については、ネット上に沢山の情報が出ている。

CPやHPの計算方法がわかったところで、
では、自分の持っているポケモンの個体値は、どうしたらわかるのだろう。
ネット上には、個体値計算機なるものがあるが、どうやって計算しているのだろう。

HPから考えてみる。
ポケモンの名前から種族値はわかる。「強化に必要な砂」からおおよそのポケモンレベルがわかる。
ならば、レベルの該当範囲内で片っ端からHP計算を行い、ゲーム画面に出ているHPと合致すれば、それが正解となる。
しかし、複数の正解が出ることも考えられるので、これらは正解候補となる。

CPも同じように計算を行うが、既にHP計算でレベル候補が絞られているので、その分だけで済む。
ゲーム画面に出ているCPと合致すれば、やはりこれも、正解候補となる。

正解候補をテーブルにして表示すれば、捕まえたポケモンの個体値の目安となる。
チームリーダーによる評価がわかれば、更に正解候補を絞り込むことができる。

ここまで来れば、後はプログラムに落とし込むだけです。

おそらく、このプログラムの難しさは、種族値テーブルの更新だと思います。
Niantic社の変更を随時、追跡しなければならないことです。

実行サンプルはこちら

Chrome, FF で動作確認しています。
以下はソースコードです。index.htmlの１本です。


<!DOCTYPE html>
<html>
<head>
<meta charset="utf-8" />
<meta name="viewport" content="width=device-width, initial-scale=1.0, user-scalable=no" />
<title>poke5 個体値計算ツール</title>
<style>
* {
	margin:0;
	padding: 0;
}
body {
	padding:5px;
	background:#ace url(kinkoi.jpg);
}
#namediv input {
	width:70px;
	height:25px;
	font-size:12px;
}
input.pk_inp {
	display:block;
	margin:2px 10px 2px 0;
	padding: 2px;
}
form.hyouka input {
	margin:2px 2px 2px 20px;
	padding: 2px;
}
span.inp_lavel {
	display:block;
	width:60px;
	float:left;
	clear:both;
}
select {
	margin:5px;
	padding: 2px;
}
table {
	margin:5px;
	padding: 2px;
}
td {
	width:100px;
	text-align:center;
}
</style>

<script type="text/javascript">

var mltbl = [
// ポケモンレベルに関するテーブル
//  ml:    ポケモンレベル
//  perc:  レベル40に対する割合
//  c:     CP補正値
//  sand:  強化に必要な砂
//  candy: 強化に必要な飴
{ml:1,perc:1.414724036,c:0.094,sand:0,candy:0},
{ml:1.5,perc:2.923932542,c:0.1351374,sand:200,candy:1},
{ml:2,perc:4.433141061,c:0.1663979,sand:200,candy:1},
{ml:2.5,perc:5.942349554,c:0.1926509,sand:200,candy:1},
{ml:3,perc:7.451557498,c:0.2157325,sand:200,candy:1},
{ml:3.5,perc:8.960765971,c:0.2365727,sand:400,candy:1},
{ml:4,perc:10.4699736,c:0.2557201,sand:400,candy:1},
{ml:4.5,perc:11.97918208,c:0.2735304,sand:400,candy:1},
{ml:5,perc:13.48839026,c:0.2902499,sand:400,candy:1},
{ml:5.5,perc:14.99759871,c:0.3060574,sand:600,candy:1},
{ml:6,perc:16.50680775,c:0.3210876,sand:600,candy:1},
{ml:6.5,perc:18.01601622,c:0.335445,sand:600,candy:1},
{ml:7,perc:19.52522443,c:0.3492127,sand:600,candy:1},
{ml:7.5,perc:21.03443291,c:0.3624578,sand:800,candy:1},
{ml:8,perc:22.54364161,c:0.3752356,sand:800,candy:1},
{ml:8.5,perc:24.05285006,c:0.3875924,sand:800,candy:1},
{ml:9,perc:25.56206057,c:0.3995673,sand:800,candy:1},
{ml:9.5,perc:27.07126902,c:0.4111936,sand:1000,candy:1},
{ml:10,perc:28.58047697,c:0.4225,sand:1000,candy:1},
{ml:10.5,perc:30.00849419,c:0.4335117,sand:1000,candy:1},
{ml:11,perc:31.43650994,c:0.4431076,sand:1000,candy:1},
{ml:11.5,perc:32.86452696,c:0.45306,sand:1300,candy:2},
{ml:12,perc:34.29254396,c:0.4627984,sand:1300,candy:2},
{ml:12.5,perc:35.72056114,c:0.4723361,sand:1300,candy:2},
{ml:13,perc:37.14857675,c:0.481685,sand:1300,candy:2},
{ml:13.5,perc:38.5765939,c:0.4908558,sand:1600,candy:2},
{ml:14,perc:40.00461398,c:0.4998584,sand:1600,candy:2},
{ml:14.5,perc:41.43263118,c:0.5087018,sand:1600,candy:2},
{ml:15,perc:42.86064648,c:0.517394,sand:1600,candy:2},
{ml:15.5,perc:44.28866364,c:0.5259425,sand:1900,candy:2},
{ml:16,perc:45.71668075,c:0.5343543,sand:1900,candy:2},
{ml:16.5,perc:47.14469795,c:0.5426358,sand:1900,candy:2},
{ml:17,perc:48.57271292,c:0.5507927,sand:1900,candy:2},
{ml:17.5,perc:50.00073006,c:0.5588306,sand:2200,candy:2},
{ml:18,perc:51.42875488,c:0.5667545,sand:2200,candy:2},
{ml:18.5,perc:52.85677208,c:0.5745692,sand:2200,candy:2},
{ml:19,perc:54.28478797,c:0.5822789,sand:2200,candy:2},
{ml:19.5,perc:55.71280516,c:0.5898879,sand:2500,candy:2},
{ml:20,perc:57.14082102,c:0.5974,sand:2500,candy:2},
{ml:20.5,perc:58.56883825,c:0.6048188,sand:2500,candy:2},
{ml:21,perc:59.99873261,c:0.6121573,sand:2500,candy:2},
{ml:21.5,perc:61.42674981,c:0.6194041,sand:3000,candy:3},
{ml:22,perc:62.8566514,c:0.6265671,sand:3000,candy:3},
{ml:22.5,perc:64.28466862,c:0.6336492,sand:3000,candy:3},
{ml:23,perc:65.7145723,c:0.640653,sand:3000,candy:3},
{ml:23.5,perc:67.14258943,c:0.647581,sand:3500,candy:3},
{ml:24,perc:68.57248641,c:0.6544356,sand:3500,candy:3},
{ml:24.5,perc:70.00050364,c:0.6612193,sand:3500,candy:3},
{ml:25,perc:71.43040741,c:0.667934,sand:3500,candy:3},
{ml:25.5,perc:72.85842452,c:0.6745819,sand:4000,candy:4},
{ml:26,perc:74.28832607,c:0.6811649,sand:4000,candy:4},
{ml:26.5,perc:75.71634322,c:0.6876849,sand:4000,candy:4},
{ml:27,perc:77.14623491,c:0.6941437,sand:4000,candy:4},
{ml:27.5,perc:78.57425209,c:0.7005429,sand:4500,candy:4},
{ml:28,perc:80.00416311,c:0.7068842,sand:4500,candy:4},
{ml:28.5,perc:81.43218023,c:0.7131691,sand:4500,candy:4},
{ml:29,perc:82.86207287,c:0.7193991,sand:4500,candy:4},
{ml:29.5,perc:84.29009019,c:0.7255756,sand:5000,candy:4},
{ml:30,perc:85.71999464,c:0.7317,sand:5000,candy:4},
{ml:30.5,perc:86.43399483,c:0.734741,sand:5000,candy:4},
{ml:31,perc:87.14799422,c:0.7377695,sand:5000,candy:4},
{ml:31.5,perc:87.86199453,c:0.7407856,sand:6000,candy:6},
{ml:32,perc:88.57599299,c:0.7437894,sand:6000,candy:6},
{ml:32.5,perc:89.28999328,c:0.7467812,sand:6000,candy:6},
{ml:33,perc:90.00399022,c:0.749761,sand:6000,candy:6},
{ml:33.5,perc:90.71799057,c:0.7527291,sand:7000,candy:8},
{ml:34,perc:91.43199898,c:0.7556855,sand:7000,candy:8},
{ml:34.5,perc:92.14599918,c:0.7586304,sand:7000,candy:8},
{ml:35,perc:92.85999484,c:0.7615638,sand:7000,candy:8},
{ml:35.5,perc:93.57399517,c:0.7644861,sand:8000,candy:10},
{ml:36,perc:94.28799738,c:0.7673972,sand:8000,candy:10},
{ml:36.5,perc:95.00199776,c:0.7702973,sand:8000,candy:10},
{ml:37,perc:95.71600203,c:0.7731865,sand:8000,candy:10},
{ml:37.5,perc:96.43000241,c:0.776065,sand:9000,candy:12},
{ml:38,perc:97.14399355,c:0.7789328,sand:9000,candy:12},
{ml:38.5,perc:97.85799386,c:0.7817901,sand:9000,candy:12},
{ml:39,perc:98.57199794,c:0.784637,sand:9000,candy:12},
{ml:39.5,perc:99.2859983,c:0.7874736,sand:10000,candy:15},
{ml:40,perc:100,c:0.7903,sand:10000,candy:15}
];

var pktbl = [
// 種族値テーブル
//  no:    ポケモン番号
//  name:  名前
//  maxcp: 最大ＣＰ(Lv=40)
//  atk:   種族値 攻撃
//  def:   種族値 防御
//  hp:    種族値 体力
{no:1,name:"フシギダネ",maxcp:981,atk:118,def:118,hp:90},
{no:2,name:"フシギソウ",maxcp:1552,atk:151,def:151,hp:120},
{no:3,name:"フシギバナ",maxcp:2568,atk:198,def:198,hp:160},
{no:4,name:"ヒトカゲ",maxcp:831,atk:116,def:96,hp:78},
{no:5,name:"リザード",maxcp:1484,atk:158,def:129,hp:116},
{no:6,name:"リザードン",maxcp:2686,atk:223,def:176,hp:156},
{no:7,name:"ゼニガメ",maxcp:808,atk:94,def:122,hp:88},
{no:8,name:"カメール",maxcp:1324,atk:126,def:155,hp:118},
{no:9,name:"カメックス",maxcp:2291,atk:171,def:210,hp:158},
{no:10,name:"キャタピー",maxcp:393,atk:55,def:62,hp:90},
{no:11,name:"トランセル",maxcp:419,atk:45,def:94,hp:100},
{no:12,name:"バタフリー",maxcp:1701,atk:167,def:151,hp:120},
{no:13,name:"ビードル",maxcp:397,atk:63,def:55,hp:80},
{no:14,name:"コクーン",maxcp:392,atk:46,def:86,hp:90},
{no:15,name:"スピアー",maxcp:1777,atk:169,def:150,hp:130},
{no:16,name:"ポッポ",maxcp:580,atk:85,def:76,hp:80},
{no:17,name:"ピジョン",maxcp:1085,atk:117,def:108,hp:126},
{no:18,name:"ピジョット",maxcp:1994,atk:166,def:157,hp:166},
{no:19,name:"コラッタ",maxcp:588,atk:103,def:70,hp:60},
{no:20,name:"ラッタ",maxcp:1549,atk:161,def:144,hp:110},
{no:21,name:"オニスズメ",maxcp:673,atk:112,def:61,hp:80},
{no:22,name:"オニドリル",maxcp:1814,atk:182,def:135,hp:130},
{no:23,name:"アーボ",maxcp:778,atk:110,def:102,hp:70},
{no:24,name:"アーボック",maxcp:1737,atk:167,def:158,hp:120},
{no:25,name:"ピカチュウ",maxcp:787,atk:112,def:101,hp:70},
{no:26,name:"ライチュウ",maxcp:2025,atk:193,def:165,hp:120},
{no:27,name:"サンド",maxcp:1194,atk:126,def:145,hp:100},
{no:28,name:"サンドパン",maxcp:2328,atk:182,def:202,hp:150},
{no:29,name:"ニドラン♀",maxcp:736,atk:86,def:94,hp:110},
{no:30,name:"ニドリーナ",maxcp:1218,atk:117,def:126,hp:140},
{no:31,name:"ニドクイン",maxcp:2338,atk:180,def:174,hp:180},
{no:32,name:"ニドラン♂",maxcp:739,atk:105,def:76,hp:92},
{no:33,name:"ニドリーノ",maxcp:1252,atk:137,def:112,hp:122},
{no:34,name:"ニドキング",maxcp:2386,atk:204,def:157,hp:162},
{no:35,name:"ピッピ",maxcp:1085,atk:107,def:116,hp:140},
{no:36,name:"ピクシー",maxcp:2353,atk:178,def:171,hp:190},
{no:37,name:"ロコン",maxcp:774,atk:96,def:122,hp:76},
{no:38,name:"キュウコン",maxcp:2157,atk:169,def:204,hp:146},
{no:39,name:"プリン",maxcp:713,atk:80,def:44,hp:230},
{no:40,name:"プクリン",maxcp:1906,atk:156,def:93,hp:280},
{no:41,name:"ズバット",maxcp:569,atk:83,def:76,hp:80},
{no:42,name:"ゴルバット",maxcp:1830,atk:161,def:153,hp:150},
{no:43,name:"ナゾノクサ",maxcp:1069,atk:131,def:116,hp:90},
{no:44,name:"クサイハナ",maxcp:1512,atk:153,def:139,hp:120},
{no:45,name:"ラフレシア",maxcp:2367,atk:202,def:170,hp:150},
{no:46,name:"パラス",maxcp:836,atk:121,def:99,hp:70},
{no:47,name:"パラセクト",maxcp:1657,atk:165,def:146,hp:120},
{no:48,name:"コンパン",maxcp:902,atk:100,def:102,hp:120},
{no:49,name:"モルフォン",maxcp:1937,atk:179,def:150,hp:140},
{no:50,name:"ディグダ",maxcp:465,atk:109,def:88,hp:20},
{no:51,name:"ダグトリオ",maxcp:1333,atk:167,def:147,hp:70},
{no:52,name:"ニャース",maxcp:638,atk:92,def:81,hp:80},
{no:53,name:"ペルシアン",maxcp:1539,atk:150,def:139,hp:130},
{no:54,name:"コダック",maxcp:966,atk:122,def:96,hp:100},
{no:55,name:"ゴルダック",maxcp:2270,atk:191,def:163,hp:160},
{no:56,name:"マンキー",maxcp:1002,atk:148,def:87,hp:80},
{no:57,name:"オコリザル",maxcp:2105,atk:207,def:144,hp:130},
{no:58,name:"ガーディ",maxcp:1110,atk:136,def:96,hp:110},
{no:59,name:"ウインディ",maxcp:2839,atk:227,def:166,hp:180},
{no:60,name:"ニョロモ",maxcp:695,atk:101,def:82,hp:80},
{no:61,name:"ニョロゾ",maxcp:1313,atk:130,def:130,hp:130},
{no:62,name:"ニョロボン",maxcp:2441,atk:182,def:187,hp:180},
{no:63,name:"ケーシィ",maxcp:1148,atk:195,def:103,hp:50},
{no:64,name:"ユンゲラー",maxcp:1859,atk:232,def:138,hp:80},
{no:65,name:"フーディン",maxcp:2887,atk:271,def:194,hp:110},
{no:66,name:"ワンリキー",maxcp:1199,atk:137,def:88,hp:140},
{no:67,name:"ゴーリキー",maxcp:1910,atk:177,def:130,hp:160},
{no:68,name:"カイリキー",maxcp:2889,atk:234,def:162,hp:180},
{no:69,name:"マダツボミ",maxcp:916,atk:139,def:64,hp:100},
{no:70,name:"ウツドン",maxcp:1475,atk:172,def:95,hp:130},
{no:71,name:"ウツボット",maxcp:2268,atk:207,def:138,hp:160},
{no:72,name:"メノクラゲ",maxcp:956,atk:97,def:182,hp:80},
{no:73,name:"ドククラゲ",maxcp:2374,atk:166,def:237,hp:160},
{no:74,name:"イシツブテ",maxcp:1193,atk:132,def:163,hp:80},
{no:75,name:"ゴローン",maxcp:1815,atk:164,def:196,hp:110},
{no:76,name:"ゴローニャ",maxcp:2916,atk:211,def:229,hp:160},
{no:77,name:"ポニータ",maxcp:1502,atk:170,def:132,hp:100},
{no:78,name:"ギャロップ",maxcp:2252,atk:207,def:167,hp:130},
{no:79,name:"ヤドン",maxcp:1204,atk:109,def:109,hp:180},
{no:80,name:"ヤドラン",maxcp:2482,atk:177,def:194,hp:190},
{no:81,name:"コイル",maxcp:1083,atk:165,def:128,hp:50},
{no:82,name:"レアコイル",maxcp:2237,atk:223,def:182,hp:100},
{no:83,name:"カモネギ",maxcp:1092,atk:124,def:118,hp:104},
{no:84,name:"ドードー",maxcp:1011,atk:158,def:88,hp:70},
{no:85,name:"ドードリオ",maxcp:2138,atk:218,def:145,hp:120},
{no:86,name:"パウワウ",maxcp:899,atk:85,def:128,hp:130},
{no:87,name:"ジュゴン",maxcp:1894,atk:139,def:184,hp:180},
{no:88,name:"ベトベター",maxcp:1269,atk:135,def:90,hp:160},
{no:89,name:"ベトベトン",maxcp:2709,atk:190,def:184,hp:210},
{no:90,name:"シェルダー",maxcp:958,atk:116,def:168,hp:60},
{no:91,name:"パルシェン",maxcp:2475,atk:186,def:323,hp:100},
{no:92,name:"ゴース",maxcp:1002,atk:186,def:70,hp:60},
{no:93,name:"ゴースト",maxcp:1716,atk:223,def:112,hp:90},
{no:94,name:"ゲンガー",maxcp:2619,atk:261,def:156,hp:120},
{no:95,name:"イワーク",maxcp:1002,atk:85,def:288,hp:70},
{no:96,name:"スリープ",maxcp:992,atk:89,def:158,hp:120},
{no:97,name:"スリーパー",maxcp:2048,atk:144,def:215,hp:170},
{no:98,name:"クラブ",maxcp:1386,atk:181,def:156,hp:60},
{no:99,name:"キングラー",maxcp:2694,atk:240,def:214,hp:110},
{no:100,name:"ビリリダマ",maxcp:857,atk:109,def:114,hp:80},
{no:101,name:"マルマイン",maxcp:1900,atk:173,def:179,hp:120},
{no:102,name:"タマタマ",maxcp:1102,atk:107,def:140,hp:120},
{no:103,name:"ナッシー",maxcp:2916,atk:233,def:158,hp:190},
{no:104,name:"カラカラ",maxcp:943,atk:90,def:165,hp:100},
{no:105,name:"ガラガラ",maxcp:1691,atk:144,def:200,hp:120},
{no:106,name:"サワムラー",maxcp:2406,atk:224,def:211,hp:100},
{no:107,name:"エビワラー",maxcp:2098,atk:193,def:212,hp:100},
{no:108,name:"ベロリンガ",maxcp:1322,atk:108,def:137,hp:180},
{no:109,name:"ドガース",maxcp:1091,atk:119,def:164,hp:80},
{no:110,name:"マタドガス",maxcp:2183,atk:174,def:221,hp:130},
{no:111,name:"サイホーン",maxcp:1679,atk:140,def:157,hp:160},
{no:112,name:"サイドン",maxcp:3300,atk:222,def:206,hp:210},
{no:113,name:"ラッキー",maxcp:1469,atk:60,def:176,hp:500},
{no:114,name:"モンジャラ",maxcp:2208,atk:183,def:205,hp:130},
{no:115,name:"ガルーラ",maxcp:2463,atk:181,def:165,hp:210},
{no:116,name:"タッツー",maxcp:921,atk:129,def:125,hp:60},
{no:117,name:"シードラ",maxcp:1979,atk:187,def:182,hp:110},
{no:118,name:"トサキント",maxcp:1006,atk:123,def:115,hp:90},
{no:119,name:"アズマオウ",maxcp:2040,atk:175,def:154,hp:160},
{no:120,name:"ヒトデマン",maxcp:926,atk:137,def:112,hp:60},
{no:121,name:"スターミー",maxcp:2303,atk:210,def:184,hp:120},
{no:122,name:"バリヤード",maxcp:1984,atk:192,def:233,hp:80},
{no:123,name:"ストライク",maxcp:2464,atk:218,def:170,hp:140},
{no:124,name:"ルージュラ",maxcp:2512,atk:223,def:182,hp:130},
{no:125,name:"エレブー",maxcp:2196,atk:198,def:173,hp:130},
{no:126,name:"ブーバー",maxcp:2254,atk:206,def:169,hp:130},
{no:127,name:"カイロス",maxcp:2770,atk:238,def:197,hp:130},
{no:128,name:"ケンタロス",maxcp:2488,atk:198,def:197,hp:150},
{no:129,name:"コイキング",maxcp:220,atk:29,def:102,hp:40},
{no:130,name:"ギャラドス",maxcp:3281,atk:237,def:197,hp:190},
{no:131,name:"ラプラス",maxcp:2603,atk:165,def:180,hp:260},
{no:132,name:"メタモン",maxcp:718,atk:91,def:91,hp:96},
{no:133,name:"イーブイ",maxcp:969,atk:104,def:121,hp:110},
{no:134,name:"シャワーズ",maxcp:3157,atk:205,def:177,hp:260},
{no:135,name:"サンダース",maxcp:2730,atk:232,def:201,hp:130},
{no:136,name:"ブースター",maxcp:2904,atk:246,def:204,hp:130},
{no:137,name:"ポリゴン",maxcp:1567,atk:153,def:139,hp:130},
{no:138,name:"オムナイト",maxcp:1345,atk:155,def:174,hp:70},
{no:139,name:"オムスター",maxcp:2685,atk:207,def:227,hp:140},
{no:140,name:"カブト",maxcp:1172,atk:148,def:162,hp:60},
{no:141,name:"カブトプス",maxcp:2517,atk:220,def:203,hp:120},
{no:142,name:"プテラ",maxcp:2608,atk:221,def:164,hp:160},
{no:143,name:"カビゴン",maxcp:3355,atk:190,def:190,hp:320},
{no:144,name:"フリーザー",maxcp:2933,atk:192,def:249,hp:180},
{no:145,name:"サンダー",maxcp:3330,atk:253,def:188,hp:180},
{no:146,name:"ファイヤー",maxcp:3272,atk:251,def:184,hp:180},
{no:147,name:"ミニリュウ",maxcp:860,atk:119,def:94,hp:82},
{no:148,name:"ハクリュー",maxcp:1609,atk:163,def:138,hp:122},
{no:149,name:"カイリュー",maxcp:3581,atk:263,def:201,hp:182},
{no:150,name:"ミュウツー",maxcp:4760,atk:330,def:200,hp:212},
{no:151,name:"ミュウ",maxcp:3090,atk:210,def:209,hp:200},
{no:152,name:"チコリータ",maxcp:801,atk:92,def:122,hp:90},
{no:153,name:"ベイリーフ",maxcp:1296,atk:122,def:155,hp:120},
{no:154,name:"メガニウム",maxcp:2227,atk:168,def:202,hp:160},
{no:155,name:"ヒノアラシ",maxcp:831,atk:116,def:96,hp:78},
{no:156,name:"マグマラシ",maxcp:1484,atk:158,def:129,hp:116},
{no:157,name:"バクフーン",maxcp:2686,atk:223,def:176,hp:156},
{no:158,name:"ワニノコ",maxcp:1011,atk:117,def:116,hp:100},
{no:159,name:"アリゲイツ",maxcp:1598,atk:150,def:151,hp:130},
{no:160,name:"オーダイル",maxcp:2721,atk:205,def:197,hp:170},
{no:161,name:"オタチ",maxcp:519,atk:79,def:77,hp:70},
{no:162,name:"オオタチ",maxcp:1667,atk:148,def:130,hp:170},
{no:163,name:"ホーホー",maxcp:640,atk:67,def:101,hp:120},
{no:164,name:"ヨルノズク",maxcp:2040,atk:145,def:179,hp:200},
{no:165,name:"レディバ",maxcp:663,atk:72,def:142,hp:80},
{no:166,name:"レディアン",maxcp:1275,atk:107,def:209,hp:110},
{no:167,name:"イトマル",maxcp:685,atk:105,def:73,hp:80},
{no:168,name:"アリアドス",maxcp:1636,atk:161,def:128,hp:140},
{no:169,name:"クロバット",maxcp:2466,atk:194,def:178,hp:170},
{no:170,name:"チョンチー",maxcp:1067,atk:106,def:106,hp:150},
{no:171,name:"ランターン",maxcp:2077,atk:146,def:146,hp:250},
{no:172,name:"ピチュー",maxcp:376,atk:77,def:63,hp:40},
{no:173,name:"ピィ",maxcp:620,atk:75,def:91,hp:100},
{no:174,name:"ププリン",maxcp:512,atk:69,def:34,hp:180},
{no:175,name:"トゲピー",maxcp:540,atk:67,def:116,hp:70},
{no:176,name:"トゲチック",maxcp:1543,atk:139,def:191,hp:110},
{no:177,name:"ネイティ",maxcp:925,atk:134,def:89,hp:80},
{no:178,name:"ネイティオ",maxcp:1975,atk:192,def:146,hp:130},
{no:179,name:"メリープ",maxcp:887,atk:114,def:82,hp:110},
{no:180,name:"モココ",maxcp:1402,atk:145,def:112,hp:140},
{no:181,name:"デンリュウ",maxcp:2695,atk:211,def:172,hp:180},
{no:182,name:"キレイハナ",maxcp:2108,atk:169,def:189,hp:150},
{no:183,name:"マリル",maxcp:420,atk:37,def:93,hp:140},
{no:184,name:"マリルリ",maxcp:1503,atk:112,def:152,hp:200},
{no:185,name:"ウソッキー",maxcp:2065,atk:167,def:198,hp:140},
{no:186,name:"ニョロトノ",maxcp:2371,atk:174,def:192,hp:180},
{no:187,name:"ハネッコ",maxcp:508,atk:67,def:101,hp:70},
{no:188,name:"ポポッコ",maxcp:882,atk:91,def:127,hp:110},
{no:189,name:"ワタッコ",maxcp:1553,atk:118,def:197,hp:150},
{no:190,name:"エイパム",maxcp:1188,atk:136,def:112,hp:110},
{no:191,name:"ヒマナッツ",maxcp:316,atk:55,def:55,hp:60},
{no:192,name:"キマワリ",maxcp:2048,atk:185,def:148,hp:150},
{no:193,name:"ヤンヤンマ",maxcp:1326,atk:154,def:94,hp:130},
{no:194,name:"ウパー",maxcp:596,atk:75,def:75,hp:110},
{no:195,name:"ヌオー",maxcp:1929,atk:152,def:152,hp:190},
{no:196,name:"エーフィ",maxcp:3000,atk:261,def:194,hp:130},
{no:197,name:"ブラッキー",maxcp:2052,atk:126,def:250,hp:190},
{no:198,name:"ヤミカラス",maxcp:1392,atk:175,def:87,hp:120},
{no:199,name:"ヤドキング",maxcp:2482,atk:177,def:194,hp:190},
{no:200,name:"ムウマ",maxcp:1781,atk:167,def:167,hp:120},
{no:201,name:"アンノーン",maxcp:1022,atk:136,def:91,hp:96},
{no:202,name:"ソーナンス",maxcp:1024,atk:60,def:106,hp:380},
{no:203,name:"キリンリキ",maxcp:1863,atk:182,def:133,hp:140},
{no:204,name:"クヌギダマ",maxcp:1045,atk:108,def:146,hp:100},
{no:205,name:"フォレトス",maxcp:2263,atk:161,def:242,hp:150},
{no:206,name:"ノコッチ",maxcp:1615,atk:131,def:131,hp:200},
{no:207,name:"グライガー",maxcp:1758,atk:143,def:204,hp:130},
{no:208,name:"ハガネール",maxcp:2439,atk:148,def:333,hp:150},
{no:209,name:"ブルー",maxcp:1124,atk:137,def:89,hp:120},
{no:210,name:"グランブル",maxcp:2440,atk:212,def:137,hp:180},
{no:211,name:"ハリーセン",maxcp:1910,atk:184,def:148,hp:130},
{no:212,name:"ハッサム",maxcp:2801,atk:236,def:191,hp:140},
{no:213,name:"ツボツボ",maxcp:300,atk:17,def:396,hp:40},
{no:214,name:"ヘラクロス",maxcp:2938,atk:234,def:189,hp:160},
{no:215,name:"ニューラ",maxcp:1868,atk:189,def:157,hp:110},
{no:216,name:"ヒメグマ",maxcp:1184,atk:142,def:93,hp:120},
{no:217,name:"リングマ",maxcp:2760,atk:236,def:144,hp:180},
{no:218,name:"マグマッグ",maxcp:750,atk:118,def:71,hp:80},
{no:219,name:"マグカルゴ",maxcp:1543,atk:139,def:209,hp:100},
{no:220,name:"ウリムー",maxcp:663,atk:90,def:74,hp:100},
{no:221,name:"イノムー",maxcp:2284,atk:181,def:147,hp:200},
{no:222,name:"サニーゴ",maxcp:1214,atk:118,def:156,hp:110},
{no:223,name:"テッポウオ",maxcp:749,atk:127,def:69,hp:70},
{no:224,name:"オクタン",maxcp:2124,atk:197,def:141,hp:150},
{no:225,name:"デリバード",maxcp:937,atk:128,def:90,hp:90},
{no:226,name:"マンタイン",maxcp:2032,atk:148,def:260,hp:130},
{no:227,name:"エアームド",maxcp:2032,atk:148,def:260,hp:130},
{no:228,name:"デルビル",maxcp:1110,atk:152,def:93,hp:90},
{no:229,name:"ヘルガー",maxcp:2529,atk:224,def:159,hp:150},
{no:230,name:"キングドラ",maxcp:2424,atk:194,def:194,hp:150},
{no:231,name:"ゴマゾウ",maxcp:1175,atk:107,def:107,hp:180},
{no:232,name:"ドンファン",maxcp:3022,atk:214,def:214,hp:180},
{no:233,name:"ポリゴン2",maxcp:2546,atk:198,def:183,hp:170},
{no:234,name:"オドシシ",maxcp:1988,atk:192,def:132,hp:146},
{no:235,name:"ドーブル",maxcp:389,atk:40,def:88,hp:110},
{no:236,name:"バルキー",maxcp:404,atk:64,def:64,hp:70},
{no:237,name:"カポエラー",maxcp:1905,atk:173,def:214,hp:100},
{no:238,name:"ムチュール",maxcp:1230,atk:153,def:116,hp:90},
{no:239,name:"エレキッド",maxcp:1073,atk:135,def:110,hp:90},
{no:240,name:"ブビィ",maxcp:1178,atk:151,def:108,hp:90},
{no:241,name:"ミルタンク",maxcp:2312,atk:157,def:211,hp:190},
{no:242,name:"ハピナス",maxcp:3219,atk:129,def:229,hp:510},
{no:243,name:"ライコウ",maxcp:3349,atk:241,def:210,hp:180},
{no:244,name:"エンテイ",maxcp:3377,atk:235,def:176,hp:230},
{no:245,name:"スイクン",maxcp:2823,atk:180,def:235,hp:200},
{no:246,name:"ヨーギラス",maxcp:904,atk:115,def:93,hp:100},
{no:247,name:"サナギラス",maxcp:1608,atk:155,def:133,hp:140},
{no:248,name:"バンギラス",maxcp:3670,atk:251,def:212,hp:200},
{no:249,name:"ルギア",maxcp:3598,atk:193,def:323,hp:212},
{no:250,name:"ホウオウ",maxcp:4650,atk:263,def:301,hp:212},
{no:251,name:"セレビィ",maxcp:3090,atk:210,def:210,hp:200},
];

// ポケモン名ボタンを、番号順または名前順に並べて表示
var changetbl = function(mode){
	var len = pktbl.length;
	var div = document.getElementById("namediv");
	while (div.firstChild){
		div.removeChild(div.firstChild);
	}

	if (mode == 0){
		pktbl.sort(function(a,b){
			if (a.no > b.no) return 1;
			if (a.no < b.no) return -1;
			return 0;
		});
	}else if (mode == 1){
		pktbl.sort(function(a,b){
			if (a.name > b.name) return 1;
			if (a.name < b.name) return -1;
			return 0;
		});
	}

	for (var i=0; i<len; i++){
		var inp = document.createElement("input");
		inp.type = "button";
		inp.value = pktbl[i].name;
		if (pktbl[i].no < 152){
			inp.style.backgroundColor = "#8ca";
		}else{
			inp.style.backgroundColor = "#c8a";
		}
		inp.onclick = function(e){setname(this.value);};
		div.appendChild(inp);
	}
};

// 入力内容のクリア
var clearinp = function(mode){
	var elm = null;
	if (mode){
		elm = document.getElementById("pkname");
		elm.value = "";
		elm = document.getElementById("pk_cp");
		elm.value = "";
		elm = document.getElementById("pk_hp");
		elm.value = "";
		elm = document.getElementById("selsand");
		elm.selectedIndex = 0;
	}

	elm = document.getElementById("pk_choice1_1");
	elm.checked = false;
	elm = document.getElementById("pk_choice1_2");
	elm.checked = false;
	elm = document.getElementById("pk_choice1_3");
	elm.checked = false;
	elm = document.getElementById("pk_choice1_4");
	elm.checked = false;
	elm = document.getElementById("pk_choice2_1");
	elm.checked = false;
	elm = document.getElementById("pk_choice2_2");
	elm.checked = false;
	elm = document.getElementById("pk_choice2_3");
	elm.checked = false;
	elm = document.getElementById("pk_choice2_4");
	elm.checked = false;
	elm = document.getElementById("pk_choice3_1");
	elm.checked = false;
	elm = document.getElementById("pk_choice3_2");
	elm.checked = false;
	elm = document.getElementById("pk_choice3_3");
	elm.checked = false;

	elm = document.getElementById("answer");
	while (elm.firstChild){
		elm.removeChild(elm.firstChild);
	}
};

// 初期処理
var init = function(){
	changetbl(0);
};

// ポケモンボタンが押されたとき呼ばれる
var setname = function(name){
	var inp = document.getElementById("pkname");
	inp.value = name;
	inp = document.getElementById("pk_cp");
	inp.focus();
};

// 入力チェック
var precheck = function(cp, hp, sand){
	var str = "該当なし";
	var n = pktbl.length;
	var inp = document.getElementById("pkname");
	var pid = -1;
	for (var i=0; i<n; i++){
		if (inp.value == pktbl[i].name){
			pid = i;
			break;
		}
	}

	if (pid < 0){
		str = "ポケモンを選択して下さい。"
	}else{
		if (!sand){
			str = "強化に必要な星の砂を選択して下さい。"
			pid = -1;
		}
		if (!hp){
			str = "HPを半角数字で入力して下さい。"
			pid = -1;
		}
		if (!cp){
			str = "CPを半角数字で入力して下さい。"
			pid = -1;
		}
	}

	return {pid:pid, str:str};
};

// CP計算
var cpcalc = function(p_atk, p_def, p_hp, m_atk, m_def, m_hp, m_c){
	return parseInt((p_atk+m_atk)*Math.sqrt(p_def+m_def)*Math.sqrt(p_hp+m_hp)*m_c*m_c/10.0);
};

// 計算処理メイン
var calcmain = function(pid, cp, hp, sand, k1, k2, k2f, ares){
	var atk_range = {min:0, max:15}; // 攻撃レンジ 枝刈
	var def_range = {min:0, max:15}; // 防御レンジ 枝刈
	var hp_range = {min:0, max:15};  // 体力レンジ 枝刈

	var kotai1 = [ // チームリーダー第一評価
		{min:37, max:45},
		{min:30, max:36},
		{min:23, max:29},
		{min:0, max:22}
	];
	var kotai2 = [ // チームリーダー第二評価
		{min:15, max:15},
		{min:13, max:14},
		{min:8, max:12},
		{min:0, max:7}
	];

	var m = mltbl.length;
	var ahp = [];
	var batk = parseInt(pktbl[pid].atk); // 種族値 攻撃
	var bdef = parseInt(pktbl[pid].def); // 種族値 防御
	var bhp = parseInt(pktbl[pid].hp);   // 種族値 体力

	if (k2){
		if (k2f & 1){ // 第二評価 攻撃ON
			atk_range = kotai2[k2-1];
			def_range.max = kotai2[k2-1].max;
			hp_range.max = kotai2[k2-1].max;
		}
		if (k2f & 2){ // 第二評価 防御ON
			def_range = kotai2[k2-1];
			atk_range.max = kotai2[k2-1].max;
			hp_range.max = kotai2[k2-1].max;
		}
		if (k2f & 4){ // 第二評価 体力ON
			hp_range = kotai2[k2-1];
			atk_range.max = kotai2[k2-1].max;
			def_range.max = kotai2[k2-1].max;
		}
	}

	// 与えられたパラメータに合致するＨＰとポケモンレベルのペア一覧を取得 枝刈
	for (var i=hp_range.min; i<=hp_range.max; i++){
		for (var j=0; j<m-1; j++){
			if (mltbl[j+1].sand == sand){
				var thp = parseInt((bhp+i) * mltbl[j].c);
				if (thp == hp){
					ahp.push({hp:i, mli:j});
				}
			}
		}
	}

	for (var i=0; i<ahp.length; i++){ // ＨＰループ
		var mlt = mltbl[ahp[i].mli];
		for (var j=atk_range.min; j<=atk_range.max; j++){ // 攻撃ループ
			for (var k=def_range.min; k<=def_range.max; k++){ // 防御ループ
				// CP計算
				var t = cpcalc(j,k,ahp[i].hp, batk,bdef,bhp, mlt.c);
				if (t == cp){ // CP一致
					var h = ahp[i].hp; // j:攻撃, k:防御, h:体力
					var tot = j+k+h;
					if (k1){
						// 合計が第一評価を満たすか
						if (tot < kotai1[k1-1].min || kotai1[k1-1].max < tot){
							continue;
						}
					}
					if (k2){
						// 攻撃が第2評価範囲を満たしているか
						if (j < atk_range.min || atk_range.max < j){
							continue;
						}
						// 防御が第2評価範囲を満たしているか
						if (k < def_range.min || def_range.max < k){
							continue;
						}
						// 体力が第2評価範囲を満たしているか
						if (h < hp_range.min || hp_range.max < h){
							continue;
						}
						// 攻撃のみチェックの時、攻撃 > 防御,体力であること
						if (k2f == 1 && (j <= k || j <= h)){
							continue;
						}
						// 防御のみチェックの時、防御 > 攻撃,体力であること
						if (k2f == 2 && (k <= j || k <= h)){
							continue;
						}
						// 体力のみチェックの時、体力 > 攻撃,防御であること
						if (k2f == 4 && (h <= j || h <= k)){
							continue;
						}
						// 攻撃/防御チェックの時、攻撃＝防御かつ、攻撃＞体力であること
						if (k2f == 3 && (j != k || j <= h)){
							continue;
						}
						// 攻撃/体力チェックの時、攻撃＝体力かつ、攻撃＞防御であること
						if (k2f == 5 && (j != h || j <= k)){
							continue;
						}
						// 防御/体力チェックの時、防御＝体力かつ、防御＞攻撃であること
						if (k2f == 6 && (k != h || k <= j)){
							continue;
						}
						// 攻撃/防御/体力チェックの時、攻撃＝防御＝体力であること
						if (k2f == 7 && (j != k || j != h)){
							continue;
						}
					}
					// 配列に追加, {％, レベル, 攻撃, 防御, ＨＰ}
					ares.push({per:tot/45, ml:mlt.ml, atk:j, def:k, hp:h});
				}
			}
		}
	}
};

// 計算ボタンが押されたとき呼ばれる
var calc = function(){
	var cp = parseInt(document.getElementById("pk_cp").value);
	var hp = parseInt(document.getElementById("pk_hp").value);
	var sel = document.getElementById("selsand");
	var sand = parseInt( sel.options[sel.selectedIndex].value);
	var pchk = precheck(cp, hp, sand); // 入力チェック
	var ares = [];
	var answerstr = pchk.str;
	var div = document.getElementById("answer");
	var k1 = parseInt(hyouka.pk_choice1.value);
	var k2 = parseInt(hyouka.pk_choice2.value);
	var k2f = 0;
	if (k2){
		if (hyouka.pk_choice3a.checked){
			k2f |= 1;
		}
		if (hyouka.pk_choice3b.checked){
			k2f |= 2;
		}
		if (hyouka.pk_choice3c.checked){
			k2f |= 4;
		}
	}

	if (pchk.pid >= 0){
		// 計算機メイン ares配列に結果が返る
		calcmain(pchk.pid, cp, hp, sand, k1, k2, k2f, ares);

		if (ares.length){
			// ％を優先して降順ソート
			ares.sort(function(a,b){
				if (a.per < b.per) return 1;
				if (a.per > b.per) return -1;
				if (a.ml < b.ml) return 1;
				if (a.ml > b.ml) return -1;
				return 0;
			});

			// 結果をテーブルの形に整形
			answerstr = "<table border='1'><tr><th>％</th><th>Lv</th><th>攻撃</th><th>防御</th><th>ＨＰ</th></tr>";
			for (var i=0; i<ares.length; i++){
				answerstr += "<tr><td>" + Math.round(ares[i].per * 10000) / 100 + "</td><td>" + ares[i].ml
 + "</td><td>" + ares[i].atk + "</td><td>" + ares[i].def + "</td><td>" + ares[i].hp + "</td></tr>";
			}
			answerstr += "</table>";
		}
	}

	div.innerHTML = answerstr;
};

</script>

</head>

<body onload="init();">
個体値計算ツール　
<input class="pk_choice1" type="radio" name="btn_choice" value="0" checked="checked" onclick="changetbl(0);">番号順　
<input class="pk_choice1" type="radio" name="btn_choice" value="1" onclick="changetbl(1);">名前順
<br />
<div id="namediv"></div>
<br />
<span class="inp_lavel">Name </span><input id="pkname" class="pk_inp" type="text" value="" readonly="true">
<span class="inp_lavel">CP </span><input id="pk_cp" class="pk_inp" type="text" value="">
<span class="inp_lavel">HP </span><input id="pk_hp" class="pk_inp" type="text" value="">
<select id="selsand">
<option value="0">強化に必要な星の砂</option>
<option value="200">200</option>
<option value="400">400</option>
<option value="600">600</option>
<option value="800">800</option>
<option value="1000">1000</option>
<option value="1300">1300</option>
<option value="1600">1600</option>
<option value="1900">1900</option>
<option value="2200">2200</option>
<option value="2500">2500</option>
<option value="3000">3000</option>
<option value="3500">3500</option>
<option value="4000">4000</option>
<option value="4500">4500</option>
<option value="5000">5000</option>
<option value="6000">6000</option>
<option value="7000">7000</option>
<option value="8000">8000</option>
<option value="9000">9000</option>
<option value="10000">10000</option>
</select>
<br />
<br />
<form class="hyouka" name="hyouka">
第一評価<br />
<input id="pk_choice1_1" class="pk_choice1" type="radio" name="pk_choice1" value="1">1st
<input id="pk_choice1_2" class="pk_choice1" type="radio" name="pk_choice1" value="2">2nd
<input id="pk_choice1_3" class="pk_choice1" type="radio" name="pk_choice1" value="3">3rd
<input id="pk_choice1_4" class="pk_choice1" type="radio" name="pk_choice1" value="4">4th
<br />
第二評価<br />
<div>
<input id="pk_choice2_1" class="pk_choice2" type="radio" name="pk_choice2" value="1">1st
<input id="pk_choice2_2" class="pk_choice2" type="radio" name="pk_choice2" value="2">2nd
<input id="pk_choice2_3" class="pk_choice2" type="radio" name="pk_choice2" value="3">3rd
<input id="pk_choice2_4" class="pk_choice2" type="radio" name="pk_choice2" value="4">4th
<br />
<input id="pk_choice3_1" class="pk_choice3" type="checkbox" name="pk_choice3a" value="1">攻撃
<input id="pk_choice3_2" class="pk_choice3" type="checkbox" name="pk_choice3b" value="2">防御
<input id="pk_choice3_3" class="pk_choice3" type="checkbox" name="pk_choice3c" value="3">ＨＰ
</div>
</form>
<br />
<input type="button" value="　計　算　" onclick="calc(); return false;">
<input type="button" value="　評価クリア　" onclick="clearinp(false); return false;">
<input type="button" value="　ALLクリア　" onclick="clearinp(true); return false;">
<br />
<br />
<hr>
<div id="answer">
</div>
</body>
</html>



Tags: プログラムメモ
author : HUNDREDSOFT | - | -

東芝のFlashAirエミュレータ

| プログラムメモ | 2016/02/19 Fri 00:33

FlashAirエミュレータ　flashair-lua-dev で少し困った。


FlashAirのエミュレータを探していて
https://flashair-developers.com/ja/documents/resources/
の下の方に表示される、flashair-lua-dev を試した。

fa オブジェクトがエミュレートできるので大変便利で助かったのですが、
ファイルを multipart/form-data でアップロードが出来なかった。

【手順】
https://flashair-developers.com/ja/documents/api/lua/reference/#request
のファイルアップロードのサンプルコードを使い、flashair.luaを呼び出すが、通信が完了しない。
(CentOS 6.7 x86)

【問題点】
「引数bodyに、<!--WLANSDFILE-->が記述された場合、置き換えで、file内容を埋め込む」
というFlashAir独自仕様が満たされていないように見えます。

【対策案】
少し冗長ですが、flashair.luaのtrequest関数を、次の様に修正する。
	local trequest = function(...)
		local param = ...
		local url      = param["url"]
		local method   = param["method"]
		local headers  = param["headers"]
		local file     = param["file"]
		local reqbody  = param["body"]
		local bufsize  = param["bufsize"]
		local redirect = param["redirect"]

		-- Hundredsoft add start
		if (file ~= nil) then
			local fp = io.open(file, "rb")
			if (fp ~= nil) then
				local r1, r2 = string.match(reqbody, "(.-)%<%!%-%-WLANSDFILE%-%-%>(.+)")
				if (r1 ~= nil and r2 ~= nil) then
					reqbody = r1 .. fp:read("*a") .. r2
				end
				fp:close()
			end
		end
		-- Hundredsoft add end

		local http = require("socket.http")
		local ltn12 = require("ltn12")

		local body = {}
		local b, c, h = http.request {
			url = url,
			sink = ltn12.sink.table(body),
			method = method,
			headers = headers,
			source = ltn12.source.string(reqbody), 
			step = nil,
			proxy = nil, 
			redirect = redirect,
			create = nil,
		}
		return table.concat(body), c, h
	end

で、動いた。

CentOS上では、メモリが豊富なので bodyに、<!--WLANSDFILE--> などと打ち込まず、
ファイル内容を、そのまま入れ込む方法でも動きますが、
実機(FlashAir)では、大きな写真を送ろうとすると、メモリ不足でコケてしまいます。
大きな画像を送るときには、<!--WLANSDFILE--> は、必須のようです。

個人的には、あとは、FTP エミュレータができると、とても嬉しい。



Tags: プログラムメモ
author : HUNDREDSOFT | - | -

DicomViewer HTML5版

| プログラムメモ | 2015/02/22 Sun 03:48

HTMLとjavascriptだけで作成したDICOM Viewerです。


こちらから、利用できます。

WindowsOSだけでなく、Androidタブレット等でも利用可能です。
ブラウザ単独で動作するので、専用のソフトウェアをインストールすることなくご利用頂けます。



どうしても心配される場合は、こちら から本アプリをダウンロードして、
外部とは切断されたネットワーク環境の中でお使い下さい。

ローカル環境で動作させることも可能ですが、IE,Chromeでは、URL指定ができません。
(複数ファイル指定やファイルドロップは可能です。)


【操作ガイド】

まず初めに、URL指定でa.dcmを表示して下さい(Readボタンをクリック)。












画像部分をクリックすると、動画の場合、上図のメッセージが表示されます。
OKを押すと、画像をメモリ上に展開して動画再生を行います。
再生速度は、DICOM情報に従います。
画像サイズや画像枚数により、非常に時間が掛る場合があります。
このようなケースでは、ブラウザから、スクリプトを停止するかどうかの警告メッセージが表示されます。

再生の停止は、再度、画面上をクリックして下さい。
停止後、再クリックすると、再生が始まります。

静止画の場合は、次のメッセージが表示されます。



IISサーバーにインストールする場合で、URL指定でDICOMファイルにアクセスする時には、
サーバー側にMIME設定を行って下さい。



DICOMファイルの拡張子が「.dcm」である場合は、
上図のように「application/x-www-form-urlencoded」を設定して下さい。
DICOMファイルに複数の拡張子がある場合は、それぞれの拡張子毎に設定が必要です。




古いブラウザをお使いの場合、上記の表示になります。
IE10以降・Firefox34以降・Chrome31以降・Safari7.1以降・Opera26以降・
Android ブラウザ4.4.3以降でお使い下さい。
Windows版Safari(5.1.7)では動作しません。

対応DICOMファイル

画像対応
Explicit VR Little Endian
Implicit VR Little Endian
Jpeg Lossy
RLE Lossless

画像未対応のフォーマットであっても、タグ等のDICOMデータは、表示します。




認識できないファイルについては、上図のように「解析不能でした。」と表示されます。


こちらから、利用できます。



Tags: プログラムメモ
author : HUNDREDSOFT | - | -

ModifiedSimpson

| プログラムメモ | 2015/02/08 Sun 00:18
Echoを使った左室容積の計算方法には幾つかあるのですが、４腔断面・２腔断面の両方が得られる場合に、良く使われるのが Modified Simpson法です。

日本医師会雑誌の「心エコーのABC」から計算部分を抜粋。


一見、簡単なプログラムで構築できそうですが、交点算出が割と面倒です。
ここでは、javascriptでの実装を目指します。



トレースラインを Q0, Q1, Q2, ...で表し、その長径をMPで表す。
MPをm分割し、分割点をP0, P1, P2, ... とする。



等間隔で分割することにして、k番目の分割線の垂線が、トレース座標のQi, Qi+1線分を横切ることを考える。
直線MPと 直線(Pk)-(Qi)のベクトル内積は負になるが、MPと (Pk)-(Qi+1)のベクトル内積は正になる。
内積符号が異なった場合に、各トレース座標間の直線(Q)と分割点(Pk)からのMPの直交ベクトルとの交点が、候補になる。
トレースラインが捻じれる場合もあるので、分割点(Pk)から最も遠い候補を解とする。
また、左右それぞれに解を求める必要があるので、MPと分割点(Pk)から解までのベクトルとで、
外積を取り、左右どちらの解であるかを確認する。

複雑そうに思えるかもしれないが、ループ回数は、トレースポイント数 x 分割数 でしかなく、
四則演算だけでコードが書ける。
この解説部分に該当する関数名は、calc() で、 85行しかない。




デモプログラムでは、実運用に近づけるため、B-DUAL画像を模して左右に、４腔断面と２腔断面を配置し、
下段に、Modified Simponによる体積を表示する。
左右画像上に表示されるSimpson値は、単画像による等分割で、回転体により体積を求めた場合の数値で、
最初に載せた文献によれば、これを計測値とすることは少ないようです。

ここでは、アルゴリズムの紹介であるので、[pix]単位の計算結果を表示している。
生体長換算値 [cm/pix]があれば、生体長換算値を３回掛ければ [mL] での結果が求まる。


Canvasは、左右と下段の3枚を使っている。
スライス算出部は、２つnewして、左右それぞれに振り、
結果が求まった時点で、callbackにより、下段のスクリプトを呼び出す仕組みにした。

実行サンプルはこちら

Chrome, FF, IE11 で動作確認しています。

以下はソースコードです。
html と js の各１本です。

[html]

<!DOCTYPE html>
<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<meta http-equiv="X-UA-Compatible" content="IE=edge" />
<meta name="viewport" content="width=device-width, initial-scale=1.0, user-scalable=no" />
<title>Modified Simpson</title>

<style type="text/css">
#div1 {
	width:800px;
	height:600px;
}
#canvas1 {
	width:400px;
	height:500px;
	background-color:#f0f0e0;
	float:left;
}
#canvas2 {
	width:400px;
	height:500px;
	background-color:#e0f0f0;
	float:left;
}
#canvas3 {
	width:800px;
	height:100px;
	background-color:#f0e0f0;
	float:left;
	clear: both;
}
</style>

<script type="text/javascript" src="DSM.js"></script>

<script type="text/javascript">
var init = function(){
	ModifiedSimpson.init(document.getElementById("canvas3"), 20);
	var dsm1 = new DSMPublic.DSM(document.getElementById("canvas1"), 
		function(hd){
			ModifiedSimpson.paint(1, hd);
		});
	var dsm2 = new DSMPublic.DSM(document.getElementById("canvas2"),
		function(hd){
			ModifiedSimpson.paint(2, hd);
		});
};
</script>
</head>


<body onload="init();">
Modified Simpson's Rule Sample<br /><br />
<div id="div1">
<canvas id="canvas1" width="400" height="500"></canvas>
<canvas id="canvas2" width="400" height="500"></canvas>
<canvas id="canvas3" width="800" height="100"></canvas>
</div>

</body>
</html>

[js]

// (c)Hundredsoft Corporation. 2015 All right reserved.
//	NAME
//		DSM.js
//
//	FUNCTION
//		Modified Simpson's rule.
//
//	REVISION
//		1.00	T.Furumoto	Jan-30-2015	Original
//
//	UTF-8で保存して下さい
//
////////////////////////////////////////////////////////////

new function() {

///////////////////////////////////
 // Private Scope
///////////////////////////////////

var m_linesz = 1;
var m_Slice = 20;
var m_Chain = 2;

var mdown = function(p, mouseX, mouseY){
	if (p.m_mode != 1){
		clear(p);
	}
};

var mup = function(p, mouseX, mouseY){
	var ret;
	if (p.m_mode == 0){
		p.m_trace = traceChain(p.m_trace);
		p.m_P = autoL(p.m_trace);
		ret = calc(p.m_trace, p.m_P);
		if (ret){
			paint(p, ret);
			p.m_mode = 1;
		}else{
			errorDisp(p);
		}
	}else if (p.m_mode == 1){
		p.m_mode = 0;
		ret = calc(p.m_trace, p.m_P);
		if (ret){
			paint(p, ret);
		}else{
			errorDisp(p);
		}
	}
};

var mmove = function(p, mouseX, mouseY){
	if (p.m_mode == 0){
		var n = p.m_trace.x.length;
		p.m_trace.x[n] = mouseX;
		p.m_trace.y[n] = mouseY;
		drawROI(p);

	}else{
		p.m_P = manualL(p.m_trace, mouseX, mouseY);
		var ret = calc(p.m_trace, p.m_P);
		if (ret){
			paint(p, ret);
		}else{
			errorDisp(p);
		}
	}
};


var intersection = function(p1x, p1y, p2x, p2y, q1x, q1y, q2x, q2y){
	var a =  p1y - p2y;
	var b = -p1x + p2x;
	var c =  q1y - q2y;
	var d = -q1x + q2x;
	var det = a*d - b*c;
	if( det * det < 0.01 ) // 平行
		return false;

	var e = a * p1x + b * p1y;
	var f = c * q1x + d * q1y;
	return {"x" : (d*e - b*f) / det, "y" : (-c*e + a*f) / det};
};

var autoL = function(trace){
	var n = trace.x.length;
	var mx = (trace.x[n-1] + trace.x[0]) / 2;
	var my = (trace.y[n-1] + trace.y[0]) / 2;
	var len2 = 0;
	var px = false;
	var py = false;
	for (var i=0; i<n; i++){
		var w = (mx - trace.x[i]) * (mx - trace.x[i])
			+ (my - trace.y[i]) * (my - trace.y[i]);
		if (w > len2){
			len2 = w;
			px = trace.x[i];
			py = trace.y[i];
		}
	}
	return {"x": px, "y": py};
};

var manualL = function(trace, mouseX, mouseY){
	var n = trace.x.length;
	// Pointから最も近い点(p)
	var px = trace.x[0];
	var py = trace.y[0];
	var lmin = (px - mouseX) * (px - mouseX) + (py - mouseY) * (py - mouseY);
	for( var i=1; i<n; i++ ){
		var l = (trace.x[i] - mouseX) * (trace.x[i] - mouseX)
			 + (trace.y[i] - mouseY) * (trace.y[i] - mouseY);
		if( l < lmin ){
			px = trace.x[i];
			py = trace.y[i];
			lmin = l;
		}
	}
	return {"x": px, "y": py};
};

var traceChain = function(o_trace){
	var n = o_trace.x.length;
	var traceX = [];
	var traceY = [];
	traceX.push(o_trace.x[0]);
	traceY.push(o_trace.y[0]);
	for (var i=1; i<n; i++){
		var dx = o_trace.x[i] - o_trace.x[i-1];
		var dy = o_trace.y[i] - o_trace.y[i-1];
		var len = dx * dx + dy * dy;
		if (len > m_Chain * m_Chain * m_Chain * m_Chain){
			var m = ~~(Math.sqrt(len) / m_Chain);
			var ex = dx / m;
			var ey = dy / m;
			var x = o_trace.x[i-1] + ex;
			var y = o_trace.y[i-1] + ey;
			for (var j=1; j<m-1; j++){
				traceX.push(x);
				traceY.push(y);
				x += ex;
				y += ey;
			}
		}
		traceX.push(o_trace.x[i]);
		traceY.push(o_trace.y[i]);
	}
	return {"x": traceX, "y": traceY};
};

var calc = function(trace, p){
	var n = trace.x.length;
	if (n < m_Slice * 2 + 2){
		return false;
	}

	var mx = (trace.x[n-1] + trace.x[0]) / 2;
	var my = (trace.y[n-1] + trace.y[0]) / 2;

	// M-Pベクトル
	var lx = p.x - mx;
	var ly = p.y - my;

	// M-Pに直行するベクトル
	var vx = ly;
	var vy = -lx;

	var out_plx = [];
	var out_ply = [];

	var rlen = [];
	var llen = [];
	for( var i=0; i<m_Slice+1; i++ ){
		rlen[i] = -9999999999;
		llen[i] =  9999999999;
	}

	var qox = mx;
	var qoy = my;
	for (var i=0; i<n+1; i++){
		var qnx = mx;
		var qny = my;
		if( i != n ){
			qnx = trace.x[i];
			qny = trace.y[i];
		}
		for (var j=0; j<m_Slice+1; j++){
			var px = (mx * (m_Slice-j) + p.x * j) / m_Slice;
			var py = (my * (m_Slice-j) + p.y * j) / m_Slice;
			// 線をまたぐか？
			if( ( (qox - px) * lx >= -(qoy - py) * ly
			   && (qnx - px) * lx <  -(qny - py) * ly )
			 || ( (qox - px) * lx <  -(qoy - py) * ly
			   && (qnx - px) * lx >= -(qny - py) * ly ) )
			{
				// 交点座標
				var t = intersection(qox, qoy, qnx, qny, px, py, px+vx, py+vy);
				if (!t){
					var d1x = qnx - px;
					var d1y = qny - py;
					var d2x = qox - px;
					var d2y = qoy - py;
					if( (d1x * d1x + d1y * d1y) > (d2x * d2x + d2y * d2y) ){
						t = {"x" : qnx, "y" : qny};
					}else{
						t = {"x" : qox, "y" : qoy};
					}
				}

				var qtlen = (t.x-px) * (t.x-px) + (t.y-py) * (t.y-py);

				// 方向確認
				if( (t.x - px) * vx + (t.y - py) * vy > 0 ){
					qtlen *= -1.;
				}
				// 最遠方点なら更新
				if( qtlen < llen[j] ){
					out_plx[j*2] = t.x;
					out_ply[j*2] = t.y;
					llen[j] = qtlen;
				}
				if( qtlen > rlen[j] ){
					out_plx[j*2+1] = t.x;
					out_ply[j*2+1] = t.y;
					rlen[j] = qtlen;
				}
			}
		}
		qox = qnx;
		qoy = qny;
	}

	return {"x": out_plx, "y": out_ply};
};

var clear = function(p){
	m_Slice = ModifiedSimpson.getSlice();
	p.m_trace.x = [];
	p.m_trace.y = [];
	p.m_ctx.font = "16px 'MS Gothic'";
	p.m_mode = 0;
	p.m_ctx.clearRect(0, 0, p.m_width, p.m_height);
	p.callback({"H":0, "D":false});
};

var drawROI = function(p){
	var ctx = p.m_ctx;
	ctx.clearRect(0, 0, p.m_width, p.m_height);

	ctx.beginPath();
	ctx.arc(p.m_trace.x[0], p.m_trace.y[0], m_linesz, 0, Math.PI*2, false);
	ctx.fill();
	ctx.closePath();
	var n = p.m_trace.x.length;
	for (var i=1; i<n+1; i++){
		ctx.beginPath();
		ctx.arc(p.m_trace.x[i], p.m_trace.y[i], m_linesz, 0, Math.PI*2, false);
		ctx.fill();
		ctx.closePath();

		ctx.lineWidth = m_linesz * 2;
		ctx.strokeStyle = "#808000";
		ctx.beginPath();
		ctx.moveTo(p.m_trace.x[i-1], p.m_trace.y[i-1]);
		ctx.lineTo(p.m_trace.x[i], p.m_trace.y[i]);
		ctx.stroke();
		ctx.closePath();
	}
};

var paint = function(p, out_pl){
	drawROI(p);

	var ctx = p.m_ctx;
	var n = p.m_trace.x.length;
	var mx = (p.m_trace.x[n-1] + p.m_trace.x[0]) / 2;
	var my = (p.m_trace.y[n-1] + p.m_trace.y[0]) / 2;

	ctx.lineWidth = m_linesz * 2;
	ctx.strokeStyle = "#808000";
	ctx.beginPath();
	ctx.moveTo(p.m_trace.x[0], p.m_trace.y[0]);
	ctx.lineTo(p.m_trace.x[n-1], p.m_trace.y[n-1]);
	ctx.stroke();

	ctx.moveTo(mx, my);
	ctx.lineTo(p.m_P.x, p.m_P.y);
	ctx.stroke();
	ctx.closePath();

	// M-Pベクトル
	var lx = p.m_P.x - mx;
	var ly = p.m_P.y - my;

	var leng = Math.sqrt(lx*lx + ly*ly);
	var area = 0;
	for (var i=0; i<n; i++){
		var j = (i == n-1) ? 0 : i+1;
		area += p.m_trace.x[i] * p.m_trace.y[j] - p.m_trace.x[j] * p.m_trace.y[i];
	}
	area = Math.abs(area);

	var v = 0;
	var h = leng / m_Slice;
	var d2 = [];
	var hpi = Math.PI * h / 4.0;
	for (var i=0; i<m_Slice+1; i++){
		if (i < m_Slice){
			var dx = out_pl.x[i*2+1] - out_pl.x[i*2  ];
			var dy = out_pl.y[i*2+1] - out_pl.y[i*2  ];
			var d = (dx*dx + dy*dy);
			d2.push(d);
			v += hpi * d;
		}
		ctx.beginPath();
		ctx.moveTo(out_pl.x[i*2  ], out_pl.y[i*2  ]);
		ctx.lineTo(out_pl.x[i*2+1], out_pl.y[i*2+1]);
		ctx.stroke();
		ctx.closePath();
	}
	ctx.fillText("Simpson: " + Math.round(v) + " [pix^3]", 20, 20);
	ctx.fillText("Area:    " + Math.round(area) + " [pix^2]", 20, 40);
	ctx.fillText("L:       " + Math.round(leng) + " [pix]", 20, 60);

	p.callback({"H":h, "D":d2});
};

var errorDisp = function(p){
	var ctx = p.m_ctx;
	var n = p.m_trace.x.length;
	ctx.beginPath();
	ctx.moveTo(p.m_trace.x[0], p.m_trace.y[0]);
	ctx.lineTo(p.m_trace.x[n-1], p.m_trace.y[n-1]);
	ctx.stroke();
	ctx.closePath();
	ctx.fillText("トレースラインが短すぎます", 20, 50);
};


///////////////////////////////////
// Public Scope
///////////////////////////////////
DSMPublic = {

DSM: function(canvas, callback){
	this.m_ctx = null;
	this.m_width = 400;
	this.m_height = 400;
	this.m_trace = {"x": [], "y": []};
	this.m_mode = 0;
	this.m_P = {"x": 0, "y": 0};
	this.callback = callback;

	var iam = this;
	var mouseDown = false;
	var mouseX = 0;
	var mouseY = 0;

	if (navigator.userAgent.indexOf('iPhone')>0 ||
		navigator.userAgent.indexOf('iPod')>0 || 
		navigator.userAgent.indexOf('iPad')>0 ||
		navigator.userAgent.indexOf('Android')>0) {

		canvas.addEventListener('touchstart',
			function(e) {
				e.preventDefault();
				var n = e.touches.length;
				if (n > 0) {
					e = e || window.event;
					mouseDown = true;
					var rect=e.target.getBoundingClientRect();
					mouseX = e.touches[n-1].pageX-rect.left;
					mouseY = e.touches[n-1].pageY-rect.top;
					mdown(iam, mouseX, mouseY);
				}
			}, false);

		canvas.addEventListener('touchend',
			function(e) {
				mouseDown = false;
				mup(iam, mouseX, mouseY);
			}, false);

		canvas.addEventListener('touchmove',
			function (e) {
				e.preventDefault();
				var rect=e.target.getBoundingClientRect();
				var n = e.touches.length;
				if (n > 0){
					mouseX = e.touches[n-1].pageX-rect.left;
					mouseY = e.touches[n-1].pageY-rect.top;
					if (mouseDown) {
						mmove(iam, mouseX, mouseY);
					}
				}
			}, false
		);

	}else{
		canvas.addEventListener('mousedown',
			function(e) {
				e.preventDefault();
				mouseDown = true;
				e = e || window.event;
				var rect=e.target.getBoundingClientRect();
				mouseX = e.clientX-rect.left;
				mouseY = e.clientY-rect.top;
				mdown(iam, mouseX, mouseY);
			}, false);


		canvas.addEventListener('mouseup',
			function(e) {
				mouseDown = false;
				mup(iam, mouseX, mouseY);
			}, false);

		canvas.addEventListener('mousemove',
			function (e) {
				e.preventDefault();
				var rect=e.target.getBoundingClientRect();
				mouseX = e.clientX-rect.left;
				mouseY = e.clientY-rect.top;
				if (mouseDown) {
					mmove(iam, mouseX, mouseY);
				}
			}, false
		);
	}

	this.m_ctx = canvas.getContext('2d');
	this.m_width = canvas.width;
	this.m_height = canvas.height;

	clear(this);
}

} // Public Scope
}; // Constructor



///////////////////////////////////
//
// Measurement Display Area
//
///////////////////////////////////

new function() {

///////////////////////////////////
 // Private Scope
///////////////////////////////////
var m_ctx = null;
var m_width = 100;
var m_height = 100;
var m_m1 = {"H": 0, "D": false};
var m_m2 = {"H": 0, "D": false};
var m_Slice = 20;

var report = function(){
	if (!m_ctx){
		return;
	}
	m_ctx.clearRect(0, 0, m_width, m_height);

	if (m_m1.D && m_m2.D){
		var h = (m_m1.H + m_m2.H) / 2;
		var v = 0;
		var hpi = Math.PI * h / 4.0;
		for (var i=0; i<m_Slice; i++){
			v += hpi * Math.sqrt( m_m1.D[i] * m_m2.D[i] );
		}
		m_ctx.fillText("Modified Simpson: " + Math.round(v) + " [pix^3]", 50, 50);
	}
};

///////////////////////////////////
// Public Scope
///////////////////////////////////
ModifiedSimpson = {
init: function(canvas, slc){
	if (slc){
		m_Slice = slc;
	}
	m_ctx = canvas.getContext('2d');
	m_width = canvas.width;
	m_height = canvas.height;
	m_ctx.font = "16px 'MS Gothic'";
	m_ctx.clearRect(0, 0, m_width, m_height);
},

paint: function(id, hd){
	if (id == 1){
		m_m1 = hd;
	}else if (id == 2){
		m_m2 = hd;
	}
	report();
},

getSlice: function(){
	return m_Slice;
}

} // Public Scope
}; // Constructor



Tags: プログラムメモ
author : HUNDREDSOFT | - | -

XPは終わりですか。

| プログラムメモ | 2014/02/06 Thu 00:00
久しぶりに、WordPressを触る機会があり、

ごちゃごちゃと「Contact Form 7」のメールフォームを、いろいろ変更していたのですが、
なぜか、XP IE8 からメールが送れない。
(XP-IE6 からも送れなかったが、IE6はさすがにもういいでしょう。)

いろいろやり過ぎてしまった感もあり、しばらく格闘したが、
やはり、XP からは送れない。

しょうがないので、CentOS(6.4),PHP(5.3.3),mysql(14.14)に 最新のWordPress(3.8.1)を入れて、
まっさらな WP環境で「Contact Form 7 (3.7)」だけ入れて試したが（WP Multibyte Patchは入れたが）、
XP IE8 からはメールが送れない。
(IE10, FF, Chrome, Safari は OK。)

気になって調べたら、3.7 の更新日はなんと、2014/2/04 だった（一昨日！）。

そう言えばあんまり考えずに、Updateやっちまったな〜と思い、
3.7を消して、3.6をダウンロードして入れ直すと、XP-IE8 からあっさりメールが送れた。
念のため、IE10,IE6,FF,Chrome,Safari からも送ってみたが問題なし。



※追記

http://ja.forums.wordpress.org/topic/69307

こちらに解決策がありました。

Tags: プログラムメモ
author : HUNDREDSOFT | - | -

Javascriptで細線化

| プログラムメモ | 2013/10/19 Sat 21:49
２値化画像のラインを細線化するアルゴリズムは沢山ありますが、
ループ回数がどうしても多くなるので、javascriptでの実装例は少ないです。

ここでは、比較的高速なZhang-Suen, NWGと、田村の方式を実装してみます。

実行サンプルはこちら

Zhang-Suenのアルゴリズムは、比較的単純であるので実装し易いです。

ラスタスキャンにより画素情報を読み込み、
調査画素を含む、3x3pixの画素値に対し、3種の条件を満たしていれば、
調査画素を白(1)から黒(0)に置き換えます。

調査画素(P1)を含む、3x3pixの画素に次のように番号を振る。


条件１．
　外周一周を眺めた時、
　(P2->P3->P4->P5->P6->P7->P8->P9->P2)
　黒→白となる並びが一つだけであること。

　(英文:A(P1)=number of 0,1 patterns(transitions from 0 to 1) in the ordered sequence of P2,P3,P4,P5,P6,P7,P8,P9,P2.
　Condition: A(P1) = 1 )

条件２．
　外周の加算(白なら+1)の結果(B)が、2<= B <=6、を満たす。

　(英文:B(P1)=P2+P3+P4+P5+P6+P7+P8+P9 (number of nonzero neighbords of P1.)
　Condition: 2 <= B(P1) <= 6 )

条件３(1)．
　P2 x P4 x P6 = 0 かつ、P4 x P6 x P8 = 0 を満たす。

条件３(2)．
　P2 x P4 x P8 = 0 かつ、P2 x P6 x P8 = 0 を満たす。

まず、条件１．２．３(1)の判定を行い、条件を満たした場合、該当画素を除去リストに追加する。
全画素調査後、除去リストに画素があれば、その画素を全て黒(0)に変える。なければ終了する。

次に、条件１．２．３(2)の判定を行い、同様の処理を続ける。
以降、条件３(1)(2)を切り替えながら、条件を満たす画素がなくなるまで、繰り返し行う。



NWG (Nagendraprasad-Wang-Gupta)について は、Zhang-Suenと非常に似た処理であるため、
説明は割愛します。
詳細は、ソースコードか、
A note on the Nagendraprasad-Wang-Gupta thinning algorithm をご覧ください。



Zhang-Suenのアルゴリズムが主に計算による判定を行うのに対し、
田村のアルゴリズムは、図形のマッチングを利用します。

大まかな処理としては、
2種類のパターンがあり、それぞれのパターンは、
3x3の除去パターンと3x3の非除去パターンに分かれる。
（ここで言う除去は、白から黒への変更を意味する。）

ラスタスキャンにより画素情報を読み込み、
各画素について3x3のパターンマッチングを行う。

パターンには次のものがある。



【処理概要】

�� パターン１の除去パターンでマッチングすれば�△悄�一致しなければ次画素へ。
�� パターン１の非除去パターンでマッチングすれば次画素へ、一致しなければ��へ。
�� 該当画素を除去リストに追加し、次画素へ。
�� �　銑�の処理を全画素に対して行う。
�� 除去リストに画素があれば、その画素を全て黒(0)に変える。なければ終了。

�� パターン２の除去パターンでマッチングすれば�Г悄�一致しなければ次画素へ。
�� パターン２の非除去パターンでマッチングすれば次画素へ、一致しなければ�┐悄�
�� 該当画素を除去リストに追加し、次画素へ。
�� �Α銑┐僚萢�を全画素に対して行う。
�� 除去リストに画素があれば、その画素を全て黒(0)に変える。なければ終了。
�� �，北瓩襦�


【マッチングコード例】

説明の為に、3x3pixで表現される１つのパターンにP0〜P8の番号を振ります。


P0〜P8を各2bitで表すものとし、
白bit:10, 黒bit:01, どちらでも良いものはbit:00　とします。
各列の先頭2bitに00を付加し、１バイト/列で表現すると、
１パターンは３バイトになります。

例えば、パターン１の除去パターンであれば、


ここで、ある画素を含む3x3pixの状態が次の様になっているとすると

これは、パターン１・非除去パターンの6番目とマッチします。

非除去パターン(6)は、P = 0x040A04で表現されるので

G AND P = P が成り立つ。

この論理式が成り立つ場合に、パターンマッチと認定できます。

田村の方式では、それほどパターンが多くないので、
割と簡単なコーディングで済みます。



<!DOCTYPE html>
<html>
<!-- 
// (c)Hundredsoft Corporation. 2013 All right reserved.
//
//    UTF-8で保存して下さい
-->
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<meta http-equiv="X-UA-Compatible" content="IE=edge" />
<meta name="viewport" content="width=device-width, initial-scale=1.0, user-scalable=no" />
<title>細線化サンプル</title>

<script type="text/javascript">
//
// Zhang-Suen Algorithm
//
var zhangsuen = function(imgdata){
    var w = imgdata.width;
    var h = imgdata.height;
    var ind = imgdata.data;

    var x, y, rAry;
    var bFlag = true;

    for (var k=0; k<100 && bFlag; k++){
        if (!(k & 1)){
            bFlag = false;
        }
        rAry = new Uint8Array(ind);
        for (y=1; y<h-1; y++){
            for (x=1; x<w-1; x++){
                var i = (y*w + x)*4;
                if (rAry[i]){
                    var a,b,p1,p2,p3,p4,p5,p6,p7,p8,p9;
                    // [p9 p2 p3]
                    // [p8 p1 p4]
                    // [p7 p6 p5]
                    p1 = 1;
                    p2 = (rAry[i-w*4  ]) ? 1 : 0;
                    p3 = (rAry[i-w*4+4]) ? 1 : 0;
                    p4 = (rAry[i    +4]) ? 1 : 0;
                    p5 = (rAry[i+w*4+4]) ? 1 : 0;
                    p6 = (rAry[i+w*4  ]) ? 1 : 0;
                    p7 = (rAry[i+w*4-4]) ? 1 : 0;
                    p8 = (rAry[i    -4]) ? 1 : 0;
                    p9 = (rAry[i-w*4-4]) ? 1 : 0;
                    a = 0;
                    if (!p2 && p3){a++;}
                    if (!p3 && p4){a++;}
                    if (!p4 && p5){a++;}
                    if (!p5 && p6){a++;}
                    if (!p6 && p7){a++;}
                    if (!p7 && p8){a++;}
                    if (!p8 && p9){a++;}
                    if (!p9 && p2){a++;}
                    b = p2+p3+p4+p5+p6+p7+p8+p9;

                    if (a == 1 && 2 <= b && b <= 6){
                        if ((!(k & 1) && p2*p4*p6 == 0 && p4*p6*p8 == 0)
                         || ( (k & 1) && p2*p4*p8 == 0 && p2*p6*p8 == 0))
                        {
                            ind[i] = ind[i+1] = ind[i+2] = 0;
                            bFlag = true;
                        }
                    }
                }
            }
        }
    }
};

//
// NWG Algorithm
//
var nwg_method = function(imgdata){
    var w = imgdata.width;
    var h = imgdata.height;
    var ind = imgdata.data;

    var x, y, rAry;
    var bFlag = true;

    for (var k=0; k<100 && bFlag; k++){
        bFlag = false;
        rAry = new Uint8Array(ind);
        for (y=1; y<h-1; y++){
            for (x=1; x<w-1; x++){
                var i = (y*w + x)*4;
                if (rAry[i]){
                    var a,b,c,e,f,p0,p1,p2,p3,p4,p5,p6,p7;
                    // [p7 p0 p1]
                    // [p6    p2]
                    // [p5 p4 p3]
                    p0 = (rAry[i-w*4  ]) ? 1 : 0;
                    p1 = (rAry[i-w*4+4]) ? 1 : 0;
                    p2 = (rAry[i    +4]) ? 1 : 0;
                    p3 = (rAry[i+w*4+4]) ? 1 : 0;
                    p4 = (rAry[i+w*4  ]) ? 1 : 0;
                    p5 = (rAry[i+w*4-4]) ? 1 : 0;
                    p6 = (rAry[i    -4]) ? 1 : 0;
                    p7 = (rAry[i-w*4-4]) ? 1 : 0;
                    a = 0;
                    if (!p0 && p1){a++;}
                    if (!p1 && p2){a++;}
                    if (!p2 && p3){a++;}
                    if (!p3 && p4){a++;}
                    if (!p4 && p5){a++;}
                    if (!p5 && p6){a++;}
                    if (!p6 && p7){a++;}
                    if (!p7 && p0){a++;}
                    b = p0+p1+p2+p3+p4+p5+p6+p7;

                    if (2 <= b && b <= 6){
                        c = 0;
                        if ((p0+p1+p2+p5 == 0 && p4+p6 == 2)
                         || (p2+p3+p4+p7 == 0 && p0+p6 == 2)){
                            c = 1;
                        }
                        if (a == 1 || c == 1){
                            e = (p2+p4) * p0 * p6;
                            f = (p0+p6) * p2 * p4;
                            if ((!(k & 1) && e == 0)
                             || ( (k & 1) && f == 0)){
                                ind[i] = ind[i+1] = ind[i+2] = 0;
                                bFlag = true;
                            }
                        }
                    }
                }
            }
        }
    }
};

//
// 田村 Algorithm
//
var tamura = function(imgdata){
    // [p0 p1 p2]
    // [p3 p4 p5]
    // [p6 p7 p8]
    //
    // P[0-8]を2bitで表す。白bit:10, 黒bit:01, その他bit:00
    // 先頭にbit00を付加し、１バイト/列で１つのパターンを３バイトで表現。
    // Bit並びは、00[p0][p1][p2]00[p3][p4][p5]00[p6][p7][p8]とする。

    // 除去するパターン(1)
    var pat1  = new Array(0x040800, 0x000900);
    // 除去しないパターン(1)
    var pat1n = new Array(0x040a09, 0x182900,
                0x001908, 0x042804, 0x081900, 0x040a04, 0x001824, 0x241800,
                0x060900, 0x000906, 0x192a11, 0x190a19, 0x112a19, 0x192819);

    // 除去するパターン(2)
    var pat2  = new Array(0x000804, 0x001800);
    // 除去しないパターン(2)
    var pat2n = new Array(0x182804, 0x001a09,
                0x001908, 0x042804, 0x081900, 0x040a04, 0x001824, 0x241800,
                0x060900, 0x000906, 0x192a11, 0x190a19, 0x112a19, 0x192819);

    var w = imgdata.width;
    var h = imgdata.height;
    var ind = imgdata.data;

    var bFlag = true;
    for (var k=0; k<100 && bFlag; k++){
        bFlag = false;
        var rAry = new Uint8Array(ind);
        var pat, patn;
        if (k & 1){
            pat = pat2;
            patn = pat2n;
        }else{
            pat = pat1;
            patn = pat1n;
        }

        var x,y;
        for (y=1; y<h-1; y++){
            for (x=1; x<w-1; x++){
                var i = (y*w + x)*4;
                if (rAry[i]){
                    var f = 0x000800;

                    if (rAry[i-w*4-4]){f |= 0x200000;}
                    else              {f |= 0x100000;}
                    if (rAry[i-w*4    ]){f |= 0x080000;}
                    else              {f |= 0x040000;}
                    if (rAry[i-w*4+4]){f |= 0x020000;}
                    else              {f |= 0x010000;}
                    if (rAry[i      -4]){f |= 0x002000;}
                    else              {f |= 0x001000;}
//                    if (rAry[i        ]){f |= 0x000800;}
//                    else              {f |= 0x000400;}
                    if (rAry[i      +4]){f |= 0x000200;}
                    else              {f |= 0x000100;}
                    if (rAry[i+w*4-4]){f |= 0x000020;}
                    else              {f |= 0x000010;}
                    if (rAry[i+w*4    ]){f |= 0x000008;}
                    else              {f |= 0x000004;}
                    if (rAry[i+w*4+4]){f |= 0x000002;}
                    else              {f |= 0x000001;}

                    // 除去するパターンに一致
                    if ((f & pat[0]) == pat[0] || (f & pat[1]) == pat[1]){
                        // 除去しないパターンに一致
                        if ((f & patn[ 0]) == patn[ 0] || (f & patn[ 1]) == patn[ 1]
                         || (f & patn[ 2]) == patn[ 2] || (f & patn[ 3]) == patn[ 3]
                         || (f & patn[ 4]) == patn[ 4] || (f & patn[ 5]) == patn[ 5]
                         || (f & patn[ 6]) == patn[ 6] || (f & patn[ 7]) == patn[ 7]
                         || (f & patn[ 8]) == patn[ 8] || (f & patn[ 9]) == patn[ 9]
                         || (f & patn[10]) == patn[10] || (f & patn[11]) == patn[11]
                         || (f & patn[12]) == patn[12] || (f & patn[13]) == patn[13]){
                            ;
                        }else{
                            ind[i] = ind[i+1] = ind[i+2] = 0;
                            bFlag = true;
                        }
                    }
                }
            }
        }
    }
};

var thinning = function(pattern){
    var chk = document.getElementById('difftm')
    var cvs = document.getElementById('IDcanvas')
    var ctx = cvs.getContext("2d");
    var imgdata = ctx.getImageData(0, 0, cvs.width, cvs.height);
    var stt = new Date;
    if (pattern == 0){
        tamura(imgdata);
    }else if (pattern == 1){
        zhangsuen(imgdata);
    }else{
        nwg_method(imgdata);
    }
    chk.innerHTML = "ConvertTime:" + ((new Date) - stt) + "ms";
    ctx.putImageData(imgdata, 0, 0);
};

var init = function(filename){
    var cvs = document.getElementById('IDcanvas')
    var ctx = cvs.getContext("2d");
    var img = new Image();
    img.onload = (function(){
        ctx.drawImage(img, 0, 0);
    });
    img.src = filename;
};
</script>
</head>

<body onload="init('test.png');">
<canvas id="IDcanvas" width="512" height="512"></canvas>
<hr />
<input type="button" value="画像１" style="width:100px; height:50px;" onclick="init('test.png'); return false;">
<input type="button" value="画像２" style="width:100px; height:50px;" onclick="init('test2.png'); return false;">
<input type="button" value="画像３" style="width:100px; height:50px;" onclick="init('test3.png'); return false;">
<br />
<input type="button" value="田村" style="width:100px; height:50px;" onclick="thinning('0'); return false;">
<input type="button" value="Zhang-Suen" style="width:100px; height:50px;" onclick="thinning('1'); return false;">
<input type="button" value="NWG" style="width:100px; height:50px;" onclick="thinning('2'); return false;">

<hr />
画像を選んでから、田村 or Zhang-Suen oe NWG をクリック！<br />
（Chrome,FireFoxだと0.5秒程度ですが、IE10だと変換に5秒程度掛ります。）<br />
<br />
<div id="difftm"></div>
</body>
</html>

説明では削除リストを使っていますが、元データをTypedArray(Uint8Array)にコピーした方が早かったので、
実装では、変えています。(内容的に変わるものではありません。)

Zhang-Suen,NWG,田村で、どの方法が良いかは、入力画像によって変わるので一概には言えません。
サンプルを動かすとわかりますが、処理時間に関しては、NWGより、Zhang-Suenの方が早いです。
NWGは、Zhang-Suenの改良版であるため若干処理が増えているためです。
人間の目で見ると僅かな違いですが、NWGの方が綺麗になっているようです。

田村に関しては、サンプルの画像１・画像２では速度が劣りますが、画像３では圧倒的に早いです。
文字系では妙なひげが出やすいですが、写真等のエッジ抽出画像に対しては良い結果が出ています。

(論文発表順では、【Hilditch(1968)】→【田村(1978)】→【Zhang-Suen(1984)】→【NWG(1989)】)

今回実装していませんが、
Hilditchの方法は細線化処理の元祖とも言えますが、十分な結果が得られない場合があるらしい。
Rosenfeldは、Zhang-Suen,NWG,田村に比べ、かなり処理が重くなるらしい。
（参考： Comparing Hilditch, Rosenfeld, Zhang-Suen,and Nagendraprasad - Wang-Gupta Thinning ）


実行サンプルはこちら


IE10/Chrome/FF/Androidで、確認しています。

IE9は、TypedArray(Uint8Array)に対応していないので、削除リストを使う方式に修正する必要があります。
IE10でも、Uint8ClampedArray（飽和演算付き8bit）には対応していないので、
このサンプルでは、Uint8Arrayを使っています。
(Chrome,FFでの、ImageData.dataの型はUint8ClampedArrayです。)

このループ数になるとChromeやFFは速いですが、IEだとかなり遅いので、
画像サイズによっては、Web Workersによるマルチスレッド化を検討する必要があります。
但し、IEではImageData.dataの型(CanvasPixelArray)が異なるので、注意が必要です。


--------------- 文献 ---------------
【田村の方式】
H. Tamura: A comparison of line thinning algorithms from digital geometry viewpoint, Proc. 4th Int. Joint Conf. on Pattern Recognition, pp. 715 - 719 (1978.11)

【Zhang-Suenの方式】
Zhang, T. Y. and Suen, Ching Y., “A Fast Parallel Algorithms ForThinning Digital Patterns”,Communication of the ACM,Vol 27, No. 3,Maret 1984, pp.236-239.

【NWGの方式】
P.S.P. Wang and Y.Y. Zhang (1989). A fast and flexible thinning algo-rithm. IEEE Transactions on Computation C-38, 741–745



Tags: プログラムメモ
author : HUNDREDSOFT | - | -

Catmull-Romスプライン

| プログラムメモ | 2013/08/14 Wed 07:37
Catmull-Romスプライン

Edwin CatmullとRaphael Romにより開発された曲線補間ですが、
Edwin Catmullはピクサー・アニメーション・スタジオの現社長だそうです。

私がCatmullのスプラインを最初に調べたのは、８年も前になりますが、
そのころは記事や論文がほとんどなくて苦労した記憶があります。

3次のB-Spline(Cardinal)に比べてやや鋭角な印象があります。
その分、Createrの意思が反映できます。
正確なデザインは、Bezierを使った方が良いですが、
ラフなデザインや、流線型に限れば簡単で重宝します。

今回も簡単にjavascriptでプログラムしてみます。
Opener（開曲線）なカーブは、検索すればいくらでもありますので、
Closed（閉曲線）なカーブを描いてみます。

とは言っても、行列式は決まっているので、
媒介変数を使って2価関数を許し、始点・終点の節点を調節するだけです。

Catmull-Splineでは４つの節点から、３次の補間式を作るので、
始点を、P[n-1],P[0],P[1],P[2]から作り、
終点は、P[n-2],P[n-1],P[0],P[1]から作ります。

関数catmullCircuitに、節点のX,Y座標列を与えると、
節点間の20個の補間座標を返します。

実行サンプルはこちら

<!DOCTYPE html>
<html>
<!-- 
// (c)Hundredsoft Corporation. 2013 All right reserved.
//
//	UTF-8で保存して下さい
-->
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<meta http-equiv="X-UA-Compatible" content="IE=edge" />
<meta name="viewport" content="width=device-width, initial-scale=1.0, user-scalable=no" />
<title>CatMull-Spline Circuit</title>

<script type="text/javascript">
new function() {

var m_x = [];
var m_y = [];

var draw = function(w, h){
	m_ctx.clearRect(0, 0, w, h);

	m_ctx.lineJoin = "round";
	m_ctx.lineCap = "round";
	m_ctx.lineWidth = 2;
	m_ctx.strokeStyle = "#000000";

	var n = m_x.length;
	var i;
	// 節点の描画
	for (i=0; i<n; i++){
		m_ctx.beginPath();
		m_ctx.arc(m_x[i], m_y[i], 4, 0, Math.PI*2, false);
		m_ctx.fill();
		m_ctx.closePath();
	}
	if (n > 2){
		var cc = catmullCircuit(m_x, m_y);
		var nc = cc.x.length;

		// スプライン描画
		m_ctx.beginPath();
		m_ctx.moveTo(cc.x[0], cc.y[0]);
		for (var i=1; i<nc; i++){
			m_ctx.lineTo(cc.x[i], cc.y[i]);
		}
		m_ctx.stroke();
		m_ctx.closePath();
	}

};

var catmullCircuit = function (xi, yi){
	var n = xi.length;
	var x=[], y=[];
	var ox=[], oy=[];

	x.push(xi[n-1]);	y.push(yi[n-1]);
	x.push(xi[0]);	y.push(yi[0]);
	for (var i=1; i<n ;i++){
		x.push(xi[i]);	y.push(yi[i]);
	}
	x.push(xi[0]);	y.push(yi[0]);
	x.push(xi[1]);	y.push(yi[1]);

	n += 1;

	for (var i=1; i<n ;i++){
		var ax = -1 * x[i-1] + 3 * x[i] - 3 * x[i+1] + 1 * x[i+2];
		var bx =  2 * x[i-1] - 5 * x[i] + 4 * x[i+1] - 1 * x[i+2];
		var cx = -1 * x[i-1] + 0 * x[i] + 1 * x[i+1] + 0 * x[i+2];
		var dx =  0 * x[i-1] + 2 * x[i] + 0 * x[i+1] + 0 * x[i+2];

		var ay = -1 * y[i-1] + 3 * y[i] - 3 * y[i+1] + 1 * y[i+2];
		var by =  2 * y[i-1] - 5 * y[i] + 4 * y[i+1] - 1 * y[i+2];
		var cy = -1 * y[i-1] + 0 * y[i] + 1 * y[i+1] + 0 * y[i+2];
		var dy =  0 * y[i-1] + 2 * y[i] + 0 * y[i+1] + 0 * y[i+2];

		// 20分割（荒い）
		for (var j=0; j<=20; j++){
			var t = j / 20;
			ox.push(((ax *t*t*t) + (bx *t*t) + (cx *t) + dx) / 2);
			oy.push(((ay *t*t*t) + (by *t*t) + (cy *t) + dy) / 2);
		}
	}
	return {"x": ox, "y":oy};
};


test2 = {
init: function(){
	var canvas = document.getElementById('c');
	m_ctx = canvas.getContext("2d");

	canvas.addEventListener('touchstart',
		function(e) {
			e.preventDefault();
			var n = e.touches.length;
			var rect=e.target.getBoundingClientRect();
			m_x.push(e.touches[n-1].pageX-rect.left);
			m_y.push(e.touches[n-1].pageY-rect.top);
			draw(canvas.width, canvas.height);
		}, false);

	canvas.addEventListener('mousedown',
		function(e) {
			e.preventDefault();
			var rect=e.target.getBoundingClientRect();
			m_x.push(e.clientX-rect.left);
			m_y.push(e.clientY-rect.top);
			draw(canvas.width, canvas.height);
		}, false);
}
}
};

</script>
</head>

<body onload="test2.init()">
<canvas id="c" width="640" height="640" style="background-color:#e0e0e0;">
</canvas>
</body>
</html>



実行サンプルはこちら


IE9以降/Chrome/Safari/FF/Androidで、確認しています。



Tags: プログラムメモ
author : HUNDREDSOFT | - | -

2D平面で加法定理を使う

| プログラムメモ | 2013/07/28 Sun 06:13
三角関数の加法定理をご存じでしょうか？

高校の数学で習うやつです。
検索しても、実際のプログラムでの使い方や有効点が見つかることは皆無に近い。
まるで、入学試験のためだけにあるようです。

図形の回転が連続するケースでは、角度を保持するプログラムがほとんどです。
しかし、こういったケースでは、sin, cos 値を保存しておいた方が、
if 文が少なくなり、不具合が少なくなったりします。

-------

ある矩形ABCDがあったとして、A点を固定にしてC点を動かす。
変化後の矩形AB’C’D’ を求める。
但し、矩形の縦横比（Aspect比）は変えないものとする。



矩形の拡大率(α)は、次式で求まる。

元の矩形の幅と高さをそれぞれ、w,hとし、変化後の幅と高さをw’,h’ とする。
また、回転角をθ とする。
A,C’ は既知であるので、B’ かD’ のいずれかがわかれば良い。

幾つかの方法が考えられるが、
単純なのは、拡大率(α)と回転角(θ) を求めて、次の行列に代入して計算する方法だろう。

この方法の欠点は、θ を求めると特定の象限において、符号の操作が必要になることだろう。
例えばatanでは、(-π/2 ≦ θ ≦ π/2)の範囲でしか、答えが求まらず、
acosでは(0 ≦ θ ≦ π)の範囲でしか求まらない。
C’ をA点を中心に回転させるとわかりやすい。
この方法は意外と間違いを犯しやすく、テストケースも増えるので、選択すべきではない。

A,C,C’ の座標がわかっているので、
ベクトルAC、AC’ の内積と外積からcosθ,sinθ を算出する。

内積の公式


外積の公式


今回のケースに当てはめ、cos,sinを解くと、


HTML5のcanvasクラスでは、Context.transformで行列演算が可能なので
行列式のパラメタ(cos,sin,-sin,cos)を指定すれば良い。

一度回転させて、後から更に回転を追加する場合を考える。
図では、C’ からC'’ に制御点が移動したとする。


まず前出の方法で、ベクトルAC’、AC'’ の内積と外積からcos(ε),sin(ε)を算出する。

既にθ 回転していて、更にε の回転が加わるとすると、加法定理により


拡大率も同様に、既に拡大率β であり、
これに拡大率α が加わるとすれば、βｘα が合成拡大率になる。


更に発展させて、拡大縮小回転が断続的に連続発生するケースを考えてみます。



サンプルを実行するとわかりやすい。

操作としては、
A点をCanvas中央点として、クリック点をPとする。
　　　　この時画面には、�，龍觀舛�表示されている。

その後、スワイプ操作でP’ まで移動し、ここで指を離す。
　　　　この時画面には、�△龍觀舛�表示されている。

再び、Q点でクリックを行う。
　　　　この時画面には、�△龍觀舛�表示されている。

スワイプでQ’ まで移動し指を離す。
　　　　この時画面には、��の矩形が表示されている。

実現方法はいろいろあるが、
スワイプOFF時に、その時点の拡大率(βxα)・回転状態(sin(θ+ε),cos(θ+ε))を保存しておき、
クリック時に、β＝βxα, sin(θ)＝sin(θ+ε), cos(θ)＝cos(θ+ε)
として再開始すれば良い。

Javascriptでのソースコードは次のようになる。
<!DOCTYPE html>
<html>
<!-- 
// (c)Hundredsoft Corporation. 2013 All right reserved.
//
//	UTF-8で保存して下さい
-->
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<meta http-equiv="X-UA-Compatible" content="IE=edge" />
<meta name="viewport" content="width=device-width, initial-scale=1.0, user-scalable=no" />
<title>拡大縮小回転</title>

<script type="text/javascript">
new function() {
var m_img = new Image();
var m_ctx = null;
var m_O = null;
var m_baseX = 0;
var m_baseY = 0;
var m_cos = 1;
var m_sin = 0;
var m_prop = 1;
var m_Lcos = 1;
var m_Lsin = 0;
var m_Lprop = 1;

var mmove = function(x, y){
	// Click座標-中心 ベクトル
	var vx0 = m_baseX - m_O.x;
	var vy0 = m_baseY - m_O.y;
	var r0 = Math.sqrt(vx0*vx0 + vy0*vy0);

	// スワイプ座標-中心 ベクトル
	var vx1 = x - m_O.x;
	var vy1 = y - m_O.y;
	var r1 = Math.sqrt(vx1*vx1 + vy1*vy1);

	var cos = (vx0*vx1 + vy0*vy1) / (r0*r1); // 内積から
	var sin = (vx0*vy1 - vy0*vx1) / (r0*r1); // 外積から
	var prop = r1 / r0;
	draw(prop, cos, sin);
};
var draw = function(prop, cos, sin){
	m_Lcos = m_cos*cos - m_sin*sin; // 加法定理(cos)
	m_Lsin = m_sin*cos + m_cos*sin; // 加法定理(sin)
	m_Lprop = m_prop * prop;     // 拡大率の加法は乗じるだけ
	m_ctx.clearRect(0, 0, m_O.x*2, m_O.y*2);
	m_ctx.save();

	 // 行列計算
	m_ctx.transform(
		m_Lcos*m_Lprop,  m_Lsin*m_Lprop,
		-m_Lsin*m_Lprop, m_Lcos*m_Lprop,  m_O.x, m_O.y);

	// 先に画像中心を原点に移動(適用順に注意)
	m_ctx.transform(1, 0, 0, 1, -m_img.width/2, -m_img.height/2);

	m_ctx.drawImage(m_img, 0, 0);
	m_ctx.restore();
};

test1 = {
init: function(){
	var mouseDown = false;
	var canvas = document.getElementById('c');

	canvas.addEventListener('touchstart',
		function(e) {
			e.preventDefault();
			mouseDown = true;
			var n = e.touches.length;
			var rect=e.target.getBoundingClientRect();
			m_baseX = e.touches[n-1].pageX-rect.left;
			m_baseY = e.touches[n-1].pageY-rect.top;
			m_cos = m_Lcos;
			m_sin = m_Lsin;
			m_prop = m_Lprop;
		}, false);

	canvas.addEventListener('mousedown',
		function(e) {
			e.preventDefault();
			mouseDown = true;
			var rect=e.target.getBoundingClientRect();
			m_baseX = e.clientX-rect.left;
			m_baseY = e.clientY-rect.top;
			m_cos = m_Lcos;
			m_sin = m_Lsin;
			m_prop = m_Lprop;
		}, false);

	canvas.addEventListener('touchend',
		function(e) {mouseDown = false;}, false);

	canvas.addEventListener('mouseup',
		function(e) {mouseDown = false;}, false);

	canvas.addEventListener('touchmove',
		function (e) {
			e.preventDefault();
			var rect = e.target.getBoundingClientRect();
			if (mouseDown) {
				var n = e.touches.length;
				mmove(e.touches[n-1].pageX-rect.left,
					e.touches[n-1].pageY-rect.top);
			}
		}, false);

	canvas.addEventListener('mousemove',
		function (e) {
			e.preventDefault();
			var rect = e.target.getBoundingClientRect();
			if (mouseDown) {
				mmove(e.clientX-rect.left, e.clientY-rect.top);
			}
		}, false);

	m_img.onload = (function(){
		m_O = {x: canvas.width/2, y: canvas.height/2};
		m_ctx = canvas.getContext("2d");
		draw(1, 1 ,0);
	});
	m_img.width = 260;
	m_img.height = 390;
	m_img.src = "test.jpg";
}
}
};
</script>
</head>
<body onload="test1.init()">
<canvas id="c" width="640" height="640" style="background-color:#e0e0e0;">
</canvas>
</body>
</html>


意外と面倒な処理ですが、コードの通り、演算処理でif文は一つもありません。
この方法の利点は、メンバ変数の使用が少ない事と、θ を計算しない（三角関数の公式は使っていますが、 関数としてのsinやcos は使っていない。使っているのは加減乗除とルートのみ。）ことで、if 文を省ける事が挙げられます。


実行サンプルはこちら


IE9以降/Chrome/Safari/FF/Androidで、確認しています。



Tags: プログラムメモ
author : HUNDREDSOFT | - | -

HTML5のCanvasサンプル

| プログラムメモ | 2013/04/28 Sun 14:12
HTML5でCanvasを使ってお絵かきをするサンプル。


実行サンプルはこちら




ある配管があったとして、配管内部の錆や汚れにより流量が少なくなるケースを想定する。
検査により内部の状態を調べ、検査資料としてまとめる。

特定箇所の断面から、配管長軸の想像図を描く。
このようなアプリがあると、わかりやすい。

ActiveXも寿命なので、タブレットで簡単に操作できるように考案したものです。


このアプリは、およそ１年前に作ったものです。
評判は良かったのですが、売り物にはなりませんでした。
やはり、IE8でcanvasに対応していないことが最大のネックです。
XPユーザーがいる限り、法人向けにHTML5は難しいです。
もう少し我慢ですね。


実際には、配管とは頸動脈のことであり、錆や汚れはポリープです。
カルテに添付するシェーマを予定していたものです。


日の目を見ないのも癪なので、一部を切り出してみました。
IE9以降/Chrome/Safari/FFで、確認しています。
コード一式はこちら



Tags: プログラムメモ
author : HUNDREDSOFT | - | -

ボタンの色を変える

| プログラムメモ | 2013/03/31 Sun 15:31
Javascriptでボタンをクリックする度に背景色を変える場合、
クロージャが利用できる。

←クリック
←クリック 



<body>
<input id="testb1" type="button" value="color" />
<br />
<input id="testb2" type="button" value="color" />
<br />

<script type='text/javascript'>
new function() {
var ctbl = new Array('#00FF00', '#00E0FF', '#FF0000');
var changecolor = function (a, b){
	var x = a;
	var y = b;
	x.style.backgroundColor = ctbl[y];
	return function(){
		y = (y + 1) % (ctbl.length);
		x.style.backgroundColor = ctbl[y];
	};
};
var test = function (){
	var btn;
	btn = document.getElementById('testb1');
	btn.onclick=changecolor(btn, 0);
	btn = document.getElementById('testb2');
	btn.onclick=changecolor(btn, 1);
}

test();

};
</script>
</body>



Tags: プログラムメモ
author : HUNDREDSOFT | - | -

asmxよ、永遠に。

| プログラムメモ | 2012/07/08 Sun 02:51
HTML5でCanvasを使ってお絵かきをして、その結果を画像ファイルとしてサーバーに保存することを考える。

方法としては何も目新しいものはなく、
canvas.toDataURL()で得られる画像(Base64でエンコードが掛っている)を
サーバーに送って、サーバー側でデコードして.pngファイルに落とすだけです。

Perl,PHP,ASP.NET等、何を使っても実現可能ですが、
開発寿命の近づいている、.asmxを使ってWEBサービスで処理してみます。
今ならWCFですが、メンテ等では見る機会があるかもしれないのでメモとして残しておきます。

下記に示す方法を何と言うのかわかりませんが、microsoftajax.jsを使っているので
ブラウザは何でも良く、Client側はMSのルールに縛られずに
（サーバーがIISであることは必須ですが、ASP用の宣言や定義とかは不要、）、
自由な .htmlや .js を並べることができます。
.asmxもBuildの必要はなく、ソースとしてアップできるので、簡単に修正できます。

出力ファイル名(相対パス)を fname、Base64でエンコードされたデータ
(先頭の'data:image/png;base64,'は削除済み)を datauri として、大概は次のようになります。

    [WebMethod]
    public void imgwrite(string fname, string datauri)
    {
        byte[] todecode_byte = null;
        try{
            System.Text.UTF8Encoding encoder = new System.Text.UTF8Encoding();
            System.Text.Decoder utf8Decode = encoder.GetDecoder();

            todecode_byte = Convert.FromBase64String(datauri);

        }catch (Exception e){
            throw new Exception("Error in imgwrite:" + e.Message);
        }

        string filename = HttpContext.Current.Request.PhysicalApplicationPath + fname.Replace("/", "\\");
        FileStream fs = null;
        BinaryWriter bw = null;
        try{
            fs = new FileStream(filename, FileMode.Append, FileAccess.Write);
            bw = new BinaryWriter(fs);
            bw.Write(todecode_byte);

        }catch (Exception e){
            throw new Exception("Error in imgwrite:" + e.Message);
        }finally{
            if (bw != null) bw.Close();
            if (fs != null) fs.Close();
        }
    }

これだけだとつまらないので、javascriptのみでサーバー側のファイル操作を可能にするサービス集合にします。
まずは、UploadService.asmx 

<%@ WebService Language="C#" CodeBehind="UploadService.asmx.cs" Class="UploadService.UploadService" %>

続いて、UploadService.asmx.cs です。 

// (c) Copyright HUNDREDSOFT Corporation 2012
//	All Rights Reserved.
//
//	NAME
//		UploadService.asmx.cs
//
using System;
using System.Data;
using System.Web;
using System.Collections;
using System.Web.Services;
using System.Web.Services.Protocols;
using System.ComponentModel;
using System.IO;
using System.Net;
using System.Text;

namespace UploadService
{
    /// <summary>
    /// imgupload()
    /// imgwrite(fname, datauri)
    /// txtwrite(fname, cname, datatxt)
    /// download(fname, uri)
    /// filecopy(fname, src)
    /// filedelete(fname)
    /// filemove(fname, src)
    /// 
    /// [fname]
    /// 出力ファイル名（相対パス）
    /// 
    /// [uri]
    /// 入力ファイル名（URLアドレス）
    /// 
    /// [src]
    /// 入力ファイル名（相対パス）
    /// 
    /// [datauri]
    /// 画像データをBase64でエンコードした文字列
    /// 
    /// [cname]
    /// 出力文字エンコード(例)
    /// "shift_jis"
    /// "EUC-JP"
    /// "UNICODE"
    /// "UTF-8"
    /// 
    /// [datatxt]
    /// 文字列
    /// 
    /// imguploadは、次のように呼び出す。
    ///  <form method='POST' enctype='multipart/form-data' action='UploadService.asmx/imgupload'>
    ///  <input type='text' name='filename'></input>
    ///  <input type='file' name='file'></input>
    ///  <input type='submit' value='imgupload'></input>
    ///  </form>
    /// 
    /// （注）imgwrite,txtwriteでは、100kを超える文字列は100k以下で分割して送信すること
    /// 
    /// </summary>
    /// 
    [WebService(Namespace = "http://localhost/")]
    [WebServiceBinding(ConformsTo = WsiProfiles.BasicProfile1_1)]
    [System.Web.Script.Services.ScriptService]
    [ToolboxItem(false)]

    public class UploadService : System.Web.Services.WebService
    {
        // function imgupload
        // Multi-formでUploadされたバイナリファイルをサーバー上に書き出します。
        //
        [WebMethod]
        public void imgupload()
        {
            FileStream fs = null;
            try{
                HttpFileCollection hfc = HttpContext.Current.Request.Files;

                if (hfc.Count > 0){
                    string fname = HttpContext.Current.Request.Form["filename"];
                    string filename = HttpContext.Current.Request.PhysicalApplicationPath + fname.Replace("/", "\\");

                    Stream his = hfc[0].InputStream;
                    byte[] bdata = new byte[his.Length];
                    his.Read(bdata, 0, (int)his.Length);

                    fs = new FileStream(filename, FileMode.Append, FileAccess.Write);
                    fs.Write(bdata, 0, bdata.Length);
                }

            }catch (Exception e){
                throw new Exception("Error in imgupload:" + e.Message);
            }finally{
                if (fs != null) fs.Close();
            }
        }

        // function imgwrite
        // Base64でエンコードされたバイナリファイルをサーバー上でデコードして書き出します。
        //
        [WebMethod]
        public void imgwrite(string fname, string datauri)
        {
            byte[] todecode_byte = null;
            try{
                System.Text.UTF8Encoding encoder = new System.Text.UTF8Encoding();
                System.Text.Decoder utf8Decode = encoder.GetDecoder();

                todecode_byte = Convert.FromBase64String(datauri);

            }catch (Exception e){
                throw new Exception("Error in imgwrite:" + e.Message);
            }

            string filename = HttpContext.Current.Request.PhysicalApplicationPath + fname.Replace("/", "\\");
            FileStream fs = null;
            BinaryWriter bw = null;
            try{
                fs = new FileStream(filename, FileMode.Append, FileAccess.Write);
                bw = new BinaryWriter(fs);
                bw.Write(todecode_byte);

            }catch (Exception e){
                throw new Exception("Error in imgwrite:" + e.Message);
            }finally{
                if (bw != null) bw.Close();
                if (fs != null) fs.Close();
            }
        }

        // function txtwrite
        // サーバー上に、様々な日本語コードでテキストファイルに出力します。
        //
        [WebMethod]
        public void txtwrite(String fname, String cname, String datatxt)
        {
            System.IO.StreamWriter sw = null;
            try{
                System.Text.Encoding enc = System.Text.Encoding.GetEncoding(cname);
                // BOMなし
                if (cname == "UTF-8"){
                    enc = new System.Text.UTF8Encoding(false);
                }else if (cname == "UNICODE"){
                    enc = new System.Text.UnicodeEncoding(false, false);
                }
                string filename = HttpContext.Current.Request.PhysicalApplicationPath + fname.Replace("/", "\\");
                sw = new System.IO.StreamWriter(filename, true, enc);
                sw.Write(datatxt);

            }catch (Exception e){
                throw new Exception("Error in txtwrite:" + e.Message);
            }finally{
                if (sw != null) sw.Close();
            }
        }

        // function download
        // URIで指定されたファイルを、サーバー上にダウンロードします。
        //
        [WebMethod]
        public void download(string fname, string uri)
        {
            BinaryReader br = null;
            FileStream fs = null;
            BinaryWriter bw = null;
            try{
                WebClient wc = new WebClient();
                br = new BinaryReader(wc.OpenRead(uri));
                byte[] read_byte = null;
                int read_count = 0;

                string filename = HttpContext.Current.Request.PhysicalApplicationPath + fname.Replace("/", "\\");
                fs = new FileStream(filename, FileMode.Create, FileAccess.Write);
                bw = new BinaryWriter(fs);

                do{
                    read_byte = br.ReadBytes(1024);
                    read_count = read_byte.Length;
                    bw.Write(read_byte);
                } while (read_count == 1024);

            }catch (Exception e){
                throw new Exception("Error in download:" + e.Message);
            }finally{
                if (br != null) br.Close();
                if (bw != null) bw.Close();
                if (fs != null) fs.Close();
            }
        }

        // function filecopy, filedelete, filemove
        // サーバー上でファイルをコピー、削除、移動します。
        //
        [WebMethod]
        public void filecopy(string fname, string src)
        {
            try{
                string psrc = HttpContext.Current.Request.MapPath(src);
                string filename = HttpContext.Current.Request.PhysicalApplicationPath + fname.Replace("/", "\\");
                File.Copy(psrc, filename, true);

            }catch (Exception e){
                throw new Exception("Error in filecopy:" + e.Message);
            }
        }

        [WebMethod]
        public void filedelete(string fname)
        {
            try{
                string filename = HttpContext.Current.Request.PhysicalApplicationPath + fname.Replace("/", "\\");
                File.Delete(filename);

            }catch (Exception e){
                throw new Exception("Error in filedelete:" + e.Message);
            }
        }

        [WebMethod]
        public void filemove(string fname, string src)
        {
            try{
                string psrc = HttpContext.Current.Request.PhysicalApplicationPath + src.Replace("/", "\\");
                string filename = HttpContext.Current.Request.PhysicalApplicationPath + fname.Replace("/", "\\");
                File.Move(psrc, filename);

            }catch (Exception e){
                throw new Exception("Error in filemove:" + e.Message);
            }
        }

        // function getfilenames
        // サーバー上のファイル一覧を取得します。
        //
        [WebMethod]
        public string getfilenames(string src)
        {
            string ret = "";
            try{
                string psrc = HttpContext.Current.Request.PhysicalApplicationPath + src.Replace("/", "\\");
                ret = getdir(psrc, ret);

            }catch (Exception e){
                throw new Exception("Error in getfilenames:" + e.Message);
            }
            return ret;
        }

        private string getdir(string dir, string ret)
        {
            string[] files = Directory.GetFiles(dir);
            foreach (string s in files){
                ret += s + ",";
            }
            string[] dirs = Directory.GetDirectories(dir);
            foreach (string s in dirs){
                ret = getdir(s, ret);
            }
            char[] tc = {','};
            return ret.TrimEnd(tc);
        }
    }
}

最後にテスト用の test.html です。（テスト用なので,IE9(canvasサポート)だけを選ばない。）


<html>
<head>
<title>test UploadService.asmx</title>
<script type="text/javascript" src="microsoftajax.js" ></script>
<script type="text/javascript" src="uploadservice.asmx/js" ></script>
<script type="text/javascript">
function test_imgwrite()
{
	var data =
	'R0lGODlhIQAkAPcAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA' +
	'AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA' +
	'AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA' +
	'AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA' +
	'AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA' +
	'AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA' +
	'AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA' +
	'AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA' +
	'AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA' +
	'AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA' +
	'AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA' +
	'AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA' +
	'/wD//////yH+DEh1bmRyZWRzb2Z0AAAsAAAAACEAJAAACM0A/wkcSLCgwYMIEypcyLDhQn8QI0qcSHEiv4IQ+2' +
	'ncyLGjx40QLw7M2M+hQY0hCZI0WRClP5ECV7Ic6BLmP5kz/9VU6U9jToE7R/b8SJRoSqEVkyo9GnNo0acgX/L0' +
	'+TNoU6o5rd4c+lNnP6Zbsc7UinPsV6lCoao9a3Op24pt38oFG5EjSZd2nZK9ezav37xx574N/Bdv1I50+SrWy5' +
	'ctT8FuCR+ebDiq5Mp9Ke+FvPRyZswfEzMeDRXsVs5wDfJbzbq169euu8qefTAgADs=';

	UploadService.UploadService.imgwrite("save/b.gif", data, 
		function(){
			alert("OK");
		}, function(){
			alert("error");
		}
	);
}

function test_txtwrite()
{
	UploadService.UploadService.txtwrite("save/sjis.txt", "shift_jis", "abcあいうえお0123456789",
		function(){}, function(){});

	UploadService.UploadService.txtwrite("save/UNICODE.txt", "UNICODE", "abcあいうえお0123456789",
		function(){}, function(){});

	UploadService.UploadService.txtwrite("save/UTF-8.txt", "UTF-8", "abcあいうえお0123456789",
		function(){}, function(){});
}

function test_download()
{
	UploadService.UploadService.download("save/a.gif", "http://www.hundredsoft.jp/hundredsoft.gif", 
		function(){
			alert("OK");
		}, function(){
			alert("error");
		}
	);
}

function test_filecopy()
{
	UploadService.UploadService.filecopy("save/b.gif", "save/a.gif", 
		function(){
			alert("OK");
		}, function(){
			alert("error");
		}
	);
}

function test_filedelete()
{
	UploadService.UploadService.filedelete("save/b.gif", 
		function(){
			alert("OK");
		}, function(){
			alert("error");
		}
	);
}

function test_filemove()
{
	UploadService.UploadService.filemove("save/b.gif", "save/a.gif", 
		function(){
			alert("OK");
		}, function(){
			alert("error");
		}
	);
}

function test_getfilenames()
{
	UploadService.UploadService.getfilenames("save",
		function(dirstr){
			//var dirs = dirstr.split(",");
			var dirs = dirstr.replace(/,/g, "\r\n");
			alert(dirs);
		}, function(){
			alert("error");
		}
	);
}
</script>
</head>
<body>
<h3>test UploadService.asmx</h3>
<hr />
 <input type="button" id="test01button" name="test01" value="test_imgwrite" onclick="test_imgwrite()"></input>
 <input type="button" id="test02button" name="test02" value="test_txtwrite" onclick="test_txtwrite()"></input>
 <input type="button" id="test03button" name="test03" value="test_download" onclick="test_download()"></input>
 <input type="button" id="test04button" name="test04" value="test_filecopy" onclick="test_filecopy()"></input>
 <input type="button" id="test05button" name="test05" value="test_filedelete" onclick="test_filedelete()"></input>
 <input type="button" id="test06button" name="test06" value="test_filemove" onclick="test_filemove()"></input>
 <input type="button" id="test07button" name="test07" value="test_getfilenames" onclick="test_getfilenames()"></input>
<br />
<hr />
test_imgupload<br />
<form method="POST" enctype="multipart/form-data" action="UploadService.asmx/imgupload" target="idmy">
保存パス　　： <input type="text" name="filename"></input><br />
画像ファイル： <input type="file" name="file"></input>
 <input type="submit" value="imgupload"></input>
</form>
<hr />
<iframe name="idmy" width="0" height="0" style="visibility:hidden"></iframe>
</body>
</html>

上記、実行環境をまとめたものを置きます。
.NET Framework V4.0が入った環境なら、App_Codeフォルダにソースだけ入れておけば良いのですが、
それ以前の環境だと、binフォルダにコンパイルしたものを入れないと、うまく動きませんでした。
　　 .NET V4.0以降 (App_Code)
　　 .NET V3.5以前 (bin)
違いはBuild済みか、そうでないかだけですので実行に違いはありません。

さて、一番問題になる環境構築について説明しますが、
作業したことがある人なら「釈迦に説法」なので、以降は適当に読み飛ばしてください。



Tags: プログラムメモ
続きを読む＞＞
author : HUNDREDSOFT | - | -

x64でも拡張倍精度を使いたい

| プログラムメモ | 2012/05/13 Sun 00:04
x64では最低でもSSE2まで使えるので、今更FPUを使うことはあまりありません。

もし、倍精度(仮数部53bit)でも精度が足りない場合、
SSE2を使って４倍精度浮動小数点ライブラリを組むか、
４倍精度に対応した x64向けのコンパイラを用意するのが普通でしょう。

しかし、計算速度を要求されるようなケースで、もし64bit精度で足りるなら、
FPUによる拡張倍精度(仮数部64bit)という選択肢もあります。
（例えば10進で、あと1桁だけ精度が欲しい場合にコンパイラまで変えないですよね）

以前にx86用の拡張倍精度クラスを書きましたが、
今度はx64用のサンプルコード。

x64では、インラインアセンブラが使えないので、組み込み関数(Compiler Intrinsics)を使うか、
MASM ( ml64.exe ) を使ってアセンブリをコーディングするかの選択になりますが、
FPUを使うので、ml64用の.asmファイルを作ることにします。
ターゲットは、x64 のWindowsで Visual Studio で開発しているものとします。
他のコンパイラだと呼び出し規約( MSのFastCall )が異なることがあるので注意。

double をオーバーラップするような Ldouble クラスを作る事にします。
メンバとして、拡張倍精度のデータを１つ持ちます。
必要なのは 10バイトですが、切りの良いところで 16バイト取っています。

まず、ヘッダーファイル「Ldouble.h」です。
基本的な演算をクラスに入れて、三角関数,指数関数等を外部関数として呼び出せるようにしています。



アセンブリで作成するのは、extern "C"の関数群です。
これが「Ldouble.asm」です。



最後に、Ldoubleクラスを使ったサンプルコード「main.cpp」です。




VS2010でBuildする際は、新規プロジェクト(VisualC++,Win32,Win32コンソールアプリケーション,空のプロジェクト)で、
ソース( Ldouble.asm, Ldouble.h, main.cpp )を追加して、
構成マネージャーでプラットフォームをx64にします。
新規でx64が選べない場合は、VC++2010の64bit対応のインストールが足りないかもしれません。
その場合は、ml64.exeも入ってないかもしれないので確認して下さい。

続いて、ソリューションエクスプローラの Ldouble.asm を右クリックして、プロパティ。
構成プロパティ・全般で項目の種類→カスタムビルドツールに変更して下さい。

全般の下に出てくるカスタムビルドツールのコマンドライン欄に、
ml64.exe /DWIN_X64 /Zi /c /Cp /Fl /Fo$(IntDir)\$(InputName).obj Ldouble.asm
出力ファイル欄に
$(IntDir)\$(InputName).obj
を指定して下さい。

これで Build が通るようになります。

例外は考えていません。ご利用の際にはFPU制御ワードを適宜設定して下さい。
(Ldouble.asmのdeffcwを変えて下さい。)
(参考) FPU制御ワードについて

 MSB
  F  E  D  C  B  A  9  8  7  6  5  4  3  2  1  0
 +-----------+-----+-----+-----+----------------+
      NC      丸め  精度   NC      例外マスク
 +-----------+-----+-----+-----+----------------+

  例外マスク bit5〜0
   0:無効操作例外
   1:デノーマル例外
   2:零除算例外
   3:オーバーフロー例外
   4:アンダーフロー例外
   5:不正確結果例外

  精度 bit9,8
   00: 24bit(単精度)
   01: reserved
   10: 53bit(倍精度)
   11: 64bit(拡張倍精度)

  丸め bit11,10
   00: 最近値（四捨五入に近い）
   01: 切り捨て(floor)
   10: 切り上げ(ceil)
   11: 零方向切り捨て


Ldoubleへの入力は、double か整数です。
出力は double型 だけです。これは、double へのキャストとして実現されます。
もちろん、LdoubleからLdoubleへのコピーで精度が落ちることはありません。

Ldoubleは、64bit精度を保持していますので、
Ldoubleを使った演算中は常に64bit精度を保ちます。
出力(キャスト) は double(53bit) 精度に丸めた値を返しますが、
Ldouble が持っている値が 53bit に丸められるわけではありません。

この方法はとてもうまく問題を解決してくれるかもしれません。
計算精度ボトムの一部だけをLdoubleによる高精度計算に置き換えて、
その後キャストするだけで、今までの計算ルーチンに戻せます。
例えば、１つの double型のスタック変数を Ldouble にしただけ(1文字修正)で問題が解決するかもしれません。
とても自然に組み込めるのは良いのですが、コードを一見しただけだと変化（影響）がわかりにくいので、
コメントは丁寧に書く必要があります。


さて、Build後にできたexeを実行すると次の結果が得られました。

as:= 1.0000000000000000e+001
bs=: 1.2300000000000000e+002

+  : 1.3300000000000000e+002: 1.3300000000000000e+002
-  : -1.1300000000000000e+002: -1.1300000000000000e+002
abs: 1.1300000000000000e+002: 1.1300000000000000e+002
if : -1.1300000000000000e+002: -1.1300000000000000e+002
*  : 1.2300000000000000e+003: 1.2300000000000000e+003
/  : 1.2300000000000001e+001: 1.2300000000000001e+001
%  : 3.0000000000000000e+000: 3.0000000000000000e+000
pi : 3.1415926535897931e+000: 3.1415926535897931e+000
sin: 4.9999999999999994e-001: 5.0000000000000000e-001
cos: 8.6602540378443871e-001: 8.6602540378443860e-001
sin: 4.9999999999999994e-001: 5.0000000000000000e-001
cos: 8.6602540378443871e-001: 8.6602540378443860e-001
tan: 5.7735026918962573e-001: 5.7735026918962573e-001
ata: 4.8234790710102493e-001: 4.8234790710102499e-001
at2: 9.8279372324732905e-001: 9.8279372324732905e-001
sqr: 1.4142135623730951e+000: 1.4142135623730951e+000
exp: 1.6439050426651380e+005: 1.6439050426651380e+005
ln : 1.2010000000000000e+001: 1.2010000000000000e+001
pow: 3.3234387002712369e+000: 3.3234387002712369e+000
int: 3.0000000000000000e+000: 3.0000000000000000e+000
dec: 3.2343870027123689e-001: 3.2343870027123672e-001

Machin's pi: 3.1415926535897936e+000: 3.1415926535897931e+000
Leibniz's pi(5000000項での真値): 3.1415924535897932384646433832795
Leibniz's pi(5000000)  double(3.1415924535897797e+000) time = 0.036[sec]
Leibniz's pi(5000000) Ldouble(3.1415924535897930e+000) time = 0.073[sec]

最後のLeibniz's piを除いては、基本演算の確認です。
左側がdouble演算による結果で、右側がLdouble(FPU)による結果です。
Leibniz's piはライプニッツの公式によるπの計算ですが、
わざと演算量を増やして、処理時間を比較しています。
5,000,000項で打ち切って、演算結果と演算に要した時間を表示します。

演算時間に関しては、double演算のおよそ倍の時間が掛っていますが、
double演算では累積誤差により 14 桁の精度しかありませんが、FPUでは 16 桁の精度があります。
三角関数などを使うと、FPUの方が10倍くらい処理時間が掛ることもあります。
doubleに比べて、11 bit多くの精度を出すため処理時間が掛かるようです。
しかし除算については、double演算とFPUで処理時間は同等でした。

逆アセンブラを見てみると、double演算ではSSE2が使われていますが並列化までは行われていません。
うまく並列化（最適化）できれば、double演算はもう少し早くなります。
(Release-Buildではdouble演算側に多少のベクトル化が行われていました。)

正確さと処理時間のトレードオフですので、どちらが良いというものではありません。
目的に即した解決策が見つけられなければ、このような方法もありかと思います。

上記、掲載ソースを圧縮したものをここに置いておきます。


Tags: プログラムメモ
author : HUNDREDSOFT | - | -

πのBBP計算でモンゴメリを使う

| プログラムメモ | 2012/05/06 Sun 09:56
現在では10兆桁を超えるπが求められているが、πの特定の桁だけを算出する方法が知られていて、
πの部分算出の記録では500兆桁を超えている。

この部分算出の元になったのが、次の公式です。


この式からd桁目を求めるために、16^dを乗じて、小数点以下を取りだすことを考える。


ここでSを次のように定義する。

すると


これが1995年にSimon Plouffe(プラウフ)により発見された、BBP(Bailey-Borwein-Plouffe)公式と呼ばれるもので、
これをプログラムしたものが、
http://www.experimentalmath.info/bbp-codes/ にあります。(piqpr8.c)
ここでは右向きバイナリ法を使って剰余(mod)を求めているが、剰余計算にモンゴメリ法を使うことを考える。

モンゴメリ法については過去にも
http://www.hundredsoft.jp/win7blog/log/eid93.html や
http://www.hundredsoft.jp/win7blog/log/eid98.html でサンプルを載せていますが、
高速なモンゴメリ法では、法は奇数であることが望ましい。

S1,S5の法は、8k+1,8k+5であるので奇数であるが、
S4,S6は偶数となるので、このままでは使えない

整数論の定理の中に次のものがある。
ならば

が成り立つ。

これをS4の第一項に当てはめると、

であるから、


同様にS6の第一項は、


これで法が奇数となり、高速なモンゴメリ法を使うことができる。
プログラムは次のようになる。 (defaultではπの100万桁(Hex)目を算出する。)


#include <stdio.h>
#include <math.h>
#include <time.h>

int expmong (int p, int ak);

int main()
{
  double pid, s1, s2, s3, s4;
  double series (int m, int n);
  void ihex (double x, int m, char c[]);
  int id = 1000000;
#define NHX 16
  char chx[NHX];

  clock_t start,end;

/*  id is the digit position.  Digits generated follow immediately after id. */
  start = clock();

  s1 = series (1, id);
  s2 = series (4, id);
  s3 = series (5, id);
  s4 = series (6, id);
  pid = 4. * s1 - 2. * s2 - s3 - s4;
  pid = pid - (int) pid + 1.;
  ihex (pid, NHX, chx);

  end = clock();
  printf (" position = %i\n fraction = %.15f \n hex digits =  %10.10s\n time = %.3f[sec]",
  id, pid, chx, (double)(end-start)/CLOCKS_PER_SEC);

  return 0;
}

void ihex (double x, int nhx, char chx[])

/*  This returns, in chx, the first nhx hex digits of the fraction of x. */

{
  int i;
  double y;
  char hx[] = "0123456789ABCDEF";

  y = fabs (x);

  for (i = 0; i < nhx; i++){
    y = 16. * (y - floor (y));
    chx[i] = hx[(int) y];
  }
}

#define eps 1e-17

double series (int m, int id)

/*  This routine evaluates the series  sum_k 16^(id-k)/(8*k+m) 
    using the modular exponentiation technique. */

{
  int k, ak, p;
  double s, t;

  s = 0.;

/*  Sum the series up to id. */

  if (m == 4){
    for (k = 0; k < id; k++){
      ak = 2 * k + 1;
      p = expmong (id - k - 1, ak);
      t = (p * (4 % ak)) % ak;
      s = s + t / ak;
      s = s - (int) s;
    }
  }else if (m == 6){
    for (k = 0; k < id; k++){
      ak = 4 * k + 3;
      p = expmong (id - k - 1, ak);
      t = (p * (8 % ak)) % ak;
      s = s + t / ak;
      s = s - (int) s;
    }
  }else{
    for (k = 0; k < id; k++){
      ak = 8 * k + m;
      t = expmong (id - k, ak);
      s = s + t / ak;
      s = s - (int) s;
    }
  }

/*  Compute a few terms where k >= id. */

  for (k = id; k <= id + 100; k++){
    ak = 8 * k + m;
    t = pow (16., (double) (id - k)) / ak;
    if (t < eps) break;
    s = s + t;
    s = s - (int) s;
  }
  return s;
}

// tのモンゴメリリダクション
int reduction(__int64 t, int n, int n2, int nb)
{
  static int masktbl[] = {
      0,
      0x00000001,0x00000003,0x00000007,0x0000000f,0x0000001f,0x0000003f,0x0000007f,0x000000ff,
      0x000001ff,0x000003ff,0x000007ff,0x00000fff,0x00001fff,0x00003fff,0x00007fff,0x0000ffff,
      0x0001ffff,0x0003ffff,0x0007ffff,0x000fffff,0x001fffff,0x003fffff,0x007fffff,0x00ffffff,
      0x01ffffff,0x03ffffff,0x07ffffff,0x0fffffff,0x1fffffff,0x3fffffff,0x7fffffff,0xffffffff,
  };
  __int64 c = t * n2;
  c &= masktbl[nb];    // mod Rの意味,(nb)bit以上を0クリア
  c = c * n + t;
  c >>= nb;            // 1/r の意味
  if (c >= n) c -= n;
  return c;
}

// モンゴメリ冪剰余計算
int expmong (int p, int ak)
{
/*  expm = 16^p mod ak. */

  int nb = 0;          // 有効Bit数
  int r;               // r: nより大きな２のべき
  for (nb=0, r=1; r<ak; nb++, r<<=1) ;

  // r2 = r^2 mod ak
  int r2 = ((__int64)r * r) % ak;

  // n*n2 ≡ -1 mod r
  int n2 = 0; /* 求めるN' */
  {
    int t=0, vi=1, ni=nb;
    while (ni--){
      if ((t & 1) == 0){
        t += ak;
        n2 += vi;
      }
      t >>= 1;
      vi <<= 1;
    }
  }

  // 冪剰余 16^p mod ak, バイナリ法の下位桁から計算
  int q = reduction(16 * r2, ak, n2, nb);
  int x = reduction(r2, ak, n2, nb);
  while (p > 0){
    if (p & 1){
      x = reduction((__int64)x * q, ak, n2, nb);
    }
    q = reduction((__int64)q * q, ak, n2, nb);
    p >>= 1;
  }
  return reduction(x, ak, n2, nb);
}

改変は、series() のm=4,m=6の処理とexpm()→expmong()です。

intel系x86コードでの64bit整数演算は早くないので、面白い結果が出ました。

BBP計算速度テスト (x86用,x64用ターゲットを変えてBuildしている)

【右向きバイナリ法(WindowsXP 32bit)】
(1回目)
position = 1000000
fraction = 0.423429797567895
hex digits = 6C65E52CB4
time = 2.763[sec]
(2回目)
position = 1000000
fraction = 0.423429797567895
hex digits = 6C65E52CB4
time = 2.804[sec]
(3回目)
position = 1000000
fraction = 0.423429797567895
hex digits = 6C65E52CB4
time = 2.763[sec]

【モンゴメリ法(WindowsXP 32bit)】
(1回目)
position = 1000000
fraction = 0.423429797567801
hex digits = 6C65E52CB4
time = 3.114[sec]
(2回目)
position = 1000000
fraction = 0.423429797567801
hex digits = 6C65E52CB4
time = 3.134[sec]
(3回目)
position = 1000000
fraction = 0.423429797567801
hex digits = 6C65E52CB4
time = 3.204[sec]

【結果】
およそ15%、モンゴメリ法が時間が掛る

【右向きバイナリ法(Windows7 64bit)】
(1回目)
position = 1000000
fraction = 0.423429797567895
hex digits = 6C65E52CB4
time = 1.788[sec]
(2回目)
position = 1000000
fraction = 0.423429797567895
hex digits = 6C65E52CB4
time = 1.801[sec]
(3回目)
position = 1000000
fraction = 0.423429797567895
hex digits = 6C65E52CB4
time = 1.786[sec]

【モンゴメリ法(Windows7 64bit)】
(1回目)
position = 1000000
fraction = 0.423429797567895
hex digits = 6C65E52CB4
time = 1.255[sec]
(2回目)
position = 1000000
fraction = 0.423429797567895
hex digits = 6C65E52CB4
time = 1.258[sec]
(3回目)
position = 1000000
fraction = 0.423429797567895
hex digits = 6C65E52CB4
time = 1.223[sec]

【結果】
およそ30%、モンゴメリ法が早い

VC2010でBuildした限りでは、
32bitではバイナリ法の方が早く、64bitではモンゴメリ法の方が早くなりました。

バイナリ法のコードも整数演算化してみたのですが、
部分的なコード変更だと、かえって成績が悪くなってしまいました。
バイナリ法のコンパイラオプションでSSE2を有効にしてみてもあまり効果はなかったのですが、
Intelコンパイラを使うと、もう少し早くなるかもしれません。

このコードでも、２〜３千万桁までは計算できますが、
expmong()の変数 r2 と series()の変数 p の型を__int64にすれば、１億桁まで計算できるようになります。
但し、x86(32bit)だと64bit整数演算が増えるとそれだけ処理が遅くなります。
x64で１億桁の計算は、私のPCでは150秒ほどでした。
拡張倍精度(仮数部64bit)と長整数(80bit程度)を組み合わせれば、このアルゴリズムでも100億桁ぐらいは行けると思いますが、これだけの計算量になると普通のPCでは難しいです。（電気代も上がるし...）



Tags: プログラムメモ
author : HUNDREDSOFT | - | -

多倍長整数でSteinのGCDを使う

| プログラムメモ | 2012/04/22 Sun 22:43
次のような、32bit系の多倍長整数クラスがあったとして、


class Bigint
{
public:
	Bigint();
	...
	... 様々な定義(四則演算、operator=,operator<<=等)
	...
private:
	int   m_sz;   // 配列数
	int*  m_num;  // 配列(リトルエンディアン)
};

16進で次のように表せる整数があれば
0x111100002222000033330000

このクラスでは次のように格納されるものとする。(注1)

(m_sz=3)
|- m_num[0]-|- m_num[1]-|- m_num[2]-|
| 3333 0000 | 2222 0000 | 1111 0000 |

これを用いて、Stein's algorithmによるGCD計算を考える。


SteinのGCDについてはWikiPediaが詳しい。
http://en.wikipedia.org/wiki/Binary_GCD_algorithm

一般に ユークリッドの互除法(拡張) の方がループ回数は少ないが、
剰余計算がネックになり、応答性が悪くなる場合が多い。

但し、実装する乗算・剰余算のアルゴリズムや数値の桁数によっては、
ループ回数の少ないユークリッドの互除法(拡張)の方が早い場合もある。

Steinのアルゴリズムで遅いと思われるのは、Shift演算の回数に尽きる。
サンプルコードをそのまま使うと、いくら負荷の少ないShift演算と言えども、
多桁演算で処理回数が多くなると大きな遅延を招きます。
WikiPediaのSteinのサンプルコードで問題となるのは、


/* From here on, u is always odd. */
do {
  while ((v & 1) == 0)  /* Loop X */
     v >>= 1;

の部分。要は右端の０を除去しているだけである。
右端に並んでいる０の数を数える関数を用意して、Shift回数を減らことを目指す。

http://en.wikipedia.org/wiki/Count_trailing_zeros にもあるが、
intel系であれば、アセンブラ(bsf)が使える。


//
// 高速化のためゼロチェックを行っていない
//  ゼロだとアクセス違反の例外。
//
int Bigint::bsf() const{
	int* p = m_num;
	__asm{
		mov edx, [p]
		mov eax, edx
bsf_l2:
		bsf ecx, [eax]
		jnz bsf_l1
		add eax, 4
		jmp bsf_l2
bsf_l1:
		sub eax, edx
		shr eax, 2
		mov edx, BITPERBLOCK
		mul edx
		add eax, ecx
	}
}

BITPERBLOCKは、配列１つに含まれる多倍長整数のBit数である。
32bitであれば、32が使われる場合が多い。(注2)
(64bitの場合、bsfqを使う。)

(m_sz=2, BITPERBLOCK=32)の次のような多倍長整数があるとき

|- m_num[0]-|- m_num[1]-|
| 0000 0000 | ABCD 0000 |
bsf()は、32+16=48 を返す。


アセンブラを使わない次の方法もある。


int Bigint::bsf() const{
	static int table[32] = {0,1,10,2,11,14,22,3,30,12,15,17,19,23,26,4,31,9,13,21,29,16,18,25,8,20,28,24,7,27,6,5};

	int count = 0;
	for (int i=0; i<m_sz; i++){
		int x = m_num[i];
		if (x == 0){
			count += BITPERBLOCK;
		}else{
			count += table[((unsigned int)(x & -x) * 0x07C4ACDDUL) >> 27];
			break;
		}
	}
	return count;
}

// 64bit/Blockの場合
int Bigint::bsf64() const{
	static int table[64] = {
		0,1,59,2,60,40,54,3,61,32,49,41,55,19,35,4,62,52,30,33,50,12,14,42,56,16,27,20,36,23,44,5,
		63,58,39,53,31,48,18,34,51,29,11,13,15,26,22,43,57,38,47,17,28,10,25,21,37,46,9,24,45,8,7,6};

	int count = 0;
	for (int i=0; i<m_sz; i++){
		__int64 x = m_num[i];
		if (x == 0){
			count += BITPERBLOCK;
		}else{
			count += table[((unsigned __int64)(x & -x) * 0x03F566ED27179461ULL) >> 58];
			break;
		}
	}
	return count;
}

http://en.wikipedia.org/wiki/Count_trailing_zeros にもありますが、
Ｍ系列を使ったマジックです。
こちらの方が詳しいです。 http://d.hatena.ne.jp/siokoshou/20090704


すると、Steinのコードは次のように書ける。


//
// *this = GCD(*this, v)
//
void Bigint::gcd(const Bigint& v) {
	if (v == 0){
		return;
	}
	if (operator==(0)){
		operator= (v);
		return;
	}

	Bigint x[2] = {*this, v};

	if (x[0] < 0)	x[0].operator-();
	if (x[1] < 0)	x[1].operator-();

	int i0 = x[0].bsf();
	int i1 = x[1].bsf();
	int shift = (i0 < i1) ? i0 : i1;
	x[0] >>= i0;
	x[1] >>= shift;

	int id = 0;
	do {
		x[id^1] >>= x[id^1].bsf();

		if (x[id] > x[id^1]){
			id ^= 1;
		}
		x[id^1] -= x[id];
	} while (x[id^1] != 0);

	operator= (x[id] << shift);
}

これで30〜40%は、処理時間を短縮できました。

これこそがSteins;Gateの選択♪

このコードを実装したサンプルの16進・２進・10進電卓です。
整数だけですがおよそ100万桁まで計算できます。


(注1) 固定長ではないBigintクラスでは、伸長・縮退等のメモリの管理も必要になるので、これほど単純なものではありません。
(注2) 32bit(Block)全てを使うことが必ずしも有利ではありません。SSEを有効に使うなら桁あふれが発生しないよう 31bit 以下にします。

Tags: プログラムメモ
author : HUNDREDSOFT | - | -

カラツバ法で似非ゴースト暗算

| プログラムメモ | 2012/04/15 Sun 17:13

ゴースト暗算より早い！...場合もある計算方法....(汗


ゴースト暗算というのを、少し前に知りました。
ゴースト暗算が筆算であるなら、乗算４回・加算３〜４回の演算になります。

例えば、48x53を考えてみると

【通常の筆算】
     4 8
   x 5 3
---------
     2 4
   1 2
   4 0
 2 0
---------
 2 5 4 4


【ゴースト暗算】
千の位と一の位を先に出して、残りを後から加算します。
通常の筆算の加算の順序を入れ替えていることがわかります。
(詳しくは書籍をご覧ください。魚の絵を使ったわかりやすいものです)


�� (4x5)x10 + 8x5 = 240          ...千の位(2)
��                   40x10 + 8x3
                  =  424         ...一の位(4)
��                   42 + 4x3    ...正確には一の位と十の位で分解して加算します
                  =  54
-----------------------------
                    2544


【カラツバ】
それならばと、ゴーストっぽいカラツバ法を考えてみました。
とは言っても単純に 基数を１０ にしただけですが。
同様に千の位と一の位を先に出して、残りを後から計算します。

�� (4x5)x100 + (8x3) = 2024  上位桁同士の乗算と下位桁同士の乗算なので繰り上がりはない。並べるだけ。
�� (4x5)     + (8x3) =  44   �，両絨棉�分と下位部分の加算。20+24+2
------------------------------
                        46
�� -(4-8)x(5-3)      =   8   カラツバ法のミソ。
------------------------------
                        54
------------------------------
                       2544

上の説明ではわかりにくいので、計算方法を簡単なパターンにします。

�‐絨矛紊伐式矛綟瓜里魍櫃韻栃造戮泙后�(→2024)
��下３桁を足します(0+2+4→6)
��上３桁を足します(2+0+2→4)
...ここで元の十の位と一の位の数「0,2」は忘れて、2464を記憶する。
�こ櫃韻訖堯�掛けられる数、それぞれ十の位と一の位の差を取って掛けます。(-4x2→-8)
�ド箙翦薪�(-8→+8)
�� ��と�△任任�た数に、�イ凌瑤鯊�します。(46+8→54)
併せて繋げると、2544の出来上がり。

�い�カラツバ法のミソで、たすき掛けの掛け算(4x3,8x5)は必要ありません。

いくつか計算サンプルを並べます。

これもカラツバの効果が良く出ています。

最も効果が出るのは、十と一の位が同じ数字の場合で、

簡単すぎて驚きます。

残念ながら全ての数字の並びで簡単ではありません。
次の例をご覧ください。

この例だと、桁上げと桁下げが発生するので、少し面倒です。

カラツバ法について詳しくは
http://ja.wikipedia.org/wiki/%E3%82%AB%E3%83%A9%E3%83%84%E3%83%90%E6%B3%95
一般には、分割統治法として有名で、乗算の他にも除算や平方根の多倍長数計算で使われます。(注)
Karatsuba博士 もこんな使い方されるとは、思わなかったでしょう。


【まとめ】
幾つかの例では、カラツバ法の方が早いような気もしますが、
ゴースト暗算の良さは、単純化された小さな加算を組み合わせることで、小学校低学年でも計算できるところにあります。
パターンによって難易度があまり上下するのは、低学年では厳しいかもしれません。

カラツバ法の場合、乗算回数が減りますが�キΔ�減算になる場合もあるのでここが暗算上の鬼門になります。
例えば、乗数や被乗数に22や77など、同じ数字が並んでいる場合は�きイ�いらなくなるので簡単ですが、最後の例などカラツバ暗算で解くのは減算による桁下げがあるので面倒です。
また、小さな桁数であれば人間にとっては、加算より乗算の方が早いので、
本来のカラツバ法の効果は、全くと言っていい程ありません。(注)

２桁２桁までしか通用しないとか「そろばん」の方が良いとか批判的な意見も散見しますが、
ゴースト暗算は、イメージ(魚パターン)で簡単にわかりやすく計算方法を解説したという点では、
今までなかった新しい試みであり、ちょっと真似してみたくなりました。


どうでも良いことですが、どこのページも「48 x 53」を参考にしているので、
48x53の答えだけは覚えちゃっている人がいるのではないだろうか。


Tags: プログラムメモ
続きを読む＞＞
author : HUNDREDSOFT | - | -

Windows7で簡単にマルチタッチを試す

| プログラムメモ | 2012/04/01 Sun 10:23
Windows7で簡単にマルチタッチを試す方法として、 "Multi-Touch Vista" があります。
名前の通り、それほど新しいソフトではないのでご存じの方もいるでしょう。

２つのマウス入力を２つのタッチポイントに見立てて、
マルチタッチ環境を実現してくれるドライバとコントロールソフトです。

Virtual環境でこれを何とか動かせないか試したのですが、
VMware, VirtualPC, VirtualBox いずれも動かせませんでした。

ドライバの認識まではできて、コントロールパネルに
「ペンとタッチ」が追加されるところまでは行くのですが
Multitouch.Configuration.WPF.exeの"Configuer device"で
Block native windows...のチェックを入れると
マウスが一切効かなくなります。
やはり実環境が必要なのかもしれません。

実環境であれば、概ね問題なく使えるのですが、
スタートメニューやタスクマネジャーを開けるとちょっと困った問題が発生します。
�．好拭璽肇瓮縫紂湿紊法▲織奪舛寮屬ぁ�のマークが表示されない。
�▲織好�マネージャー上で操作できない。

�△砲弔い討蓮�
Multitouch.Service.Console.exe,Multitouch.Driver.Console.exeを
管理者として実行すれば可能ですが�，量簑蠅�あるので、
管理者権限で実行しているWindow上にタッチマークを表示することができずに、
管理者アプリのWindow背面にマークが表示されてしまいます。

"Multi-Touch Vista"は、コード公開されているので
multitouchvista-30935.zipをダウンロードして
Multitouch.InputProviders.slnをVS2010で開くと、
DebugCursor.csがマークを表示していることがわかります。

Z-Orderに関する問題ですが、常に最前面にユーザーWindowを表示する方法は、
WindowsＯＳのバージョンによっても変わります。

・Formの指定やCreateWindowでTopMostを指定する
・SetForegroundWindowsやAttachThreadInputの呼び出し
・BringWindowToTopの呼び出し

等がありますが、

Windows7でC#なら、Control.BringToFront()を使うのが簡単です。

DebugCursorを呼び出しているのは、DeviceStatusクラスですので
SyncUpdateLocation()に一行追加します。


void SyncUpdateLocation(object state)
{
   debugCursor.BringToFront(); //<-------------- この行を追加
   debugCursor.Location = new Point(Location.X - debugCursor.Width / 2, Location.Y - debugCursor.Height / 2);
}


これをBuildして、Outputフォルダに出力されたファイルを上書きすれば、
�．好拭璽肇瓮縫紂湿紊法▲織奪舛寮屬ぁ�のマークが表示されない。
�▲織好�マネージャー上で操作できない。
の問題が解決します。



Tags: プログラムメモ　VirtualPC関連　Windows7(x64)関連
author : HUNDREDSOFT | - | -

Android着信音のループ再生

| プログラムメモ | 2012/02/26 Sun 08:38
昔の着メロ屋さんだと、着信音のループ部分にかなりのこだわりがあり、
繰り返し部分がスムーズになるよう調節を行っています。
いわゆる、ループ素材というやつです。

既に様々なサイトで紹介されていますが、
Android携帯では、MP3,wav等、様々な音楽ファイルを着信音に設定できます。
しかし、自作の着信音を設定すると、１回しか再生されないため、
短い音声で着信音を作成すると、１回だけ再生した後、無音になってしまいます。

但し、Oggで作成した音声ファイルでメタタグ「ANDROID_LOOP=true」を設定した場合はループ再生が可能になります。

もちろんAudacity等でOggを作成していれば、メタタグ編集で入れられますが、
音楽編集に他のソフトを使っていたりすると、面倒です。
開発を止めていたwav2mldでも簡単に追加できそうなので、方法を考えてみました。

PCでのプログラムとしては、OggVorbis標準の http://xiph.org のものを使ったとして
libvorbis-1.x.x\lib\info.c の vorbis_comment_add_tagをCallすれば良い。

例えば、encoder_example.c であれば、

vorbis_comment_init(&vc);
vorbis_comment_add_tag(&vc,"ENCODER","encoder_example.c");
の後に、
vorbis_comment_add_tag(&vc, "ANDROID_LOOP", "true");
を入れれば、Android用ループ音声が作成できます。

作成したoggファイルは、どこかにUpLoadしてダウンロードしても良いし、
Wi-Fi が使える環境なら、Es explorerなどでPCから直接コピーすることもできます。
コピー出力場所は、
/sdcard/media/audio/ringtones
です。

wav2mld Ver.2.35 で機能追加しました。
これで、自作着信音がAndroid携帯上でループ再生できるようになります。

Tags: プログラムメモ
author : HUNDREDSOFT | - | -

ニュートン法を使った２進多桁の整数除算(その２)

| プログラムメモ | 2011/02/20 Sun 08:56
前回、除算 x = p / q　をNewton法を使って解くプログラムを紹介しましたが、
単独の初期値を使った方法では、あまり早くありません。
数万bitを超える多桁演算の場合、64bit精度の初期値を使ったとしても、あまり速度に差が出ません。
そこで大きな桁(数万bitを超える)で計算を開始するのではなく、除数(q)の先頭部分を使った小さな桁で計算を開始して、これをNewton法の初期値とし、もう少し大きな桁でNewton法を解く、
これを繰り返して最終的な答えを得る方法を紹介します。
少しコード修正しました(2011.2.11) 

Bigint newtondiv(const Bigint& q, const int& n, const Bigint& init){
  Bigint x(init), m(0), c2(2);
  c2 <<= n;
  while(m != x){
    m = x;
    x *= c2 - q * x;
    x >>= n;
  }
  return x;
}

Bigint iDiv(Bigint p, Bigint q){
  int pb = p.length(); // 被除数のbit数
  int qb = q.length(); // 除数のbit数
  int n;
  int b = 64;
  int step = qb / 4;

  Bigint d = 1;
  d <<= b;      // Newton法の初期値として, 1 << b を与える。

  d = newtondiv(q >> (qb-b), b*2, d);  // 1回目
  b += step;
  d <<= step;

  d = newtondiv(q >> (qb-b), b*2, d);  // 2回目
  b += step;
  d <<= step;

  d = newtondiv(q >> (qb-b), b*2, d);  // 3回目
  d <<= (pb-b);
  n = pb + qb;

  d = newtondiv(q, n, d);  // 4回目

  d *= p;        // d = p / q
  d >>= n;

  // 誤差修正
  if (p >= (d+1)*q){
      d += 1;
  }
  return d;
}

newtondivは、Newton法で （2n / q） を求める関数です。
初回は先頭64bit程度で計算します。
2回目は、初回の計算結果を初期値として、64+(除数のBit数/4)での計算。
3回目は、２回目の計算結果を初期値として、64+(除数のBit数/2)での計算。
4回目は、３回目の計算結果を初期値として、全桁での計算。
と計算桁を増やしていきます。

除数・被除数のとり方にもよりますが、計算時間は明らかに短くなります。
この方法は、Bigintの演算回数としては、効果はないのですが、
可変長のFFT乗算を用いる場合、演算量は２の冪数桁単位に増加するので、
Totalの計算量としては、こちらの方が有利です。

4分割としていますが、利用するBigint(多倍長整数)クラスで最速となるよう調節する必要があります。
これについては、「高橋 大介, 金田 康正. "多倍長平方根の高速計算法". 情報処理学会研究報告 95-HPC-58, pp.51-56.」が詳しいです。
文献に従えば、初期値に倍精度結果を使う・３〜４分割程度・FFT再利用、を行えば最大限の高速化を望めますが、このサンプルコードではそこまで対応できていません。

ここでは説明を簡単にするため初期値に（2^64）を与えています。
より正確な初期値のとり方や、詳細アルゴリズムはこちら
整数除算の詳細解説

除数・非除数の桁差が小さければ、このような方法もあります。
回復法を使った２進多桁の除算

同様の計算方法を行っている、整数平方根解説はこちら
整数平方根での詳細解説

このアルゴリズムを利用した16進・２進・10進電卓です。
整数だけですがおよそ100万桁まで計算できます。


(注)ここでのBigintクラスは、多倍長整数を格納する仮想のクラスであり、実在するものではありません。

Tags: プログラムメモ
続きを読む＞＞
author : HUNDREDSOFT | - | -

ニュートン法を使った２進多桁の整数平方根(その２)

| プログラムメモ | 2011/02/06 Sun 10:27
前回、平方根 √s　をNewton法を使って解くプログラムを紹介しましたが、
単独の初期値を使った方法では、あまり早くありません。
数万bitを超える多桁演算の場合、64bit精度の初期値を使ったとしても、あまり速度に差が出ません。
そこで大きな桁(数万bitを超える)で計算を開始するのではなく、(s)の先頭部分を使った小さな桁で計算を開始して、これをNewton法の初期値とし、もう少し大きな桁でNewton法を解く、
これを繰り返して最終的な答えを得る方法を紹介します。


Bigint newtonsqrt(const Bigint& q, const int& n, const Bigint& init){
  Bigint x(init), m(0), c3(3);
  c3 <<= n;
  while(m != x){
    m = x;
    x *= c3 - ((m * m * q) >> n);
    x >>= n+1;
  }
  return x;
}

Bigint iSqrt(Bigint s){
  int n = s.length(); // bit数
  int b;
  if (n & 1)    b = 65; // 初期演算Bit数(n:奇数)
  else          b = 64; // 初期演算Bit数(n:偶数)

  Bigint x = 1;
  x <<= (b / 2);  // Newton法の初期値として, 1 << (b/2) を与える。

  Bigint w;
  int m = n / 8;  // 4分割程度

  while(1){
    w = s >> (n - b);        // (s)先頭から指定bit数(n-b)切り取り
    x = newtonsqrt(w, b, x);

    if (b + 2 * m > n)
      break;

    b += 2 * m;
    x <<= m;
  }
  x <<= (n - b) / 2;         // 必ず割り切れる
  x = newtonsqrt(s, n, x);

  x *= s;         // x = s * (2^n) / √s = (2^n) * √s
  x >>= n;        // x = √s

  // 誤差修正
  w = x + 1;
  if (s >= w * w){
     x = w;
  }
  return x;
}

newtonsqrtは、Newton法で （2n / √q） を求める関数です。
初回は先頭64bit程度で計算します。
2回目は、初回の計算結果を初期値として、(n/8)桁増やして計算。
3回目は、２回目の計算結果を初期値として、更に(n/8)桁増やして計算。
と次々と桁を増やしていきます。

計算時間は明らかに短くなります。
この方法は、Bigintの演算回数としては、効果はないのですが、
可変長のFFT乗算を用いる場合、演算量は２の冪数桁単位に増加するので、
Totalの計算量としては、こちらの方が有利です。

分割数は、利用するBigint(多倍長整数)クラスで最速の性能が出るように調節する必要があります。
これについては、「高橋 大介, 金田 康正. "多倍長平方根の高速計算法". 情報処理学会研究報告 95-HPC-58, pp.51-56.」が詳しいです。
文献に従えば、初期値に倍精度結果を使う・３〜４分割程度・FFT再利用、を行えば最大限の高速化を望めますが、このサンプルコードではそこまで対応できていません。

また、この方法は平方根に限らず、Newton法を使った様々な計算に応用できます。

整数平方根での詳細解説

整数除算での詳細解説

このアルゴリズムを利用した16進・２進・10進電卓です。
整数だけですがおよそ100万桁まで計算できます。


(注)ここでのBigintクラスは、多倍長整数を格納する仮想のクラスであり、実在するものではありません。

Tags: プログラムメモ
続きを読む＞＞
author : HUNDREDSOFT | - | -

モンゴメリ冪剰余 サンプルコード

| プログラムメモ | 2010/12/26 Sun 11:11
前回、モンゴメリ乗算 サンプルコードを書きましたが、よく使われる冪剰余についてもC++サンプルを載せておきます。 


// モンゴメリ冪剰余演算 (a^b) mod nを求める
// gcd(n,r) == 1 となること。(nは奇数であれば良い)
//
#ifndef MONTGOMERY_INCLUDED
#define MONTGOMERY_INCLUDED

#include "Bigint.h"

class montgomery{
public:
    montgomery(const Bigint& n) : m_n(), m_r2(), m_nb(0), m_n2(){
        m_n = n;
        // n を法とする
        // Bit長を調べる(nb)
        m_nb = n.length();

        // r: nより大きな２のべき
        Bigint r(1);
        r <<= m_nb;

        // r2 = r^2 mod n
        m_r2 = r * r;
        m_r2 %= n;

        // n*n2 ≡ -1 mod r
        m_n2 = 0; /* 求めるN' */
        Bigint t = 0;
        Bigint vi = 1;
        int ni = m_nb;
        while (ni--){ /* Rのトップビットを除いたビット数分繰り返す */
            if ((t & 1) == 0){
                /* ゼロになっているビットがあったら、N'のその部分を1にする（NはRと互いに素なので必ず奇数）*/
                t += n; /* 掛け算だが、二進数一桁の掛け算なので実質は足し算 */
                m_n2 += vi; /* N'のその部分を1にする */
            }
            t >>= 1; /* 必ず端数が出るが切り捨てる */
            vi <<= 1; /* Rは2の冪なので、絶対端数は出ない */
        }
    }

    ~montgomery(){}

    Bigint reduction(const Bigint& t){
        // tのモンゴメリリダクション
        Bigint c = t * m_n2;
        c.hibitmask(m_nb); // mod Rの意味,(m_nb)bit以上を0クリア
        c *= m_n;
        c += t;
        c >>= m_nb; // 1/r の意味
        if (c >= m_n) c -= m_n;
        return c;
    }

    // 冪剰余 a^b mod n, バイナリ法の下位桁から計算
    Bigint exp(const Bigint& a, const Bigint& b){
        Bigint p = reduction(a * m_r2);
        Bigint x = reduction(m_r2);
        Bigint y(b);
        while (y > 0){
            if (y & 1){
                x = reduction(x * p);
            }
            p = reduction(p * p);
            y >>= 1;
        }
        return reduction(x);
    }

protected:
    montgomery() : m_n(), m_r2(), m_nb(0), m_n2(){}

private:
    Bigint m_n;   // 法(m_n)
    Bigint m_r2;  // rを(m_n)より大きな２のべきとして、(m_r2)=(r)^2 mod n
    Bigint m_n2;  // (m_n)*(m_n2) ≡ -1 mod r を満たす(m_n2)
    int    m_nb;  // (m_n)の有効Bit数, 例えば(m_n)=5 のとき(m_nb)=3
};

#endif


wikipediaの説明にあるようにバイナリ法を用いています。
初めに法(n)とするコンストラクタを呼び出します。
冪剰余は、exp()メソッドで行います。


このアルゴリズムを利用した16進・２進・10進電卓です。
整数だけですがおよそ100万桁まで計算できます。


(注)Bigintクラスは、多倍長整数を格納する仮想のクラスであり、実在するものではありません。

Tags: プログラムメモ
author : HUNDREDSOFT | - | -

VC++で拡張倍精度をどうしても使いたい

| プログラムメモ | 2010/12/19 Sun 08:25
VC++で拡張倍精度浮動小数点形式をどうしても使いたい

Visual Stdioでは拡張倍精度浮動小数(80bit)が使えません。
long double と書いても、64bit精度になります。
通常、仮数部53bit(暗黙Bit含め)で不足することはないでしょうが、
どうしても仮数部64bitが必要という場合のサンプルクラスです。
様々なコンストラクタ、四則演算・比較の他、FPUでサポートしている三角関数や平方根を入れました。


/* (c) Hundredsoft Corporation 2010 All rights reserved.

	Ldouble.h -  VC 拡張倍精度クラス

*/

#ifndef LDOUBLE_INCLUDED
#define LDOUBLE_INCLUDED

#include <memory.h>

typedef struct {
	unsigned char d[12];
} LONG_DOUBLE_STR, *LPLONG_DOUBLE_STR;

void Ldouble_init()
{
	unsigned short fcw=0x033f;	// round to nearest, 64 bits, exceptions disabled
	__asm {
		fldcw     [fcw]
	}
}

class Ldouble
{
public:
	Ldouble() {
		memset(&m_data, 0, sizeof(m_data));
	}
	virtual ~Ldouble(){}

	Ldouble(const Ldouble& v) {
		if (this != &v){
			memcpy(&m_data, &v.m_data, sizeof(m_data));
		}
	}
	Ldouble(const __int64& v) {
		LPLONG_DOUBLE_STR m = &m_data;
		__asm {
			mov  eax, [v]
			fild qword ptr [eax]
			mov  eax, [m]
			fstp tbyte ptr [eax]
		}
	}
	Ldouble(const unsigned __int64& v) {
		__int64 t = 0x8000000000000000;
		int i;
		for (i=63; i>=0; i--, t>>=1){
			if (v & t) break;
		}
		if (i < 0){
			memset(&m_data, 0, sizeof(m_data));
		}else{
			t = v;
			t <<= 63-i; // 仮数部
			int ex = i + 16383; // 指数部
			memcpy(&m_data, &t, sizeof(__int64));
			*(int*)&m_data.d[8] = ex;
		}
	}
	Ldouble(const int& v) {
		LPLONG_DOUBLE_STR m = &m_data;
		__asm {
			mov  eax, [v]
			fild dword ptr [eax]
			mov  eax, [m]
			fstp tbyte ptr [eax]
		}
	}
	Ldouble(const unsigned int& v) {
		__int64 i = v;
		LPLONG_DOUBLE_STR m = &m_data;
		__asm {
			lea  eax, [i]
			fild qword ptr [eax]
			mov  eax, [m]
			fstp tbyte ptr [eax]
		}
	}
	Ldouble(const short& v) {
		LPLONG_DOUBLE_STR m = &m_data;
		int i = v;
		__asm {
			lea  eax, [i]
			fild dword ptr [eax]
			mov  eax, [m]
			fstp tbyte ptr [eax]
		}
	}
	Ldouble(const unsigned short& v) {
		LPLONG_DOUBLE_STR m = &m_data;
		int i = v;
		__asm {
			lea  eax, [i]
			fild dword ptr [eax]
			mov  eax, [m]
			fstp tbyte ptr [eax]
		}
	}
	Ldouble(const char& v) {
		LPLONG_DOUBLE_STR m = &m_data;
		int i = v;
		__asm {
			lea  eax, [i]
			fild dword ptr [eax]
			mov  eax, [m]
			fstp tbyte ptr [eax]
		}
	}
	Ldouble(const unsigned char& v) {
		LPLONG_DOUBLE_STR m = &m_data;
		int i = v;
		__asm {
			lea  eax, [i]
			fild dword ptr [eax]
			mov  eax, [m]
			fstp tbyte ptr [eax]
		}
	}
	Ldouble(const double& v) {
		LPLONG_DOUBLE_STR m = &m_data;
		__asm {
			mov  eax, [v]
			fld  qword ptr [eax]
			mov  eax, [m]
			fstp tbyte ptr [eax]
		}
	}
	Ldouble& operator=(const Ldouble& v){
		if (this != &v){
			memcpy(&m_data, &v.m_data, sizeof(m_data));
		}
		return *this;
	}

	// Unary -								
	Ldouble operator-(){
		Ldouble v = *this;
		v.m_data.d[9] ^= 0x80;
		return v;
	}
	Ldouble operator++(){
		LPLONG_DOUBLE_STR m = &m_data;
		__asm {
			mov  eax, [m]
			fld  tbyte ptr [eax]
			fld1
			fadd
			mov  eax, [m]
			fstp tbyte ptr [eax]
		}
		return *this;
	}
	Ldouble operator++(int){
		Ldouble v = *this;
		LPLONG_DOUBLE_STR m = &m_data;
		__asm {
			mov  eax, [m]
			fld  tbyte ptr [eax]
			fld1
			fadd
			mov  eax, [m]
			fstp tbyte ptr [eax]
		}
		return v;
	}
	Ldouble operator--(){
		LPLONG_DOUBLE_STR m = &m_data;
		__asm {
			mov  eax, [m]
			fld  tbyte ptr [eax]
			fld1
			fsub
			mov  eax, [m]
			fstp tbyte ptr [eax]
		}
		return *this;
	}
	Ldouble operator--(int){
		Ldouble v = *this;
		LPLONG_DOUBLE_STR m = &m_data;
		__asm {
			mov  eax, [m]
			fld  tbyte ptr [eax]
			fld1
			fsub
			mov  eax, [m]
			fstp tbyte ptr [eax]
		}
		return v;
	}

	// Binary operator +=, -=, *=, /=
	Ldouble& operator+=(const Ldouble& v){
		LPLONG_DOUBLE_STR m = &m_data;
		LPLONG_DOUBLE_STR mv = (LPLONG_DOUBLE_STR)&v.m_data;
		__asm {
			mov  eax, [m]
			fld  tbyte ptr [eax]
			mov  eax, [mv]
			fld  tbyte ptr [eax]
			fadd
			mov  eax, [m]
			fstp tbyte ptr [eax]
		}
		return *this;
	}
	Ldouble& operator-=(const Ldouble& v){
		LPLONG_DOUBLE_STR m = &m_data;
		LPLONG_DOUBLE_STR mv = (LPLONG_DOUBLE_STR)&v.m_data;
		__asm {
			mov  eax, [m]
			fld  tbyte ptr [eax]
			mov  eax, [mv]
			fld  tbyte ptr [eax]
			fsub
			mov  eax, [m]
			fstp tbyte ptr [eax]
		}
		return *this;
	}
	Ldouble& operator*=(const Ldouble& v){
		LPLONG_DOUBLE_STR m = &m_data;
		LPLONG_DOUBLE_STR mv = (LPLONG_DOUBLE_STR)&v.m_data;
		__asm {
			mov  eax, [m]
			fld  tbyte ptr [eax]
			mov  eax, [mv]
			fld  tbyte ptr [eax]
			fmul
			mov  eax, [m]
			fstp tbyte ptr [eax]
		}
		return *this;
	}
	Ldouble& operator/=(const Ldouble& v){
		LPLONG_DOUBLE_STR m = &m_data;
		LPLONG_DOUBLE_STR mv = (LPLONG_DOUBLE_STR)&v.m_data;
		__asm {
			mov  eax, [m]
			fld  tbyte ptr [eax]
			mov  eax, [mv]
			fld  tbyte ptr [eax]
			fdiv
			mov  eax, [m]
			fstp tbyte ptr [eax]
		}
		return *this;
	}

	// compare
	// this <  v: -
	// this == v: 0
	// this >  v: +
	int comp(const Ldouble& v) const{
		LPLONG_DOUBLE_STR m = (LPLONG_DOUBLE_STR)&m_data;
		LPLONG_DOUBLE_STR mv = (LPLONG_DOUBLE_STR)&v.m_data;
		short sts;
		__asm {
			mov  eax, [m]
			fld  tbyte ptr [eax]
			mov  eax, [mv]
			fld  tbyte ptr [eax]
			lea  eax, [sts]
			fcompp
			fnstsw word ptr [eax]
		}
		return ((sts & 0x4100) - 1) ^ 0x3fff;
	}
	bool operator==(const Ldouble& v) const{
		return (comp(v) == 0);
	}
	bool operator!=(const Ldouble& v) const{
		return (comp(v) != 0);
	}
	bool operator<(const Ldouble& v) const{
		return (comp(v) < 0);
	}
	bool operator<=(const Ldouble& v) const{
		return (comp(v) <= 0);
	}
	bool operator>(const Ldouble& v) const{
		return (comp(v) > 0);
	}
	bool operator>=(const Ldouble& v) const{
		return (comp(v) >= 0);
	}

	// cast
	operator double() {
		LPLONG_DOUBLE_STR m = (LPLONG_DOUBLE_STR)&m_data;
		double t;
		__asm {
			mov  eax, [m]
			fld  tbyte ptr [eax]
			fstp qword ptr [t]
		}
		return t;
	}
	operator int() {
		LPLONG_DOUBLE_STR m = (LPLONG_DOUBLE_STR)&m_data;
		double t;
		__asm {
			mov  eax, [m]
			fld  tbyte ptr [eax]
			fstp qword ptr [t]
		}
		return (int)t;
	}

	//friend Binary operator	 +,-,*,/
	friend  const Ldouble operator+(const  Ldouble&a, const Ldouble&b){
		Ldouble p(a);
		p += b;
		return p;
	}
	friend  const Ldouble operator-(const  Ldouble&a, const Ldouble&b){
		Ldouble p(a);
		p -= b;
		return p;
	}
	friend  const Ldouble operator*(const  Ldouble&a, const Ldouble&b){
		Ldouble p(a);
		p *= b;
		return p;
	}
	friend  const Ldouble operator/(const  Ldouble&a, const Ldouble&b){
		Ldouble p(a);
		p /= b;
		return p;
	}


	LONG_DOUBLE_STR	m_data;
};

namespace LD {
	Ldouble abs(const Ldouble& s){
		Ldouble v(s);
		v.m_data.d[9] &= 0x7f;
		return v;
	}

	Ldouble sin(const Ldouble& s){
		Ldouble v(s);
		LPLONG_DOUBLE_STR m = &v.m_data;
		__asm {
			mov  eax, [m]
			fld  tbyte ptr [eax]
			fsin
			fstp tbyte ptr [eax]
		}
		return v;
	}

	Ldouble cos(const Ldouble& s){
		Ldouble v(s);
		LPLONG_DOUBLE_STR m = &v.m_data;
		__asm {
			mov  eax, [m]
			fld  tbyte ptr [eax]
			fcos
			fstp tbyte ptr [eax]
		}
		return v;
	}

	Ldouble tan(const Ldouble& s){
		Ldouble v(s);
		double t;
		LPLONG_DOUBLE_STR m = &v.m_data;
		__asm {
			mov  eax, [m]
			fld  tbyte ptr [eax]
			fptan
			lea  ebx, [t]
			fstp qword ptr [ebx]
			fstp tbyte ptr [eax]
		}
		return v;
	}

	Ldouble atan(const Ldouble& s){
		Ldouble v(s);
		LPLONG_DOUBLE_STR m = &v.m_data;
		__asm {
			mov  eax, [m]
			fld  tbyte ptr [eax]
			fld1
			fpatan
			fstp tbyte ptr [eax]
		}
		return v;
	}

	Ldouble atan2(const Ldouble& x, const Ldouble& y){
		Ldouble v(y);
		LPLONG_DOUBLE_STR xm = (LPLONG_DOUBLE_STR)&x.m_data;
		LPLONG_DOUBLE_STR vm = &v.m_data;
		__asm {
			mov  eax, [xm]
			fld  tbyte ptr [eax]
			mov  eax, [vm]
			fld  tbyte ptr [eax]
			fpatan
			fstp tbyte ptr [eax]
		}
		return v;
	}

	Ldouble sqrt(const Ldouble& s){
		Ldouble v(s);
		LPLONG_DOUBLE_STR m = &v.m_data;
		__asm {
			mov  eax, [m]
			fld  tbyte ptr [eax]
			fsqrt
			fstp tbyte ptr [eax]
		}
		return v;
	}

	Ldouble pi(){
		Ldouble v;
		LPLONG_DOUBLE_STR m = &v.m_data;
		__asm {
			mov  eax, [m]
			fldpi
			fstp tbyte ptr [eax]
		}
		return v;
	}

	Ldouble log(const Ldouble& s){
		Ldouble v(s);
		LPLONG_DOUBLE_STR m = &v.m_data;
		__asm {
			fldln2
			mov  eax, [m]
			fld  tbyte ptr [eax]
			fyl2x
			fstp tbyte ptr [eax]
		}
		return v;
	}

	Ldouble log10(const Ldouble& s){
		Ldouble v(s);
		LPLONG_DOUBLE_STR m = &v.m_data;
		__asm {
			fldlg2
			mov  eax, [m]
			fld  tbyte ptr [eax]
			fyl2x
			fstp tbyte ptr [eax]
		}
		return v;
	}

	Ldouble exp(const Ldouble& s){
		Ldouble v(s);
		LPLONG_DOUBLE_STR m = &v.m_data;
		double half = 0.5;
		bool sign = ((m->d[9] & 0x80) != 0);
		m->d[9] &= 0x7f;

		__asm {
			fldl2e                 // log2e
			mov  eax, [m]
			fld  tbyte ptr [eax]   // log2e, v
			fmul                   // log2e*v
			fld  st(0)             // log2e*v, log2e*v
			lea  ebx, [half]
			fld  qword ptr [ebx]   // log2e*v, log2e*v, 0.5
			fsub                   // log2e*v, log2e*v-0.5
			frndint                // log2e*v, (int)(log2e*v-0.5)
			fxch                   // (int)(log2e*v-0.5), log2e*v
			fsub  st(0), st(1)     // (int)(log2e*v-0.5), log2e*v - (int)(log2e*v-0.5)
			f2xm1                  // (int)(log2e*v-0.5), 2^(log2e*v - (int)(log2e*v-0.5))-1
			fld1                   // (int)(log2e*v-0.5), 2^(log2e*v - (int)(log2e*v-0.5))-1, 1
			fadd                   // (int)(log2e*v-0.5), 2^(log2e*v - (int)(log2e*v-0.5))
			fscale                 // (int)(log2e*v-0.5), 2^(log2e*v - (int)(log2e*v-0.5)) * 2^(int)(log2e*v-0.5)
			fstp tbyte ptr [eax]   // (int)(log2e*v-0.5)
			fstp qword ptr [ebx]   // 
		}

		if (sign){
			__asm {
				fld1
				mov  eax, [m]
				fld  tbyte ptr [eax]
				fdiv
				fstp tbyte ptr [eax]
			}
		}
		return v;
	}

	Ldouble pow(const Ldouble& x, const Ldouble& y){
		Ldouble u(x),v(y);
		LPLONG_DOUBLE_STR mx = &u.m_data;
		LPLONG_DOUBLE_STR m = &v.m_data;
		double half = 0.5;
		bool xsign = ((mx->d[9] & 0x80) != 0);
		bool ysign = ((m->d[9] & 0x80) != 0);
		if (xsign){
			int ny = (int)v;
			Ldouble t = ny;
			if (t != y){
				throw "Invalid value.";
			}
			xsign = (ny & 1) ? true : false;
		}
		mx->d[9] &= 0x7f;
		m->d[9] &= 0x7f;

		__asm {
			fld1                   // 1
			mov  eax, [mx]
			fld  tbyte ptr [eax]   // 1, x
			fyl2x                  // log2x
			mov  eax, [m]
			fld  tbyte ptr [eax]   // log2x, v
			fmul                   // log2x*v
			fld  st(0)             // log2x*v, log2x*v
			lea  ebx, [half]
			fld  qword ptr [ebx]   // log2x*v, log2x*v, 0.5
			fsub                   // log2x*v, log2x*v-0.5
			frndint                // log2x*v, (int)(log2x*v-0.5)
			fxch                   // (int)(log2x*v-0.5), log2x*v
			fsub  st(0), st(1)     // (int)(log2x*v-0.5), log2x*v - (int)(log2x*v-0.5)
			f2xm1                  // (int)(log2x*v-0.5), 2^(log2x*v - (int)(log2x*v-0.5))-1
			fld1                   // (int)(log2x*v-0.5), 2^(log2x*v - (int)(log2x*v-0.5))-1, 1
			fadd                   // (int)(log2x*v-0.5), 2^(log2x*v - (int)(log2x*v-0.5))
			fscale                 // (int)(log2x*v-0.5), 2^(log2x*v - (int)(log2x*v-0.5)) * 2^(int)(log2x*v-0.5)
			fstp tbyte ptr [eax]   // (int)(log2x*v-0.5)
			fstp qword ptr [ebx]   // 
		}

		if (ysign){
			__asm {
				fld1
				mov  eax, [m]
				fld  tbyte ptr [eax]
				fdiv
				fstp tbyte ptr [eax]
			}
		}
		if (xsign){
			m->d[9] |= 0x80;
		}
		return v;
	}
} // end of namespace

#endif

Ldoubleクラスを使う前に、Ldouble_init()を呼び出して、FPUに拡張倍精度を指示する必要があります。
このコードは動作しますが、あまり速くありません。
実行速度についてはコンパイラが拡張倍精度をサポートしてくれないと、どうにもなりません。
実行速度をそれほど求めていないが、仮数部64bitの精度が必要であるか、
最大値が10進で300桁以上になる場合（拡張倍精度ではおよそ5000桁まで可能）にのみ有効です。
Win32アプリで拡張倍精度計算を行い、同時に実行速度も求めるなら、VC以外のコンパイラ(C++ Builder等)を利用した方が良いでしょう。



Tags: プログラムメモ
author : HUNDREDSOFT | - | -

多倍長サンプルアプリ

| プログラムメモ | 2010/12/05 Sun 18:33
多桁の除算や平方根計算を実装したサンプルの16進・２進・10進電卓です。
整数だけですがおよそ100万桁まで計算できます。




乗算は、実数FFTとKaratsuba法を組み合わせています。
この内、実数FFTには、京大の大浦博士のDFTコード(fftsg.c)の一部をFPU化して使っています。
VC++だと拡張倍精度の long double が使えないため、できるだけFPUを使って精度を高めています。
精度が上がると一度に乗算できる桁数が上がり、応答性が向上します。
このときの精度については 「(*1)FFT多倍長乗算器のVLSI設計」 を参考に、最大値どうしの乗算が最大の誤差を与えると仮定し、１度に行えるFFT乗算桁数を決定しています。
また、２進→１０進変換にもFFT乗算が必要になります。
これには神奈川工科大学 平山博士の「(*2)分割統治法による多倍長演算の高速化」を参考にしました。


(*1) 矢崎俊志, 阿部公輝: 「FFT 多倍長乗算器のVLSI 設計」, 日本応用数理学会論文誌, Vol.15, No.3, pp.385-401 (2005-9)

(*2) 平山弘: 「分割統治法による多倍長演算の高速化 (数式処理における理論と応用の研究)」, 数理解析研究所講究録 (2000), 1138: 247-255




Tags: プログラムメモ
author : HUNDREDSOFT | - | -
1/2 >>
| TOP |