#include <stdio.h>
#include <vector>
#include <vector>
#include <iostream>
#include <iomanip>
#include <stdarg.h>

// modexp3 0857bde36f680cda9535784bc4aec5f4344131071b419f732ac9c74d0e61db49dd958c7344236e0279df009c6e66aec6ba574c2820d4aeb0c4d814c8e184c6ea7e6d8aa3e15d1c251c78c5364ea2b3edb3c19e90739afa765506242e78fcdc71a87efdfe2df6ce6039fc62cb3b360cb77cd5574292282df352886cbc3fcfbff2 10001 F765A3A0C9C291D81A56FE73794A746B8DA23DBE155D0D495B49D581B5C6545F449A10FDF1C26A92FBD1F43A0687044927A6A21B69A73999E6083D03ACDAFFA6409F1BC71D810628F6E18F76231ED6E22D54ED2502E66F8A33D0D5F07B3EB605F7418110E2EF9A5EE77B070F4EADFCF3D70C53E870F29C9D4F229F2CB6C25383
//   = 0001FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF003021300906052B0E03021A050004147020B30D42ED7CDF0E849F30E8B2E137941E9691
// see also 
// c:\local\cvsprj\secphone\trunk\machine\echocancel\fir_fxp\nr_class.cpp
// c:\local\cvsprj\secphone\trunk\machine\echocancel\fir_fxp\nr_class.h

typedef long long carry_t;
typedef unsigned long base_t;

#define BITSINBASE (8*sizeof(base_t))

typedef std::vector<base_t> number_t;
class Bitmask : public std::vector<base_t> {
public:
    Bitmask(size_t bits)
    {
        for (size_t i=0 ; i<bits ; i++)
        {
            push_back(static_cast<base_t>(1<<i));
        }
    }
};
class Bitsmask : public std::vector<base_t> {
public:
    Bitsmask(size_t bits)
    {
        for (size_t i=0 ; i<bits ; i++)
        {
            push_back(static_cast<base_t>((1<<(i+1)) - 1));
        }
    }
};
Bitmask bitmask(BITSINBASE);    // 1, 2, 4, 8, 16, ...
Bitsmask bitsmask(BITSINBASE);  // 1, 3, 7, 15, 31, ...

class Number {
public:
    Number() : _negative(false) {}
    Number(int val)
    {
        *this = val;
    }
    Number(const char* hexstr)
    {
        clear();
        int j=0;
        for (int i=strlen(hexstr)-1 ; i>=0 ; i--, j++) {
            int shift= j & (sizeof(base_t)*2-1);
            if (shift) {
                _number.back() |= digit2val(hexstr[i])<<(4*shift);
            }
            else {
                _number.push_back( static_cast<base_t>(digit2val(hexstr[i])) );
            }
        }
    }
    void swap(Number &a)
    {
        _number.swap(a._number);

        bool tmp= _negative;
        _negative= a._negative;
        a._negative= tmp;
    }

	// specify number, in normal reading order ( msb first )
    static Number CreateNumber(int n, ...)
    {
        Number val;
        va_list ap;
        va_start(ap, n);
        val.setvalue(n, ap);
        va_end(ap);

        return val;
    }
    void setvalue(int n, va_list ap)
    {
        _negative= false;
        _number.clear();
		_number.resize(n);
        if (n--) {
            int val= va_arg(ap, int);
            if (val<0) {
                _number[n]= -val;
                _negative= true;
            }
            else if (val>0) {
                _number[n]= val;
            }
        }
        while (n--)
        {
            _number[n]= va_arg(ap, base_t);
        }
    }
    base_t& operator[](size_t n)
    {
        return _number[n];
    }
    base_t operator[](size_t n) const
    {
        return _number[n];
    }
    void resize(size_t n)
    {
        _number.resize(n);
    }

    Number& operator=(const Number& num)
    {
        _negative= num._negative;
        _number= num._number;
        return *this;
    }
    int digit2val(char c)
    {
        return c<='9' ? c-'0'
                : c<='F' ? c-'A'+10
                : c<='f' ? c-'a'+10
                :0;
    }
    Number& operator=(int val)
    {
        _negative= val<0;
        _number.clear();
        if (val)
            _number.push_back(static_cast<base_t>(_negative?-val:val));
        return *this;
    }

    bool is_negative() const
    {
        return _negative;
    }
    static bool is_even(base_t a)
    {
        return (a&1)==0;
    }

