fork(2) download
copy 
  1. #include <iostream>
  2. #include <cstdlib>
  3. #include <cstdint>
  4. #include <cassert>
  5. #include <string>
  6.  
  7. typedef int64_t CONTAINER;
  8. unsigned int BIT = 62; /* BIT <= 64bit - 2bit = 62bit */
  9.  
  10. CONTAINER CMASK = 1ULL << BIT;
  11. CONTAINER BMASK = CMASK - 1;
  12.  
  13. namespace QZ {
  14.   class mpz_base_class {
  15.  
  16.   private:
  17.     unsigned int n;
  18.     CONTAINER *b;
  19.  
  20.   private:
  21.     static void LeftShift(mpz_base_class &a);
  22.     static void RightShift(mpz_base_class &a, int &lsb_out);
  23.     static void LeftShift2(int &msb, mpz_base_class &a);
  24.     static void LeftShift3(mpz_base_class &a, int msb);
  25.  
  26.   protected:
  27.     static bool iszero(mpz_base_class const &n);
  28.  
  29.   public:
  30.     mpz_base_class(); //
  31.     virtual ~mpz_base_class();
  32.     mpz_base_class(const signed long c);
  33.     mpz_base_class(mpz_base_class const &ob);
  34.     mpz_base_class &operator=(mpz_base_class const &ob);
  35.  
  36.     friend bool operator==(mpz_base_class const &a, mpz_base_class const &b);
  37.     friend bool operator!=(mpz_base_class const &a, mpz_base_class const &b) { return !(a == b); }
  38.  
  39.     friend std::ostream &operator<<(std::ostream &stream, mpz_base_class c);
  40.  
  41.     friend mpz_base_class operator+(mpz_base_class const &a, mpz_base_class const &b);
  42.     friend mpz_base_class operator-(mpz_base_class const &a, mpz_base_class const &b);
  43.  
  44.     friend bool operator<(mpz_base_class const &a, mpz_base_class const &b);
  45.     friend bool operator<=(mpz_base_class const &a, mpz_base_class const &b) { return (a < b) ? true : (a == b) ? true : false; }
  46.     friend bool operator>(mpz_base_class const &a, mpz_base_class const &b) { return (a < b) ? false : (a == b) ? false : true; }
  47.     friend bool operator>=(mpz_base_class const &a, mpz_base_class const &b) { return (a < b) ? false : (a == b) ? true : true; }
  48.  
  49.     friend mpz_base_class operator*(mpz_base_class const &a, mpz_base_class const &b);
  50.     friend void div2(mpz_base_class const &a, mpz_base_class const &b, mpz_base_class &q, mpz_base_class &r);
  51.     friend mpz_base_class operator/(mpz_base_class const &a, mpz_base_class const &b) { mpz_base_class q, r; div2(a, b, q, r); return q; }
  52.     friend mpz_base_class operator%(mpz_base_class const &a, mpz_base_class const &b) { mpz_base_class q, r; div2(a, b, q, r); return r; }
  53.  
  54.     bool testBit(int n);
  55.     mpz_base_class operator>>(int n) { int dmy; mpz_base_class r = *this; for (int i = 0; i < n; i++) RightShift(r, dmy); return r; }
  56.     mpz_base_class operator<<(int n) { mpz_base_class r = *this; for (int i = 0; i < n; i++) LeftShift(r); return r; }
  57.     int bitLength();
  58.  
  59.     unsigned int get_mpzclass_n() { return this->n; }
  60.     CONTAINER *get_mpzclass_b() { return this->b; }
  61.  
  62.     virtual void dump() const {
  63.       std::cout << "n = " << this->n << ": ";
  64.       for (unsigned int i = 0; i < this->n; i++)
  65.         std::cout << i << ":" << this->b[i] << ",";
  66.       std::cout << std::endl;
  67.     }
  68.  
  69.     void test42() {
  70.       // *this == 4 + 256;
  71.       this->b[2] = 0;
  72.     }
  73.   };
  74.  
  75.   class mpz_class : public mpz_base_class {
  76.   private:
  77.     bool negative;
  78.  
  79.   public:
  80.     mpz_class();
  81.     ~mpz_class();
  82.     mpz_class(const signed long c);
  83.     mpz_class(mpz_class const &ob);
  84.     mpz_class &operator=(mpz_class const &ob);
  85.     mpz_class &operator=(mpz_base_class const &ob);
  86.  
  87.     friend bool operator==(mpz_class const &a, mpz_class const &b);
  88.     friend bool operator!=(mpz_class const &a, mpz_class const &b) { return !(a == b); }
  89.  
  90.     friend std::ostream &operator<<(std::ostream &stream, mpz_class c);
  91.  
  92.     friend mpz_class operator-(mpz_class const &a);
  93.  
  94.     friend bool operator<(mpz_class const &a, mpz_class const &b);
  95.     friend bool operator<=(mpz_class const &a, mpz_class const &b) { return (a < b) ? true : (a == b) ? true : false; }
  96.     friend bool operator>(mpz_class const &a, mpz_class const &b) { return (a < b) ? false : (a == b) ? false : true; }
  97.     friend bool operator>=(mpz_class const &a, mpz_class const &b) { return (a < b) ? false : (a == b) ? true : true; }
  98.  
  99.     friend mpz_class operator+(mpz_class const &a, mpz_class const &b);
  100.     friend mpz_class operator-(mpz_class const &a, mpz_class const &b);
  101.  
  102.     friend mpz_class operator*(mpz_class const &a, mpz_class const &b);
  103.     friend mpz_class operator/(mpz_class const &a, mpz_class const &b);
  104.     friend mpz_class operator%(mpz_class const &a, mpz_class const &b);
  105.  
  106.     mpz_class &operator+=(mpz_class const &n) { *this = *this + n; return *this; }
  107.     mpz_class &operator-=(mpz_class const &n) { *this = *this - n; return *this; }
  108.     mpz_class &operator*=(mpz_class const &n) { *this = *this * n; return *this; }
  109.     mpz_class &operator/=(mpz_class const &n) { *this = *this / n; return *this; }
  110.     mpz_class &operator%=(mpz_class const &n) { *this = *this % n; return *this; }
  111.  
  112.     mpz_class operator>>(int n) { mpz_class r; r = this->mpz_base_class::operator>>(n); return r; }
  113.     mpz_class operator<<(int n) { mpz_class r; r = this->mpz_base_class::operator<<(n); return r; }
  114.  
  115.     bool get_mpzclass_negative() { return this->negative; }
  116.  
  117.     virtual void dump() const {
  118.       mpz_base_class::dump();
  119.       std::cout << "sign = " << this->negative << std::endl;
  120.     }
  121.   };
  122.  
  123.   mpz_base_class::mpz_base_class() {
  124. //    std::cout << "mpz_base_class constructor1" << std::endl;
  125.     this->n = 0;
  126.     this->b = 0;
  127.   }  
  128.   mpz_base_class::~mpz_base_class() {
  129. //    std::cout << "mpz_base_class destructor" << std::endl;
  130.     if (this->b) delete [] this->b;
  131.     this->n = 0;
  132.     this->b = 0;
  133.   }
  134.   mpz_base_class::mpz_base_class(const signed long c) {
  135. //    std::cout << "mpz_base_class constructor2" << std::endl;
  136.     signed long cc;
  137.     cc = c;
  138.     if (c == 0) {
  139.       this->n = 1;
  140.       this->b = new CONTAINER [1];
  141.       this->b[0] = 0;
  142.       return;
  143.     }
  144.     if (c < 0) {
  145.       cc = -c; /* for mpz_class not for mpz_base_class */
  146.     }
  147.     assert(cc > 0);
  148.  
  149.     CONTAINER *newbody;
  150.     unsigned int i, j, r, pos;
  151.  
  152.     this->b = 0;
  153.     this->n = 0;
  154.  
  155.     while (cc > 0) {
  156.       r = 0;
  157.       pos = 1;
  158.       for (i = 0; i < BIT; i++) {
  159.         if ((cc & 1) != 0) {
  160.           r = (r | pos);
  161.         } else {
  162.         }
  163.         pos = pos * 2;
  164.         cc = cc / 2;
  165.       }
  166.       newbody = new CONTAINER [this->n + 1];
  167.       for (j = 0; j < this->n; j++)
  168.         newbody[j] = this->b[j];
  169.       newbody[j] = r;
  170.       delete [] this->b;
  171.       this->b = newbody;
  172.       this->n++;
  173.     }
  174.     assert(this->n != 0);
  175.   }
  176.  
  177.   mpz_base_class::mpz_base_class(const mpz_base_class &ob) {
  178. //    std::cout << "mpz_base_class constructor3" << std::endl;
  179.     this->n = ob.n;
  180.     if (ob.n == 0) {
  181.       this->b = 0;
  182.     } else {
  183.       this->b = new CONTAINER [ob.n];
  184.       for (unsigned int i = 0; i < ob.n; i++)
  185.         this->b[i] = ob.b[i];
  186.     }
  187.     assert((this->n == 0 && this->b == 0) || (this->n > 0 && this->b != 0));
  188.   }
  189.  
  190.   mpz_base_class &mpz_base_class::operator=(const mpz_base_class &ob) {
  191. //    std::cout << "mpz_base_class operator=" << std::endl;
  192.     CONTAINER *org_b = this->b; // for checking self-assignment
  193.     this->n = ob.n;
  194.     if (ob.n == 0) {
  195.       this->b = 0;
  196.       delete [] org_b;
  197.     } else {
  198.       this->b = new CONTAINER [ob.n];
  199.       for (unsigned int i = 0; i < ob.n; i++)
  200.         this->b[i] = ob.b[i];
  201.       delete [] org_b;
  202.     }
  203.     assert((this->n == 0 && this->b == 0) || (this->n > 0 && this->b != 0));
  204.     return *this;
  205.   }
  206.  
  207.   bool operator==(mpz_base_class const &a, mpz_base_class const &b) {
  208. //    std::cout << "mpz_base_class operator==" << std::endl;
  209.     unsigned int i, less;
  210.     bool isEqual = true;
  211.  
  212.     if (a.n < b.n)
  213.       less = a.n;
  214.     else
  215.       less = b.n;
  216.     for (i = 0; i < less; i++)
  217.       if (a.b[i] != b.b[i]) {
  218.         isEqual = false;
  219.         goto label;
  220.       }
  221.     if (a.n < b.n) {
  222.       for (; i < b.n; i++)
  223.         if (b.b[i] != 0) {
  224.           isEqual = false;
  225.           goto label;
  226.         }
  227.     } else {
  228.       for (; i < a.n; i++)
  229.         if (a.b[i] != 0) { // be carefull! not b.b[] but a.b[]
  230.           isEqual = false;
  231.           goto label;
  232.         }
  233.     }
  234.   label:
  235.     return isEqual;
  236.   }
  237.  
  238.  
  239.   std::ostream &operator<<(std::ostream &stream, mpz_base_class c) {
  240. //    std::cout << "mpz_base_class operator<<" << std::endl;
  241.     std::basic_string<char> mod10;
  242.     mpz_base_class q, n;
  243.     int mod;
  244.     bool fzero;
  245.  
  246.     fzero = true;
  247.     n = c;
  248.     for(;;) {
  249.       /* div10(n, q, mod); */
  250.       {
  251.         q = n;
  252.         mod = 0;
  253.         for (unsigned int i = 0; i < n.n * BIT; i++) {
  254.           /* shift_div(mod, q); */
  255.           {
  256.             unsigned int cy, i;
  257.  
  258.             cy = 0;
  259.             for (i = 0; i < n.n; i++) {
  260.               q.b[i] = q.b[i] << 1;
  261.               if (cy)
  262.                 q.b[i] = q.b[i] | 0x01;
  263.               cy = (q.b[i] >> BIT);
  264.               q.b[i] = q.b[i] & BMASK;
  265.             }
  266.             mod = mod << 1;
  267.             if (cy)
  268.               mod = mod | 0x01;
  269.  
  270.           }
  271.           if (mod >= 10) {
  272.             q.b[0] = q.b[0] | 0x01;
  273.             mod = mod - 10;
  274.           }
  275.         }
  276.       }
  277.       if (mpz_class::iszero(q) && mod == 0)
  278.         break;
  279.       mod10 = (char)('0' + mod) + mod10;
  280.       fzero = false;
  281.       n = q;
  282.     }
  283.  
  284.     if (fzero)
  285.       stream << '0';
  286.     else
  287.       stream << mod10;
  288.     return stream;
  289.   }
  290.  
  291.   mpz_base_class operator+(mpz_base_class const &a, mpz_base_class const &b) {
  292.     mpz_base_class r;
  293.     mpz_base_class const * more_ab;
  294.     unsigned int less_n;
  295.     if (a.n > b.n) {
  296.       r.n = a.n;
  297.       more_ab = &a;
  298.       less_n = b.n;
  299.     } else {
  300.       r.n = b.n;
  301.       more_ab = &b;
  302.       less_n = a.n;
  303.     }
  304.     r.b = new CONTAINER [r.n];
  305.     unsigned int i, carry = 0;
  306.     for (i = 0; i < less_n; i++) {
  307.       CONTAINER s = a.b[i] + b.b[i] + carry;
  308.       if (s >= CMASK) {
  309.         carry = 1;
  310.         r.b[i] = (s  & BMASK);
  311.       } else {
  312.         carry = 0;
  313.         r.b[i] = s; /* need */
  314.       }
  315.     }
  316.     for (;i < r.n; i++) {
  317.       CONTAINER s = more_ab->b[i] + carry;
  318.       if (s >= CMASK) {
  319.         carry = 1;
  320.         r.b[i] = (s  & BMASK); /* need */
  321.       } else {
  322.         carry = 0;
  323.         r.b[i] = (s  & BMASK); /* need */ /* do not to insert break keyword.*/
  324.       }
  325.     }
  326.     if (carry) {
  327.       CONTAINER *new_body = new CONTAINER [r.n + 1];
  328.       for (i = 0; i < r.n; i++)
  329.         new_body[i] = r.b[i];
  330.       new_body[i] = carry;
  331.       delete [] r.b;
  332.       r.b = new_body;
  333.       // r.b[r.n] = 1; /* NOT need */
  334.       r.n++;
  335.     }
  336.     return r;
  337.   }
  338.  
  339.  
  340.   bool mpz_base_class::iszero(mpz_base_class const &n) {
  341.     bool isZero = true;
  342.     for (unsigned int i = 0; i < n.n; i++)
  343.       if (n.b[i] != 0) {
  344.         isZero = false;
  345.         break;
  346.       }
  347.     return isZero;
  348.   }
  349.  
  350.   mpz_base_class operator-(mpz_base_class const &a, mpz_base_class const &b) {
  351.     int vCount, moreCount, i;
  352.     CONTAINER diff;
  353.     int borrowR = 0;
  354.     bool flagMore, flagLess;
  355.     mpz_base_class r;
  356.     flagMore = false;
  357.     flagLess = false;
  358.     moreCount = a.n;
  359.     vCount = a.n;
  360.     r = a;
  361.     if (a.n > b.n) {
  362.       flagMore = true;
  363.       vCount = b.n;
  364.       moreCount = a.n;
  365.       r = a;
  366.     }
  367.     if (a.n < b.n) {
  368.       flagLess = true;
  369.       vCount = a.n;  
  370.       moreCount = b.n;
  371.       r = b;
  372.     }
  373.  
  374.     borrowR = 0;
  375.     for (i = 0; i < vCount; i++) {
  376.       diff = a.b[i] - b.b[i] - borrowR;
  377.       if (diff < 0) {
  378.         borrowR = 1;
  379.         r.b[i] = diff + CMASK;
  380.       } else {
  381.         borrowR = 0;
  382.         r.b[i] = diff;
  383.       }
  384.     }
  385.     if (flagMore) {
  386.       for (;i < moreCount; i++) {
  387.         diff = a.b[i] - borrowR;
  388.         if (diff < 0) {
  389.           borrowR = 1;
  390.           r.b[i] = diff + CMASK;
  391.         } else {
  392.           borrowR = 0;
  393.           r.b[i] = diff;
  394.         }
  395.       }
  396.     }
  397.     if (flagLess) {
  398.       for (;i < moreCount; i++) {
  399.         if (b.b[i] > 0) {
  400.           borrowR = 1;
  401.         }
  402.       }
  403.     }
  404.     if (borrowR > 0) throw std::underflow_error("negative value.");
  405.     return r;
  406.   }
  407.  
  408.   bool operator<(mpz_base_class const &a, mpz_base_class const &b) { // if a < b then true
  409.     try {
  410.       mpz_base_class r = a - b;
  411.     } catch (std::underflow_error &e) {
  412.       return true;
  413.     }
  414.     return false;
  415.   }
  416.  
  417.   mpz_base_class operator*(mpz_base_class const &a, mpz_base_class const &b) {
  418.     mpz_base_class Rshift = b;
  419.     mpz_base_class Lshift = a;
  420.     mpz_base_class r = 0;
  421.     int lsb;
  422.     while (!mpz_base_class::iszero(Rshift)) {
  423.       mpz_base_class::RightShift(Rshift, lsb);
  424.       if (lsb) {
  425.         r = r + Lshift;
  426.       }
  427.       mpz_base_class::LeftShift(Lshift);
  428.     }
  429.     return r;
  430.   }
  431.  
  432.   void mpz_base_class::LeftShift(mpz_base_class &a) {
  433.     unsigned int i, msb;
  434.     msb = 0;
  435.     for (i = 0; i < a.n; i++) {
  436.       a.b[i] *= 2;
  437.       if (msb > 0)
  438.         a.b[i] = (a.b[i] | 0x01);
  439.  
  440.       if (a.b[i] > BMASK) {
  441.         msb = 1;
  442.         a.b[i] -= CMASK;
  443.       } else {
  444.         msb = 0;
  445.       }
  446.     }
  447.     if (msb > 0) {
  448.       CONTAINER *newbody;
  449.       newbody = new CONTAINER [a.n + 1];
  450.       newbody[a.n] = 1;
  451.       for (unsigned int i = 0; i < a.n; i++) {
  452.         newbody[i] = a.b[i];
  453.       }
  454.       a.n++;
  455.       delete [] a.b;
  456.       a.b = newbody;
  457.     }
  458.   }
  459.  
  460.   void mpz_base_class::RightShift(mpz_base_class &a, int &lsb_out) {
  461.     int i, lsb;
  462.     lsb = 0;
  463.     for (i = a.n - 1; i >= 0; --i) {
  464.       if (lsb > 0)
  465.         a.b[i] = (a.b[i] | CMASK);
  466.  
  467.       if ((a.b[i] & 0x01) > 0)
  468.         lsb = 1;
  469.       else
  470.         lsb = 0;
  471.       a.b[i] /= 2;
  472.     }
  473.     lsb_out = lsb;
  474.     if (a.b[a.n - 1] == 0) {
  475.       CONTAINER *newbody;
  476.       newbody = new CONTAINER [a.n - 1];
  477.       for (unsigned int i = 0; i < a.n - 1; i++)
  478.         newbody[i] = a.b[i];
  479.       a.n--;
  480.       delete [] a.b;
  481.       a.b = newbody;
  482.     }
  483.   }
  484.  
  485.   void mpz_base_class::LeftShift2(int &msb, mpz_base_class &a) {
  486.     unsigned int i, msb_v;
  487.     msb_v = 0;
  488.     for (i = 0; i < a.n; i++) {
  489.       a.b[i] *= 2;
  490.       if (msb_v > 0)
  491.         a.b[i] = (a.b[i] | 0x01);
  492.  
  493.       if (a.b[i] > BMASK) {
  494.         msb_v = 1;
  495.         a.b[i] -= CMASK;
  496.       } else {
  497.         msb_v = 0;
  498.       }
  499.     }
  500.     msb = msb_v;
  501.   }
  502.  
  503.   void mpz_base_class::LeftShift3(mpz_base_class &a, int msb) {
  504.     unsigned int i, msb_v;
  505.     msb_v = msb;
  506.     for (i = 0; i < a.n; i++) {
  507.       a.b[i] *= 2;
  508.       if (msb_v > 0)
  509.         a.b[i] = (a.b[i] | 0x01);
  510.  
  511.       if (a.b[i] > BMASK) {
  512.         msb_v = 1;
  513.         a.b[i] -= CMASK;
  514.       } else {
  515.         msb_v = 0;
  516.       }
  517.     }
  518.     if (msb_v > 0) {
  519.       CONTAINER *newbody;
  520.       newbody = new CONTAINER [a.n + 1];
  521.       newbody[a.n] = 1;
  522.       for (unsigned int i = 0; i < a.n; i++) {
  523.         newbody[i] = a.b[i];
  524.       }
  525.       a.n++;
  526.       delete [] a.b;
  527.       a.b = newbody;
  528.     }
  529.   }
  530.   void div2(mpz_base_class const &a, mpz_base_class const &b, mpz_base_class &q, mpz_base_class &r) {
  531.     if (b == mpz_base_class(0)) {
  532.       std::cout << "null devided." << std::endl;
  533.       return;
  534.     }
  535.  
  536.     /* LeftRest <- LeftShiftBuff */
  537.     mpz_base_class LeftRest(0);
  538.     mpz_base_class LeftShiftBuff = a;
  539.     mpz_base_class Quotient(0);
  540.     int msb, lsb;
  541.     msb = 0;
  542.     int n = a.n * BIT;
  543.     while (n > 0) {
  544.       //    std::cout << "LeftShiftBuff=" << LeftShiftBuff << std::endl;
  545.       //    std::cout << "LeftRest=" << LeftRest << std::endl;
  546.       mpz_base_class::LeftShift2(msb, LeftShiftBuff);
  547.       mpz_base_class::LeftShift3(LeftRest, msb);
  548.       lsb = 0;
  549.       if (b < LeftRest || b == LeftRest) {
  550.         LeftRest = LeftRest - b;
  551.         lsb = 1;
  552.       }
  553.       mpz_base_class::LeftShift3(Quotient, lsb);
  554.       //    std::cout << "Quotient=" << Quotient << std::endl;
  555.       n--;
  556.     }
  557.     q = Quotient;
  558.     r = LeftRest;
  559.   }
  560.  
  561.   bool mpz_base_class::testBit(int n) {
  562.     unsigned int test_byte = n / BIT;
  563.     int test_mask = 1 << (n % BIT);
  564.     if (test_byte >= this->n) return false;
  565.     if ((this->b[test_byte] & test_mask) > 0) return true;
  566.     return false;
  567.   }
  568.  
  569. /*----------------------------------------------*/
  570.   mpz_class::mpz_class() : mpz_base_class() {
  571. //    std::cout << "mpz_class constructor1" << std::endl;
  572.     this->negative = false;
  573.   }
  574.  
  575.   mpz_class::~mpz_class() {
  576. //    std::cout << "mpz_class destructor1" << std::endl;
  577.     this->negative = false;
  578.   }
  579.  
  580.   mpz_class::mpz_class(const signed long c) : mpz_base_class(c) {
  581. //    std::cout << "mpz_class costructor2" << std::endl;
  582.     if (c >= 0) {
  583.       this->negative = false;
  584.     } else {
  585.       this->negative = true;
  586.     }
  587.   }
  588.  
  589.   mpz_class::mpz_class(const mpz_class &ob) : mpz_base_class(ob) {
  590. //    std::cout << "mpz_class constructor3" << std::endl;
  591.     this->negative = ob.negative;
  592.   }
  593.  
  594.   mpz_class &mpz_class::operator=(mpz_class const &ob) {
  595. //    std::cout << "mpz_class operator=(mpz_class&)" << std::endl;
  596.     this->mpz_base_class::operator=(ob);
  597.     this->negative = ob.negative;
  598.     return *this;
  599.   }
  600.  
  601.   mpz_class &mpz_class::operator=(mpz_base_class const &ob) {
  602. //   std::cout << "mpz_class operator=(mpz_base_class&)" << std::endl;
  603.     this->mpz_base_class::operator=(ob);
  604.     this->negative = false;
  605.     return *this;
  606.   }
  607.  
  608.   bool operator==(mpz_class const &a, mpz_class const &b) {
  609. //    std::cout << "mpz_class::operator==" << std::endl;
  610.     if (a.negative == b.negative)
  611.       return operator==((mpz_base_class)a, (mpz_base_class)b);
  612.     return false;
  613.   }
  614.  
  615.   std::ostream &operator<<(std::ostream &stream, mpz_class c) {
  616.     if (c.negative)
  617.       stream << '-';
  618.     return operator<<(stream, (mpz_base_class)c);
  619.   }
  620.  
  621.   mpz_class operator-(mpz_class const &a) {
  622.     mpz_class r = a;
  623.     r.negative = !r.negative;
  624.     return r;
  625.   }
  626.  
  627.   mpz_class operator+(mpz_class const &a, mpz_class const &b) {
  628.     mpz_class r;
  629.     if (!a.negative && !b.negative) {
  630.       r = (mpz_base_class)a + (mpz_base_class)b; /* operator=(mpz_class, mpz_base_class) */
  631.       r.negative = false;
  632.     } else if (a.negative && !b.negative) {
  633.       mpz_class t = -a; assert(t.negative == false);
  634.       if (operator<((mpz_base_class)t, (mpz_base_class)b) == true) { /* t < b */
  635.         r = (mpz_base_class)b - (mpz_base_class)t;
  636.         r.negative = false;
  637.       } else { /* b < t */
  638.         r = (mpz_base_class)t - (mpz_base_class)b;
  639.         r.negative = true;
  640.       }
  641.       if (mpz_base_class::iszero(r))
  642.         r.negative = false;
  643.     } else if (!a.negative && b.negative) {
  644.       mpz_class t = -b; assert(t.negative == false);
  645.       if (operator<((mpz_base_class)t, (mpz_base_class)a) == true) { /* t < a */
  646.         r = (mpz_base_class)a - (mpz_base_class)t;
  647.         r.negative = false;
  648.       } else { /* a < t */
  649.         r = (mpz_base_class)t - (mpz_base_class)a;
  650.         r.negative = true;
  651.       }
  652.       if (mpz_base_class::iszero(r))
  653.         r.negative = false;
  654.     } else if (a.negative && b.negative) {      
  655.       QZ::mpz_class aa = a, bb = b;
  656.       aa = -aa; bb = -bb;
  657.       r = (mpz_base_class)aa + (mpz_base_class)bb;
  658.       r.negative = true;
  659.       if (mpz_base_class::iszero(r))
  660.         r.negative = false;
  661.     }
  662.     return r;
  663.   }
  664.  
  665.   mpz_class operator-(mpz_class const &a, mpz_class const &b) {
  666.     QZ::mpz_class r = -b;
  667.     return a + r;
  668.   }
  669.  
  670.   bool operator<(mpz_class const &a, mpz_class const &b) {
  671.     QZ::mpz_class r = a - b;
  672.     if (r.negative)
  673.       return true;
  674.     return false;
  675.   }
  676.  
  677.   mpz_class operator*(mpz_class const &a, mpz_class const &b) {
  678.     mpz_class x = (a > 0) ? a : -a;
  679.     mpz_class y = (b > 0) ? b : -b;
  680.     assert(x >= 0 && y >= 0);
  681.     mpz_class r;
  682.     r = (mpz_base_class)x * (mpz_base_class)y;
  683.  
  684.     if (a >= 0 && b >= 0)  {
  685.     } else if (a < 0 && b >= 0) {
  686.       r = -r;
  687.     } else if (a >= 0 && b < 0) {
  688.       r = -r;
  689.     } if (a < 0 && b < 0) {
  690.     }
  691.     if (mpz_base_class::iszero(r)) r.negative = false;
  692.     return r;
  693.   }
  694.  
  695.   mpz_class operator/(mpz_class const &a, mpz_class const &b) {
  696.     mpz_class x = (a >= 0) ? a : -a;
  697.     mpz_class y = (b >= 0) ? b : -b;
  698.     mpz_class q, r;
  699.     div2(x, y, q, r);
  700.     if (a >= 0 && b >= 0) {
  701.     } else if (a < 0 && b >= 0) {
  702.       q = -q; r = -r;
  703.       if (mpz_base_class::iszero(q)) q.negative = false;
  704.       if (mpz_base_class::iszero(r)) r.negative = false;
  705.     } else if (a >= 0 && b < 0) {
  706.       q = -q;
  707.       if (mpz_base_class::iszero(q)) q.negative = false;
  708.       if (mpz_base_class::iszero(r)) r.negative = false;
  709.     } else if (a < 0 && b < 0) {
  710.       r = -r;
  711.       if (mpz_base_class::iszero(q)) q.negative = false;
  712.       if (mpz_base_class::iszero(r)) q.negative = false;
  713.     }
  714.     return q;
  715.   }
  716.  
  717.   mpz_class operator%(mpz_class const &a, mpz_class const &b) {
  718.     mpz_class x = (a >= 0) ? a : -a;
  719.     mpz_class y = (b >= 0) ? b : -b;
  720.     mpz_class q, r;
  721.     div2(x, y, q, r);
  722.     if (a >= 0 && b >= 0) {
  723.     } else if (a < 0 && b >= 0) {
  724.       q = -q; r = -r;
  725.       if (mpz_base_class::iszero(q)) q.negative = false;
  726.       if (mpz_base_class::iszero(r)) r.negative = false;
  727.     } else if (a >= 0 && b < 0) {
  728.       q = -q;
  729.       if (mpz_base_class::iszero(q)) q.negative = false;
  730.       if (mpz_base_class::iszero(r)) r.negative = false;
  731.     } else if (a < 0 && b < 0) {
  732.       r = -r;
  733.       if (mpz_base_class::iszero(q)) q.negative = false;
  734.       if (mpz_base_class::iszero(r)) r.negative = false;
  735.     }
  736.     return r;
  737.   }
  738.  
  739.   int mpz_base_class::bitLength() {
  740.     if (this->n == 0)
  741.       return 0;
  742.  
  743.     int n = this->n;
  744.     int i;
  745.     for (i = n - 1; !(b[i] > 0); --i) {}
  746.  
  747.     CONTAINER msb = this->b[i];
  748.  
  749.     int j = 0;
  750.     while (msb > 0) {
  751.       j++;
  752.       msb = msb >> 1;
  753.     }
  754.     return BIT * i + j;
  755.   }
  756.  
  757. } /* namespace QZ */
  758. /*-----------------------------------------------------------------*/
  759.  
  760. QZ::mpz_class euclid(QZ::mpz_class a, QZ::mpz_class b) {
  761.   if (b == 0) return a;
  762.   return euclid(b, a % b);
  763. }
  764.  
  765. class Ratio {
  766.   QZ::mpz_class a;
  767.   QZ::mpz_class b;
  768. public:
  769.   Ratio() : a(0), b(0) {}
  770.   Ratio(QZ::mpz_class x, QZ::mpz_class y) : a(x), b(y) {}
  771.   friend std::ostream &operator<<(std::ostream &s, Ratio r) {
  772.     if (r.b == 1)
  773.       s << r.a;
  774.     else
  775.       s << "(" << r.a << "/" << r.b << ")";
  776.     return s;
  777.   }
  778.   friend Ratio operator+(Ratio const &x, Ratio const &y) {
  779.     Ratio r;
  780.     QZ::mpz_class c = euclid(x.b, y.b);
  781.     QZ::mpz_class zy = x.b / c;
  782.     QZ::mpz_class zx = y.b / c;
  783.     r.b = y.b * zy;
  784.     assert(r.b == x.b * zx);
  785.     r.a = x.a * zx + y.a * zy;
  786.     c = euclid(r.a, r.b);
  787.     r.a = r.a / c;
  788.     r.b = r.b / c;
  789.     return r;
  790.   }
  791.   friend Ratio operator-(Ratio const &x, Ratio const &y) {
  792.     Ratio r;
  793.     QZ::mpz_class c = euclid(x.b, y.b);
  794.     QZ::mpz_class zy = x.b / c;
  795.     QZ::mpz_class zx = y.b / c;
  796.     r.b = y.b * zy;
  797.     assert(r.b == x.b * zx);
  798.     r.a = x.a * zx - y.a * zy;
  799.     c = euclid(r.a, r.b);
  800.     r.a = r.a / c;
  801.     r.b = r.b / c;
  802.     return r;
  803.   }
  804.   friend Ratio operator*(Ratio const &x, Ratio const &y) {
  805.     Ratio r;
  806.     QZ::mpz_class xa, xb, ya, yb;
  807.     QZ::mpz_class c1 = euclid(x.a, y.b);
  808.     QZ::mpz_class c2 = euclid(x.b, y.a);
  809.     xa = x.a / c1; yb = y.b / c1; 
  810.     xb = x.b / c2; ya = y.a / c2; 
  811.     r.a = xa * ya;
  812.     r.b = xb * yb;
  813.     return r;
  814.   }
  815. };
  816.  
  817. class Field5 {
  818.   Ratio a;
  819.   Ratio b;
  820. public:
  821.   Field5() : a(Ratio()), b(Ratio()) {}
  822.   Field5(Ratio x, Ratio y) : a(x), b(y) {}
  823.   friend std::ostream &operator<<(std::ostream &s, Field5 r) {
  824. //    s << "(" << r.a << "," << r.b << ")";
  825.     s << r.b;
  826.     return s;
  827.   }
  828.   friend Field5 operator-(Field5 const &x, Field5 const &y) {
  829.     Field5 r;
  830.     r.a = x.a - y.a; r.b = x.b - y.b;
  831.     return r;
  832.   }
  833.   friend Field5 operator*(Field5 const &x, Field5 const &y) {
  834.     Field5 r;
  835.     r.a = x.a * y.a + Ratio(5, 1) * x.b * y.b;
  836.     r.b = x.a * y.b + x.b * y.a;
  837.     return r;
  838.   }  
  839. };
  840.  
  841. QZ::mpz_class fib(int n) {
  842.   static QZ::mpz_class t1 = 1, t2 = 1;
  843.   if (n == 1 || n == 2)
  844.     return 1;
  845.   QZ::mpz_class s = t1 + t2;
  846.   t1 = t2; t2 = s;
  847.   return s;
  848. }
  849.  
  850. Field5 Field5Power(Field5 a, int n) { /* =a^n */
  851.   unsigned int mask = 1 << (sizeof(int) * 8 - 1);
  852.   Field5 r(Ratio(1, 1), Ratio(0, 1));
  853.   for (;;) {
  854.     if ((n & mask) > 0)
  855.       r = r * a;
  856.     mask = (mask >> 1);
  857.     /* out of loop */
  858.     if (mask == 0)
  859.       break;
  860.     r = r * r;
  861.   }
  862.   return r;
  863. }
  864.  
  865. int const M = 200;
  866. int main() {
  867.   Field5 phi1(Ratio(1, 2), Ratio(1, 2));
  868.   Field5 phi2(Ratio(1, 2), Ratio(-1, 2));
  869.  
  870.   for (int n = 1; n <= M; n++) {
  871.     QZ::mpz_class s = fib(n);
  872.  
  873.     Field5 r1, r2, r3;
  874.     r1 = Field5Power(phi1, n);
  875.     r2 = Field5Power(phi2, n);
  876.     r3 = r1 - r2;
  877.  
  878.     std::cout << n << ":: " << s << " : " << r3 << std::endl;
  879.   }
  880.   return 0;
  881. }
  882. /* end */
  883.  
