2014年11月27日 (木)

LX用64ビット整数の整備

■LX用64ビット整数の整備
通常では32ビットまでしか扱えないHP200LXのプログラミング環境で,64
ビット整数を扱えるようにするため,独自の64ビット整数を定義し,自作
の演算関数を整備した.
 
 int64:64ビットデータを格納する構造体
 
主な関数(仮の名称)
 int64 lxshr(int64 a, int n);    //右シフト a>>n
 int64 lxshl(int64 b, int n);    //左シフト b<<n
 int  lxcmp(int64 x, int64 y);  //大小判定
 int64 lxadd(int64 x, int64 y);  //加算 x+y
 int64 lxsub(int64 x, int64 y);  //減算 x-y
 int64 lxmul(int64 x, int64 y);  //乗算 x*y
 int64 lxdiv(int64 n, int64 d, int64 *r); //除算 n/d, r=余り
 int64 lxdiv0(int64 n, int64 d, int64 *r); //除算 n/d, r=余り
 
演算時間
 lxmul()  :32ビット乗算の約 10倍
 lxdiv()  :32ビット乗算の約 30倍  
 lxdiv0() :32ビット乗算の約300倍 回復型除算
 ※速度を要求される場合には不向き 
 
[HP200LX/LXDIC]

|

2014年11月25日 (火)

Sunrise::仕様確定

■Sunrise::仕様確定
Ver 9.9
 
変更点
・浮動小数点数による計算を中止する.
・事前計算データ表参照による補正計算方式を廃止する.
・計算には,小数部19ビットの固定小数点数を用いる.
 
仕様概要
・対応地域
  東経: -180°~ 180°
  北緯:  -90°~ 90°
  時差:  -12 ~ +12  (日本は+9:JST=UTC+9)
・計算誤差:日出日入時刻:1分以内
      方位,高度 :1度以内
・表示項目
  場所(緯度,経度)
  日付,時刻
  日出時刻,日出方位
  南中時刻,南中高度
  日入時刻,日入方位
  太陽高度,太陽方位(指定した時刻の値)
  日出日入一覧表画面の表示
  グラフ表示(高度の時間変化,高度・方位関係)
・設定ファイルの位置情報登録数:10地点
 
計算
・プログラム内蔵の太陽赤緯/均時差テーブルから得たデータに,暦と
 季節のずれに基づく補正を加える.
・この赤緯/均時差,及び,位置情報を計算式に代入し,時角を求める.
 時角から,日出日入南中時刻を計算する.
・得られた時刻における方位角および高度を計算する.
・三角関数の計算は,小数部19ビットの固定小数点数で行う.
 三角関数値は,内蔵の三角関数値表から,内挿計算により算出する.
・固定小数点数の計算は,全て32ビット整数の演算で行う.
・グラフ用データの計算は,日出から日入までの時間を50等分して,
 各時刻における方位と高度を計算する.
・計算は現行の暦に基づいて行う.
 
コンパイラー
・Borland C++ 3.1(PC上でコンパイル)
 
[HP200LX/LXDIC]

|

2014年11月19日 (水)

Sunrise::固定小数点数の乗算

■固定小数点数の乗算
Sunrise
 
・固定小数点数の乗算も,小数部19ビットの場合は,計算の途中でオーバーフ
 ローが発生する.
   a*b ==> (a*b)>>19
 a=0x80000 , b=0x80000 の場合(実数では a=1.0, b=1.0)
  a*b = 0x40 0000 0000  : 39ビット(32ビットを越える)
  (a*b)>>19  : 20ビット(結果は20ビット以下になる)
 三角関数値の演算は,1.0以下の数値の演算なので,必ず 1.0以下の結果が得
 られる.すなわち,20ビット以下になる.
 
 したがって,64ビット整数をつかわなくても,工夫すれば,32ビット整数で
 高速に計算することができる.
 
 実際の計算は,32ビット整数を16ビット整数に分割して行う.
 //符号なし32ビット整数 a ==> (a1,a2): a1, a2 は32ビット整数型
  a1 = a>>16;         //a を16+16ビットに分割
  a2 = a&0x0000ffff; 
  b1 = b>>16;         //b を16+16ビットに分割
  b2 = b&0x0000ffff; 
 //積(a1, a2)*(b1,b2) = (a1*b1, a2*b1+a1*b2, a2*b2)
  c1 = (a1*b1);       //32--64ビット:c1,c2,c3は32ビット整数型
  c2 = (a2*b1+a1*b2); //16--48ビット:c1,c2,c3はオーバーフローしない※
  c3 = (a2*b2);       // 0--32ビット
 //c1,c2,c3を19ビット右シフト(c1,c2,c3はベースが16ビットずつ異なる)
  c  = (c1<<(32-19))+(c2>>(19-16))+(c3>>19);
 (演算回数) 乗算 4回,加算 3回,シフト 5回
 