    bool is_even() const
    {
        return (basesize()==0) || Number::is_even(_number[0]);
    }
    bool is_zero() const
    {
        if (basesize()==0)
            return true;
        for (number_t::const_iterator i= _number.begin() ; i!=_number.end() ; ++i)
            if (*i)
                return false;
        return true;
    }
    void clear()
    {
        _number.clear();
        _negative= false;
    }
    friend std::ostream& operator<< (std::ostream& os, const Number& num)
    {
        number_t::const_reverse_iterator i= num._number.rbegin();

        while (i!=num._number.rend() && (*i)==0)
            ++i;

        if (i==num._number.rend()) {
            os << "0";
            return os;
        }
        if (num._negative)
            os << '-';

        // if (os.flags()&ios_base::hex)
        os << "0x";

        // the first part
        if (i!=num._number.rend())
            os << std::hex << std::uppercase << *i++;
        while (i!=num._number.rend())
        {
            os << std::hex << std::uppercase << std::setfill('0') << std::setw(2*sizeof(base_t)) << *i++;
        }
        return os;
    }

/*

x = a+b*G+c*G^2+d*G^3

x     = {           a,           b,           c,           d }                      {     d,           c,           b,           a}

x<<01 = {       a<<01, b<<01|a>>31, c<<01|b>>31, d<<01|c>>31,       d>>31 }         { d>>31, d<<01|c>>31, c<<01|b>>31, b<<01|a>>31, a<<01 }
x<<16 = {       a<<16, b<<16|a>>16, c<<16|b>>16, d<<16|c>>16,       d>>16 }         { d>>16, d<<16|c>>16, c<<16|b>>16, b<<16|a>>16, a<<16 }
x<<31 = {       a<<31, b<<31|a>>01, c<<31|b>>01, d<<31|c>>01,       d>>01 }         { d>>01, d<<31|c>>01, c<<31|b>>01, b<<31|a>>01, a<<31 }
x<<32 = {           0,           a,           b,           c,           d }         {     d,           c,           b,           a,     0 }
x<<33 = {           0,       a<<01, b<<01|a>>31, c<<01|b>>31, d<<01|c>>31,  d>>31 } { d>>31, d<<01|c>>31, c<<01|b>>31, b<<01|a>>31, a<<01, 0 }
x<<48 = {           0,       a<<16, b<<16|a>>16, c<<16|b>>16, d<<16|c>>16,  d>>16 } { d>>16, d<<16|c>>16, c<<16|b>>16, b<<16|a>>16, a<<16, 0 }
x<<63 = {           0,       a<<31, b<<31|a>>01, c<<31|b>>01, d<<31|c>>01,  d>>01 } { d>>01, d<<31|c>>01, c<<31|b>>01, b<<31|a>>01, a<<31, 0 }

*/
//      4 << 3 << 2 << 1 << 0
//   6 << 5 << 4 << 3 << 2 << 1 << 0
    Number& operator<<=(unsigned int sh)
    {
        unsigned int bitshift= numbits(sh);
        unsigned int byteshift= numbases(sh);
		if (byteshift)
			_number.resize(basesize()+byteshift);
		if (bitshift==0) {
			for (number_t::reverse_iterator i= _number.rbegin()+byteshift ; i != _number.rend() ; ++i)
			{
				*(i-byteshift) = *i;
			}
			for (number_t::iterator i= _number.begin() ; i != _number.begin()+byteshift ; ++i)
				*i = 0;
			return *this;
		}
		_number.resize(basesize()+1);
		byteshift++;
        carry_t carry= 0;
        for (number_t::reverse_iterator i= _number.rbegin() ; i != _number.rend() ; ++i)
        {
			carry <<= BITSINBASE;
			if (i<_number.rend()-byteshift)
				carry |= (carry_t)*(i+byteshift) << bitshift;
			base_t value= static_cast<base_t>(carry>>BITSINBASE);
			if (i<_number.rend()-byteshift)
				value |= *(i+byteshift)>>(BITSINBASE-bitshift);
			*i =  value;
        }
		if (_number.back()==0)
			_number.resize(basesize()-1);
        return *this;
    }
/*
x = a+b*G+c*G^2+d*G^3

x     = {           a,           b,           c,           d }                      {     d,           c,           b,           a}

x>>01 = { b<<31|a>>01, c<<31|b>>01, d<<31|c>>01,       d>>01 }                      { d>>01, d<<31|c>>01, c<<31|b>>01, b<<31|a>>01 }
x>>16 = { b<<16|a>>16, c<<16|b>>16, d<<16|c>>16,       d>>16 }                      { d>>16, d<<16|c>>16, c<<16|b>>16, b<<16|a>>16 }
x>>31 = { b<<01|a>>31, c<<01|b>>31, d<<01|c>>31,       d>>31 }                      { d>>31, d<<01|c>>31, c<<01|b>>31, b<<01|a>>31 }
x>>32 = {           b,           c,           d }                                   {     d,           c,           b }
x>>33 = { c<<31|b>>01, d<<31|c>>01,       d>>01 }                                   { d>>01, d<<31|c>>01, c<<31|b>>01 }
x>>48 = { c<<16|b>>16, d<<16|c>>16,       d>>16 }                                   { d>>16, d<<16|c>>16, c<<16|b>>16 }
x>>63 = { c<<01|b>>31, d<<01|c>>31,       d>>31 }                                   { d>>31, d<<01|c>>31, c<<01|b>>31 }

*/
//        5 >> 4 >> 3 >> 2 >> 1 >> 0
//      4 >> 3 >> 2 >> 1 >> 0
    Number& operator>>=(unsigned int sh)
    {
        if (_number.empty())
            return *this;
        unsigned int bitshift= numbits(sh);
        unsigned int byteshift= numbases(sh);
        if (byteshift>basesize()) {
            clear();
            return *this;
        }
		if (bitshift==0) {
			for (number_t::iterator i= _number.begin()+byteshift ; i != _number.end() ; ++i)
			{
				*(i-byteshift) = *i;
			}
			_number.resize(basesize()-byteshift);
			return *this;
		}
        carry_t carry= 0;
        for (number_t::iterator i= _number.begin() ; i != _number.end() ; ++i)
        {
			carry >>= BITSINBASE;
			if (i<_number.end()-byteshift)
				carry |= *(i+byteshift);
			base_t value= static_cast<base_t>(carry >> bitshift);
			if (i<_number.end()-byteshift-1)
				value |= *(i+byteshift+1) << (BITSINBASE-bitshift);
			*i = value;
        }

        _number.resize(basesize()-byteshift);
		if (_number.back()==0)
			_number.resize(basesize()-1);
        return *this;
    }
    Number& operator+=(const Number& val)
    {
        if (_negative==val._negative) {
            add_abs(val._number);
        }
        else switch (compare_abs(val._number)) {
            case -1: {
                         number_t x= _number;
                         _number= val._number;
                         sub_abs(x);
                         _negative= !_negative;
                     }
                     break;
            case 0: clear(); break;
            case 1: {
                        sub_abs(val._number);
                    }
        }
        return *this;
    }
    Number& operator-=(const Number& val)
    {
        if (_negative!=val._negative) {
            add_abs(val._number);
        }
        else switch (compare_abs(val._number)) {
            case -1: {
                         number_t x= _number;
                         _number= val._number;
                         sub_abs(x);
                         _negative= !_negative;
                     }
                     break;
            case 0: clear(); break;
            case 1: {
                        sub_abs(val._number);
                    }
        }
        return *this;
    }
    Number& operator*=(base_t a)
    {
        carry_t carry= 0;
        for (number_t::iterator i= _number.begin() ; i != _number.end() ; ++i)
        {
            carry >>= BITSINBASE;
            carry += (carry_t)(*i) * a;
            (*i) = carry;
        }
		carry >>= BITSINBASE;
		if (carry)
			_number.push_back(carry);
        return *this;
    }
    Number& operator*=(const Number& a)
    {
        Number b= *this;
        clear();
        for (size_t i=0 ; i<a.basesize() ; i++) {

            *this += b*a._number[i];

            b<<=BITSINBASE;
        }
        _negative= (_negative!=a._negative);
        return *this;
    }

