bignum.h 7.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326
  1. /***
  2. * This software was written by Nils Maier. No copyright is claimed, and the
  3. * software is hereby placed in the public domain.
  4. *
  5. * In case this attempt to disclaim copyright and place the software in the
  6. * public domain is deemed null and void in your jurisdiction, then the
  7. * software is Copyright 2004,2013 Nils Maier and it is hereby released to the
  8. * general public under the following terms:
  9. * Redistribution and use in source and binary forms, with or without
  10. * modification, are permitted.
  11. * There's ABSOLUTELY NO WARRANTY, express or implied.
  12. */
  13. #ifndef BIGNUM_H
  14. #define BIGNUM_H
  15. #include <cstring>
  16. #include <algorithm>
  17. #include <memory>
  18. #include <stdint.h>
  19. namespace bignum {
  20. template <size_t dim> class ulong {
  21. public:
  22. typedef char char_t;
  23. typedef std::make_unsigned<char_t>::type uchar_t;
  24. private:
  25. std::unique_ptr<char_t[]> buf_;
  26. public:
  27. inline ulong() : buf_(new char_t[dim]()) {}
  28. inline ulong(size_t t) : buf_(new char_t[dim]())
  29. {
  30. memcpy(buf_.get(), (char_t*)&t, sizeof(t));
  31. }
  32. inline ulong(const ulong<dim>& rhs) : buf_(new char_t[dim]())
  33. {
  34. memcpy(buf_.get(), rhs.buf_.get(), dim);
  35. }
  36. explicit inline ulong(const char_t* data, size_t size)
  37. : buf_(new char_t[dim]())
  38. {
  39. if (size > dim) {
  40. throw std::bad_alloc();
  41. }
  42. memcpy(buf_.get(), data, size);
  43. }
  44. virtual ~ulong() {}
  45. ulong<dim>& operator=(const ulong<dim>& rhs)
  46. {
  47. memcpy(buf_.get(), rhs.buf_.get(), dim);
  48. return *this;
  49. }
  50. bool operator==(const ulong<dim>& rhs) const
  51. {
  52. return memcmp(buf_.get(), rhs.buf_.get(), dim) == 0;
  53. }
  54. bool operator!=(const ulong<dim>& rhs) const
  55. {
  56. return memcmp(buf_.get(), rhs.buf_.get(), dim) != 0;
  57. }
  58. bool operator>(const ulong<dim>& rhs) const
  59. {
  60. const auto b1 = buf_.get();
  61. const auto b2 = rhs.buf_.get();
  62. for (ssize_t i = dim - 1; i >= 0; --i) {
  63. for (ssize_t j = 1; j >= 0; --j) {
  64. char_t t = ((uchar_t)(b1[i] << 4 * (1 - j))) >> 4;
  65. char_t r = ((uchar_t)(b2[i] << 4 * (1 - j))) >> 4;
  66. if (t != r) {
  67. return t > r;
  68. }
  69. }
  70. }
  71. return false;
  72. }
  73. bool operator>=(const ulong<dim>& rhs) const
  74. {
  75. return *this == rhs || *this > rhs;
  76. }
  77. bool operator<(const ulong<dim>& rhs) const { return !(*this >= rhs); }
  78. bool operator<=(const ulong<dim>& rhs) const
  79. {
  80. return *this == rhs || *this < rhs;
  81. }
  82. ulong<dim> operator+(const ulong<dim>& rhs) const
  83. {
  84. ulong<dim> rv;
  85. const auto b1 = buf_.get();
  86. const auto b2 = rhs.buf_.get();
  87. const auto rb = rv.buf_.get();
  88. bool base = false;
  89. for (size_t i = 0; i < dim; ++i) {
  90. for (ssize_t j = 0; j < 2; ++j) {
  91. char_t t = ((uchar_t)(b1[i] << 4 * (1 - j))) >> 4;
  92. char_t r = ((uchar_t)(b2[i] << 4 * (1 - j))) >> 4;
  93. if (base) {
  94. t++;
  95. }
  96. if (r + t >= 16) {
  97. rb[i] += (t + r - 16) << j * 4;
  98. base = true;
  99. }
  100. else {
  101. rb[i] += (t + r) << j * 4;
  102. base = false;
  103. }
  104. }
  105. }
  106. return rv;
  107. }
  108. ulong<dim>& operator+=(const ulong<dim>& rhs)
  109. {
  110. *this = *this + rhs;
  111. return *this;
  112. }
  113. ulong<dim>& operator++()
  114. {
  115. *this = *this + 1;
  116. return *this;
  117. }
  118. ulong<dim> operator++(int)
  119. {
  120. ulong<dim> tmp = *this;
  121. *this = *this + 1;
  122. return tmp;
  123. }
  124. ulong<dim> operator-(const ulong<dim>& rhs) const
  125. {
  126. ulong<dim> rv;
  127. const auto b1 = buf_.get();
  128. const auto b2 = rhs.buf_.get();
  129. const auto rb = rv.buf_.get();
  130. bool base = false;
  131. for (size_t i = 0; i < dim; ++i) {
  132. for (ssize_t j = 0; j < 2; ++j) {
  133. char_t t = ((uchar_t)(b1[i] << 4 * (1 - j))) >> 4;
  134. char_t r = ((uchar_t)(b2[i] << 4 * (1 - j))) >> 4;
  135. if (base) {
  136. t--;
  137. }
  138. if (t >= r) {
  139. rb[i] += (t - r) << j * 4;
  140. base = false;
  141. }
  142. else {
  143. rb[i] += (t + 16 - r) << j * 4;
  144. base = true;
  145. }
  146. }
  147. }
  148. return rv;
  149. }
  150. ulong<dim>& operator-=(const ulong<dim>& rhs)
  151. {
  152. *this = *this - rhs;
  153. return *this;
  154. }
  155. ulong<dim>& operator--()
  156. {
  157. *this = *this - 1;
  158. return *this;
  159. }
  160. ulong<dim> operator--(int)
  161. {
  162. ulong<dim> tmp = *this;
  163. *this = *this - 1;
  164. return tmp;
  165. }
  166. ulong<dim> operator*(const ulong<dim>& rhs) const
  167. {
  168. ulong<dim> c = rhs, rv;
  169. const ulong<dim> null;
  170. size_t cap = c.capacity();
  171. while (c != null) {
  172. ulong<dim> tmp = *this;
  173. tmp.mul(cap - 1);
  174. rv += tmp;
  175. ulong<dim> diff(1);
  176. diff.mul(cap - 1);
  177. c -= diff;
  178. cap = c.capacity();
  179. }
  180. return rv;
  181. }
  182. ulong<dim>& operator*=(const ulong<dim>& rhs)
  183. {
  184. *this = *this* rhs;
  185. return *this;
  186. }
  187. ulong<dim> operator/(const ulong<dim>& rhs) const
  188. {
  189. ulong<dim> quotient, remainder;
  190. div(rhs, quotient, remainder);
  191. return quotient;
  192. }
  193. ulong<dim>& operator/=(const ulong<dim>& rhs)
  194. {
  195. *this = *this / rhs;
  196. return *this;
  197. }
  198. ulong<dim> operator%(const ulong<dim>& rhs) const
  199. {
  200. ulong<dim> quotient, remainder;
  201. div(rhs, quotient, remainder);
  202. return remainder;
  203. }
  204. ulong<dim>& operator%=(const ulong<dim>& rhs)
  205. {
  206. *this = *this % rhs;
  207. return *this;
  208. }
  209. ulong<dim> mul_mod(const ulong<dim>& mul, const ulong<dim>& mod) const
  210. {
  211. if (capacity() + mul.capacity() <= dim) {
  212. return (*this * mul) % mod;
  213. }
  214. ulong<dim * 2> et(buf_.get(), dim), emul(mul.buf_.get(), dim),
  215. emod(mod.buf_.get(), dim), erv = (et * emul) % emod;
  216. ulong<dim> rv;
  217. erv.binary(rv.buf_.get(), dim);
  218. return rv;
  219. }
  220. std::unique_ptr<char_t[]> binary() const
  221. {
  222. ulong<dim> c = *this;
  223. std::unique_ptr<char_t[]> rv;
  224. rv.swap(c.buf_);
  225. return rv;
  226. }
  227. void binary(char_t* buf, size_t len) const
  228. {
  229. memcpy(buf, buf_.get(), std::min(dim, len));
  230. }
  231. size_t length() const { return dim; }
  232. private:
  233. size_t capacity() const
  234. {
  235. size_t rv = dim * 2;
  236. const auto b = buf_.get();
  237. for (ssize_t i = dim - 1; i >= 0; --i) {
  238. char_t f = b[i] >> 4;
  239. char_t s = (b[i] << 4) >> 4;
  240. if (!f && !s) {
  241. rv -= 2;
  242. continue;
  243. }
  244. if (!f) {
  245. --rv;
  246. }
  247. return rv;
  248. }
  249. return rv;
  250. }
  251. void mul(size_t digits)
  252. {
  253. ulong<dim> tmp = *this;
  254. auto bt = tmp.buf_.get();
  255. auto b = buf_.get();
  256. memset(b, 0, dim);
  257. const size_t npar = digits % 2;
  258. const size_t d2 = digits / 2;
  259. for (size_t i = d2; i < dim; ++i) {
  260. for (size_t j = 0; j < 2; ++j) {
  261. char_t c = ((uchar_t)(bt[(dim - 1) - i] << 4 * (1 - j))) >> 4;
  262. char_t r = c << (npar * (1 - j) * 4 + (1 - npar) * j * 4);
  263. ssize_t idx = i - d2 - npar * j;
  264. if (idx >= 0) {
  265. b[(dim - 1) - idx] += r;
  266. }
  267. }
  268. }
  269. }
  270. void div(const ulong<dim>& d, ulong<dim>& q, ulong<dim>& r) const
  271. {
  272. ulong<dim> tmp = d;
  273. r = *this;
  274. q = 0;
  275. size_t cr = r.capacity();
  276. const size_t cd = d.capacity();
  277. while (cr > cd) {
  278. tmp = d;
  279. tmp.mul(cr - cd - 1);
  280. ulong<dim> qt(1);
  281. qt.mul(cr - cd - 1);
  282. ulong<dim> t = tmp;
  283. t.mul(1);
  284. if (r >= t) {
  285. tmp = t;
  286. qt.mul(1);
  287. }
  288. while (r >= tmp) {
  289. r -= tmp;
  290. q += qt;
  291. }
  292. cr = r.capacity();
  293. }
  294. while (r >= d) {
  295. r -= d;
  296. ++q;
  297. }
  298. }
  299. };
  300. } // namespace bignum
  301. #endif