※通常は,c2で発生する桁上りの処理が必要になるが,Sunriseの計算では必
 要ない.
 最大値をとる場合 a=0x80000, b=0x80000(==>1.0,1.0)を考える.
   a1=0x0008, a2=0x0000
   b1=0x0008, b2=0x0000
 a1, b1が常に小さいので,c2=a1*b2+a2*b1 が32ビット以上になる事はない.
 a2, b2が最大値0xffffをとる場合でも,a1*b2, a2*b1は0x7fff8:19ビット
 
[HP200LX/LXDIC]
 

|

2014年11月18日 (火)

Sunrise::固定小数点数の除算

■固定小数点数の除算
Sunrise
固定小数点数の計算では,小数部のビット数を大きくするほど,計算誤差を
小さくできる.(Sunriseの計算では19ビット)
しかし,32ビット整数の計算では,小数部を16ビット以上にすると,計算の
途中でオーバーフローが発生する.
 小数部19ビットの場合の除算:
   a/b ==> (a<<19)/b
 a=0x80000 の場合(実数では a=1.0 に相当)
  a ==> a<<19   20ビット ==> 39ビット(32ビットを越える)
  1000 0000 0000 0000 0000 = 0x80000
  0100 0000 0000 0000 0000 0000 0000 0000 0000 0000 = 0x4000000000
したがって,正しい計算結果を得るためには,特別の工夫が必要になる.
 
試した主な方法
1.分子と分母を同じ数で割って,小さくしてから計算する.
  (a<<19)/b ===>  ((a<<19)/128)/(b/128) = (a<<12)/(b>>7)
  計算は高速だが,分母で切捨てたビット数分だけ誤差が拡大する.
  
2.64ビット整数で計算する.
  HP200LXの環境では,64ビット整数は扱えないので,独自の演算関数を自
  作する必要があるが,単純な低速アルゴリズムを用いると,演算速度が
  浮動小数点数よりも数倍遅くなってしまい,まったく使いものにならない.
 
3.1の計算結果を補正する.(採用した方法)
  商の概算値を q0 とする.q0=(a<<12)/(b>>7)
  概算値と真の値の差を b倍して考える.正しい商 q=(a<<19)/b
     d = b*q0-(a<<19) (d>0)
  ここで,b*q0 ,(a<<19) は共に32ビットを越えるが,d が32ビットを越
  えることはない.これは,概算値と真の値の差が小さいことによる.
  したがって,d の値は途中のオーバーフローにもかかわらず,正しい値
  になっている.(計算例参照)
  商の真の値 q は,
    q = q0-d/b    となる.(q<q0)
  (これは,正確な値で,64ビット整数による計算値と一致する)
  
 (まとめ)固定小数点数の除算  a/b (小数部19ビット)
    q0=(a<<12)/(b>>7)
    d = b*q0-(a<<19)  (d>0)
    q = q0-d/b
 (32ビット整数の演算回数)乗算1回,除算2回,減算2回,シフト3回
  