    Number& operator/=(const Number& val)
    {
        Number q, r;
        divmod(val, q, r);
        *this= q;
        return *this;
    }
    Number& operator%=(const Number& val)
    {
        Number q, r;
        divmod(val, q, r);
        *this= r;
        return *this;
    }

    Number operator-(const Number& val) const
    {
        Number result= *this;
        result -= val;
        return result;
    }
    Number operator+(const Number& val) const
    {
        Number result= *this;
        result += val;
        return result;
    }
    Number operator/(const Number& val) const
    {
        Number result= *this;
        result /= val;
        return result;
    }
    Number operator%(const Number& val) const
    {
        Number result= *this;
        result %= val;
        return result;
    }
    Number operator*(base_t a) const
    {
        Number result= *this;
        result *= a;
        return result;
    }
    Number operator*(const Number& a) const
    {
        Number result= *this;
        result *= a;
        return result;
    }


    bool operator>(const Number &val) const
    {
        Number x= *this;
        x -= val;
        return !x._negative && !x.is_zero();
    }
    virtual bool operator==(const Number &val) const
    {
        Number x= *this;
        x -= val;
        return x.is_zero();
    }
    bool operator<(const Number &val) const
    {
        Number x= *this;
        x -= val;
        return x._negative && !x.is_zero();
    }
    bool operator>=(const Number &val) const
    {
        return !(*this < val);
    }
    bool operator!=(const Number &val) const
    {
        return !(*this == val);
    }
    bool operator<=(const Number &val) const
    {
        return !(*this > val);
    }
    size_t basesize() const
    {
        return _number.size();
    }
    size_t bitsize() const
    {
        return basesize()*BITSINBASE;
    }
    void negate()
    {
        _negative = ! _negative;
    }
    bool bittest(unsigned int bit) const
    {
        return (_number[numbases(bit)] & bitmask[numbits(bit)])!=0;
    }
    int highestbit() const
    {
        int byte;
        for (byte=basesize()-1 ; _number[byte]==0 && byte>=0 ; byte--)
            ;
        if (byte<0)
            return -1;
        int bit;
        for (bit=BITSINBASE-1 ; (_number[byte]&bitmask[bit])==0 && bit>=0 ; bit--)
            ;
        return byte*BITSINBASE+bit;
    }
    int highestbase() const
    {
        int byte;
        for (byte=basesize()-1 ; _number[byte]==0 && byte>=0 ; byte--)
            ;
        return byte;
    }


private:
    bool _negative;
    number_t _number;

    unsigned int numbases(unsigned int bitcount) const { return bitcount/BITSINBASE; }
    unsigned int numbits(unsigned int bitcount) const { return bitcount&(BITSINBASE-1); }

