// 誤解法(法Bでの零割り)チェック #pragma GCC optimize ( "O3" ) #pragma GCC optimize( "unroll-loops" ) #pragma GCC target ( "sse4.2,fma,avx2,popcnt,lzcnt,bmi2" ) #include #include #include #include using namespace std; using ll = long long; #define MAIN main #define TYPE_OF( VAR ) remove_const::type >::type #define UNTIE ios_base::sync_with_stdio( false ); cin.tie( nullptr ) #define CEXPR( LL , BOUND , VALUE ) constexpr const LL BOUND = VALUE #define CIN( LL , A ) LL A; cin >> A #define ASSERT( A , MIN , MAX ) assert( MIN <= A && A <= MAX ) #define CIN_ASSERT( A , MIN , MAX ) CIN( TYPE_OF( MAX ) , A ); ASSERT( A , MIN , MAX ) #define FOR( VAR , INITIAL , FINAL_PLUS_ONE ) for( TYPE_OF( FINAL_PLUS_ONE ) VAR = INITIAL ; VAR < FINAL_PLUS_ONE ; VAR ++ ) #define REPEAT( HOW_MANY_TIMES ) FOR( VARIABLE_FOR_REPEAT , 0 , HOW_MANY_TIMES ) #define QUIT return 0 #define RETURN( ANSWER ) cout << ( ANSWER ) << "\n"; QUIT template inline T Residue( const T& a , const T& p ){ return a >= 0 ? a % p : p - ( - a - 1 ) % p - 1; } #define POWER( ANSWER , ARGUMENT , EXPONENT ) \ TYPE_OF( ARGUMENT ) ANSWER{ 1 }; \ { \ TYPE_OF( ARGUMENT ) ARGUMENT_FOR_SQUARE_FOR_POWER = ( ARGUMENT ); \ TYPE_OF( EXPONENT ) EXPONENT_FOR_SQUARE_FOR_POWER = ( EXPONENT ); \ while( EXPONENT_FOR_SQUARE_FOR_POWER != 0 ){ \ if( EXPONENT_FOR_SQUARE_FOR_POWER % 2 == 1 ){ \ ANSWER *= ARGUMENT_FOR_SQUARE_FOR_POWER; \ } \ ARGUMENT_FOR_SQUARE_FOR_POWER *= ARGUMENT_FOR_SQUARE_FOR_POWER; \ EXPONENT_FOR_SQUARE_FOR_POWER /= 2; \ } \ } \ #define POWER_MOD( ANSWER , ARGUMENT , EXPONENT , MODULO ) \ TYPE_OF( ARGUMENT ) ANSWER{ 1 }; \ { \ TYPE_OF( ARGUMENT ) ARGUMENT_FOR_SQUARE_FOR_POWER = ( MODULO + ( ARGUMENT ) % MODULO ) % MODULO; \ TYPE_OF( EXPONENT ) EXPONENT_FOR_SQUARE_FOR_POWER = ( EXPONENT ); \ while( EXPONENT_FOR_SQUARE_FOR_POWER != 0 ){ \ if( EXPONENT_FOR_SQUARE_FOR_POWER % 2 == 1 ){ \ ANSWER = ( ANSWER * ARGUMENT_FOR_SQUARE_FOR_POWER ) % MODULO; \ } \ ARGUMENT_FOR_SQUARE_FOR_POWER = ( ARGUMENT_FOR_SQUARE_FOR_POWER * ARGUMENT_FOR_SQUARE_FOR_POWER ) % MODULO; \ EXPONENT_FOR_SQUARE_FOR_POWER /= 2; \ } \ } \ template class TwoByTwoMatrix { public: // private: T m_M00; T m_M01; T m_M10; T m_M11; public: inline TwoByTwoMatrix( const T& M00 , const T& M01 , const T& M10 , const T& M11 ) noexcept; inline TwoByTwoMatrix( const T& n = T() ) noexcept; inline TwoByTwoMatrix& operator=( const TwoByTwoMatrix& mat ) noexcept; inline TwoByTwoMatrix& operator*=( const TwoByTwoMatrix& mat ) noexcept; inline TwoByTwoMatrix operator*( const TwoByTwoMatrix& mat ); }; template inline TwoByTwoMatrix::TwoByTwoMatrix( const T& M00 , const T& M01 , const T& M10 , const T& M11 ) noexcept : m_M00( M00 ) , m_M01( M01 ) , m_M10( M10 ) , m_M11( M11 ) {} template inline TwoByTwoMatrix::TwoByTwoMatrix( const T& n ) noexcept : m_M00( n ) , m_M01() , m_M10() , m_M11( n ) {} template inline TwoByTwoMatrix& TwoByTwoMatrix::operator=( const TwoByTwoMatrix& mat ) noexcept { if( &mat != this ){ m_M00 = mat.m_M00; m_M01 = mat.m_M01; m_M10 = mat.m_M10; m_M11 = mat.m_M11; } return *this; } template inline TwoByTwoMatrix& TwoByTwoMatrix::operator*=( const TwoByTwoMatrix& mat ) noexcept { return operator=( *this * mat ); } template inline TwoByTwoMatrix TwoByTwoMatrix::operator*( const TwoByTwoMatrix& mat ) { return TwoByTwoMatrix( m_M00 * mat.m_M00 + m_M01 * mat.m_M10 , m_M00 * mat.m_M01 + m_M01 * mat.m_M11 , m_M10 * mat.m_M00 + m_M11 * mat.m_M10 , m_M10 * mat.m_M01 + m_M11 * mat.m_M11 ); } template class QuotientRing { protected: INT m_n; const INT* m_p_M; public: inline QuotientRing() noexcept; inline QuotientRing( const INT& n , const INT* const & p_M = nullptr ) noexcept; inline QuotientRing( const QuotientRing& n ) noexcept; inline QuotientRing& operator+=( const QuotientRing& n ) noexcept; inline QuotientRing& operator+=( const INT& n ) noexcept; // operator<が定義されていても負の数は正に直さず剰余を取ることに注意。 inline QuotientRing& operator-=( const QuotientRing& n ) noexcept; inline QuotientRing& operator-=( const INT& n ) noexcept; inline QuotientRing& operator*=( const QuotientRing& n ) noexcept; inline QuotientRing& operator*=( const INT& n ) noexcept; inline const INT& Represent() const noexcept; inline const INT& GetModulo() const noexcept; // m_nの正負やm_p_Mの一致込みの等号。 static inline bool Equal( const QuotientRing& n0 , const QuotientRing& n1 ) noexcept; template static QuotientRing Power( const QuotientRing& n , const T& exponent ); }; template inline bool operator==( const QuotientRing& n0 , const QuotientRing& n1 ) noexcept; template inline bool operator!=( const QuotientRing& n0 , const QuotientRing& n1 ) noexcept; template inline QuotientRing operator+( const QuotientRing& n0 , const T& n1 ) noexcept; template inline QuotientRing operator-( const QuotientRing& n0 , const T& n1 ) noexcept; template inline QuotientRing operator*( const QuotientRing& n0 , const T& n1 ) noexcept; template inline QuotientRing Power( const QuotientRing& n , const T& exponent ); template inline QuotientRing::QuotientRing() noexcept : m_n() , m_p_M( nullptr ) {} template inline QuotientRing::QuotientRing( const INT& n , const INT* const & p_M ) noexcept : m_n( p_M == nullptr ? n : n % *p_M ) , m_p_M( p_M ) {} template inline QuotientRing::QuotientRing( const QuotientRing& n ) noexcept : m_n( n.m_n ) , m_p_M( n.m_p_M ) {} template inline QuotientRing& QuotientRing::operator+=( const QuotientRing& n ) noexcept { if( m_p_M == nullptr ){ m_p_M = n.m_p_M; } m_n += n.m_n; if( m_p_M != nullptr ){ m_n %= *m_p_M; } return *this; } template inline QuotientRing& QuotientRing::operator+=( const INT& n ) noexcept { m_n += n; if( m_p_M != nullptr ){ m_n %= *m_p_M; } return *this; } template inline QuotientRing& QuotientRing::operator-=( const QuotientRing& n ) noexcept { if( m_p_M == nullptr ){ m_p_M = n.m_p_M; } m_n -= n.m_n; if( m_p_M != nullptr ){ m_n %= *m_p_M; } return *this; } template inline QuotientRing& QuotientRing::operator-=( const INT& n ) noexcept { m_n -= n; if( m_p_M != nullptr ){ m_n %= *m_p_M; } return *this; } template inline QuotientRing& QuotientRing::operator*=( const QuotientRing& n ) noexcept { if( m_p_M == nullptr ){ m_p_M = n.m_p_M; } m_n *= n.m_n; if( m_p_M != nullptr ){ m_n %= *m_p_M; } return *this; } template inline QuotientRing& QuotientRing::operator*=( const INT& n ) noexcept { m_n *= n; if( m_p_M != nullptr ){ m_n %= *m_p_M; } return *this; } template inline const INT& QuotientRing::Represent() const noexcept { return m_n; } template inline const INT& QuotientRing::GetModulo() const noexcept { static const INT zero{ 0 }; return m_p_M == nullptr ? zero : *m_p_M; } template inline bool QuotientRing::Equal( const QuotientRing& n0 , const QuotientRing& n1 ) noexcept { return n0.m_n == n1.m_n && n0.m_p_M == n1.m_p_M; } template template QuotientRing QuotientRing::Power( const QuotientRing& n , const T& exponent ) { QuotientRing answer{ 1 , n.m_p_M }; QuotientRing power{ n }; while( exponent != 0 ){ if( exponent % 2 == 1 ){ answer *= power; } power *= power; exponent /= 2; } return answer; } template inline bool operator==( const QuotientRing& n0 , const QuotientRing& n1 ) noexcept { return QuotientRing::Equal( n0 , n1 ); } template inline bool operator!=( const QuotientRing& n0 , const QuotientRing& n1 ) noexcept { return ! QuotientRing::Equal( n0 , n1 ); } template inline QuotientRing operator+( const QuotientRing& n0 , const T& n1 ) noexcept { return QuotientRing( n0 ).operator+=( n1 ); } template inline QuotientRing operator-( const QuotientRing& n0 , const T& n1 ) noexcept { return QuotientRing( n0 ).operator-=( n1 ); } template inline QuotientRing operator*( const QuotientRing& n0 , const T& n1 ) noexcept { return QuotientRing( n0 ).operator*=( n1 ); } template inline QuotientRing Power( const QuotientRing& n , const T& exponent ) { return QuotientRing::template Power( n , exponent ); } inline ll sgn( const ll& n ) { return n == 0 ? 0 : n > 0 ? 1 : -1; } inline int Solve() { CEXPR( ll , bound_N , 1000000000000000000 ); CIN_ASSERT( N , 1 , bound_N ); CEXPR( ll , bound_B , 1000000000 ); CIN_ASSERT( B , 1 , bound_B ); CEXPR( ll , bound_Aij , 1000000000 ); CEXPR( int , size , 2 ); ll A00 , A01 , A10 , A11; ll* p_A[size][size] = { { &A00 , &A01 } , { &A10 , & A11 } }; FOR( i , 0 , size ){ FOR( j , 0 , size ){ CIN_ASSERT( Aij , -bound_Aij , bound_Aij ); *( p_A[i][j] ) = Aij; } } ll det_A = A00 * A11 - A01 * A10; ll det_A_minus_E = ( A00 - 1 ) * ( A11 - 1 ) - A01 * A10; ll Delta , sgn_Delta; if( det_A_minus_E != 0 ){ ll tr_A = A00 + A11; ll D = tr_A * tr_A - 4 * det_A; if( D >= 0 ){ ll det_A_plus_E = ( A00 + 1 ) * ( A11 + 1 ) - A01 * A10; sgn_Delta = N % 2 == 0 ? sgn( det_A ) * sgn( det_A_plus_E ) : sgn( det_A ); } else { if( tr_A == 0 && det_A == 1 && N % 4 == 0 ){ sgn_Delta = 0; } else if( tr_A == -1 && det_A == 1 && N % 3 == 0 ){ sgn_Delta = 0; } else { sgn_Delta = sgn( det_A_minus_E ); } } POWER_MOD( det_A_minus_E_inv , det_A_minus_E , B - 2 , B ); QuotientRing A00mod{ A00 , &B }; QuotientRing A01mod{ A01 , &B }; QuotientRing A10mod{ A10 , &B }; QuotientRing A11mod{ A11 , &B }; TwoByTwoMatrix > A{ A00mod , A01mod , A10mod , A11mod }; POWER( power_A , A , N ); QuotientRing one{ 1 , &B }; QuotientRing zero{ 0 , &B }; ll det_power_A_minus_E = ( ( power_A.m_M00 - one ) * ( power_A.m_M11 - one ) - power_A.m_M01 * power_A.m_M10 ).Represent(); Delta = Residue( ( ( det_A * det_power_A_minus_E ) % B ) * det_A_minus_E_inv , B ); } else { if( det_A == 1 ){ N %= B; RETURN( ( N * N ) % B ); } else if( det_A == 0 ){ RETURN( 0 ); } else if( N % 2 == 0 ){ sgn_Delta = det_A == -1 ? 0 : 1; } else { sgn_Delta = sgn( det_A ); } POWER_MOD( det_A_minus_1_inv , det_A - 1 , B - 2 , B ); POWER_MOD( power_det_A , det_A , N , B ); Delta = Residue( ( ( ( ( ( N % B ) * det_A ) % B ) * ( power_det_A - 1 ) ) % B ) * det_A_minus_1_inv , B ); } RETURN( sgn_Delta >= 0 ? Delta : B - Delta ); } int MAIN() { UNTIE; CEXPR( int , bound_T , 1000 ); CIN_ASSERT( T , 1 , bound_T ); REPEAT( T ){ Solve(); } QUIT; }