●計算例 0x80000=524288
  a=524012=0x7FEEC,  b=524280=0x7FFF8
  真値  a/b ==>   0x7FEEC*0x80000/0x7FFF8
            = 0x7FEF3 = 524019 ※
  概算値 a/b ==> q0 = (0x7FEEC<<12)/(0X7FFF8>>7)
            = 0x7FEEC000/0x0FFF
            = 0x7FF6B = 524139 
            (120の誤差は無視できない)
  計算値 d=q0*b-(a<<19)
            = 0x7FF6B*0x7FFF8-(0x7FEEC<<19)
            = 0x3FFB1804AB-0x3FF7600000(上位3桁が同じ)
            (オーバーフロー部分は,両方とも3Fで等しい)
            = 0x03B804AB
       d/b    = 0x03B804AB/0x7FFF8
            = 0x77 = 119
      q = q0-d/b = 0x7FF6B-0x0077
            = 0x7FEF4 = 524020 ※
            (誤差:120/524288 ==> 1/524288)
  ※実数の計算では,524019.99591 になるので,どちらも正しい.
   (切上げか切捨ての違い)
         
[HP200LX/LXDIC]

|

2014年11月 7日 (金)

Sunrise::三角関数の整数演算化

■Sunrise::三角関数の整数演算化
三角関数とその演算を完全に整数で行う.
 
・三角関数の計算
三角関数表の参照と内挿計算により行う.
 sin 表:900+2点(long)3608バイト
  定義域:0~9000(角度の100倍)
  値域 :0~524288(関数値の0x80000倍)
 asin表:1024+2点(int)2052バイト
 (補助表)256+2点(int)516バイト※
  定義域:0~524288(関数値の0x80000倍)
  値域 :0~9000(角度の100倍)
 ※変化の大きい部分(1.0付近)は,刻みを小さくする.
 
・四則演算
固定小数を用いて計算する.
このプログラムで扱うのは,基本的に三角関数同士の演算であり,
1.0以下の数値になる.精度を確保するため,小数部のビット数を大
きくとる.
 小数部:19ビット
  1.0 ==> 0x80000,  0.5 ==> 0x40000
 演算の概略(小数部8ビットの場合)
  和:a+b   ==>  a+b
  差:a-b   ==>  a-b
  積:a*b   ==>  (a*b)>>8
  商:a/b   ==>  (a<<8)/b
 小数部のビット数を大きくとる場合,計算途中でオーバーフロー
 が発生するため,実際の計算はもう少し複雑になる.
 
・計算誤差
積と商は,計算のたびに誤差が累積する.
これを避けるため,三角関数の積は,和に変換して計算する.
  con(x)cos(y) = (cos(x+y)+cos(x-y))/2 など
テスト版の誤差(浮動小数版との比較)
 低緯度・中緯度の場合:0.08分(角度換算0.02度)
 高緯度の場合    :0.32分(角度換算0.08度)※
 北極南極近傍       :1.28分(角度換算0.32度)※
 
※高緯度地域の特殊性
 高緯度地域には,特徴的な3つの時期がある.
1.1日中昼の時期(日没なし)
2.1日で昼と夜が交替する時期(日出日入あり)
3.1日中夜の時期(日出なし)
それぞれの時期の境界では,昼の長さが1日で1時間程度変化する
ことがある.この様な場合は,計算誤差が大きくなる.
(1日のずれで,数十分の差がある)
 
[HP200LX/LXDIC]

|

2014年11月 4日 (火)

sprintf()

■sprintf()
 
シスマネ上のプログラムを作成する場合,標準の入出力関数は使えない.
(stdio.hをインクルードすること自体ができない)
しかし,stdio.h で宣言されている関数の内,文字や数字の出力に不可欠
な sprintf() は使える.(sprintfはバッファへの出力のみ)
 
使用手順
・stdio.h をインクルードせず,代替のプロトタイプ宣言をする.
 int sprintf(char *, const char *, ...);
 
・エラー(undefined symbol)回避のダミー関数を置く.
 int _terminate(void){Done=TRUE;return 0;}
 int _checknull(void){return 0;}
 int _restorezero(void){return 0;}
 int _cleanup(void){return 0;}
 int errno = 0;
 _terminate()は Done=TRUE がないと,タスクスイッチに失敗し,ハン
 グアップすることがある.
 
[HP200LX/LXDIC]

|

2014年11月 1日 (土)

三角関数計算の整数演算化

■三角関数計算の整数演算化
 Sunrise