    void add_abs(const number_t& val)
    {
        carry_t carry= 0;

        if (basesize() < val.size())
            _number.resize(val.size());

        number_t::const_iterator i_val= val.begin();
        number_t::iterator i_sum= _number.begin();
        while (i_sum!=_number.end())
        {
            carry += (carry_t)(*i_sum) + (i_val!=val.end() ? (*i_val):0);
            (*i_sum) = static_cast<base_t>(carry);
            carry >>= BITSINBASE;
            ++i_sum;
            if (i_val!=val.end()) ++i_val;
        }
        if (carry)
            _number.push_back(static_cast<base_t>(carry));
		if (_number.back()==0)
			_number.resize(basesize()-1);
	}
    void sub_abs(const number_t& val)
    {
        carry_t carry= 0;

        if (basesize() < val.size())
            _number.resize(val.size());
        number_t::const_iterator i_val= val.begin();
        number_t::iterator i_sum= _number.begin();
        while (i_sum!=_number.end())
        {
            carry += (carry_t)(*i_sum) -  (i_val!=val.end() ? (*i_val):0);
            (*i_sum) = static_cast<base_t>(carry);
            carry >>= BITSINBASE;
            ++i_sum;
            if (i_val!=val.end()) ++i_val;
        }
        if (carry)
            _number.push_back(static_cast<base_t>(carry));
		if (_number.back()==0)
			_number.resize(basesize()-1);
    }
    int compare_abs(const number_t& val) const
    {
        if (basesize() < val.size())
            return -1;
        if (basesize() > val.size())
            return 1;

        number_t::const_reverse_iterator i_val= val.rbegin();
        number_t::const_reverse_iterator i_num= _number.rbegin();
        while (i_num!=_number.rend() && i_val!=val.rend())
        {
            if ((*i_num)<(*i_val))
                return -1;
            if ((*i_num)>(*i_val))
                return 1;
            ++i_num;
            ++i_val;
        }
        return 0;
    }
    // handbook of applied cryptography - algorithm 14.20
    //
    // input
    //    x=sum(x[i]*b^i, i=0..n)
    //    y=sum(y[i]*b^i, i=0..t)   y[t]!=0
    //    1<=t<=n
    // output
    //    q=sum(x[i]*b^i, i=0..n-t)
    //    r=sum(x[i]*b^i, i=0..t)
    //    such that  x=q*y+r,  0<=r<y
    // 
    // q[j]=0  for j=0 to n-t
    // while (x>=y*b^(n-t))
    //    q[n-t]++; 
    //    x -= y*b^(n-t)
    // for i=n downto t+1
    //    if (x[i]==y[i])
    //       q[i-t-1]=b-1
    //    else
    //       q[i-t-1]=trunc((x[i]*b+x[i-1])/y[t])
    //    while (q[i-t-1]*(y[t]*b+y[t-1]) > x[i]*b^2+x[i-1]*b+x[i-2])
    //       q[i-t-1]--;
    //    x -= q[i-t-1]*y*b^(i-t-1)
    //    if (x<0) 
    //       x += y*b^(i-t-1)
    //       q[i-t-1] --;
    // r=x
    // return (q,r)
    void divmod(const Number& y, Number& q, Number& r) const
    {
        Number x= *this;
        int t= y.highestbase(); 
		if (t<1) {
			divmod(y._number[0], q, r);
            // todo: correct sign of result.
			return;
		}
        int n= x.basesize()-1;
		if (n<t) {
			r= *this;
			q.clear();
			return;
		}
        q.resize(n-t+1);
        r.resize(t+1);
        for (int j=0 ; j <= n-t ; j++) {
            q._number[j] = 0;
        }
        Number yy= y;
        yy <<= BITSINBASE*(n-t);
        while (x >= yy) {
            q._number[n-t]++;
            x -= yy;
        }
        for (int i=n ; i>=t+1 ; i--) {
            yy >>= BITSINBASE;
            if (x._number[i]==y._number[t]) {
                q._number[i-t-1]= 0xFFFFFFFF;
            }
            else {
                q._number[i-t-1]= ((carry_t)x._number[i]<<32+x._number[i-1])/y._number[t];

            }
            while(((carry_t)y._number[t]<<32+y._number[t-1])*q._number[i-t-1] > x._number[i]<<64+x._number[i-1]<<32+x._number[i-2]) {
                q._number[i-t-1]--;
                x -= yy*q._number[i-t-1];
            }
            if (x._negative) {
                x += yy;
                q._number[i-t-1]--;
            }
        }

        // x=q*y+r    : x= q*y+r
        // x=q*-y+r   : x=-q*y+r
        // -x=q*y+r   : x=-q*y-r
        // -x=q*-y+r  : x= q*y-r

        r= x;
        q._negative=x.is_negative()!=y.is_negative();
    }
	void divmod(base_t y, Number& q, Number& r) const
	{
		carry_t carry= 0;
		// sum(x[i]*B^i) = sum(q[i]*B^i)*y+r
		// -> q[m] = x[i]/y
		q.clear();
		q.resize(basesize());

		for (int i=basesize()-1 ; i>=0 ; i--)
		{
			q._number[i]= ((carry<<32)+_number[i])/y;
			carry= ((carry<<32)+_number[i])%y;
		}
		r= carry;
	}

    void shrink()
    {
        while (!_number.empty() && _number.back()==0)
            _number.pop_back();
    }


};
Number operator<<(const Number&num, unsigned int sh)
{
    Number x= num;
    x <<= sh;
    return x;
}

Number operator>>(const Number&num, unsigned int sh)
{
    Number x= num;
    x >>= sh;
    return x;
}

class Modular : public Number {
public:
    Modular(int n, const Number& m)
        : Number(n), _modulus(m)
    {
    }
    Modular(const Number& n, const Number& m)
        : Number(n), _modulus(m)
    {
    }
    Modular& operator+=(const Number& num)
    {
        Number::operator+=(num);
        modtrunc();
        return *this;
    }
    Modular& operator<<=(unsigned int sh)
    {
        Number::operator<<=(sh);
        modtrunc();
        return *this;
    }
    Modular& operator*=(const Number& mult)
    {
        Number bits= *this;
        Modular lshifter= Modular(mult, _modulus);
        clear();

        for (size_t i_bits= 0 ; i_bits<bits.basesize() ; ++i_bits)
        {
            base_t bitval= bits[i_bits];
            for (size_t i=0 ; i<BITSINBASE ; i++)
            {
                if (bitval & bitmask[i]) {
                    *this += lshifter;
                }
                lshifter <<= 1;
            }
        }
        return *this;
    }
	virtual bool operator==(const Number& val) const
	{
		Number diff= (*this - val) % _modulus;
		return diff.is_zero();
	}
    Modular modexp(const Number& exp)
    {
        Modular num2= Modular(1, _modulus);

        num2 *= *this;

        Modular result= Modular(1, _modulus);

        for (size_t expbit=0 ; expbit<exp.bitsize() ; expbit++) {
            if (exp.bittest(expbit)) {
                  result *= num2;
            }
            num2 *= num2;
        }
        return result;
    }



protected:
    const Number _modulus;