Time limit exceeded #stdin #stdout 5s 4388KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
1:: 1 : 1
2:: 1 : 1
3:: 2 : 2
4:: 3 : 3
5:: 5 : 5
6:: 8 : 8
7:: 13 : 13
8:: 21 : 21
9:: 34 : 34
10:: 55 : 55
11:: 89 : 89
12:: 144 : 144
13:: 233 : 233
14:: 377 : 377
15:: 610 : 610
16:: 987 : 987
17:: 1597 : 1597
18:: 2584 : 2584
19:: 4181 : 4181
20:: 6765 : 6765
21:: 10946 : 10946
22:: 17711 : 17711
23:: 28657 : 28657
24:: 46368 : 46368
25:: 75025 : 75025
26:: 121393 : 121393
27:: 196418 : 196418
28:: 317811 : 317811
29:: 514229 : 514229
30:: 832040 : 832040
31:: 1346269 : 1346269
32:: 2178309 : 2178309
33:: 3524578 : 3524578
34:: 5702887 : 5702887
35:: 9227465 : 9227465
36:: 14930352 : 14930352
37:: 24157817 : 24157817
38:: 39088169 : 39088169
39:: 63245986 : 63245986
40:: 102334155 : 102334155
41:: 165580141 : 165580141
42:: 267914296 : 267914296
43:: 433494437 : 433494437
44:: 701408733 : 701408733
45:: 1134903170 : 1134903170
46:: 1836311903 : 1836311903
47:: 2971215073 : 2971215073
48:: 4807526976 : 4807526976
49:: 7778742049 : 7778742049
50:: 12586269025 : 12586269025
51:: 20365011074 : 20365011074
52:: 32951280099 : 32951280099
53:: 53316291173 : 53316291173
54:: 86267571272 : 86267571272
55:: 139583862445 : 139583862445
56:: 225851433717 : 225851433717
57:: 365435296162 : 365435296162
58:: 591286729879 : 591286729879
59:: 956722026041 : 956722026041
60:: 1548008755920 : 1548008755920
61:: 2504730781961 : 2504730781961
62:: 4052739537881 : 4052739537881
63:: 6557470319842 : 6557470319842
64:: 10610209857723 : 10610209857723
65:: 17167680177565 : 17167680177565
66:: 27777890035288 : 27777890035288
https://ideone.com/IpZR5Y
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!