これまでの演算の整数化は,東京の事前計算データをもとに補正を加える
もので,緯度が北緯20~54度に限定されていた.
さらに,この方法では,任意時刻の太陽高度と方位角の計算ができないた
め,グラフ表示の高速化にはならない.
 
Sunrise計算方式の変遷
・浮動小数全域版
 すべて浮動小数点数による三角関数の計算.
 サイズが大きく,動作も遅いが,例外なく全地点の計算ができる.
 (グラフ表示に4秒程度かかる)
 
・国内限定版
 データ表参照と整数演算による補正計算.
 動作は軽快だが,地域が限定される上に,グラフ表示ができない.
 
・ハイブリッド版
 北緯20~54度の計算は表参照により,その他は,浮動小数点数の計算
 により実行する.
 サイズが大きく,グラフ表示が遅い.(2秒程度)
 
結局のところ,全域にわたって高速に計算し,グラフ表示も高速化にする
ことはできなかった.
今後は,根本的な考え方を見直す.すなわち,三角関数の計算そのものを
整数演算により高速に実行する事を考える.
 
・整数化全域版(テスト中)
 全計算を整数化して,一覧表のスクロールやグラフ表示も高速化する.
 目標は,グラフ表示時間:0.5秒, 計算誤差:0.3分,0.3度以下
 
●三角関数計算の整数演算化
・三角関数のデータ表(整数)を作成する.
  sin, asin, acos (引数の刻みは,求める精度に依存)
  cos ,tan, atan2などは,上記3つから計算できる.
・三角関数値の計算は,表参照と内挿計算.(整数演算)
・小数を含む全ての数値は,整数化して計算を実行する.
 (簡易的な固定小数点を使う)
  
[HP200LX/LXDIC]

|

2014年10月31日 (金)

三角関数の計算回数

■三角関数の計算回数
 Sunrise
 
1.計算式どおりに計算する場合
(赤緯δ,均時差e,北緯φ,東経θ, 時差td,補正定数Ds)
 時刻計算              //三角関数計算回数:6回
  a =  -sin(Ds)-sin(δ)*sin(φ);     //3
  b =  cos(δ)*cos(φ);              //2
  時角 t = 180/M_PI*acos(a/b);       //1
 角度計算               //三角関数計算回数:15回
  高度 h = asin(sin(φ)*sin(δ)+cos(φ)*cos(δ)*cos(t));//6
  方位 A
    sinA = cos(δ)*sin(t)/cos(h);         //3
    cosA = (sin(h)*sin(φ)-sin(d))/cos(h)/cos(φ); //5
    A = atan2(sinA, cosA) + M_PI;          //1
(合計計算回数): 36回
  1地点につき,時刻1回,方位2回(6+15*2=36)
  
2.1度計算した値を保持する場合
 時刻計算              //三角関数計算回数:5回
  (同一地点では固定値のもの:日出日入時刻の値)
  δのsin,cos  ==>  変数sind, cosd に保持  //2回
  φのsin,cos  ==>  変数sinp, cosp に保持  //2回
  時角 t = 180/M_PI*acos(a/b);             //1回
 角度計算              //三角関数計算回数:6回
  (任意時刻の値)
  時角 t のsin,cos  ==>  変数sind, cosd に保持  //2回
  高度 h = asin(sinp*sind+cosp*cosd*cost);      //1回
  高度 h のsin,cos  ==>  変数sinh, cosh に保持  //2回
  方位 A = atan2(sinA, cosA) + M_PI;          //1回
     (sinA = cosd*sint/cosh;      )
     (cosA = (sinh*sinp-sind)/cosh/cosp;)
(合計計算回数):17回
  1地点につき,時刻1回,方位2回(5+6*2=17)
  (南中高度は太陽赤緯と緯度の加減算で計算できる)
  
・一覧表画面の1ページ分の三角関数計算回数
 360回 ==> 170回
 
・グラフ計算の場合の三角関数計算回数
 グラフ計算は,25点の時刻における方位と高度を計算する.
  1の場合:6+15*25 = 381回/1グラフ
  2の場合:5+6*25  = 155回/1グラフ
 381回 ==> 155回
 
いずれの場合も,実行時間が,半分程度に短縮される.
 