	// this is not sufficient, for doing m % R ( from the montgomery modexp!!
    void modtrunc()
    {
		/*
        number_t::reverse_iterator i_val= _number.rbegin();
        while (i_val!=_number.rend() && (*i_val)==0)
            ++i_val;
        number_t::const_reverse_iterator i_mod= _modulus._number.rbegin();
        while (i_mod!=_modulus._number.rend() && (*i_mod)==0)
            ++i_mod;
        int diff= distance(i_val, _number.rend()) - distance(i_mod, _modulus._number.rend());
        if (diff<0)
            return;
        if (diff>1) {
            printf("!!!diff = %d\n", diff);
            exit(1);
        }
        if (diff==1) {
            sub_abs(_modulus._number);
            return;
        }
        while (i_val!=_number.rend() && (*i_val)==(*i_mod)) {
            ++i_val;
            ++i_mod;
        }
        if (i_val!=_number.rend() && (*i_val)>=(*i_mod)) {
            sub_abs(_modulus._number);
        }
		*/
		*this %= _modulus;
    }
};



// for x,y - calculates a,b,v such that a*x+b*y=v = gcd(x,y)
void gcd_calc(const Number& x, Number& a, const Number& y, Number& b, Number& v)
{
    // handbook of applied cryptography - algorithm 14.61

    int nshifts= 0;
    Number xx= x;
    Number yy= y;
    while (xx.is_even() && yy.is_even()) {
        xx>>=1;
        yy>>=1;
        nshifts++;
    }
    Number u = xx;
    v = yy;
    Number A= 1;
    Number B;
    Number C;
    Number D= 1;
//    std::cout << "(a,b,c,d,u,v) : " << A << ',' << B << ',' << C << ',' << D << ',' << u << ',' << v << '\n';
    while (!u.is_zero()) {
        while (u.is_even()) {
            u>>=1;
            if (!(A.is_even() && B.is_even())) {
                A+=yy; 
                B-=xx;
            }
            A>>=1;
            B>>=1;
        }
        while (v.is_even()) {
            v>>=1;
            if (!(C.is_even() && D.is_even())) {
                C+=yy; 
                D-=xx;
//                std::cout << "c,d+-=xy : " << C << ',' << D << '\n';
            }
            C>>=1;
            D>>=1;
//            std::cout << "c,d : " << C << ',' << D << '\n';
        }
        if (u>=v) {
            u-=v;
            A-=C;
            B-=D;
//            std::cout << "if(a,b,u) : " << A << ',' << B << ',' << u << '\n';
        }
        else {
            v-=u;
            C-=A;
            D-=B;
//            std::cout << "else(c,d,v) : " << C << ',' << D << ',' << v << '\n';
        }
    }
    a= C;
    b= D;
    v<<=nshifts;
}
/*
void gcd_calc(const Number& x, Number& a, const Number& y, Number& b, Number& v)
{
    Number m[2][3];
    int n,i;
    long q;
    m[0][0]= 1; m[0][1]= 0; m[0][2]= x;
    m[1][0]= 0; m[1][1]= 1; m[1][2]= y;
    i=0; n=0;
    while (!m[1-i][2].is_zero())
    {
        q= m[i][2]/m[1-i][2];
        m[i][2] -= q*m[1-i][2];
        m[i][1] -= q*m[1-i][1];
        m[i][0] -= q*m[1-i][0];
        i=1-i;
        n++;
    }
	if (m[i][0]<0) { m[i][0]+=y; m[i][1]-=x; }
	a=m[i][0]; b=m[i][1];  v=m[i][2];
}
*/
/*

   gcd=x*a+b*y

 */
