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

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


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

例えば、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


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

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

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

‐絨矛紊伐式矛綟瓜里魍櫃韻栃造戮泙后(→2024)
下3桁を足します(0+2+4→6)
上3桁を足します(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など、同じ数字が並んでいる場合はきイいらなくなるので簡単ですが、最後の例などカラツバ暗算で解くのは減算による桁下げがあるので面倒です。
また、小さな桁数であれば人間にとっては、加算より乗算の方が早いので、
本来のカラツバ法の効果は、全くと言っていい程ありません。(注)

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


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



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

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

Windows7で簡単にマルチタッチを試す方法として、 "Multi-Touch Vista" があります。
名前の通り、それほど新しいソフトではないのでご存じの方もいるでしょう。

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

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を表示する方法は、
WindowsOSのバージョンによっても変わります。

・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 | - | -

FaceBookのメモ、サンプルアプリ

タイムラインから戻す。(ウォールに)

2012/3/24からFaceBookのウォールが廃止になり、タイムラインになる。
厄介なのは、一見、試験運用中と思わせるタイムラインのプレビューをつついてしまうと、2度と、ウォールに戻れなくなることだ。
これは、IE8以降が該当するらしい。

現状の回避策としては、IE起動後、F12でIEの開発Viewを表示して、
メニューの「ブラウザモード」をIE7に変えるしかないらしい。

3/24にははっきりしますが
新インタフェースはIE8以降しか対応しないのだろうか。
MSも、3月からのWindowsUpdateで無条件にIE8化をするらしいので、
これに併せたということだろうか。


-------
Win7+IE9で信頼済みサイトに設置した「いいねボタン」がうまく機能しない。

一般的には
http://developers.facebook.com/docs/reference/plugins/like/
で、Likeボタンのひな型を作成して、HPに張り付けます。

ネットワーク条件は種々様々であるので一概に言えませんが、
Win7+IE9の組み合わせで、信頼済みサイトに設置した「いいねボタン」を押すと
http://www.facebook.com/connect/connect_to_external_page_widget_loggedin.php?social_plugin=like&external_page_url=<指定のurl>#_=_
のブランクページが開いてしまう。
通常なら、コメント入力が現れるはずだ。

駄目な場合には上で書いた、F12でIE7互換に変えても失敗する。
http://developers.facebook.com/tools/lint/
でも特に問題は検出されない。

ベータ版のWin8CP+IE10だと少し複雑で
・メトロ版・・・ブランクページが開く。
・DeskTop版・・・ブランクは出ず、いいねも押せるが、コメントが開かない。

少なくとも、イントラ環境下では「いいね」を使わない方が良さそうだ。


-------
FaceBookサンプルアプリ javaScript SDKでuser情報の表示

(2012.8)現在、この方法でuser情報は表示できなくなりました。

ログオンしていることが条件になるが、
ログオンしていなくても、「いいね」を押すとログオンが要求される。
この時に再度、FB.init()を呼び出してuser情報を取得できた。

いずれにしても、ユーザーに気づかれずに、
FaceBookのuser情報が取得できてしまうことは、とても危険だ。

Deveropper登録してアプリIDを取得すれば、FaceBookのSDKを使ったページが好きなHP上に構築できる。
FaceBookキャンバスは、https:のアドレスが必要だが、
FaceBookキャンバスを利用しないのであれば、どこにでも置ける。
つまり、野良アプリである。
忘れそうなのでメモ。

一般的な「いいねボタン」を配置する、次のようなhtmlを用意する。


<!DOCTYPE html>

<!--- xmlnsは、IE不具合のおまじない、ほとんど効果はない --->
<html xmlns:fb="http://www.facebook.com/2008/fbml">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8">

<!--- ここのメタ部は、いいねボタンでのみ必要なものです。 --->
<meta property="og:title" content="xxxxxx" />
<meta property="og:type" content="xxxxx" />
<meta property="og:url" content="http://www.xxxxxxxxxx.xx.xx/" />
<meta property="og:image" content="http://www.xxxxxxxxx.xx.xx/xxxx.gif" />
<meta property="og:site_name" content="xxxxxx" />
<meta property="fb:admins" content="nnnnnnnnnnnnnnn" />
<!--- いいねメタ ここまで。 --->

</head>
<body>

<!--- FaceBookアプリ --->
<div id="fb-root"></div>
<script script type="text/javascript" src="http://www.xxxxxxxxxxx.xx.xx/core.js"></script>
<!--- FaceBookアプリ ここまで。 --->

<!--- ようこそ:xxさん --->
<div id="username"></div>

<!--- ユーザー情報表示部 --->
<div id="userinf"></div>

<!--- 友人情報表示部 --->
<div id="friendinfo"></div>
<br />

<!--- いいねボタン --->
<div id="iinebutton"></div>

</body>
</html>



channel.htmlは、
https://developers.facebook.com/docs/reference/javascript/
にある通り、

<script src="//connect.facebook.net/en_US/all.js"></script>
の一行だけ。

続いて、FaceBookアプリのcore.jsを次のように書く。


// Copyright(C) Hundredsoft Corporation 2012 All rights reserved. 
//
//  FaceBookアクセスサンプル (javascript SDK)  core.js
//
//   2012.3 時点で意味があって、かつBasic Permissionで引けそうな項目を並べた
//

var m_login = false;

window.fbAsyncInit = function() {
   FB.init({
      appId      : 'xxxxxxxxxxxxxxx', // App ID
      channelUrl : 'channel.html', // Channel File
      status     : true, // check login status
      cookie     : true, // enable cookies to allow the server to access the session
      xfbml      : true  // parse XFBML
   });

   // Additional initialization code here

   // Asyncの場合、いいねはSDK読み込み後に作らないと日本語にならない
   document.getElementById("iinebutton").innerHTML = '<div class="fb-like" data-href="http://www.xxxxxxxxxxx.xx.xx/" data-send="true" data-width="450" data-show-faces="true"></div>';

   FB.getLoginStatus(function(response){
      if (response.status === 'connected') {
         // Logged in
         m_login = true;
         showUser();
      }
      else{
         //本来はここでFB.loginでログインWindowを出させる。
         //今回はユーザの別操作でログインさせるため、無処理。
      }
   });

   // いいねを押したとき呼ばれる
   FB.Event.subscribe('edge.create', function (response) {
      if (!m_login){
         //ちょっと無理させます。
         FB.init({
            appId      : 'xxxxxxxxxxxxxxx', // App ID
            channelUrl : 'channel.html', // Channel File
            status     : true, // check login status
            cookie     : true, // enable cookies to allow the server to access the session
            xfbml      : true  // parse XFBML
         });

         FB.getLoginStatus(function(response){
            if (response.status === 'connected') {
               // Logged in
               m_login = true;
               showUser();
            }
         });
      }
   });
};

function showUser(){
   FB.api('/me', function(response){
      //日本語名称は,userテーブルではなく,profile。しかし権限追加が必要。
      document.getElementById("username").innerHTML = '<h2>ようこそ:' + response.name + ' さん</h2>';

      var query = FB.Data.query('select uid, name, profile_url, sex, email, current_location, hometown_location, work_history, education_history, hs_info, pic_small from user where uid={0}', response.id);
      query.wait(function(rows) {
         var name_e = (rows[0].name!= undefined ?        rows[0].name : "");
         var url =    (rows[0].profile_url!= undefined ? rows[0].profile_url : "");
         var email =  (rows[0].email!= undefined ?       rows[0].email : "");
         var pic_e =  (rows[0].pic_small!= undefined ?   rows[0].pic_small : "");

         var current_loc = "";
         if (rows[0].current_location!= undefined){
             current_loc = (rows[0].current_location.city!= undefined ? rows[0].current_location.city : "") + "," +
                           (rows[0].current_location.state!= undefined ? rows[0].current_location.state : "") + "," +
                           (rows[0].current_location.country!= undefined ? rows[0].current_location.country : "");
         }
         var hometown = "";
         if (rows[0].hometown_location!= undefined){
             hometown = (rows[0].hometown_location.city!= undefined ? rows[0].hometown_location.city : "") + "," +
                        (rows[0].hometown_location.state!= undefined ? rows[0].hometown_location.state : "") + "," +
                        (rows[0].hometown_location.country!= undefined ? rows[0].hometown_location.country : "");
         }
         var highsc = "";
         if (rows[0].hs_info!= undefined){
             highsc = (rows[0].hs_info.hs1_name!= undefined ? rows[0].hs_info.hs1_name : "");
         }
         var edname = "";
         if (rows[0].education_history!= undefined){
             if (rows[0].education_history[0]!= undefined){
                edname = (rows[0].education_history[0].name!= undefined ? rows[0].education_history[0].name : "");
             }
         }
         var wkname = "";
         if (rows[0].work_history!= undefined){
             if (rows[0].work_history[0]!= undefined){
                wkname = (rows[0].work_history[0].company_name!= undefined ? rows[0].work_history[0].company_name : "");
             }
         }

         var jsex = "";
         if (rows[0].sex!= undefined){
            if (rows[0].sex == "female"){
               jsex = "女性";
            }else{
               jsex = "男性";
            }
         }
         document.getElementById('userinf').innerHTML =  
            '名前: '   +  name_e                                                    + "<br />" +
            '性別: '   +  jsex                                                      + "<br />" +
            'URL : '   +  '<a href="' + url + '" target="_blank">' + url + '</a>'   + "<br />" +
            'Mail: '   +  (email!= undefined ? email : "")                          + "<br />" +
            '住所: '   +  current_loc                                               + "<br />" +
            '田舎: '   +  hometown                                                  + "<br />" +
            '会社: '   +  wkname                                                    + "<br />" +
            '大学: '   +  edname                                                    + "<br />" +
            '高校: '   +  highsc                                                    + "<br />" +
               '<img src="' +  pic_e + '" width="50" />' + "<br /><br /><br />";
      });

      showFriend();
   });
}


function showFriend(){
   FB.api('/me/friends', function(response){
      var users = "";
      var i;
      for(i=0;i<response.data.length;i++){
         var f_id = (response.data[i].id!= undefined ?   response.data[i].id : "");
         var f_name = (response.data[i].name!= undefined ?   response.data[i].name : "");
         users = users + '<a href="http://www.facebook.com/profile.php?id=' + f_id + '" target="_blank"><img src="http://graph.facebook.com/' + f_id + '/picture" width="50" />' + f_name + '</a>&nbsp;&nbsp;&nbsp;&nbsp;';
      }
      document.getElementById("friendinfo").innerHTML = '現在の友人は、' + i + ' 人です。<br />' + users;
   });
}

// Load the SDK Asynchronously
// ja-JPに変えないと,いいねが日本語にならない
//
(function(d){
   var js, id = 'facebook-jssdk', ref = d.getElementsByTagName('script')[0];
   if (d.getElementById(id)) {return;}
   js = d.createElement('script'); js.id = id; js.async = true;
   js.src = "//connect.facebook.net/ja_JP/all.js";
   ref.parentNode.insertBefore(js, ref);
}(document));



この3つのファイルを適当なHPに置けば、FaceBook情報を表示してくれる。

こちらにサンプルを置いておきます。
IEだと「いいね」を押してログオンしても情報が表示されないときがあるので、
その時はリロード(再読み込み)してください。




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

Android着信音のループ再生

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

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

但し、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 | - | -

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

前回、除算 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回目は、2回目の計算結果を初期値として、64+(除数のBit数/2)での計算。
4回目は、3回目の計算結果を初期値として、全桁での計算。
と計算桁を増やしていきます。

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

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

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

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

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

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


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


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

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

前回、平方根 √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回目は、2回目の計算結果を初期値として、更に(n/8)桁増やして計算。
と次々と桁を増やしていきます。

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

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

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

整数平方根での詳細解説

整数除算での詳細解説

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


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


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

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

前回、モンゴメリ乗算 サンプルコードを書きましたが、よく使われる冪剰余についても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より大きな2のべき
        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)より大きな2のべきとして、(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進・2進・10進電卓です。
整数だけですがおよそ100万桁まで計算できます。


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


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

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

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 | - | -

多倍長サンプルアプリ

多桁の除算平方根計算を実装したサンプルの16進・2進・10進電卓です。
整数だけですがおよそ100万桁まで計算できます。




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


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

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





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

ニュートン法を使った2進多桁の整数平方根

多桁の整数値(s)の平方根 √s があったとして
ニュートン法(2次漸化式)を使って(s)の平方根の整数部分を求める方法について考えます。
与える引数(s)に基数のべき乗を掛けることで任意の精度に拡張できます。 __前回同様__ 高次漸化式にも応用できますが、ここでは2次漸化式を使っています。

√s を展開すると除算が必要になるので, (1/√s)を求め、これに(s)を乗じることとします。
(1/√s) の2次漸化式は、

X(i+1) = X(i) * (3 - s * X(i)^2) / 2

と表せます。
これは、

Bigint iSqrt(Bigint s){
  int n = s.length(); // bit数

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

  while(m != x){
    m = x;
    x *= t3 - ((s * x * x) >> n);
    x >>= (n + 1);
  }
  x *= s;         // x = s / √s = √s
  x >>= n;

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


と書けます。
考え方としては、計算各項を (* 2n) することで整数を一時的に固定小数点化するものです。
演算量を減らす変形を重ねているのでわかりにくいかもしれませんが、非常に短いコードです。
初期値については、平方根は (n/2)桁程度になるので、その逆数もまた 2n倍することを考えれば、
(n/2)桁程度であることから決定していますが、初期値が甘いのでループ回数はやや多いです。
ニュートン法を使った除算については、__こちら__をご覧ください。

「Newton法初期値の改良バージョン」はこちら
ニュートン法を使った2進多桁の整数平方根(その2)

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


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

ニュートン法を使った2進多桁の整数除算

除算 x = p / q があったとして
ニュートン法(2次漸化式)を使って(q)の逆数を求め、それに(p)を乗算する方法は多桁演算の定石です。
この場合、多桁の浮動小数点型に変換して演算することが多いでしょう。
ここでは、ニュートン法を整数演算のみで実装する方法について考えます。
高次漸化式にも応用できますが、ここでは2次漸化式を使っています。

1/q の2次漸化式は、

X(i+1) = X(i) * (2 - q * X(i))

と表せます。
ここで被除数のBit数をpb, 除数のBit数をqbとして、n = pb + qb とすれば


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

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

  while(m != x){
    m = x;
    x *= c2 - q * x;
    x >>= n;
  }
  x *= p;        // x = p / q
  x >>= n;

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


と書けます。
考え方としては、計算各項を (* 2n) することで整数を一時的に固定小数点化するものです。
演算量を減らす変形を重ねているのでわかりにくいかもしれませんが、非常に短いコードです。
但し、初期値が甘いのでループ回数はやや多いです。
被除数と除数の桁差が少ない場合は、回復法を使った2進多桁の除算が有利です。


「Newton法初期値の改良バージョン」はこちら
ニュートン法を使った2進多桁の整数除算(その2)

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


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

モンゴメリ乗算 サンプルコード

ウィキぺディアでモンゴメリ乗算を検索すると、わかり易く詳しい説明が載っている。
残念ながらこの説明全体をプログラミングしたものが見つからなかったのでC++ソースを書いてみた。


// モンゴメリ乗算剰余演算 c = (a*b) mod nを求める
// gcd(n,r) == 1 となること。(nは奇数であれば良い)
//
Bigint Montgomery(const Bigint& a, const Bigint& b, const Bigint& n) const{
  // n を法とする
  // Bit長を調べる(nb)
  int nb = n.length();

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

  // r2 = r^2 mod n
  Bigint r2(r*r);
  r2 %= n;

  // n*n2 ≡ -1 mod r
  Bigint n2 = 0; /* 求めるN' */
  Bigint t = 0;
  Bigint vi = 1;
  int ni = nb;
  while (ni--){ /* Rのトップビットを除いたビット数分繰り返す */
    if ((t & 1) == 0){
      /* ゼロになっているビットがあったら、N'のその部分を1にする(NはRと互いに素なので必ず奇数)*/
      t += n; /* 掛け算だが、二進数一桁の掛け算なので実質は足し算 */
      n2 += vi; /* N'のその部分を1にする */
    }
    t >>= 1; /* 必ず端数が出るが切り捨てる */
    vi <<= 1; /* Rは2の冪なので、絶対端数は出ない */
  }
  // ここまでで、今後計算に必要になる、r2,n,nb,n2が得られる。
  // つまりnを法とするモンゴメリクラスを作成するなら
  // 引数nをコンストラクタとするクラスを作成し、
  // r2,n,nb,n2をメンバ変数とする。


  // aのモンゴメリ表現をam, bのモンゴメリ表現をbmとする
  // t は作業領域
  t = a * r2;
  Bigint am = t * n2;
  am.hibitmask(nb); // mod Rの意味,(nb)bit以上を0クリア
  am *= n;
  am += t;
  am >>= nb; // 1/Rの意味
  if (am >= n) am -= n;

  t = b * r2;
  Bigint bm = t * n2;
  bm.hibitmask(nb); // mod Rの意味,(nb)bit以上を0クリア
  bm *= n;
  bm += t;
  bm >>= nb; // 1/Rの意味
  if (bm >= n) bm -= n;

  // cm: am * bm のモンゴメリリダクション
  t = (am * bm);
  Bigint cm = t * n2;
  cm.hibitmask(nb); // mod Rの意味,(nb)bit以上を0クリア
  cm *= n;
  cm += t;
  cm >>= nb; // 1/Rの意味
  if (cm >= n) cm -= n;

  // cmのモンゴメリリダクション
  t = cm;
  Bigint c = t * n2;
  c.hibitmask(nb); // mod Rの意味,(nb)bit以上を0クリア
  c *= n;
  c += t;
  c >>= nb; // 1/Rの意味
  if (c >= n) c -= n;

  return c;
}



ここでは、できるだけwikipediaの説明と合致するように書きましたが、
実利用の際には、r2 = r^2 mod n の計算が重いので、
(これが早ければ、そもそもモンゴメリ演算の必要がない)
nを法とするモンゴメリクラスを作成するのが良いです。
つまりコード中にも書いてあるように、引数nをコンストラクタとしたモンゴメリクラスを作成し、r2,n,nb,n2はメンバ変数にして、モンゴメリリダクション, 乗算剰余演算, 加減算剰余演算, 冪剰余演算等は、モンゴメリクラスのMethodにします。

※ここで書いた乗算は「乗算剰余演算」の項の上段にある方法です。

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


Tags: モンゴメリ冪剰余サンプルコード     プログラムメモ
author : HUNDREDSOFT | - | -

回復法を使った2進多桁の除算(その2)

前回同様、除算 a = x / y があったとして
式を変形して、

a = (x * (2^(2k)-1)/y + 2^k) >> 2k  (k:=被除数のbit数)

とすれば、逆数を求める部分((2^(2k)-1)/y)の回復法もまた特殊な被除数のため、簡単なものに変形できる。

Bigint iDiv(Bigint x, Bigint y){
  int xb = x.length(); // 被除数のbit数
  int yb = y.length(); // 除数のbit数
  int n = xb * 2 - 1; // ループ回数

  Bigint a(0), m(0);
  while (n--){
    m |= 1;
    a <<= 1;
    if (m >= y){
      m -= y;
      a |= 1;
    }
    m <<= 1;
  }
  a *= x;
  m = 1;
  m <<= xb;
  a += m;
  return a >> (xb*2-1);
}

この方法は、被除数bit数の倍のループとなり、これが最大の減算回数になります。
しかし誤差が出ないので乗算は1回だけで済みます。
前回の方法はbit差に依存したループ回数なので多桁にも応用できますが、この方法だと被除数が大きいとループ回数が多くなってしまいます。

桁差が小さい時に有利な方法

ニュートン法を使った2進多桁の整数除算

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


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

回復法を使った2進多桁の除算

除算 a = x / y があったとして
回復法を使うと被除数のbit数分の加算・減算を繰り返し行うことになる
式を変形して、

a = (x * ((2^k)/y) >> k  (k:=被除数のbit数)

とすれば、逆数を求める部分((2^k)/y)の回復法は特殊な被除数のため、簡単なものに変形できる。

Bigint iDiv(Bigint x, Bigint y){
  Bigint a(0), m(1);
  int xb = x.length(); // 被除数のbit数
  int yb = y.length(); // 除数のbit数
  int n = xb - yb + 1; // ループ回数(bit差+1)

  m <<= yb-1;
  while (n--){
    a <<= 1;
    if (m >= y){
      m -= y;
      a |= 1;
    }
    m <<= 1;
  }
  a *= x;
  a >>= (xb - 1);

  // 誤差修正(ここでのaは,正解より1or2小さい場合がある)
  m = y * a;
  while(x > m){
    m += y;
    if (x > m)
      a += 1;
  }
  return a;
}


と書ける。
計算量は最大でもbit差程度の加減算と2回の乗算なので
bit差が数100のオーダーならNewton法を使うより有利かも。

ニュートン法を使った2進多桁の整数除算

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


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

PerformanceCounterを使ったマイクロ秒時計

PerformanceCounterを使ったマイクロ秒時計

Linuxには、gettimeofdayがありますが、
Windowsには、マイクロ秒まで計測可能なAPIは用意されていません。

PerformanceCounterを使えばかなり正確な時間間隔を知ることはできますが、
時計として利用することは簡単ではありません。
ある程度観測をしていると、ansi標準の_ftimeとはズレていきます。
このため、_ftimeをベースとしてミリ秒以下の部分をPerformanceCounterから得ることを考えましたが、秒の変わり目で時計が逆に進む場合があります。
そこで、ftimeとPerformanceCounterから得られる時間差を調節し、
時計が逆に進むことがないようなマイクロ秒時計を作ってみました。

usec.h

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

#ifndef USEC_INCLUDED
#define USEC_INCLUDED
class usec
{
public:
	usec();
	~usec(){}
	__int64 get();
private:
	__int64 microfunc(__int64 a, __int64 b);
	__int64 usec::getstm();
	LARGE_INTEGER m_ti;
	LARGE_INTEGER m_tfi;
	__int64 m_msi;
	__int64 m_lastus;
	bool m_isSuported;
};
#endif



usec.cpp

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

#include <windows.h>
#include <time.h>
#include <sys/timeb.h>
#include "usec.h"

usec::usec()
{
	m_isSuported = false;
	if (QueryPerformanceFrequency(&m_tfi)){
		if (QueryPerformanceCounter(&m_ti)){
			m_isSuported = true;
		}
	}
	m_msi = getstm();
}

__int64 usec::microfunc(__int64 a, __int64 b)
{
	__int64 ah, al, a0, m = 0;
	int i;
	ah = (a >> 44) & 0xffff;
	al = a << 20;
	for (i=0; i<128; i++){
		m <<= 1;
		if (ah < 0) m |= 1;
		ah <<= 1;
		if (al < 0) ah |= 1;
		al <<= 1;
		if (m >= b){
			m -= b;
			al |= 1;
		}
	}
	a0 = al;
	ah = (a0 >> 44) & 0xffff;
	al = a0 << 20;
	m = 0;
	for (i=0; i<128; i++){
		m <<= 1;
		if (ah < 0) m |= 1;
		ah <<= 1;
		if (al < 0) ah |= 1;
		al <<= 1;
		if (m >= 22634874){
			m -= 22634874;
			al |= 1;
		}
	}
	return a0 - al;
}

__int64 usec::getstm()
{
	_timeb  tb;
	_ftime(&tb);
	struct tm* stm = localtime(&tb.time);
	__int64 stms = (stm->tm_hour * 3600000 + stm->tm_min * 60000 + stm->tm_sec * 1000 + tb.millitm);
	return stms * 1000;
}

__int64 usec::get()
{
	__int64 utm;
	if (m_isSuported){
		LARGE_INTEGER t, tf;
		QueryPerformanceFrequency(&tf);
		QueryPerformanceCounter(&t);
		__int64 stms = getstm();
		__int64 diff = t.QuadPart-m_ti.QuadPart;
		if (diff > 0){
			utm = m_msi + microfunc(diff, tf.QuadPart);
			diff = utm - stms;
		}else{
			utm = stms;
			diff = 999999999;
		}
		__int64 absdiff = (diff < 0) ? -diff : diff;
		if (m_tfi.QuadPart != tf.QuadPart || absdiff > 10000){
			m_tfi = tf;
			m_ti = t;
			m_msi = stms;
			if (diff > 0){
				if (utm - m_lastus > 10 )
					utm = (m_lastus + utm) / 2;
				else
					utm = m_lastus + 2;
			}else{
				if (absdiff > 3600000000)
					utm = m_msi;
				else
					utm = (utm + m_msi) / 2;
			}
		}
		m_lastus = utm;
	}else{
		utm = getstm();
	}
	return utm;
}





テスト用コード
main.cpp

#include <tchar.h>
#include <windows.h>
#include <stdio.h>
#include "usec.h"

#include <time.h>
#include <sys/timeb.h>

int main(int argc, char**argv)
{
	usec ut;
	while(1){
		char work[256];
		__int64 utm = ut.get();

		_timeb  tb;
		_ftime(&tb);
		struct tm* stm = localtime(&tb.time);

		sprintf(work, "%02d:%02d:%02d.%06d [%02d.%03d]",
				(int)((utm / 3600000000) % 24),
				(int)((utm / 60000000) % 60),
				(int)((utm / 1000000) % 60),
				(int)(utm % 1000000),
				stm->tm_sec,
				tb.millitm
				);
		printf("%s\n", work);
		Sleep(300);
	}
}
usec::microfuncは、PerformanceCounterから得られる値(LARGE_INTEGER)から
マイクロ秒を得る128bit除算を行っています。

ある時間差(t0〜t1)があったとして
t0,t1はQueryPerformanceCounterによって得られる値とします。
t = t1 - t0、この時間間隔は、周波数:f(QueryPerformanceFrequencyによって得られる)として
t / f で得られますが、
整数演算でマイクロ秒を得たいので
t * 1000000 / f
となります。
このままdoubleで演算しても良いのですが整数演算に拘って、
t * (1024*1024 - 48576) / f
= 1024*1024 * t / f * (1 - 1 / 21.58629776021)
= 1024*1024 * t / f * (1 - 1024*1024 / 22634874)
A = (t << 20) / f
と置けば
≒A - (A << 20) / 22634874 [μSec]となります。
これは、64bitを超える演算になりますのでforで回して計算しています。

テスト用コードでは、usecクラスから得られたマイクロ秒時間と共に、
_ftimeによる「秒、ミリ秒」も表示しています。

usec::getでは、ftimeとPerformanceCounterから得られる時間差を調節しているので、
それほど精度の高い時計ではありませんが、ミリ秒単位の時計だとループ間隔を短くSleep(10)とかにすると、同じ時間が続いてしまいます。
この時計の場合は、ループ間隔(Sleep)を短くしても、同じ時間が並ぶようなことはありません。

※高分解能パフォーマンスカウンタがないPCでは、ミリ秒単位の出力となります。



usec.zipソースはこちら



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

素材8

雪結晶



08_snow1.swf





.


Tags: 素材
author : HUNDREDSOFT | - | -

スパイダーソリティア

とうとう負けた。1764連勝で途切れ。

SpiderSolitaire.gif
author : HUNDREDSOFT | - | -

素材7

リンゴ

07_リンゴ4.pdf





.


Tags: 素材
author : HUNDREDSOFT | - | -

素材5

バラとチョコレート
06_choco5.pdf




.


Tags: 素材
author : HUNDREDSOFT | - | -

素材4

グラスとか

05_グラス6.pdf




.


Tags: 素材
author : HUNDREDSOFT | - | -

素材3

新年

04a_新年.pdf



.


Tags: 素材
author : HUNDREDSOFT | - | -

素材2

たき火。

03_たき火.pdf




.


Tags: 素材
author : HUNDREDSOFT | - | -