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より大きな２のべき
        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: プログラムメモ

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: プログラムメモ

多桁の除算や平方根計算を実装したサンプルの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: プログラムメモ