void gcd_calc2(const Number& x, const Number& y, Number& g)
{
    Number u= x;
    Number v= y;
    int shift;

    // GCD(0,x) := x
    if (u.is_zero()) {
      g= v;
      return;
    }
    if (v.is_zero()) {
      g= u;
      return;
    }

    // Let shift := lg K, where K is the greatest power of 2
    // dividing both u and v. 
    for (shift = 0; u.is_even()  && v.is_even() ; ++shift) {
        u >>= 1;
        v >>= 1;
    }

    while (u.is_even())
      u >>= 1;

    // From here on, u is always odd.
    do {
        while (v.is_even())
          v >>= 1;

        // Now u and v are both odd, so diff(u, v) is even.
        // Let u = min(u, v), v = diff(u, v)/2. 
        if (u < v) {
            v -= u;
        } else {
            Number diff = u; diff -= v;
            u = v;
            v = diff;
        }
        v >>= 1;
    } while (!v.is_zero());

    g = u;
    g <<= shift;

}
class Montgomery {
public:
    Montgomery(const Number& modulus, const Number& base)
        : _modulus(modulus), _base(base), _modexpvars(false)
    {
        Number g;
        gcd_calc(modulus, _inversemodulus, base, _inversebase, g);
        _inversemodulus.negate();
		if (_inversemodulus.is_negative()) {
			_inversemodulus += base;
			_inversebase += modulus;
		}

        Modular x(modulus, base);
        x *= _inversemodulus;
		if (x!=Number(-1)) {
			std::cout << "montgomery,b="<<base<<", m("<<_modulus<<")*m^-1 ("<<_inversemodulus<<")= " << x << '\n';
		}
    }
    // handbook of applied cryptography - algorithm 14.36
    //
    // input
    //    m = sum(m[i]*b^i, i=0..n-1)
    //    x = sum(m[i]*b^i, i=0..n-1)
    //    y = sum(m[i]*b^i, i=0..n-1)
    //    0<=x,y<m
    //    R = b^n    gcd(m,b)=1
    //    m' = -m^-1(mod b)
    // output
    //    A=sum(a[i]*b^i,i=0..n) = x*y*R^-1 (mod m)
    //
    // A=0
    // for i=0 to n-1
    //    u= (a[0]+x[i]*y[0])*m' (mod b)
    //    A = (A+x[i]*y+u*m)/b
    // A -= m  if (A>=m)
    // return(A)
    void multiply(const Number& x, const Number& y, Number& a)
    {
        a.clear();
        a.resize(_modulus.basesize()+1);
        for (size_t i=0 ; i<x.basesize() ; i++) {
            base_t u = (a[0]+x[i]*y[0]) * _inversemodulus[0];
            
			a += y*x[i] + _modulus*u;
			a >>= BITSINBASE;
        }
        if (a > _modulus)
            a-=_modulus;
    }

    // returns value * base^-1 ( mod modulus )
    // handbook of applied cryptography - algorithm 14.32
    //
    // input
    //    m = sum(m[i]*b^i, i=0..n-1)
    //    R = b^n    gcd(m,b)=1
    //    m' = -m^-1(mod b)
    //    T = sum(t[i]*b^i, i=0..2*n-1)
    //    T < n*R
    // output
    //    T*R^-1 mod m
    //
    // A=T
    // for i=0 to n-1
    //    u = a[i]*m' mod b
    //    A += u*m*b^i
    // A /= b^n
    // A -= m  if A>=m
    // return A
    //
    //
    // T  = val
    // A  = a
    // m  = _modulus
    // b  = _base
    // n  = _modulus.size()
    // m' = _inversemodulus
    // R  = b^n = 1<<(BITSINBASE*n)
    void reduce(const Number& val, Number& a)
    {
        a= val;
        for (size_t i=0 ; i<_modulus.basesize() ; i++) {
            base_t u= a[i]*_inversemodulus[0];

            Number um= _modulus;
            um *= u;
            um <<= i*BITSINBASE;

            a+=um;
        }
        a>>=_modulus.bitsize();
        if (a>=_modulus)
            a-=_modulus;
    }
    // handbook of applied cryptography - algorithm 14.94
    //
    // input
    //    m = sum(m[i]*b^i, i=0..l-1)
    //    R = b^l
    //    m' = -m^-1(mod b)
    //    e = sum(e[i]*2^i, i=0..t)  e[t]=1
    //    x   1<=x<m
    // output
    //    x^e(mod m)
    //
    // xt=Montmul(x,R^2 mod m)
    // A =R mod m
    //
    // for i=t downto 0
    //    A=Montmul(A,A)
    //    A=Montmul(A,xt) if (e[i])
    // A=Montmul(A,1)
    // return A
    //
    void modexp1(const Number& num, const Number& exp, Number& result)
    {
        if (!_modexpvars) {
            _r_mod_m= _base % _modulus;
            _rr_mod_m= (_base * _base) % _modulus;
            _modexpvars= true;
        }
        Number a= _r_mod_m;                 // a = R mod m
        Number b;
        multiply(num, _rr_mod_m, b);        // b = num*R mod m

		// num : x   = input value
        // a : A       = R mod m
        // b : x-tilde = Mont(x,R^2 mod m)
        result= 1;

        for (int bit= exp.highestbit() ; bit>=0 ; bit--)
        {
            Number a2;
            multiply(a, a, a2);             // a := a * a * R^-1 mod m = R mod m
            a= a2;
            if (exp.bittest(bit))
            {
                Number ab;
                multiply(a, b, ab);         // a := a * b * R^-1 mod m = R * num  mod m
                a= ab;
            }
            bit--;
        }
        multiply(a, Number(1), result);
    }
    void modexp(const Number& num, const Number& exp, Number& result)
    {
        if (!_modexpvars) {
            _r_mod_m= _base % _modulus;
            _rr_mod_m= (_base * _base) % _modulus;
            _modexpvars= true;
        }

        Number a= _r_mod_m;
        Number b= (num * _base) % _modulus;


        for (int bit= exp.highestbit() ; bit>=0 ; bit--)
	    {
			if (exp.bittest(bit)) {
				Number ab;
				multiply(a, b, ab);
				a= ab;
			}
			Number bb;
			multiply(b, b, bb);
			b= bb;
		}
        multiply(a, Number(1), result);
    }

	const Number& base() const
    {
        return _base;
    }
    const Number& modulus() const
    {
        return _modulus;
    }

private:
    Number _modulus;
    Number _base;
    Number _inversebase;
    Number _inversemodulus;

    Number _r_mod_m;
    Number _rr_mod_m;
    bool _modexpvars;
};


