fork download
copy 
  1. #include <iostream>
  2. #include <string>
  3. #include <cctype>
  4. #include <cassert>
  5. #include <exception>
  6.  
  7. enum {
  8.   BIT = 4,   // for test
  9. //  BIT = 30,
  10.   CMASK = 1 << BIT,
  11.   BMASK = CMASK - 1
  12. };
  13.  
  14. namespace QZ {
  15. class mpz_class {
  16. private:
  17.   int n;
  18.   int *b;
  19. private:
  20.   friend mpz_class sub(mpz_class const &a, mpz_class const &b);
  21.   friend mpz_class mul2(mpz_class const &a, mpz_class const &b);
  22.   friend void div2(mpz_class const &a, mpz_class const &b, mpz_class &q, mpz_class &r);
  23.   static bool iszero(mpz_class &n);
  24.   static void LeftShift(mpz_class &a);
  25.   static void RightShift(mpz_class &a, int &lsb);
  26.   static void LeftShift2(int &msb, mpz_class &a);
  27.   static void LeftShift3(mpz_class &a, int msb);
  28.  
  29. public:
  30.   mpz_class();
  31.   mpz_class(const mpz_class &ob);
  32.   mpz_class(const unsigned long n);
  33.   ~mpz_class();
  34.   mpz_class &operator=(mpz_class const n);
  35.   friend bool operator==(mpz_class const &a, mpz_class const &b);
  36.   friend bool operator!=(mpz_class const &a, mpz_class const &b) { return !(a == b); }
  37.   friend mpz_class operator+(mpz_class const &a, mpz_class const &b);
  38.   friend mpz_class operator-(mpz_class const &a, mpz_class const &b) { return sub(a, b); }
  39.  
  40.   friend mpz_class operator*(mpz_class const &a, mpz_class const &b) { return mul2(a, b); }
  41.   friend mpz_class operator/(mpz_class const &a, mpz_class const &b) { mpz_class q, r; div2(a, b, q, r); return q; }
  42.   friend mpz_class operator%(mpz_class const &a, mpz_class const &b) { mpz_class q, r; div2(a, b, q, r); return r; }
  43.  
  44.   friend bool operator<(mpz_class const &a, mpz_class const &b);
  45.   friend bool operator<=(mpz_class const &a, mpz_class const &b) { return (a < b) ? true : (a == b) ? true : false; }
  46.   friend bool operator>(mpz_class const &a, mpz_class const &b) { return (a < b) ? false : (a == b) ? false : true; }
  47.   friend bool operator>=(mpz_class const &a, mpz_class const &b) { return (a < b) ? false : (a == b) ? true : true; }
  48.  
  49.  
  50.   mpz_class &operator+=(mpz_class const &n) { *this = *this + n; return *this; }
  51.   mpz_class &operator-=(mpz_class const &n) { *this = *this - n; return *this; }
  52.   mpz_class &operator*=(mpz_class const &n) { *this = *this * n; return *this; }
  53.   mpz_class &operator/=(mpz_class const &n) { *this = *this / n; return *this; }
  54.   mpz_class &operator%=(mpz_class const &n) { *this = *this % n; return *this; }
  55.  
  56.  
  57.   friend std::ostream &operator<<(std::ostream &stream, mpz_class c);
  58.   void dump() const {
  59.     std::cout << "n = " << this->n << ": ";
  60.     for (int i = 0; i < this->n; i++)
  61.       std::cout << i << ":" << this->b[i] << ",";
  62.     std::cout << std::endl;
  63.   }
  64. /*--------------------------------------------------------*/
  65. public:
  66.   void print_tmp1(int prec);
  67.   void bit_output(int i, int K); /* i: 現在着目している桁, k:一つの桁でビットを表示する個数 */
  68. }; /* class QZ::mpz_class */
  69.  
  70. void mpz_class::bit_output(int i, int K) { /* i: 現在着目している桁, k:一つの桁でビットを表示する個数 */
  71.   int pbuf = 0;
  72.   int mask = 1;
  73.   for (int j = 0; j < K; j++) {
  74.     if (((this->b[i]) & mask) > 0) {
  75.       pbuf = (pbuf | 0x01);
  76.     } else {
  77.     }
  78.     pbuf = (pbuf << 1);    
  79.     mask = (mask << 1);
  80.   }
  81.   pbuf = (pbuf >> 1);    
  82.   for (int j = 0; j < K; j++) {
  83.     if ((pbuf & 0x01) > 0)
  84.       std::cout << '1';
  85.     else
  86.       std::cout << '0';
  87.     pbuf = (pbuf >> 1);
  88.   }
  89. }
  90.  
  91. int f(int n, int BIT) { /* f(1)=1, f(2)=2, f(3)=3, f(4)=4, f(5)=1, f(6)=2, f(7)=3, f(8)=4,f(9)=1,f(10)=2 etc */
  92.   int s;
  93.   return ((s = (n % BIT)) == 0) ? BIT : s;
  94. }
  95.  
  96. void mpz_class::print_tmp1(int prec) {
  97.   int max_need_buff_position = (prec - 1) / 4;
  98.   for (int i = max_need_buff_position; i >= 0 /* i is signed  */; --i) {
  99.  
  100.     if (i == max_need_buff_position && i > this->n - 1) { /* 要求される精度におけるMSB 桁で、メモリ上に存在しない */
  101.       for (int j = 0; j < f(prec, BIT); j++)
  102.         std::cout << "0";
  103.  
  104.     } else if (i == max_need_buff_position && i <= this->n - 1) { /* 要求される精度におけるMSBを含む桁で、メモリ上に存在する */
  105.       this->QZ::mpz_class::bit_output(i, f(prec, BIT));
  106.  
  107.     } else if (i < max_need_buff_position && i > this->n - 1) { /* 要求される精度におけるMSB 桁未満で、メモリ上に存在しない */
  108.       for (int j = 0; j < BIT; j++)
  109.         std::cout << "0";
  110.  
  111.     } else if (i < max_need_buff_position && i <= this->n - 1) { /* 要求される精度におけるMSB 桁未満で、メモリ上に存在する */
  112.       this->QZ::mpz_class::bit_output(i, BIT);
  113.  
  114.     } else { /* for debug */
  115.       assert(0); /* not reach */
  116.     }
  117.   }
  118.   std::cout << std::endl;
  119. }
  120.  
  121. /*--------------------------------------------------------*/
  122. mpz_class::mpz_class() {
  123.   this->n = 1;
  124.   this->b = new int[1];
  125.   this->b[0] = 0;
  126. }
  127.  
  128. mpz_class::mpz_class(const mpz_class &ob) {
  129.   this->n = ob.n;
  130.   if (ob.n == 0) {
  131.     this->b = 0;
  132.   } else {
  133.     this->b = new int[ob.n];
  134.     for (int i = 0; i < ob.n; i++)
  135.       this->b[i] = ob.b[i];
  136.   }
  137. }
  138.  
  139. mpz_class::mpz_class(const unsigned long c) {
  140.   if (c <= 0) {
  141.     this->n = 1;
  142.     this->b = new int[1];
  143.     this->b[0] = 0;
  144.   } else {
  145.     int *newbody;
  146.     int i, j, r, cc, pos;
  147.     this->b = 0;
  148.     this->n = 0;
  149.     cc = c;
  150.     while (cc > 0) {
  151.       r = 0;
  152.       pos = 1;
  153.       for (i = 0; i < BIT; i++) {
  154.         if ((cc & 1) != 0) {
  155.           r = (r | pos);
  156.         } else {
  157.         }
  158.         pos = pos * 2;
  159.         cc = cc / 2;
  160.       }
  161.       newbody = new int[this->n + 1];
  162.       for (j = 0; j < this->n; j++)
  163.         newbody[j] = this->b[j];
  164.       newbody[j] = r; /* j == this->n */
  165.       delete this->b;
  166.       this->b = newbody;
  167.       this->n++;
  168.     }
  169.   }
  170.   assert(this->n != 0);
  171. }
  172.  
  173. mpz_class::~mpz_class() {
  174.   delete [] this->b;
  175.   this->b = 0;
  176.   this->n = 0;
  177. }
  178.  
  179. mpz_class &mpz_class::operator=(mpz_class ob) {
  180.   this->n = ob.n;
  181.   if (ob.n == 0) {
  182.     delete [] this->b;
  183.     this->b = 0;
  184.   } else {
  185.     delete [] this->b;
  186.     this->b = new int[ob.n];
  187.     for (int i = 0; i < ob.n; i++) {
  188.       this->b[i] = ob.b[i];
  189.     }
  190.   }
  191.   return *this;
  192. }
  193.  
  194. bool operator==(mpz_class const &a, mpz_class const &b) {
  195.   int i, less;
  196.   bool isEqual = true;
  197.   if (a.n < b.n) {
  198.     less = a.n;
  199.   } else {
  200.     less = b.n;
  201.   }
  202.   for (i = 0; i < less; i++) {
  203.     if ((a.b[i]) != (b.b[i])) {
  204.       isEqual = false;
  205.     }      
  206.   }
  207.   if (a.n < b.n) {
  208.     for (; i < b.n; i++)
  209.       if (b.b[i] != 0) {
  210.         isEqual = false;
  211.       }
  212.   } else {
  213.     for (;i < a.n; i++)
  214.       if (a.b[i] != 0) {
  215.         isEqual = false;
  216.       }
  217.   }
  218.   return isEqual;
  219. }
  220.  
  221. mpz_class operator+(mpz_class const &a, mpz_class const &b) {
  222.   mpz_class r;
  223.   mpz_class const * more_ab;
  224.   int less_n;
  225.   if (a.n > b.n) {
  226.     r.n = a.n;
  227.     more_ab = &a;
  228.     less_n = b.n;
  229.   } else {
  230.     r.n = b.n;
  231.     more_ab = &b;
  232.     less_n = a.n;
  233.   }
  234.   r.b = new int [r.n];
  235.   int i, carry = 0;
  236.   for (i = 0; i < less_n; i++) {
  237.     int s = a.b[i] + b.b[i] + carry;
  238.     if (s >= CMASK) {
  239.       carry = 1;
  240.       r.b[i] = (s  & BMASK);
  241.     } else {
  242.       carry = 0;
  243.       r.b[i] = s; /* need */
  244.     }
  245.   }
  246.   for (;i < r.n; i++) {
  247.     int s = more_ab->b[i] + carry;
  248.     if (s >= CMASK) {
  249.       carry = 1;
  250.       r.b[i] = (s  & BMASK); /* need */
  251.     } else {
  252.       carry = 0;
  253.       r.b[i] = (s  & BMASK); /* need */ /* not to insert break keyword.*/
  254.     }
  255.   }
  256.   if (carry) {
  257.     int *new_body = new int [r.n + 1];
  258.     for (i = 0; i < r.n; i++)
  259.       new_body[i] = r.b[i];
  260.     new_body[i] = carry;
  261.     delete [] r.b;
  262.     r.b = new_body;
  263.     r.b[r.n] = 1; /* need */
  264.     r.n++;
  265.   }
  266.   return r;
  267. }
  268.  
  269. mpz_class sub(mpz_class const &a, mpz_class const &b) {
  270.   int vCount, moreCount, i;
  271.   int diff;
  272.   int borrowR = 0;
  273.   bool flagMore, flagLess;
  274.   mpz_class r;
  275.   flagMore = false;
  276.   flagLess = false;
  277.   moreCount = a.n;
  278.   vCount = a.n;
  279.   r = a;
  280.   if (a.n > b.n) {
  281.     flagMore = true;
  282.     vCount = b.n;
  283.     moreCount = a.n;
  284.     r = a;
  285.   }
  286.   if (a.n < b.n) {
  287.     flagLess = true;
  288.     vCount = a.n;  
  289.     moreCount = b.n;
  290.     r = b;
  291.   }
  292.  
  293.   borrowR = 0;
  294.   for (i = 0; i < vCount; i++) {
  295.     diff = a.b[i] - b.b[i] - borrowR;
  296.     if (diff < 0) {
  297.       borrowR = 1;
  298.       r.b[i] = diff + CMASK;
  299.     } else {
  300.       borrowR = 0;
  301.       r.b[i] = diff;
  302.     }
  303.   }
  304.   if (flagMore) {
  305.     for (;i < moreCount; i++) {
  306.       diff = a.b[i] - borrowR;
  307.       if (diff < 0) {
  308.         borrowR = 1;
  309.         r.b[i] = diff + CMASK;
  310.       } else {
  311.         borrowR = 0;
  312.         r.b[i] = diff;
  313.       }
  314.     }
  315.   }
  316.   if (flagLess) {
  317.     for (;i < moreCount; i++) {
  318.       if (b.b[i] > 0) {
  319.         borrowR = 1;
  320.       }
  321.     }
  322.   }
  323.   if (borrowR > 0) throw std::underflow_error("negative value.");
  324.   return r;
  325. }
  326.  
  327. mpz_class mul2(mpz_class const &a, mpz_class const &b) {
  328.   mpz_class Rshift = b;
  329.   mpz_class Lshift = a;
  330.   mpz_class r = 0;
  331.   int lsb;
  332.   while (!mpz_class::iszero(Rshift)) {
  333.     mpz_class::RightShift(Rshift, lsb);
  334.     if (lsb) {
  335.       r = r + Lshift;
  336.     }
  337.     mpz_class::LeftShift(Lshift);
  338.   }
  339.   return r;
  340. }
  341.  
  342. void mpz_class::LeftShift(mpz_class &a) {
  343.   int i, msb;
  344.   msb = 0;
  345.   for (i = 0; i < a.n; i++) {
  346.     a.b[i] *= 2;
  347.     if (msb > 0)
  348.       a.b[i] = (a.b[i] | 0x01);
  349.  
  350.     if (a.b[i] > BMASK) {
  351.       msb = 1;
  352.       a.b[i] -= CMASK;
  353.     } else {
  354.       msb = 0;
  355.     }
  356.   }
  357.   if (msb > 0) {
  358.     int *newbody;
  359.     newbody = new int [a.n + 1];
  360.     newbody[a.n] = 1;
  361.     for (int i = 0; i < a.n; i++) {
  362.       newbody[i] = a.b[i];
  363.     }
  364.     a.n++;
  365.     delete a.b;
  366.     a.b = newbody;
  367.   }
  368. }
  369.  
  370. void mpz_class::RightShift(mpz_class &a, int &lsb_out) {
  371.   int i, lsb;
  372.   lsb = 0;
  373.   for (i = a.n - 1; i >= 0; --i) {
  374.     if (lsb > 0)
  375.       a.b[i] = (a.b[i] | CMASK);
  376.  
  377.     if ((a.b[i] & 0x01) > 0)
  378.       lsb = 1;
  379.     else
  380.       lsb = 0;
  381.     a.b[i] /= 2;
  382.   }
  383.   lsb_out = lsb;
  384.   if (a.b[a.n - 1] == 0) {
  385.     int *newbody;
  386.     newbody = new int [a.n - 1];
  387.     for (int i = 0; i < a.n - 1; i++)
  388.       newbody[i] = a.b[i];
  389.     a.n--;
  390.     delete a.b;
  391.     a.b = newbody;
  392.   }
  393. }
  394.  
  395. void mpz_class::LeftShift2(int &msb, mpz_class &a) {
  396.   int i, msb_v;
  397.   msb_v = 0;
  398.   for (i = 0; i < a.n; i++) {
  399.     a.b[i] *= 2;
  400.     if (msb_v > 0)
  401.       a.b[i] = (a.b[i] | 0x01);
  402.  
  403.     if (a.b[i] > BMASK) {
  404.       msb_v = 1;
  405.       a.b[i] -= CMASK;
  406.     } else {
  407.       msb_v = 0;
  408.     }
  409.   }
  410.   msb = msb_v;
  411. }
  412.  
  413. void mpz_class::LeftShift3(mpz_class &a, int msb) {
  414.   int i, msb_v;
  415.   msb_v = msb;
  416.   for (i = 0; i < a.n; i++) {
  417.     a.b[i] *= 2;
  418.     if (msb_v > 0)
  419.       a.b[i] = (a.b[i] | 0x01);
  420.  
  421.     if (a.b[i] > BMASK) {
  422.       msb_v = 1;
  423.       a.b[i] -= CMASK;
  424.     } else {
  425.       msb_v = 0;
  426.     }
  427.   }
  428.   if (msb_v > 0) {
  429.     int *newbody;
  430.     newbody = new int [a.n + 1];
  431.     newbody[a.n] = 1;
  432.     for (int i = 0; i < a.n; i++) {
  433.       newbody[i] = a.b[i];
  434.     }
  435.     a.n++;
  436.     delete a.b;
  437.     a.b = newbody;
  438.   }
  439. }
  440.  
  441. void div2(mpz_class const &a, mpz_class const &b, mpz_class &q, mpz_class &r) {
  442.   if (b == mpz_class(0)) {
  443.     std::cout << "null devided." << std::endl;
  444.     return;
  445.   }
  446.  
  447.   /* LeftRest <- LeftShiftBuff */
  448.   mpz_class LeftRest(0);
  449.   mpz_class LeftShiftBuff = a;
  450.   mpz_class Quotient(0);
  451.   int msb, lsb;
  452.   msb = 0;
  453.   int n = a.n * BIT;
  454.   while (n > 0) {
  455. //    std::cout << "LeftShiftBuff=" << LeftShiftBuff << std::endl;
  456. //    std::cout << "LeftRest=" << LeftRest << std::endl;
  457.     mpz_class::LeftShift2(msb, LeftShiftBuff);
  458.     mpz_class::LeftShift3(LeftRest, msb);
  459.     lsb = 0;
  460.     if (b <= LeftRest) {
  461.       LeftRest = LeftRest - b;
  462.       lsb = 1;
  463.     }
  464.     mpz_class::LeftShift3(Quotient, lsb);
  465. //    std::cout << "Quotient=" << Quotient << std::endl;
  466.     n--;
  467.   }
  468.   q = Quotient;
  469.   r = LeftRest;
  470. }
  471.  
  472. bool operator<(mpz_class const &a, mpz_class const &b) {
  473.   try {
  474.     sub(a, b);
  475.   } catch (std::underflow_error &e) {
  476.     return true;
  477.   }
  478.   return false;
  479. }
  480.  
  481. bool mpz_class::iszero(mpz_class &n) {
  482.   bool isZero = true;
  483.   for (int i = 0; i < n.n; i++)
  484.     if (n.b[i] != 0) {
  485.       isZero = false;
  486.       break;
  487.     }
  488.   return isZero;
  489. }
  490.  
  491. std::ostream &operator<<(std::ostream &stream, mpz_class c) {
  492.   std::basic_string<char> mod10;
  493.   mpz_class q, n;
  494.   int mod;
  495.   bool fzero;
  496.  
  497.   fzero = true;
  498.   n = c;
  499.   for(;;) {
  500.     /* div10(n, q, mod); */
  501.     {
  502.       q = n;
  503.       mod = 0;
  504.       for (int i = 0; i < n.n * BIT; i++) {
  505.         /* shift_div(mod, q); */
  506.         {
  507.           int cy, i;
  508.  
  509.           cy = 0;
  510.           for (i = 0; i < n.n; i++) {
  511.             q.b[i] = q.b[i] << 1;
  512.             if (cy)
  513.               q.b[i] = q.b[i] | 0x01;
  514.             cy = (q.b[i] >> BIT);
  515.             q.b[i] = q.b[i] & BMASK;
  516.           }
  517.           mod = mod << 1;
  518.           if (cy)
  519.             mod = mod | 0x01;
  520.  
  521.         }
  522.         if (mod >= 10) {
  523.           q.b[0] = q.b[0] | 0x01;
  524.           mod = mod - 10;
  525.         }
  526.       }
  527.     }
  528.     if (mpz_class::iszero(q) && mod == 0)
  529.       break;
  530.     mod10 = (char)('0' + mod) + mod10;
  531.     fzero = false;
  532.     n = q;
  533.   }
  534.   if (fzero)
  535.     stream << '0';
  536.   else
  537.     stream << mod10;
  538.   return stream;
  539. }
  540. } /* namespace QZ */
  541.  
  542. int main() {
  543.   QZ::mpz_class s;
  544.   QZ::mpz_class t;
  545.   QZ::mpz_class const1;
  546.   int prec = 16;
  547.  
  548.   const1 = 1;
  549.   for (int n = 0; n < prec; n++)
  550.     const1 *= 2; /* const1 << prec */
  551.   s = const1;
  552.   int n;
  553.   for (n = 2; ; n++) {
  554.     t = const1 / n;
  555.     if (t == 0) break;
  556.     if (n % 2 == 0)
  557.       s -= t;
  558.     else
  559.       s += t;
  560.   }
  561.   std::cout << "n = " << n << std::endl;
  562.   s.QZ::mpz_class::print_tmp1(prec);
  563.   return 0;
  564. }
  565.  
  566. /* end */
  567.  
Success #stdin #stdout 0.76s 17088KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
n = 65537
1011000101011100
https://ideone.com/jAQCcF
language:
C++ (gcc 8.3)
created:
7 years ago
visibility:
public

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!