fork download
copy 
  1. #pragma GCC optimize ("O3")
  2. // #pragma GCC target ("avx")
  3. #include<stdio.h>
  4. #include<iostream>
  5. #include<vector>
  6. #include<algorithm>
  7. #include<string>
  8. #include<string.h>
  9.  
  10. #ifdef LOCAL
  11. #define eprintf(...) fprintf(stderr, __VA_ARGS__)
  12. #else
  13. #define NDEBUG
  14. #define eprintf(...) do {} while (0)
  15. #endif
  16. #include<cassert>
  17. #include<random>
  18.  
  19. using namespace std;
  20.  
  21. typedef long long LL;
  22. typedef vector<int> VI;
  23. using ULL = unsigned long long;
  24.  
  25. #define REP(i,n) for(int i=0, i##_len=(n); i<i##_len; ++i)
  26. #define EACH(i,c) for(__typeof((c).begin()) i=(c).begin(),i##_end=(c).end();i!=i##_end;++i)
  27.  
  28. template<class T> inline void amin(T &x, const T &y) { if (y<x) x=y; }
  29. template<class T> inline void amax(T &x, const T &y) { if (x<y) x=y; }
  30. template<class Iter> void rprintf(const char *fmt, Iter begin, Iter end) {
  31.     for (bool sp=0; begin!=end; ++begin) { if (sp) putchar(' '); else sp = true; printf(fmt, *begin); }
  32.     putchar('\n');
  33. }
  34.  
  35. struct MontgomeryReduction {
  36.     ULL MOD;
  37.     // MOD * NEG_INV % (2^32) == (2^32) - 1;
  38.     ULL NEG_INV;
  39.     // R2 == (2^64) % MOD;
  40.     ULL R2;
  41.     MontgomeryReduction(ULL MOD_): MOD(MOD_) {
  42. 	NEG_INV = 0;
  43. 	ULL s = 1, t = 0;
  44. 	REP (i, 32) {
  45. 	    if (~t & 1) {
  46. 		t += MOD;
  47. 		NEG_INV += s;
  48. 	    }
  49. 	    t /= 2;
  50. 	    s *= 2;
  51. 	}
  52.  
  53. 	R2 = (1ULL<<32) % MOD;
  54. 	R2 = R2 * R2 % MOD;
  55.     }
  56.  
  57.     // return x * R % MOD;
  58.     ULL generate(ULL x) const {
  59. 	assert(x < MOD);
  60. 	return reduce(x * R2);
  61.     }
  62.  
  63.     // return x / R % MOD;
  64.     ULL reduce(ULL x) const {
  65. 	assert(x < MOD * MOD);
  66. 	x = (x + ((unsigned)x * (unsigned)NEG_INV) * MOD) >> 32;
  67. 	return x < MOD? x: x-MOD;
  68.     }
  69. };
  70.  
  71.  
  72. struct BarrettReduction {
  73.     ULL MOD;
  74.     __int128 m;
  75.  
  76.     BarrettReduction(ULL MOD_): MOD(MOD_) {
  77. 	m = (__int128(1)<<64) / MOD;
  78.     }
  79.  
  80.     ULL reduce(ULL x) const {
  81. 	ULL q = (x*m)>>64;
  82. 	x -= q * MOD;
  83. 	return x < MOD? x: x-MOD;
  84.     }
  85. };
  86.  
  87. struct BarrettReduction64 {
  88.     ULL MOD;
  89.     ULL mh, ml;
  90.  
  91.     BarrettReduction64(ULL MOD_): MOD(MOD_) {
  92. 	ULL m = ~0ULL / MOD;
  93. 	if (m*MOD+MOD == 0ULL) m++;
  94. 	mh = m>>32;
  95. 	ml = m & ~0u;
  96.     }
  97.  
  98.     ULL reduce(ULL x) const {
  99. 	ULL z = (x & ~0u) * ml;
  100. 	z = (x & ~0u) * mh + (x>>32) * ml + (z>>32);
  101. 	z = (x>>32) * mh + (z>>32);
  102. 	x -= z * MOD;
  103. 	return x < MOD? x: x-MOD;
  104.     }
  105. };
  106.  
  107. const int SIZE = 1e8;
  108.  
  109. void test1(ULL MOD) {
  110.     ULL p = 1;
  111.     for (int i=1; i<SIZE; i++) {
  112. 	p = p * i % MOD;
  113.     }
  114.     printf("%llu\n", p);
  115. }
  116.  
  117. void test2(ULL MOD) {
  118.     MontgomeryReduction MR(MOD);
  119.     ULL p = MR.generate(1);
  120.     for (int i=1; i<SIZE; i++) {
  121. 	p = MR.reduce(p * MR.generate(i));
  122.     }
  123.     p = MR.reduce(p);
  124.     printf("%llu\n", p);
  125. }
  126.  
  127. void test3(ULL MOD) {
  128.     ULL p = 1;
  129.     BarrettReduction BR(MOD);
  130.     for (int i=1; i<SIZE; i++) {
  131. 	p = BR.reduce(p * i);
  132.     }
  133.     printf("%llu\n", p);
  134. }
  135.  
  136. void test4(ULL MOD) {
  137.     ULL p = 1;
  138.     BarrettReduction64 BR(MOD);
  139.     for (int i=1; i<SIZE; i++) {
  140. 	p = BR.reduce(p * i);
  141.     }
  142.     printf("%llu\n", p);
  143. }
  144.  
  145. int prime[] = {
  146.      998244353,
  147.     1000000007,
  148.     1000000009,
  149.     1000000021,
  150.     1000000033,
  151.     1000000087,
  152.     1000000093,
  153.     1000000097,
  154.     1000000103,
  155.     1000000123,
  156. };
  157.  
  158. void MAIN() {
  159.     REP (i, 10) {
  160. 	test3(prime[i]);
  161.     }
  162. }
  163.  
  164. int main() {
  165.     int TC = 1;
  166. //    scanf("%d", &TC);
  167.     REP (tc, TC) MAIN();
  168.     return 0;
  169. }
  170.  
Success #stdin #stdout 4.67s 4184KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
464016596
388742192
727351778
459413696
731607535
889548704
423123421
601786560
597154944
512094170
https://ideone.com/sNBQTz
language:
C++14 (gcc 8.3)
created:
4 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!