[HP200LX/LXDIC]

|

2014年10月20日 (月)

太陽方位角の変化

■太陽方位角の変化
太陽の方位角は,時間の変化に対して一定の割合で変化するものではない.
北半球の夏場は,昼の12時前後に大きく変化し,日出日入前後では変化
が小さくなる.そして,真夜中の変化率も大きくなるが,昼には及ばない.
冬場には,真夜中の変化が大きくなり,昼の変化を上まわる.
南半球では,夏場と冬場が逆になる.
 
したがって,0時の方位角と12時の方位角の差が,180度になるとは限ら
ない.(南中時刻が12時の場合は180度になる)
 
例)東京 12/20
Sr6902
1時間ごとの変化
Sr6901
この例では,12時の方位は,真南(180度)から少し進んでいる.ここから
12時間たつと,夜中の方位角の変化が大きいため,0時の方位角は,昼間よ
少し余分に進む.結果として,0時と12時の間に方位角は180度以上進む.
 
  
東京 6/20
Sr6911
 
明石(南中時刻がほぼ12時) 方位角の差は約180度
Sr6912
 
沖縄  方位差180度からのずれが大きい
Sr6913
(太陽がほぼ真上を通過するため,方位は南中時刻の前後でほぼ180度変化する)
 
これらの結果は,実生活上の経験と乖離しているように感じるかもしれないが,
これは,方位角が太陽の方向を地平面に投影した角度であり,太陽の位置情報
(高度と方位)の一側面でしかないことによる.
実際,太陽が真上近くにある時は,方位角にはあまり意味がない.
 
[HP200LX/LXDIC]

|

2014年10月11日 (土)

整数演算化

■整数演算化
前回の結果により,基準年のデータ表参照により,すべての年の日出日入
データを導出できることが分かった.(整数演算のみ使用)
この表から整数演算のみで,隣接する他の地点のデータを導出でるかを検
討する.
 
(整数演算化)
 想定する範囲:東経 -180°~ 180°
        北緯 20°~ 54°
 対象となる項目:日出時刻,日入時刻,南中時刻,
         日出方位,日入方位,南中高度
         
・日出時刻,日入時刻
 経度の変化(M)は,経度15度=1時間で補正する.
 同一の経度のもとで,緯度を変化させると,日出日入時刻が変化する.
 理科年表の補正法(北緯20~50)では,この変化量を2つの因子に分解し
 ている.
   日出 T = T0 + M - n*N (N,n は表から求める)
   日入 T = T0 + M + n*N 
 ここで,n は昼長の1/2(東京)のみに依存し,N は緯度にのみ依存する.
 ただし,ある程度以上の精度を求めると,nもNと同様に,緯度の両方に依
 存してしまい,2つの因子に分解したことが無意味になる.
 
 コンピュータプログラムとしては,最初から単純に1つの補正値 nN を計
 算する方が効率的で精度が高い.(nNは緯度と昼長に依存する)
 北緯(20~54°)は1度刻み,昼長/2(290~440分)5分刻みに,nN-値を
 計算し参照テーブルを作成する.(35×31×2=2170バイト)
 予備的な計算では,誤差は 0.2 分程度に抑えられた.
 
・南中時刻
 経度の補正のみ.経度15度=1時間.
 
・南中高度
 東京における南中高度を,目的地点の緯度と東京の緯度の差により補正
 する.
 
・日出方位,日入方位,
 昼長が長くなると,その分,日入の方位角は増加し,日出の方位角は減
 少する.(この変化量は,緯度に依存する.)
 北緯(20~54°)は1度刻み,昼長/2(290~440分)5分刻みに,角度
 補正値を計算し,参照テーブルを作成する.(35×31×2=2170バイト)
 予備的な計算では,誤差は 0.2°程度に抑えられた.
 (PCでの計算:浮動小数点を用いた計算数値との差)
 
※整数演算化できた項目
 日出時刻,日入時刻,南中時刻,昼の長さ,
 日出方位,日入方位,南中高度
※整数演算化できていない項目
 任意時刻の太陽方位,太陽高度
 
[HP200LX/LXDIC]

|

«Sunrise::データ表参照による計算