#include <windows.h>
#include <stdlib.h>
#include <string.h>


void run_tests()
{

    {
    // test left+right shift
    Number a1= 0xabcdef12;
    a1 <<= BITSINBASE*3;
    a1 >>= BITSINBASE*3;
    if (a1!=0xabcdef12) {
        std::cout << "shift test#1 failed\n";
        std::cout << " >> " << (BITSINBASE*3) << " =" << a1 << "\n";
        exit(1);
    }
    a1 <<= BITSINBASE+8;
    a1 >>= BITSINBASE+8;
    if (a1!=0xabcdef12) {
        std::cout << "shift test#2a failed\n";
        std::cout << " >> " << (BITSINBASE+8) << " =" << a1 << "\n";
        exit(1);
    }
    a1 <<= 2*BITSINBASE-8;
    a1 >>= 2*BITSINBASE-8;
    if (a1!=0xabcdef12) {
        std::cout << "shift test#2b failed\n";
        std::cout << " >> " << (BITSINBASE+8) << " =" << a1 << "\n";
        exit(1);
    }
    Number x= 1;
	x <<=32;
	x <<=32;
	x <<=32;
	x >>=32;
	x >>=32;
	x >>=32;
	if (x!=1) {
        std::cout << "shift test#3a failed\n";
        std::cout << " >> " << (BITSINBASE*16) << " =" << x << "\n";
        exit(1);
	}
	x <<=63;
    for (size_t i=0 ; i<BITSINBASE*16 ; i++)
        x<<=1;
    for (size_t i=0 ; i<BITSINBASE*16 ; i++)
        x>>=1;
	x >>=63;
    if (x!=Number(1)) {
        std::cout << "shift test#3b failed\n";
        std::cout << " >> " << (BITSINBASE*16) << " =" << x << "\n";
        exit(1);
    }
    x>>=1;
    if (!x.is_zero()) {
        std::cout << "shift test failed\n";
        std::cout << " >> 1 = " << x << "\n";
        exit(1);
    }
    std::cout << "shift test done\n";
    }

    // test add/sub
    {
    if ((Number(1)-Number(2))!=Number(-1))
        std::cout << "ERROR: 1 - 2 = " << (Number(1)-Number(2)) << '\n';
    Number x,y,z;
    srand(0);
	x.resize(8);
	y.resize(8);
	z.resize(8);
    for (int i=0 ; i<8 ; i++)
        x[i]= static_cast<base_t>(rand());
    for (int i=0 ; i<8 ; i++)
        y[i]= static_cast<base_t>(rand());
    for (int i=0 ; i<8 ; i++)
        z[i]= static_cast<base_t>(rand());

    if (x > y) {
        x.swap(y);
    }
    Number x0= x;
    x-= y;
    x-= z;
    x+= y;
    x+= z;

    if (x0!=x) {
        std::cout << "ERROR: x-y-z+y+z!=x : " << x << " !!!\n";
        std::cout << " x = " << x0 << '\n';
        std::cout << " y = " << y << '\n';
        std::cout << " z = " << z << '\n';
    }

    std::cout << "add/sub test done\n";
    }

    // testing scalar mult
    {
        Number a;

        for (int i=0 ; i<10 ; i++)
        {
            a+=1;
            a*=609;
        }
        if (a!=Number::CreateNumber(3, 0x16B628FC, 0x703EEFC4, 0xE75A36AA)) {
            std::cout << "ERROR: (a+1)*609 .. 10 = " << a << '\n';
        }

        std::cout << "scalar mult test done\n";
    }
    // test mul
    {
    Number x, y;

    x=123456789;
    y=987654321;

    x*=y;
    x*=y;
    x*=y;

    if (x!=Number::CreateNumber(4, 0x16E838, 0x96AC8881, 0x2D61A1AF, 0x0CDF1765)) {
        std::cout << "ERROR: x * y**3 = " << x << '\n';
    }

    std::cout << "mul test done\n";
    }

    // testing %
    {
    Number x,y,z;

    x= Number("0857BDE36F680CDA9535784BC4AEC5F4344131071B419F732AC9C74D0E61DB49DD958C7344236E0279DF009C6E66AEC6BA574C2820D4AEB0C4D814C8E184C6EA7E6D8AA3E15D1C251C78C5364EA2B3EDB3C19E90739AFA765506242E78FCDC71A87EFDFE2DF6CE6039FC62CB3B360CB77CD5574292282DF352886CBC3FCFBFF2");
    y= Number("10001");
    z= Number("F765A3A0C9C291D81A56FE73794A746B8DA23DBE155D0D495B49D581B5C6545F449A10FDF1C26A92FBD1F43A0687044927A6A21B69A73999E6083D03ACDAFFA6409F1BC71D810628F6E18F76231ED6E22D54ED2502E66F8A33D0D5F07B3EB605F7418110E2EF9A5EE77B070F4EADFCF3D70C53E870F29C9D4F229F2CB6C25383");

    z %= x;
    if (z!=Number::CreateNumber(32, 0x57520DD, 0x2AF91D15, 0x33485DDE, 0x317E07C1, 0xA23FAFEF, 0xFEEDFD3D, 0x826E41C7, 0x14B07D01, 0x2AA927EF, 0x39BEF44B, 0x2D8EE281, 0x84E537C6, 0x0BC3018F, 0xB18F6F93, 0x998DE242, 0x20D07715, 0xEE366736, 0x95F4D5F4, 0xBD33384F, 0x3AB074F4, 0xD065F6C7, 0xEA581022, 0x921EBCAC, 0xC699BD25, 0xE0DEBB45, 0xADFA3978, 0x55E3D609, 0x998E8C2A, 0xB2E1715D, 0xE265680C, 0xF5AE4DD9, 0x7C399519)) {
        std::cout << "ERROR: z %% x = " << z << '\n';
    }
    std::cout << "modulus test done\n";
    }

    // testing Modular
    {
/*
    Number x,y,z;

    x= Number("0857BDE36F680CDA9535784BC4AEC5F4344131071B419F732AC9C74D0E61DB49DD958C7344236E0279DF009C6E66AEC6BA574C2820D4AEB0C4D814C8E184C6EA7E6D8AA3E15D1C251C78C5364EA2B3EDB3C19E90739AFA765506242E78FCDC71A87EFDFE2DF6CE6039FC62CB3B360CB77CD5574292282DF352886CBC3FCFBFF2");
    y= Number("10001");
    z= Number("F765A3A0C9C291D81A56FE73794A746B8DA23DBE155D0D495B49D581B5C6545F449A10FDF1C26A92FBD1F43A0687044927A6A21B69A73999E6083D03ACDAFFA6409F1BC71D810628F6E18F76231ED6E22D54ED2502E66F8A33D0D5F07B3EB605F7418110E2EF9A5EE77B070F4EADFCF3D70C53E870F29C9D4F229F2CB6C25383");
    Number q;
    DWORD t0= GetTickCount();
    Modular xx= Modular(x, z);
    q= xx.modexp(y);

    std::cout << (GetTickCount()-t0) << " ticks\n";

    Number r= Number("1FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF003021300906052B0E03021A050004147020B30D42ED7CDF0E849F30E8B2E137941E9691");
    if (r!=q) {
        std::cout << "got:      " << q << '\n';
        std::cout << "expected: " << r << '\n';
    }
    std::cout << "Modular calc test done\n";
*/
    }

    // testing gcd
    {
    Number a;
    Number b;
    Number g;

    // (-7-29*i)*693 + (8+33*i)*609 = 21

    gcd_calc2(Number(693), Number(609), g);
    if (g!=0x15) {
        std::cout << "ERROR: gcd=" << g << '\n';
    }
    gcd_calc(Number(693), a, Number(609), b, g);
    if (a!=-0xB5) {
        std::cout << "ERROR: a="<< a << '\n';
    }
    if (b!=0xCE) {
        std::cout << "ERROR: b="<< b << '\n';
    }
    if (g!=0x15) {
        std::cout << "ERROR: g="<< g << '\n';
    }

    if (a*693+b*609!=g) {
        std::cout << "error in gcd calc\n";
    }
    std::cout << "gcd test done\n";
    }

    // testing montgomery
    {
    Montgomery mg(1234567, Number(1)<<32);
    Number ra;
    mg.reduce(Number(123456789)<<32, ra);
    if (((ra*mg.base())%mg.modulus())!=((Number(123456789)<<32)%mg.modulus()))
        std::cout << "ERROR: t1: ra= " << ra << '\n';

    Number r123;
    mg.reduce(Number(123), r123);
    if (r123!=0xEF4BA) {
        std::cout << "ERROR: 123 -> " << r123 << '\n';
    }
    Number r456;
    mg.reduce(Number(456), r456);
    if (r456!=0x1B075)
        std::cout << "ERROR: 456 -> " << r456 << '\n';
    Number rprod;
    mg.reduce(Number(123*456), rprod);
    if (rprod!=0x906A)
        std::cout << "ERROR: 123*456 -> " << rprod << '\n';

    if (((r123*r456*mg.base())%mg.modulus())!=(rprod%mg.modulus())) {
        std::cout << "ERROR: r(123)*r(456) != r(123*456) (mod m)\n";
    }

	Number rmult;
	mg.multiply(Number(123), Number(456), rmult);
	if (rmult!=Number(0x906A))
		std::cout << "ERROR: mult : " << rmult << '\n';

	mg.modexp(Number(123), Number(31), ra);
	if (ra != Number(0x3EBD8))
		std::cout << "ERROR: 123^31 -> " << ra << '\n';

    std::cout << "montgomery test done\n";
    }
}
int main(int argc, char **argv)
{
    run_tests();
    if (argc!=4) {
        printf("Usage: modexp num exp mod\n");
        return 1;
    }
    char *numstr= argv[1];
    char *expstr= argv[2];
    char *modstr= argv[3];

    // strip leading zeros
    while (*numstr=='0') numstr++;
    while (*expstr=='0') expstr++;
    while (*modstr=='0') modstr++;

    Number num= Number(numstr);
    Number exp= Number(expstr);
    Number mod= Number(modstr);

    if (num.bitsize() > mod.bitsize()) {
        printf("number cannot be larger than mod\n");
        return 1;
    }

    // std::cout << "num=" << num << "\n";
    // std::cout << "mod=" << mod << "\n";
    // std::cout << "exp=" << exp << "\n";

    Number result;
/*
    try {
        Modular mm= Modular(num, mod);
        result= mm.modexp(exp);
    } catch(char *msg) {
        std::cout << "EXCEPTION: " << msg << '\n';
    }
    std::cout << "modexp: " << result << '\n';
*/
    Montgomery m(mod, Number(1)<<32);
    m.modexp(num, exp, result);
    std::cout << "montgomery: " << result << '\n';

}