fork download
copy 
  1. #pragma GCC target("sse,sse2,sse3,ssse3,sse4,popcnt,abm,mmx,avx,avx2")
  2. #include <algorithm>
  3. #include <array>
  4. #include <bitset>
  5. #include <cassert>
  6. #include <climits>
  7. #include <cmath>
  8. #include <functional>
  9. #include <iostream>
  10. #include <limits>
  11. #include <list>
  12. #include <map>
  13. #include <memory>
  14. #include <numeric>
  15. #include <queue>
  16. #include <random>
  17. #include <set>
  18. #include <sstream>
  19. #include <stack>
  20. #include <string.h>
  21. #include <type_traits>
  22. #include <unordered_map>
  23. #if __cplusplus < 201703L
  24. namespace std {
  25. template <typename _Tp>
  26.   constexpr bool is_void_v = is_void<_Tp>::value;
  27. template <typename _Tp>
  28.   constexpr bool is_null_pointer_v = is_null_pointer<_Tp>::value;
  29. template <typename _Tp>
  30.   constexpr bool is_integral_v = is_integral<_Tp>::value;
  31. template <typename _Tp>
  32.   constexpr bool is_floating_point_v = is_floating_point<_Tp>::value;
  33. template <typename _Tp>
  34.   constexpr bool is_array_v = is_array<_Tp>::value;
  35. template <typename _Tp>
  36.   constexpr bool is_pointer_v = is_pointer<_Tp>::value;
  37. template <typename _Tp>
  38.   constexpr bool is_lvalue_reference_v =
  39.     is_lvalue_reference<_Tp>::value;
  40. template <typename _Tp>
  41.   constexpr bool is_rvalue_reference_v =
  42.     is_rvalue_reference<_Tp>::value;
  43. template <typename _Tp>
  44.   constexpr bool is_member_object_pointer_v =
  45.     is_member_object_pointer<_Tp>::value;
  46. template <typename _Tp>
  47.   constexpr bool is_member_function_pointer_v =
  48.     is_member_function_pointer<_Tp>::value;
  49. template <typename _Tp>
  50.   constexpr bool is_enum_v = is_enum<_Tp>::value;
  51. template <typename _Tp>
  52.   constexpr bool is_union_v = is_union<_Tp>::value;
  53. template <typename _Tp>
  54.   constexpr bool is_class_v = is_class<_Tp>::value;
  55. template <typename _Tp>
  56.   constexpr bool is_function_v = is_function<_Tp>::value;
  57. template <typename _Tp>
  58.   constexpr bool is_reference_v = is_reference<_Tp>::value;
  59. template <typename _Tp>
  60.   constexpr bool is_arithmetic_v = is_arithmetic<_Tp>::value;
  61. template <typename _Tp>
  62.   constexpr bool is_fundamental_v = is_fundamental<_Tp>::value;
  63. template <typename _Tp>
  64.   constexpr bool is_object_v = is_object<_Tp>::value;
  65. template <typename _Tp>
  66.   constexpr bool is_scalar_v = is_scalar<_Tp>::value;
  67. template <typename _Tp>
  68.   constexpr bool is_compound_v = is_compound<_Tp>::value;
  69. template <typename _Tp>
  70.   constexpr bool is_member_pointer_v = is_member_pointer<_Tp>::value;
  71. template <typename _Tp>
  72.   constexpr bool is_const_v = is_const<_Tp>::value;
  73. template <typename _Tp>
  74.   constexpr bool is_volatile_v = is_volatile<_Tp>::value;
  75. template <typename _Tp>
  76.   constexpr bool is_trivial_v = is_trivial<_Tp>::value;
  77. template <typename _Tp>
  78.   constexpr bool is_trivially_copyable_v =
  79.     is_trivially_copyable<_Tp>::value;
  80. template <typename _Tp>
  81.   constexpr bool is_standard_layout_v = is_standard_layout<_Tp>::value;
  82. #pragma GCC diagnostic push
  83. #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
  84. template <typename _Tp>
  85.   [[deprecated("use is_standard_layout_v && is_trivial_v instead")]]
  86.   constexpr bool is_pod_v = is_pod<_Tp>::value;
  87. #pragma GCC diagnostic pop
  88. template <typename _Tp>
  89.   constexpr bool is_literal_type_v = is_literal_type<_Tp>::value;
  90. template <typename _Tp>
  91.   constexpr bool is_empty_v = is_empty<_Tp>::value;
  92. template <typename _Tp>
  93.   constexpr bool is_polymorphic_v = is_polymorphic<_Tp>::value;
  94. template <typename _Tp>
  95.   constexpr bool is_abstract_v = is_abstract<_Tp>::value;
  96. template <typename _Tp>
  97.   constexpr bool is_final_v = is_final<_Tp>::value;
  98. template <typename _Tp>
  99.   constexpr bool is_signed_v = is_signed<_Tp>::value;
  100. template <typename _Tp>
  101.   constexpr bool is_unsigned_v = is_unsigned<_Tp>::value;
  102. template <typename _Tp, typename... _Args>
  103.   constexpr bool is_constructible_v =
  104.     is_constructible<_Tp, _Args...>::value;
  105. template <typename _Tp>
  106.   constexpr bool is_default_constructible_v =
  107.     is_default_constructible<_Tp>::value;
  108. template <typename _Tp>
  109.   constexpr bool is_copy_constructible_v =
  110.     is_copy_constructible<_Tp>::value;
  111. template <typename _Tp>
  112.   constexpr bool is_move_constructible_v =
  113.     is_move_constructible<_Tp>::value;
  114. template <typename _Tp, typename _Up>
  115.   constexpr bool is_assignable_v = is_assignable<_Tp, _Up>::value;
  116. template <typename _Tp>
  117.   constexpr bool is_copy_assignable_v = is_copy_assignable<_Tp>::value;
  118. template <typename _Tp>
  119.   constexpr bool is_move_assignable_v = is_move_assignable<_Tp>::value;
  120. template <typename _Tp>
  121.   constexpr bool is_destructible_v = is_destructible<_Tp>::value;
  122. template <typename _Tp, typename... _Args>
  123.   constexpr bool is_trivially_constructible_v =
  124.     is_trivially_constructible<_Tp, _Args...>::value;
  125. template <typename _Tp>
  126.   constexpr bool is_trivially_default_constructible_v =
  127.     is_trivially_default_constructible<_Tp>::value;
  128. template <typename _Tp>
  129.   constexpr bool is_trivially_copy_constructible_v =
  130.     is_trivially_copy_constructible<_Tp>::value;
  131. template <typename _Tp>
  132.   constexpr bool is_trivially_move_constructible_v =
  133.     is_trivially_move_constructible<_Tp>::value;
  134. template <typename _Tp, typename _Up>
  135.   constexpr bool is_trivially_assignable_v =
  136.     is_trivially_assignable<_Tp, _Up>::value;
  137. template <typename _Tp>
  138.   constexpr bool is_trivially_copy_assignable_v =
  139.     is_trivially_copy_assignable<_Tp>::value;
  140. template <typename _Tp>
  141.   constexpr bool is_trivially_move_assignable_v =
  142.     is_trivially_move_assignable<_Tp>::value;
  143. template <typename _Tp>
  144.   constexpr bool is_trivially_destructible_v =
  145.     is_trivially_destructible<_Tp>::value;
  146. template <typename _Tp, typename... _Args>
  147.   constexpr bool is_nothrow_constructible_v =
  148.     is_nothrow_constructible<_Tp, _Args...>::value;
  149. template <typename _Tp>
  150.   constexpr bool is_nothrow_default_constructible_v =
  151.     is_nothrow_default_constructible<_Tp>::value;
  152. template <typename _Tp>
  153.   constexpr bool is_nothrow_copy_constructible_v =
  154.     is_nothrow_copy_constructible<_Tp>::value;
  155. template <typename _Tp>
  156.   constexpr bool is_nothrow_move_constructible_v =
  157.     is_nothrow_move_constructible<_Tp>::value;
  158. template <typename _Tp, typename _Up>
  159.   constexpr bool is_nothrow_assignable_v =
  160.     is_nothrow_assignable<_Tp, _Up>::value;
  161. template <typename _Tp>
  162.   constexpr bool is_nothrow_copy_assignable_v =
  163.     is_nothrow_copy_assignable<_Tp>::value;
  164. template <typename _Tp>
  165.   constexpr bool is_nothrow_move_assignable_v =
  166.     is_nothrow_move_assignable<_Tp>::value;
  167. template <typename _Tp>
  168.   constexpr bool is_nothrow_destructible_v =
  169.     is_nothrow_destructible<_Tp>::value;
  170. template <typename _Tp>
  171.   constexpr bool has_virtual_destructor_v =
  172.     has_virtual_destructor<_Tp>::value;
  173. template <typename _Tp>
  174.   constexpr size_t alignment_of_v = alignment_of<_Tp>::value;
  175. template <typename _Tp>
  176.   constexpr size_t rank_v = rank<_Tp>::value;
  177. template <typename _Tp, unsigned _Idx = 0>
  178.   constexpr size_t extent_v = extent<_Tp, _Idx>::value;
  179. #ifdef _GLIBCXX_BUILTIN_IS_SAME_AS
  180. template <typename _Tp, typename _Up>
  181.   constexpr bool is_same_v = _GLIBCXX_BUILTIN_IS_SAME_AS(_Tp, _Up);
  182. #else
  183. template <typename _Tp, typename _Up>
  184.   constexpr bool is_same_v = std::is_same<_Tp, _Up>::value;
  185. #endif
  186. template <typename _Base, typename _Derived>
  187.   constexpr bool is_base_of_v = is_base_of<_Base, _Derived>::value;
  188. template <typename _From, typename _To>
  189.   constexpr bool is_convertible_v = is_convertible<_From, _To>::value;
  190. }
  191. #endif
  192. namespace OY {
  193. #define cin OY::inputHelper<1024, 20>::getInstance()
  194. #define getchar() ({char c=cin.getChar_Checked();cin.next();c;})
  195. #define cout OY::outputHelper<1048576>::getInstance()
  196. #define putchar cout.putChar
  197. #define endl '\n'
  198. #define putlog(...) OY::printLog(", ", __VA_ARGS__)
  199.     template <uint64_t _BufferSize = 1024, uint64_t _BlockSize = 20>
  200.     class inputHelper {
  201.         FILE *m_filePtr;
  202.         char m_buf[_BufferSize], *m_end, *m_cursor;
  203.         bool m_ok;
  204.         void flush() {
  205.             uint64_t a = m_end - m_cursor;
  206.             if (a >= _BlockSize) return;
  207.             memmove(m_buf, m_cursor, a);
  208.             uint64_t b = fread(m_buf + a, 1, _BufferSize - a, m_filePtr);
  209.             m_cursor = m_buf;
  210.             if (a + b < _BufferSize) {
  211.                 m_end = m_buf + a + b;
  212.                 *m_end = EOF;
  213.             }
  214.         }
  215.  
  216.     public:
  217.         explicit inputHelper(const char *inputFileName) : m_ok(true) {
  218.             if (!*inputFileName)
  219.                 m_filePtr = stdin;
  220.             else
  221.                 m_filePtr = fopen(inputFileName, "rt");
  222.             m_end = m_cursor = m_buf + _BufferSize;
  223.         }
  224.         ~inputHelper() { fclose(m_filePtr); }
  225.         static inputHelper<_BufferSize, _BlockSize> &getInstance() {
  226. #ifdef OY_LOCAL
  227.             static inputHelper<_BufferSize, _BlockSize> s_obj("in.txt");
  228. #else
  229.             static inputHelper<_BufferSize, _BlockSize> s_obj("");
  230. #endif
  231.             return s_obj;
  232.         }
  233.         static constexpr bool isBlank(char c) { return c == ' ' || c == '\t' || c == '\n' || c == '\r'; }
  234.         static constexpr bool isEndline(char c) { return c == '\n' || c == EOF; }
  235.         const char &getChar_Checked() {
  236.             if (m_cursor < m_end) return *m_cursor;
  237.             uint64_t b = fread(m_buf, 1, _BufferSize, m_filePtr);
  238.             m_cursor = m_buf;
  239.             if (b < _BufferSize) {
  240.                 m_end = m_buf + b;
  241.                 *m_end = EOF;
  242.             }
  243.             return *m_cursor;
  244.         }
  245.         const char &getChar_Unchecked() const { return *m_cursor; }
  246.         void next() { ++m_cursor; }
  247.         void setState(bool _ok) { m_ok = _ok; }
  248.         template <typename _Tp, std::enable_if_t<std::is_signed_v<_Tp> & std::is_integral_v<_Tp>> * = nullptr>
  249.         inputHelper<_BufferSize, _BlockSize> &operator>>(_Tp &ret) {
  250.             while (isBlank(getChar_Checked())) next();
  251.             flush();
  252.             if (getChar_Unchecked() == '-') {
  253.                 next();
  254.                 if (isdigit(getChar_Unchecked())) {
  255.                     ret = -(getChar_Unchecked() - '0');
  256.                     while (next(), isdigit(getChar_Unchecked())) ret = ret * 10 - (getChar_Unchecked() - '0');
  257.                 } else
  258.                     m_ok = false;
  259.             } else {
  260.                 if (isdigit(getChar_Unchecked())) {
  261.                     ret = getChar_Unchecked() - '0';
  262.                     while (next(), isdigit(getChar_Unchecked())) ret = ret * 10 + (getChar_Unchecked() - '0');
  263.                 } else
  264.                     m_ok = false;
  265.             }
  266.             return *this;
  267.         }
  268.         template <typename _Tp, std::enable_if_t<std::is_unsigned_v<_Tp> & std::is_integral_v<_Tp>> * = nullptr>
  269.         inputHelper<_BufferSize, _BlockSize> &operator>>(_Tp &ret) {
  270.             while (isBlank(getChar_Checked())) next();
  271.             flush();
  272.             if (isdigit(getChar_Unchecked())) {
  273.                 ret = getChar_Unchecked() - '0';
  274.                 while (next(), isdigit(getChar_Unchecked())) ret = ret * 10 + (getChar_Unchecked() - '0');
  275.             } else
  276.                 m_ok = false;
  277.             return *this;
  278.         }
  279.         template <typename _Tp, std::enable_if_t<std::is_floating_point_v<_Tp>> * = nullptr>
  280.         inputHelper<_BufferSize, _BlockSize> &operator>>(_Tp &ret) {
  281.             bool neg = false, integer = false, decimal = false;
  282.             while (isBlank(getChar_Checked())) next();
  283.             flush();
  284.             if (getChar_Unchecked() == '-') {
  285.                 neg = true;
  286.                 next();
  287.             }
  288.             if (!isdigit(getChar_Unchecked()) && getChar_Unchecked() != '.') {
  289.                 m_ok = false;
  290.                 return *this;
  291.             }
  292.             if (isdigit(getChar_Unchecked())) {
  293.                 integer = true;
  294.                 ret = getChar_Unchecked() - '0';
  295.                 while (next(), isdigit(getChar_Unchecked())) ret = ret * 10 + (getChar_Unchecked() - '0');
  296.             }
  297.             if (getChar_Unchecked() == '.') {
  298.                 next();
  299.                 if (isdigit(getChar_Unchecked())) {
  300.                     if (!integer) ret = 0;
  301.                     decimal = true;
  302.                     _Tp unit = 0.1;
  303.                     ret += unit * (getChar_Unchecked() - '0');
  304.                     while (next(), isdigit(getChar_Unchecked())) {
  305.                         unit *= 0.1;
  306.                         ret += unit * (getChar_Unchecked() - '0');
  307.                     }
  308.                 }
  309.             }
  310.             if (!integer && !decimal)
  311.                 m_ok = false;
  312.             else if (neg)
  313.                 ret = -ret;
  314.             return *this;
  315.         }
  316.         inputHelper<_BufferSize, _BlockSize> &operator>>(char &ret) {
  317.             while (isBlank(getChar_Checked())) next();
  318.             ret = getChar_Checked();
  319.             next();
  320.             return *this;
  321.         }
  322.         inputHelper<_BufferSize, _BlockSize> &operator>>(std::string &ret) {
  323.             while (isBlank(getChar_Checked())) next();
  324.             if (getChar_Checked() != EOF) {
  325.                 ret.clear();
  326.                 do {
  327.                     ret += getChar_Checked();
  328.                     next();
  329.                 } while (!isBlank(getChar_Checked()) && getChar_Unchecked() != EOF);
  330.             } else
  331.                 m_ok = false;
  332.             return *this;
  333.         }
  334.         explicit operator bool() { return m_ok; }
  335.     };
  336.     template <uint64_t _BufferSize = 1048576>
  337.     class outputHelper {
  338.         FILE *m_filePtr = nullptr;
  339.         char m_buf[_BufferSize], *m_end, *m_cursor;
  340.         char m_tempBuf[50], *m_tempBufCursor, *m_tempBufDot;
  341.         uint64_t m_floatReserve, m_floatRatio;
  342.  
  343.     public:
  344.         outputHelper(const char *outputFileName, int prec = 6) : m_end(m_buf + _BufferSize) {
  345.             if (!*outputFileName)
  346.                 m_filePtr = stdout;
  347.             else
  348.                 m_filePtr = fopen(outputFileName, "wt");
  349.             m_cursor = m_buf;
  350.             m_tempBufCursor = m_tempBuf;
  351.             precision(prec);
  352.         }
  353.         static outputHelper<_BufferSize> &getInstance() {
  354. #ifdef OY_LOCAL
  355.             static outputHelper<_BufferSize> s_obj("out.txt");
  356. #else
  357.             static outputHelper<_BufferSize> s_obj("");
  358. #endif
  359.             return s_obj;
  360.         }
  361.         ~outputHelper() {
  362.             flush();
  363.             fclose(m_filePtr);
  364.         }
  365.         void precision(int prec) {
  366.             m_floatReserve = prec;
  367.             m_floatRatio = pow(10, prec);
  368.             m_tempBufDot = m_tempBuf + prec;
  369.         }
  370.         outputHelper<_BufferSize> &flush() {
  371.             fwrite(m_buf, 1, m_cursor - m_buf, m_filePtr);
  372.             fflush(m_filePtr);
  373.             m_cursor = m_buf;
  374.             return *this;
  375.         }
  376.         void putChar(const char &c) {
  377.             if (m_cursor == m_end) flush();
  378.             *m_cursor++ = c;
  379.         }
  380.         void putS(const char *c) {
  381.             while (*c) putChar(*c++);
  382.         }
  383.         template <typename _Tp, std::enable_if_t<std::is_signed_v<_Tp> & std::is_integral_v<_Tp>> * = nullptr>
  384.         outputHelper<_BufferSize> &operator<<(const _Tp &ret) {
  385.             _Tp _ret = _Tp(ret);
  386.             if (_ret >= 0) {
  387.                 do {
  388.                     *m_tempBufCursor++ = '0' + _ret % 10;
  389.                     _ret /= 10;
  390.                 } while (_ret);
  391.                 do putChar(*--m_tempBufCursor);
  392.                 while (m_tempBufCursor > m_tempBuf);
  393.             } else {
  394.                 putChar('-');
  395.                 do {
  396.                     *m_tempBufCursor++ = '0' - _ret % 10;
  397.                     _ret /= 10;
  398.                 } while (_ret);
  399.                 do putChar(*--m_tempBufCursor);
  400.                 while (m_tempBufCursor > m_tempBuf);
  401.             }
  402.             return *this;
  403.         }
  404.         template <typename _Tp, std::enable_if_t<std::is_unsigned_v<_Tp> & std::is_integral_v<_Tp>> * = nullptr>
  405.         outputHelper<_BufferSize> &operator<<(const _Tp &ret) {
  406.             _Tp _ret = _Tp(ret);
  407.             do {
  408.                 *m_tempBufCursor++ = '0' + _ret % 10;
  409.                 _ret /= 10;
  410.             } while (_ret);
  411.             do putChar(*--m_tempBufCursor);
  412.             while (m_tempBufCursor > m_tempBuf);
  413.             return *this;
  414.         }
  415.         template <typename _Tp, std::enable_if_t<std::is_floating_point_v<_Tp>> * = nullptr>
  416.         outputHelper<_BufferSize> &operator<<(const _Tp &ret) {
  417.             if (ret < 0) {
  418.                 putChar('-');
  419.                 return *this << -ret;
  420.             }
  421.             _Tp _ret = ret * m_floatRatio;
  422.             uint64_t integer = _ret;
  423.             if (_ret - integer >= 0.4999999999) integer++;
  424.             do {
  425.                 *m_tempBufCursor++ = '0' + integer % 10;
  426.                 integer /= 10;
  427.             } while (integer);
  428.             if (m_tempBufCursor > m_tempBufDot) {
  429.                 do putChar(*--m_tempBufCursor);
  430.                 while (m_tempBufCursor > m_tempBufDot);
  431.                 putChar('.');
  432.             } else {
  433.                 putS("0.");
  434.                 for (int i = m_tempBufDot - m_tempBufCursor; i--;) putChar('0');
  435.             }
  436.             do putChar(*--m_tempBufCursor);
  437.             while (m_tempBufCursor > m_tempBuf);
  438.             return *this;
  439.         }
  440.         outputHelper<_BufferSize> &operator<<(const char &ret) {
  441.             putChar(ret);
  442.             return *this;
  443.         }
  444.         outputHelper<_BufferSize> &operator<<(const std::string &ret) {
  445.             putS(ret.data());
  446.             return *this;
  447.         }
  448.     };
  449.     template <uint64_t _BufferSize, uint64_t _BlockSize>
  450.     inputHelper<_BufferSize, _BlockSize> &getline(inputHelper<_BufferSize, _BlockSize> &ih, std::string &ret) {
  451.         ret.clear();
  452.         if (ih.getChar_Checked() == EOF)
  453.             ih.setState(0);
  454.         else {
  455.             while (!inputHelper<_BufferSize, _BlockSize>::isEndline(ih.getChar_Checked())) {
  456.                 ret += ih.getChar_Unchecked();
  457.                 ih.next();
  458.             }
  459.             ih.next();
  460.         }
  461.         return ih;
  462.     }
  463.     template <typename D, typename T, typename... S>
  464.     void printLog(D delim, const T &x, S... rest) {
  465.         cout << x;
  466.         if constexpr (sizeof...(rest) > 0) {
  467.             cout << delim;
  468.             printLog(delim, rest...);
  469.         }
  470.     }
  471. }
  472. using OY::getline;
  473. using std::vector;
  474. #define int long long
  475. const int MAXN = 524291;
  476. const long double pi = 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989380952572010654858632788659361533818279682303019520353018529689957736225994138912497217752834791315155748572424541506959508295331168617278558890750983817546374649393192550604009277016711390098488240128583616035637076601047101819429555961989467678374494482553797747268471040475346462080466842590694912933136770289891521047521620569660240580381501935112533824300355876402474964732639141992726042699227967823547816360093417216412199245863150302861829745557067498385054945885869269956909272107975093029553211653449872027559602364806654991198818347977535663698074265425278625518184175746728909777727938000816470600161452491921732172147723501414419735685481613611573525521334757418494684385233239073941433345477624168625189835694855620992192221842725502542568876717904946016534668049886272327917860857843838279679766814541009538837863609506800642251252051173929848960841284886269456042419652850222106611863067442786220391949450471237137869609563643719172874677646575739624138908658326459958133904780275900994657640789512694683983525957098258226205224894077267194782684826014769909026401363944374553050682034962524517493996514314298091906592509372216964615157098583874105978859597729754989301617539284681382686838689427741559918559252459539594310499725246808459872736446958486538367362226260991246080512438843904512441365497627807977156914359977001296160894416948685558484063534220722258284886481584560285060168427394522674676788952521385225499546667278239864565961163548862305774564980355936345681743241125150760694794510965960940252288797108931456691368672287489405601015033086179286809208747609178249385890097149096759852613655497818931297848216829989487226588048575640142704775551323796414515237462343645428584447952658678210511413547357395231134271661021359695362314429524849371871101457654035902799344037420073105785390621983874478084784896833214457138687519435064302184531910484810053706146806749192781911979399520614196634287544406437451237181921799983910159195618146751426912397489409071864942319615679452080951465502252316038819301420937621378559566389377870830390697920773467221825625996615014215030680384477345492026054146659252014974428507325186660021324340881907104863317346496514539057962685610055081066587969981635747363840525714591028970641401109712062804390397595156771577004203378699360072305587631763594218731251471205329281918261861258673215791984148488291644706095752706957220917567116722910981690915280173506712748583222871835209353965725121083579151369882091444210067510334671103141267111369908658516398315019701651511685171437657618351556508849099898599823873455283316355076479185358932261854896321329330898570642046752590709154814165498594616371802709819943099244889575712828905923233260972997120844335732654893823911932597463667305836041428138830320382490375898524374417029132765618093773444030707469211201913020330380197621101100449293215160842444859637669838952286847831235526582131449576857262433441893039686426243410773226978028073189154411010446823252716201052652272111660396665573092547110557853763466820653109896526918620564769312570586356620185581007293606598764861179104533488503461136576867532494416680396265797877185560845529654126654085306143444318586769751456614068007002378776591344017127494704205622305389945613140711270004078547332699390814546646458807972708266830634328587856983052358089330657574067954571637752542021149557615814002501262285941302164715509792592309907965473761255176567513575178296664547791745011299614890304639947132962107340437518957359614589019389713111790429782856475032031986915140287080859904801094121472213179476477726224142548545403321571853061422881375850430633217518297986622371721591607716692547487389866549494501146540628433663937900397692656721463853067360965712091807638327166416274888800786925602902284721040317211860820419000422966171196377921337575114959501566049631862947265473642523081770367515906735023507283540567040386743513622224771589150495309844489333096340878076932599397805419341447377441842631298608099888687413260472156951623965864573021631598193195167353812974167729478672422924654366800980676928238280689964004824354037014163149658979409243237896907069779422362508221688957383798623001593776471651228935786015881617557829735233446042815126272037343146531977774160319906655418763979293344195215413418994854447345673831624993419131814809277771038638773431772075456545322077709212019051660962804909263601975988281613323166636528619326686336062735676303544776280350450777235547105859548702790814356240145171806246436267945612753181340783303362542327839449753824372058353114771199260638133467768796959703098339130771098704085913374641442822772634659470474587847787201927715280731767907707157213444730605700733492436931138350493163128404251219256517980694113528013147013047816437885185290928545201165839341965621349143415956258658655705526904965209858033850722426482939728584783163057777560688876446248246857926039535277348030480290058760758251047470916439613626760449256274204208320856611906254543372131535958450687724602901618766795240616342522577195429162991930645537799140373404328752628889639958794757291746426357455254079091451357111369410911939325191076020825202618798531887705842972591677813149699009019211697173727847684726860849003377024242916513005005168323364350389517029893922334517220138128069650117844087451960121228599371623130171144484640903890644954440061986907548516026327505298349187407866808818338510228334508504860825039302133219715518430635455007668282949304137765527939751754613953984683393638304746119966538581538420568533862186725233402830871123282789212507712629463229563989898935821167456270102183564622013496715188190973038119800497340723961036854066431939509790190699639552453005450580685501956730229219139339185680344903982059551002263535361920419947455385938102343955449597783779023742161727111723643435439478221818528624085140066604433258885698670543154706965747458550332323342107301545940516553790686627333799585115625784322988273723198987571415957811196358330059408730681216028764962867446047746491599505497374256269010490377819868359381465741268049256487985561453723478673303904688383436346553794986419270563872931748723320837601123029911367938627089438799362016295154133714248928307220126901475466847653576164773794675200490757155527819653621323926406160136358155907422020203187277605277219005561484255518792530343513984425322341576233610642506390497500865627109535919465897514131034822769306247435363256916078154781811528436679570611086153315044521274739245449454236828860613408414863776700961207151249140430272538607648236341433462351897576645216413767969031495019108575984423919862916421939949072362346468441173940326591840443780513338945257423995082965912285085558215725031071257012668302402929525220118726767562204154205161841634847565169998116141010029960783869092916030288400269104140792886215078424516709087000699282120660418371806535567252532567532861291042487761825829765157959847035622262934860034158722980534989650226291748788202734209222245339856264766914905562842503912757710284027998066365825488926488025456610172967026640765590429099456815065265305371829412703369313785178609040708667114965583434347693385781711386455873678123014587687126603489139095620099393610310291616152881384379099042317473363948045759314931405297634757481193567091101377517210080315590248530906692037671922033229094334676851422144773793937517034436619910403375111735471918550464490263655128162288244625759163330391072253837421821408835086573917715096828874782656995995744906617583441375223970968340800535598491754173818839994469748676265516582765848358845314277568790029095170283529716344562129640435231176006651012412006597558512761785838292041974844236080071930457618932349229279650198751872127267507981255470958904556357921221033346697499235630254947802490114195212382815309114079073860251522742995818072471625916685451333123948049470791191532673430282441860414263639548000448002670496248201792896476697583183271314251702969234889627668440323260927524960357996469256504936818360900323809293459588970695365349406034021665443755890045632882250545255640564482465151875471196218443965825337543885690941130315095261793780029741207665147939425902989695946995565761218656196733786236256125216320862869222103274889218654364802296780705765615144632046927906821207388377814233562823608963208068222468012248261177185896381409183903673672220888321513755600372798394004152970028783076670944474560134556417254370906979396122571429894671543578468788614445812314593571984922528471605049221242470141214780573455105008019086996033027634787081081754501193071412233908663938339529425786905076431006383519834389341596131854347546495569781038293097164651438407007073604112373599843452251610507027056235266012764848308407611830130527932054274628654036036745328651057065874882256981579367897669742205750596834408697350201410206723585020072452256326513410559240190274216248439140359989535394590944070469120914093870012645600162374288021092764579310657922955249887275846101264836999892256959688159205600101655256375678;
  477. #if 0
  478. using cp = std::complex<long double>;
  479. #else
  480. template<typename T>
  481. struct cp_base {
  482.     T m_real, m_imag;
  483.     constexpr cp_base(T r, T i) : m_real(r), m_imag(i) {}
  484.     constexpr cp_base() {}
  485.     constexpr cp_base(T r) : m_real(r), m_imag(0) {}
  486.     constexpr cp_base operator+(const cp_base &o) const {
  487.         return cp_base(m_real + o.m_real, m_imag + o.m_imag);
  488.     }
  489.     constexpr cp_base operator-(const cp_base &o) const {
  490.         return cp_base(m_real - o.m_real, m_imag - o.m_imag);
  491.     }
  492.     constexpr cp_base operator*(const cp_base &o) const {
  493.         return cp_base(m_real * o.m_real - m_imag * o.m_imag, 
  494.                   m_real * o.m_imag + m_imag * o.m_real);
  495.     }
  496.     constexpr cp_base operator+=(const cp_base &o) {
  497.         return *this = *this + o;
  498.     }
  499.     constexpr cp_base operator-=(const cp_base &o) {
  500.         return *this = *this - o;
  501.     }
  502.     constexpr cp_base operator*=(const cp_base &o) {
  503.         return *this = *this * o;
  504.     }
  505.     constexpr cp_base conj() const {return cp_base(m_real, -m_imag); }
  506.     constexpr T real() const {return m_real; }
  507.     constexpr T imag() const {return m_imag; }
  508.     constexpr T& real() {return m_real; }
  509.     constexpr T& imag() {return m_imag; }
  510. };
  511. using cp = cp_base<long double>;
  512. #endif
  513. int n, len;
  514. vector<int> rev;
  515. inline constexpr int fpow(int base, int power, const int& mod) {
  516.     int ans = 1;
  517.     for (; power; base = 1ll * base * base % mod, power >>= 1) 
  518.         if (power & 1) ans = 1ll * ans * base % mod;
  519.     return ans;
  520. }
  521. inline void init() {
  522.     int x = 0;
  523.     for (len = 1; (len <<= 1) < n; x++);
  524.     len <<= 1;
  525.     rev.clear();
  526.     rev.resize(len);
  527.     for (int i = 1; i < len; i++)
  528.         rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << x);
  529. }
  530. inline void init(int l) {
  531.     rev.clear();
  532.     rev.resize(l);
  533.     for (int i = 1; i < l; i++) {
  534.         rev[i] = (rev[i >> 1] >> 1) | ((i & 1) * (l >> 1));
  535.     }
  536. }
  537. inline void FFT(vector<cp>& a, int n, int type)
  538. {
  539.     for (int i = 0; i < n; ++i)
  540.         if (i < rev[i])
  541.             std::swap(a[i], a[rev[i]]);
  542.     for (int i = 1; i < n; i <<= 1)
  543.     {
  544.         cp wn(std::cos(pi / i), type * std::sin(pi / i));
  545.         for (int j = 0; j < n; j += (i << 1))
  546.         {
  547.             cp w0(1);
  548.             for (int k = 0; k < i; ++k, w0 *= wn)
  549.             {
  550.                 cp x = a[j + k], y = w0 * a[i + j + k];
  551.                 a[j + k] = x + y;
  552.                 a[i + j + k] = x - y;
  553.             }
  554.         }
  555.     }
  556. }
  557. inline void mul(int n, int m, const int& p, vector<int> a, vector<int> b, vector<int>& c) {
  558.     static cp da, db, dc, dd, rel(0.5, 0), img(0, -0.5);
  559.     static const int up = 1 << 15;
  560.     static int v1, v2, v3;
  561.     int len = 1;
  562.     while (len < n + m) len <<= 1;
  563.     a.resize(len);
  564.     b.resize(len);
  565.     rev.clear();
  566.     rev.resize(len);
  567.     for (int i = 1; i < len; i++) {
  568.     	rev[i] = (rev[i >> 1] >> 1) | ((i & 1) * (len >> 1));
  569. 	}
  570.     for (int i = 0; i < len; i++) {
  571.         a[i] = i < n ? (a[i] % p + p) % p : 0;
  572.         b[i] = i < m ? (b[i] % p + p) % p : 0;
  573.     }
  574.     vector<cp> P(len, cp(0, 0)), Q(len, cp(0, 0)), R(len, cp(0, 0));
  575.     for (int i = 0; i < len; i++) {
  576.         P[i] = cp(a[i] >> 15, a[i] & 32767);
  577.         Q[i] = cp(b[i] >> 15, b[i] & 32767);
  578.     }
  579.     FFT(P, len, 1);
  580.     FFT(Q, len, 1);
  581.     for (int i = 0; i < len; i++) {
  582.         R[i] = P[(len - i) & (len - 1)].conj();
  583.     }
  584.     for (int i = 0; i < len; i++) {
  585.         auto x = P[i];
  586.         auto y = R[i];
  587.         P[i] = (x + y) * rel;
  588.         R[i] = (x - y) * img;
  589.     }
  590.     for (int i = 0; i < len; i++) {
  591.         P[i] = P[i] * Q[i];
  592.         Q[i] = R[i] * Q[i];
  593.     }
  594.     FFT(P, len, -1);
  595.     FFT(Q, len, -1);
  596.     c.clear();
  597.     c.resize(len);
  598.     for (int i = 0; i < len; i++) {
  599.         v1 = (int)(P[i].real() / len + 0.5);
  600.         v2 = (int)(P[i].imag() / len + 0.5) + (int)(Q[i].real() / len + 0.5);
  601.         v3 = (int)(Q[i].imag() / len + 0.5);
  602.         c[i] = 1ll * v1 * up % p * up % p + 1ll * v2 * up % p + v3 % p;
  603.         c[i] = (c[i] % p + p) % p;
  604.     }
  605. }
  606.  
  607. static vector<int> fac, mfac, invfac, minvfac, inv, minv;
  608.  
  609. inline void init(int n, int m, const int &p) {
  610.     int tmp = (n << 1) | 1;
  611.     fac.resize(tmp + 1);
  612.     mfac.resize(tmp + 1);
  613.     invfac.resize(tmp + 1);
  614.     minvfac.resize(tmp + 1);
  615.     inv.resize(tmp + 1);
  616.     minv.resize(tmp + 1);
  617.     fac[0] = mfac[0] = 1;
  618.  
  619.     for (int i = 1; i <= tmp; i++) {
  620.         fac[i] = 1ll * fac[i - 1] * i % p;
  621.         mfac[i] = 1ll * mfac[i - 1] * (m - n + i - 1 + p) % p;
  622.     }
  623.  
  624.     invfac[tmp] = fpow(fac[tmp], p - 2, p);
  625.     minvfac[tmp] = fpow(mfac[tmp], p - 2, p);
  626.  
  627.     for (int i = tmp; i; i--) {
  628.         invfac[i - 1] = 1ll * invfac[i] * i % p;
  629.         minvfac[i - 1] = 1ll * minvfac[i] * (m - n + i - 1 + p) % p;
  630.         inv[i] = 1ll * invfac[i] * fac[i - 1] % p;
  631.         minv[i] = 1ll * minvfac[i] * mfac[i - 1] % p;
  632.     }
  633.  
  634.     minv[0] = 1;
  635. }
  636.  
  637. inline void LagrangeInterpolation_ex(int n, int m, const int &p, vector<int> a, vector<int> &b) {
  638.     vector<int> f, g;
  639.     int len = 1;
  640.     while (len < n * 3ll + 2) len <<= 1;
  641.     rev.clear();
  642.     rev.resize(len);
  643.     a.resize(len);
  644.     f.resize(len);
  645.     g.resize(len);
  646.     init(n, m, p);
  647.     rev.clear();
  648.     rev.resize(len);
  649.     for (int i = 1; i < len; i++) {
  650.     	rev[i] = (rev[i >> 1] >> 1) | ((i & 1) * (len >> 1));
  651. 	}
  652.  
  653.     for (int i = 0; i <= n; i++) {
  654.         f[i] = 1ll * invfac[i] * invfac[n - i] % p * a[i] % p;
  655.  
  656.         if ((n - i) & 1) {
  657.         	f[i] = p - f[i];
  658. 		}
  659.     }
  660.  
  661.     for (int i = 0; i <= (n << 1); i++) {
  662.         g[i] = minv[i + 1];
  663.     }
  664. // n + 1, (n << 1) + 1
  665.     mul(n + 1, (n << 1) + 1, p, f, g, f);
  666. 	b.clear();
  667. 	b.resize(n + 1);
  668.     for (int i = n; i <= (n << 1); i++) {
  669.         b[i - n] = 1ll * mfac[i + 1] * minvfac[i - n] % p * f[i] % p;
  670.     }
  671. }
  672. /*
  673.  
  674. signed main() {
  675. 	int n, m;
  676. 	vector<int> a, b;
  677. 	scanf("%lld%lld", &n, &m);
  678. 	a.resize(n + 1);
  679. 	for (int i = 0; i <= n; i++) scanf("%lld", &a[i]);
  680. 	LagrangeInterpolation_ex(n, m, 998244353, a, b);
  681. 	for (int i = 0; i <= n; i++) printf("%lld ", b[i]);
  682. 	return 0;
  683. }
  684.  
  685. */
  686.  
  687. inline int solve(int n, const int &p) {
  688. 	static vector<int> f, fd, st;
  689. 	st.resize(log2(n) + 5);
  690. 	f.resize(2);
  691. 	int top = 0;
  692. 	int s = n;
  693. 	while (n) {
  694. 		st[++top] = n;
  695. 		n >>= 1;
  696. 	}
  697. 	n = st[top--];
  698. 	f[0] = 1;
  699. 	f[1] = s + 1;
  700. 	while (top--) {
  701. 		f.resize(f.size() << 1);
  702. 		LagrangeInterpolation_ex(n, n + 1, p, f, fd);
  703. 		std::copy(fd.begin(), fd.end(), f.begin() + n + 1);
  704. 		f[n << 1 | 1] = 0;
  705. 		int tmp = 1ll * n * fpow(s, p - 2, p) % p;
  706. 		n <<= 1;
  707. 		LagrangeInterpolation_ex(n, tmp, p, f, fd);
  708. 		for (int i = 0; i <= n; i++) {
  709. 			f[i] = 1ll * f[i] * fd[i] % p;
  710. 		}
  711. 		if (!(st[top + 1] & 1)) continue;
  712. 		for (int i = 0; i <= n; i++) {
  713. 			f[i] = 1ll * f[i] * (1ll * s * i % p + n + 1) % p;
  714. 		}
  715. 		f[++n] = 1;
  716. 		for (int i = 1; i <= n; i++) {
  717. 			f[n] = 1ll * f[n] * (1ll * s * n % p + i) % p;
  718. 		}
  719. 	}
  720. 	int res = f[0];
  721. 	for (int i = 1; i < s; i++) {
  722. 		res = 1ll * res * f[i] % p;
  723. 	}
  724. 	return res;
  725. }
  726. inline int factorial(int n, const int &p) {
  727. 	int tn = p - 1 - n;
  728. 	int res = 0;
  729. 	if (tn <= n) {
  730. 		res = fpow(factorial(tn, p), p - 2, p);
  731. 		return (tn & 1) ? res : p - res;
  732. 	}
  733. 	int k = sqrt(n);
  734. 	res = solve(k, p);
  735. 	for (int i = k * k + 1; i <= n; i++) {
  736. 		res = 1ll * res * i % p;
  737. 	}
  738. 	return res % p;
  739. }
  740. signed main() {
  741. 	int T, n, p;
  742. 	cin >> T;
  743. 	while (T--) {
  744. 		cin >> n >> p;
  745. 		cout << factorial(n, p) << endl;
  746. 	}
  747. }
Runtime error #stdin #stdout #stderr 0.5s 27000KB
comments (?)
stdin
copy 
3
1073741823 2147483629
1073741823 2147481811
2147483628 2147483629
stdout
copy 
Standard output is empty
stderr
copy 
terminate called after throwing an instance of 'std::length_error'
  what():  vector::_M_default_append
https://ideone.com/kNhga9
hack for FACTMODP
language:
C++14 (gcc 8.3)
created:
3 years ago
visibility:
secret

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!