Advertisement
cyberjab

Untitled

Sep 16th, 2024
13
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 20.68 KB | None | 0 0
  1. class Num {
  2. public:
  3. typedef uint64_t word;
  4.  
  5. std::vector<word> words;
  6. bool neg;
  7.  
  8. static word word_mask() {
  9. return (word)-1;
  10. }
  11.  
  12. static size_t word_bits() {
  13. return sizeof(word) * CHAR_BIT;
  14. }
  15.  
  16. static word word_half_mask() {
  17. return word_mask() >> word_bits() / 2;
  18. }
  19.  
  20. static word char_to_word(char c) {
  21. switch (c) {
  22. case '0': return 0;
  23. case '1': return 1;
  24. case '2': return 2;
  25. case '3': return 3;
  26. case '4': return 4;
  27. case '5': return 5;
  28. case '6': return 6;
  29. case '7': return 7;
  30. case '8': return 8;
  31. case '9': return 9;
  32. case 'a': case 'A': return 10;
  33. case 'b': case 'B': return 11;
  34. case 'c': case 'C': return 12;
  35. case 'd': case 'D': return 13;
  36. case 'e': case 'E': return 14;
  37. case 'f': case 'F': return 15;
  38. case 'g': case 'G': return 16;
  39. case 'h': case 'H': return 17;
  40. case 'i': case 'I': return 18;
  41. case 'j': case 'J': return 19;
  42. case 'k': case 'K': return 20;
  43. case 'l': case 'L': return 21;
  44. case 'm': case 'M': return 22;
  45. case 'n': case 'N': return 23;
  46. case 'o': case 'O': return 24;
  47. case 'p': case 'P': return 25;
  48. case 'q': case 'Q': return 26;
  49. case 'r': case 'R': return 27;
  50. case 's': case 'S': return 28;
  51. case 't': case 'T': return 29;
  52. case 'u': case 'U': return 30;
  53. case 'v': case 'V': return 31;
  54. case 'w': case 'W': return 32;
  55. case 'x': case 'X': return 33;
  56. case 'y': case 'Y': return 34;
  57. case 'z': case 'Z': return 35;
  58. default: return word_mask();
  59. }
  60. }
  61.  
  62. static word word_gcd(word a, word b) {
  63. while (1) {
  64. if (a == 0) return b;
  65. b %= a;
  66. if (b == 0) return a;
  67. a %= b;
  68. }
  69. }
  70.  
  71. Num() : neg(false) {}
  72.  
  73. Num(size_t n, word w, bool neg = false) : words(n, w), neg(neg) {}
  74.  
  75. Num(const word* a, const word* b, bool neg = false) : words(a, b), neg(neg) {
  76. truncate();
  77. }
  78.  
  79. Num(const Num& a) {
  80. words = a.words;
  81. neg = a.neg;
  82. }
  83.  
  84. Num& operator = (const Num& a) {
  85. words = a.words;
  86. neg = a.neg;
  87. return *this;
  88. }
  89.  
  90. Num(int i) : neg(i < 0) {
  91. unsigned u = (i < 0) ? (unsigned)i : (unsigned)i;
  92. if (sizeof(u) <= word_bits()) {
  93. if (u > 0) push_back(u);
  94. }
  95. else {
  96. for (; u; u >>= word_bits()) push_back(u);
  97. }
  98. }
  99.  
  100. Num(const char* c, word base = 10, char** endptr = NULL) : neg(false) {
  101. // read sign
  102. if (*c == '-') {
  103. c++;
  104. neg = true;
  105. }
  106. // read digits
  107. for (; *c; c++) {
  108. mul_word(base);
  109. word b = char_to_word(*c);
  110. if (b >= base) break;
  111. add_word(b);
  112. }
  113. if (endptr) *endptr = (char*)c;
  114. }
  115.  
  116. void resize(size_t n) {
  117. words.resize(n);
  118. }
  119.  
  120. void pop_back() {
  121. words.pop_back();
  122. }
  123.  
  124. void push_back(word b) {
  125. words.push_back(b);
  126. }
  127.  
  128. word& back() {
  129. return words.back();
  130. }
  131.  
  132. const word& back() const {
  133. return words.back();
  134. }
  135.  
  136. size_t size() const {
  137. return words.size();
  138. }
  139.  
  140. word& operator [] (size_t i) {
  141. return words[i];
  142. }
  143.  
  144. const word& operator [] (size_t i) const {
  145. return words[i];
  146. }
  147.  
  148. Num& set_neg(bool neg) {
  149. this->neg = neg;
  150. return *this;
  151. }
  152.  
  153. Num& truncate() {
  154. while (size() > 0 && words.back() == 0) pop_back();
  155. return *this;
  156. }
  157.  
  158. size_t bitlength() const {
  159. if (size() == 0) return 0;
  160. size_t last = size() - 1;
  161. size_t result = word_bitlength((*this)[last]) + last * word_bits();
  162. return result;
  163. }
  164.  
  165. size_t count_trailing_zeros() const {
  166. for (size_t i = 0; i < size(); i++) {
  167. word w = (*this)[i];
  168. if (w) return word_count_trailing_zeros(w) + i * word_bits();
  169. }
  170. return 0;
  171. }
  172.  
  173. static int cmp_abs(const Num& a, const Num& b) {
  174. size_t na = a.size(), nb = b.size();
  175. if (na != nb) {
  176. return na < nb ? -1 : +1;
  177. }
  178. for (size_t i = na; i-- > 0;) {
  179. word wa = a[i], wb = b[i];
  180. if (wa != wb) {
  181. return wa < wb ? -1 : +1;
  182. }
  183. }
  184. return 0;
  185. }
  186.  
  187. static int cmp(const Num& a, const Num& b) {
  188. if (a.size() == 0 && b.size() == 0) return 0;
  189. if (!a.neg && !b.neg) return +cmp_abs(a, b);
  190. if (a.neg && b.neg) return -cmp_abs(a, b);
  191. return a.neg && !b.neg ? -1 : +1;
  192. }
  193.  
  194. static size_t word_bitlength(word a) {
  195. for (int i = word_bits() - 1; i >= 0; i--) if ((a >> i) & 1) return i + 1;
  196. return 0;
  197. }
  198.  
  199. static size_t word_count_trailing_zeros(word a) {
  200. for (int i = 0; i < (int)word_bits(); i++) if ((a >> i) & 1) return i;
  201. return word_bits();
  202. }
  203.  
  204. static word add_carry(word* a, word b) {
  205. return (*a += b) < b;
  206. }
  207.  
  208. static word sub_carry(word* a, word b) {
  209. word tmp = *a;
  210. return (*a = tmp - b) > tmp;
  211. }
  212.  
  213. static word word_mul_hi(word a, word b) {
  214. size_t n = word_bits() / 2;
  215. word a_hi = a >> n;
  216. word a_lo = a & word_half_mask();
  217. word b_hi = b >> n;
  218. word b_lo = b & word_half_mask();
  219. word tmp = ((a_lo * b_lo) >> n) + a_hi * b_lo;
  220. tmp = (tmp >> n) + ((a_lo * b_hi + (tmp & word_half_mask())) >> n);
  221. return tmp + a_hi * b_hi;
  222. }
  223.  
  224. static Num& add_unsigned_overwrite(Num& a, const Num& b) {
  225. size_t i, na = a.size(), nb = b.size(), n = std::max(na, nb);
  226. a.resize(n);
  227. word carry = 0;
  228. for (i = 0; i < nb; i++) {
  229. carry = add_carry(&a[i], carry);
  230. carry += add_carry(&a[i], b[i]);
  231. }
  232. for (; i < n && carry; i++) carry = add_carry(&a[i], carry);
  233. if (carry) a.push_back(carry);
  234. return a.truncate();
  235. }
  236.  
  237. static Num& sub_unsigned_overwrite(Num& a, const Num& b) {
  238. //assert(cmp_abs(a, b) >= 0);
  239. size_t i, na = a.size(), nb = b.size();
  240. word carry = 0;
  241. for (i = 0; i < nb; i++) {
  242. carry = sub_carry(&a[i], carry);
  243. carry += sub_carry(&a[i], b[i]);
  244. }
  245. for (; i < na && carry; i++) carry = sub_carry(&a[i], carry);
  246. //assert(!carry);
  247. return a.truncate();
  248. }
  249.  
  250. static Num mul_long(const Num& a, const Num& b) {
  251. size_t na = a.size(), nb = b.size(), nc = na + nb + 1;
  252. Num c(nc, 0, a.neg ^ b.neg), carries(nc, 0);
  253. for (size_t ia = 0; ia < na; ia++) {
  254. for (size_t ib = 0; ib < nb; ib++) {
  255. size_t i = ia + ib, j = i + 1;
  256. // WARNING: Might overflow if word size is chosen too small
  257. carries[i + 1] += add_carry(&c[i], a[ia] * b[ib]);
  258. carries[j + 1] += add_carry(&c[j], word_mul_hi(a[ia], b[ib]));
  259. }
  260. }
  261. return add_unsigned_overwrite(c, carries).truncate();
  262. }
  263.  
  264. static Num mul_karatsuba(const Num& a, const Num& b) {
  265. size_t na = a.size(), nb = b.size(), n = std::max(na, nb), m2 = n / 2 + (n & 1);
  266. Num a_parts[2], b_parts[2];
  267. split(a, a_parts, 2, m2);
  268. split(b, b_parts, 2, m2);
  269. m2 *= word_bits();
  270. Num z0 = a_parts[0] * b_parts[0];
  271. Num z1 = (a_parts[0] + a_parts[1]) * (b_parts[0] + b_parts[1]);
  272. Num z2 = a_parts[1] * b_parts[1];
  273. Num result = z2;
  274. result <<= m2;
  275. result += z1 - z2 - z0;
  276. result <<= m2;
  277. result += z0;
  278. return result;
  279. }
  280.  
  281. static Num mul(const Num& a, const Num& b) {
  282. size_t karatsuba_threshold = 20;
  283. if (a.size() > karatsuba_threshold && b.size() > karatsuba_threshold) {
  284. return mul_karatsuba(a, b);
  285. }
  286. return mul_long(a, b);
  287. }
  288.  
  289. static Num add_signed(const Num& a, bool a_neg, const Num& b, bool b_neg) {
  290. if (a_neg == b_neg) return add_unsigned(a, b).set_neg(a_neg);
  291. if (cmp_abs(a, b) >= 0) return sub_unsigned(a, b).set_neg(a_neg);
  292. return sub_unsigned(b, a).set_neg(b_neg);
  293. }
  294.  
  295. Num& operator >>= (size_t n_bits) {
  296. if (n_bits == 0) return *this;
  297. size_t n_words = n_bits / word_bits();
  298. if (n_words >= size()) {
  299. resize(0);
  300. return *this;
  301. }
  302. n_bits %= word_bits();
  303. if (n_bits == 0) {
  304. for (size_t i = 0; i < size() - n_words; i++) {
  305. (*this)[i] = (*this)[i + n_words];
  306. }
  307. }
  308. else {
  309. word hi, lo = (*this)[n_words];
  310. for (size_t i = 0; i < size() - n_words - 1; i++) {
  311. hi = (*this)[i + n_words + 1];
  312. (*this)[i] = (hi << (word_bits() - n_bits)) | (lo >> n_bits);
  313. lo = hi;
  314. }
  315. (*this)[size() - n_words - 1] = lo >> n_bits;
  316. }
  317. resize(size() - n_words);
  318. return truncate();
  319. }
  320.  
  321. Num& operator <<= (size_t n_bits) {
  322. if (n_bits == 0) return *this;
  323. size_t n_words = n_bits / word_bits();
  324. n_bits %= word_bits();
  325. size_t old_size = size();
  326. size_t n = old_size + n_words + (n_bits != 0);
  327. resize(n);
  328. if (n_bits == 0) {
  329. for (size_t i = n; i-- > n_words;) {
  330. (*this)[i] = (*this)[i - n_words];
  331. }
  332. }
  333. else {
  334. word lo, hi = 0;
  335. for (size_t i = n - 1; i > n_words; i--) {
  336. lo = (*this)[i - n_words - 1];
  337. (*this)[i] = (hi << n_bits) | (lo >> (word_bits() - n_bits));
  338. hi = lo;
  339. }
  340. (*this)[n_words] = hi << n_bits;
  341. }
  342. for (size_t i = 0; i < n_words; i++) (*this)[i] = 0;
  343. return truncate();
  344. }
  345.  
  346. static void div_mod(const Num& numerator, Num denominator, Num& quotient, Num& remainder) {
  347. quotient = 0;
  348. remainder = numerator;
  349. if (cmp_abs(remainder, denominator) >= 0) {
  350. int n = numerator.bitlength() - denominator.bitlength();
  351. denominator <<= n;
  352. for (; n >= 0; n--) {
  353. if (cmp_abs(remainder, denominator) >= 0) {
  354. sub_unsigned_overwrite(remainder, denominator);
  355. quotient.set_bit(n);
  356. }
  357. denominator >>= 1;
  358. }
  359. }
  360. quotient.set_neg(numerator.neg ^ denominator.neg);
  361. remainder.set_neg(numerator.neg);
  362. }
  363.  
  364. static void div_mod_half_word(const Num& numerator, word denominator, Num& quotient, word& remainder) {
  365. remainder = 0;
  366. Num dst(numerator.size(), 0);
  367.  
  368. for (size_t i = numerator.size(); i-- > 0;) {
  369. word dst_word = 0;
  370. word src_word = numerator[i];
  371. word parts[2];
  372. parts[0] = src_word >> word_bits() / 2;
  373. parts[1] = src_word & word_half_mask();
  374.  
  375. for (size_t j = 0; j < 2; j++) {
  376. remainder <<= word_bits() / 2;
  377. remainder |= parts[j];
  378.  
  379. word div_word = remainder / denominator;
  380. word mod_word = remainder % denominator;
  381. remainder = mod_word;
  382.  
  383. dst_word <<= word_bits() / 2;
  384. dst_word |= div_word;
  385. }
  386.  
  387. dst[i] = dst_word;
  388. }
  389.  
  390. quotient = dst.truncate().set_neg(numerator.neg);
  391. }
  392.  
  393. static void split(const Num& a, Num* parts, size_t n_parts, size_t n) {
  394. size_t i = 0;
  395. for (size_t k = 0; k < n_parts; k++) {
  396. Num& part = parts[k];
  397. part.resize(n);
  398. for (size_t j = 0; j < n && i < a.size(); j++) part[j] = a[i++];
  399. part = part.truncate();
  400. }
  401. }
  402.  
  403. static Num div(const Num& numerator, const Num& denominator) {
  404. Num quotient, remainder;
  405. div_mod(numerator, denominator, quotient, remainder);
  406. return quotient;
  407. }
  408.  
  409. static Num mod(const Num& numerator, const Num& denominator) {
  410. Num quotient, remainder;
  411. div_mod(numerator, denominator, quotient, remainder);
  412. return remainder;
  413. }
  414.  
  415. static Num add_unsigned(const Num& a, const Num& b) {
  416. Num result(a);
  417. return add_unsigned_overwrite(result, b);
  418. }
  419.  
  420. static Num sub_unsigned(const Num& a, const Num& b) {
  421. Num result(a);
  422. return sub_unsigned_overwrite(result, b);
  423. }
  424.  
  425. static Num add(const Num& a, const Num& b) {
  426. Num result = add_signed(a, a.neg, b, b.neg);
  427. return result;
  428. }
  429.  
  430. static Num sub(const Num& a, const Num& b) {
  431. Num result = add_signed(a, a.neg, b, !b.neg);
  432. return result;
  433. }
  434.  
  435. static Num gcd(const Num& a0, const Num& b0) {
  436.  
  437. if (a0.size() == 1 && b0.size() == 1) {
  438. return Num(1, word_gcd(a0[0], b0[0]));
  439. }
  440.  
  441. Num a(a0), b(b0);
  442. a.neg = b.neg = false;
  443.  
  444. if (a.size() == 0) return b0;
  445. if (b.size() == 0) return a0;
  446.  
  447. size_t n = a.count_trailing_zeros();
  448. size_t m = b.count_trailing_zeros();
  449. if (n > m) {
  450. std::swap(n, m);
  451. a.words.swap(b.words);
  452. }
  453.  
  454. a >>= n;
  455. b >>= n;
  456.  
  457. do {
  458. b >>= b.count_trailing_zeros();
  459. if (cmp_abs(a, b) > 0) a.words.swap(b.words);
  460. sub_unsigned_overwrite(b, a);
  461. } while (b.size() > 0);
  462.  
  463. a <<= n;
  464.  
  465. return a;
  466. }
  467.  
  468. typedef void (*random_func)(uint8_t* bytes, size_t n_bytes);
  469.  
  470. static Num random_bits(size_t n_bits, random_func func) {
  471. if (n_bits == 0) return 0;
  472. size_t partial_bits = n_bits % word_bits();
  473. size_t n_words = n_bits / word_bits() + (partial_bits > 0);
  474. size_t n_bytes = n_words * sizeof(word);
  475. Num result(n_words, 0);
  476. uint8_t* bytes = (uint8_t*)&result[0];
  477. func(bytes, n_bytes);
  478. if (partial_bits) {
  479. size_t too_many_bits = word_bits() - partial_bits;
  480. result.back() >>= too_many_bits;
  481. }
  482. return result;
  483. }
  484.  
  485. static Num random_inclusive(const Num& inclusive, random_func func) {
  486. size_t n_bits = inclusive.bitlength();
  487. while (true) {
  488. Num result = random_bits(n_bits, func);
  489. if (result <= inclusive) return result;
  490. }
  491. }
  492.  
  493. static Num random_exclusive(const Num& exclusive, random_func func) {
  494. size_t n_bits = exclusive.bitlength();
  495. while (true) {
  496. Num result = random_bits(n_bits, func);
  497. if (result < exclusive) return result;
  498. }
  499. }
  500.  
  501. static Num random_second_exclusive(const Num& inclusive_min_val, const Num& exclusive_max_val, random_func func) {
  502. return inclusive_min_val + random_exclusive(exclusive_max_val - inclusive_min_val, func);
  503. }
  504.  
  505. static Num random_both_inclusive(const Num& inclusive_min_val, const Num& inclusive_max_val, random_func func) {
  506. return inclusive_min_val + random_inclusive(inclusive_max_val - inclusive_min_val, func);
  507. }
  508.  
  509. Num& set_bit(size_t i) {
  510. size_t i_word = i / word_bits();
  511. size_t i_bit = i % word_bits();
  512. if (size() <= i_word) resize(i_word + 1);
  513. (*this)[i_word] |= ((word)1) << i_bit;
  514. return *this;
  515. }
  516.  
  517. word get_bit(size_t i) const {
  518. size_t i_word = i / word_bits();
  519. size_t i_bit = i % word_bits();
  520. if (i_word >= size()) return 0;
  521. return ((*this)[i_word] >> i_bit) & 1;
  522. }
  523.  
  524. void clr_bit(size_t i) {
  525. size_t i_word = i / word_bits();
  526. size_t i_bit = i % word_bits();
  527. if (i_word >= size()) return;
  528. word mask = 1;
  529. mask <<= i_bit;
  530. (*this)[i_word] &= ~mask;
  531. }
  532.  
  533. Num& mul_word(word b) {
  534. word carry = 0;
  535. for (size_t i = 0; i < size(); i++) {
  536. word a = (*this)[i];
  537. word tmp = a * b;
  538. carry = add_carry(&tmp, carry);
  539. carry += word_mul_hi(a, b);
  540. (*this)[i] = tmp;
  541. }
  542. if (carry) push_back(carry);
  543. return truncate();
  544. }
  545.  
  546. Num& add_word(word carry, size_t i0 = 0) {
  547. if (i0 >= size()) resize(i0 + 1);
  548. for (size_t i = i0; i < size() && carry; i++) {
  549. carry = add_carry(&(*this)[i], carry);
  550. }
  551. if (carry) push_back(carry);
  552. return truncate();
  553. }
  554.  
  555. void print(
  556. std::vector<char>& text,
  557. word base = 10,
  558. const char* alphabet = "0123456789abcdefghijklmnopqrstuvwxyz"
  559. ) const {
  560. if (size() == 0) {
  561. text.push_back('0');
  562. }
  563. else {
  564. Num tmp(*this);
  565. while (tmp.size() > 0) {
  566. word remainder;
  567. div_mod_half_word(tmp, base, tmp, remainder);
  568. text.push_back(alphabet[remainder]);
  569. }
  570. if (neg) text.push_back('-');
  571. std::reverse(text.begin(), text.end());
  572. }
  573. text.push_back('\0');
  574. }
  575.  
  576. friend std::ostream& operator << (std::ostream& out, const Num& num) {
  577. std::vector<char> tmp;
  578. num.print(tmp);
  579. out << &tmp[0];
  580. return out;
  581. }
  582.  
  583. double to_double() const {
  584. if (size() == 0) return 0.0;
  585. double d = 0.0, base = ::pow(2.0, word_bits());
  586. for (size_t i = size(); i-- > 0;) d = d * base + (*this)[i];
  587. return neg ? -d : d;
  588. }
  589.  
  590. bool can_convert_to_int(int* result) {
  591. if (*this < Num(INT_MIN) || *this > Num(INT_MAX)) return false;
  592.  
  593. unsigned temp = 0;
  594.  
  595. if (word_bits() >= sizeof(temp) * CHAR_BIT) {
  596. if (words.size() > 0) {
  597. temp = (*this)[0];
  598. }
  599. }
  600. else {
  601. for (size_t i = words.size(); i-- > 0;) {
  602. temp <<= word_bits();
  603. temp += (*this)[i];
  604. }
  605. }
  606.  
  607. *result = neg ? temp : temp;
  608.  
  609. return true;
  610. }
  611.  
  612. Num pow(size_t exponent) const {
  613. Num result(1), p(*this);
  614. for (; exponent; exponent >>= 1) {
  615. if (exponent & 1) {
  616. result = result * p;
  617. exponent--;
  618. }
  619. p = p * p;
  620. }
  621. return result;
  622. }
  623.  
  624. Num mod_pow(Num exponent, const Num& modulus) const {
  625. Num result(1), base = (*this) % modulus;
  626. for (; exponent.size() > 0; exponent >>= 1) {
  627. if (exponent.get_bit(0)) {
  628. result = (result * base) % modulus;
  629. }
  630. base = (base * base) % modulus;
  631. }
  632. return result;
  633. }
  634.  
  635. Num sqrt() const {
  636. Num n = *this;
  637. int bit = bitlength();
  638. if (bit & 1) bit ^= 1;
  639. Num result = 0;
  640. for (; bit >= 0; bit -= 2) {
  641. Num tmp = result;
  642. tmp.set_bit(bit);
  643. if (n >= tmp) {
  644. n -= tmp;
  645. result.set_bit(bit + 1);
  646. }
  647. result >>= 1;
  648. }
  649. return result;
  650. }
  651.  
  652. Num& operator ++() {
  653. add_word(1);
  654. return *this;
  655. }
  656.  
  657. Num& operator += (const Num& b) { return *this = add(*this, b); }
  658. Num& operator -= (const Num& b) { return *this = sub(*this, b); }
  659. Num& operator *= (const Num& b) { return *this = mul(*this, b); }
  660. Num& operator /= (const Num& b) { return *this = div(*this, b); }
  661. Num& operator %= (const Num& b) { return *this = mod(*this, b); }
  662.  
  663. bool operator == (const Num& b) const { return cmp(*this, b) == 0; }
  664. bool operator != (const Num& b) const { return cmp(*this, b) != 0; }
  665. bool operator <= (const Num& b) const { return cmp(*this, b) <= 0; }
  666. bool operator >= (const Num& b) const { return cmp(*this, b) >= 0; }
  667. bool operator < (const Num& b) const { return cmp(*this, b) < 0; }
  668. bool operator > (const Num& b) const { return cmp(*this, b) > 0; }
  669.  
  670. Num operator + (const Num& b) const { return add(*this, b); }
  671. Num operator - (const Num& b) const { return sub(*this, b); }
  672. Num operator * (const Num& b) const { return mul(*this, b); }
  673. Num operator / (const Num& b) const { return div(*this, b); }
  674. Num operator % (const Num& b) const { return mod(*this, b); }
  675. Num operator - () const { return Num(*this).set_neg(!neg); }
  676.  
  677. Num operator >> (size_t n_bits) const { return Num(*this) >>= n_bits; }
  678. Num operator << (size_t n_bits) const { return Num(*this) <<= n_bits; }
  679. };
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement