oct-inttypes.cc

Go to the documentation of this file.
00001 /*
00002 
00003 Copyright (C) 2004-2012 John W. Eaton
00004 Copyright (C) 2008-2009 Jaroslav Hajek
00005 
00006 This file is part of Octave.
00007 
00008 Octave is free software; you can redistribute it and/or modify it
00009 under the terms of the GNU General Public License as published by the
00010 Free Software Foundation; either version 3 of the License, or (at your
00011 option) any later version.
00012 
00013 Octave is distributed in the hope that it will be useful, but WITHOUT
00014 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00015 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00016 for more details.
00017 
00018 You should have received a copy of the GNU General Public License
00019 along with Octave; see the file COPYING.  If not, see
00020 <http://www.gnu.org/licenses/>.
00021 
00022 */
00023 
00024 #ifdef HAVE_CONFIG_H
00025 #include <config.h>
00026 #endif
00027 
00028 #include "lo-error.h"
00029 
00030 #include "oct-inttypes.h"
00031 
00032 template<class T>
00033 const octave_int<T> octave_int<T>::zero (static_cast<T> (0));
00034 
00035 template<class T>
00036 const octave_int<T> octave_int<T>::one (static_cast<T> (1));
00037 
00038 // define type names.
00039 #define DECLARE_OCTAVE_INT_TYPENAME(TYPE, TYPENAME) \
00040   template <> \
00041   OCTAVE_API const char * \
00042   octave_int<TYPE>::type_name () { return TYPENAME; }
00043 
00044 DECLARE_OCTAVE_INT_TYPENAME (int8_t, "int8")
00045 DECLARE_OCTAVE_INT_TYPENAME (int16_t, "int16")
00046 DECLARE_OCTAVE_INT_TYPENAME (int32_t, "int32")
00047 DECLARE_OCTAVE_INT_TYPENAME (int64_t, "int64")
00048 DECLARE_OCTAVE_INT_TYPENAME (uint8_t, "uint8")
00049 DECLARE_OCTAVE_INT_TYPENAME (uint16_t, "uint16")
00050 DECLARE_OCTAVE_INT_TYPENAME (uint32_t, "uint32")
00051 DECLARE_OCTAVE_INT_TYPENAME (uint64_t, "uint64")
00052 
00053 #ifndef OCTAVE_INT_USE_LONG_DOUBLE
00054 
00055 // Define comparison operators
00056 
00057 template <class xop>
00058 bool
00059 octave_int_cmp_op::emulate_mop (uint64_t x, double y)
00060 {
00061   static const double xxup = std::numeric_limits<uint64_t>::max ();
00062   // This converts to the nearest double. Unless there's an equality, the
00063   // result is clear.
00064   double xx = x;
00065   if (xx != y)
00066     return xop::op (xx, y);
00067   else
00068     {
00069       // If equality occured we compare as integers.
00070       if (xx == xxup)
00071         return xop::gtval;
00072       else
00073         return xop::op (x, static_cast<uint64_t> (xx));
00074     }
00075 }
00076 
00077 template <class xop>
00078 bool
00079 octave_int_cmp_op::emulate_mop (int64_t x, double y)
00080 {
00081   static const double xxup = std::numeric_limits<int64_t>::max ();
00082   static const double xxlo = std::numeric_limits<int64_t>::min ();
00083   // This converts to the nearest double. Unless there's an equality, the
00084   // result is clear.
00085   double xx = x;
00086   if (xx != y)
00087     return xop::op (xx, y);
00088   else
00089     {
00090       // If equality occured we compare as integers.
00091       if (xx == xxup)
00092         return xop::gtval;
00093       else if (xx == xxlo)
00094         return xop::ltval;
00095       else
00096         return xop::op (x, static_cast<int64_t> (xx));
00097     }
00098 
00099 }
00100 
00101 // We define double-int operations by reverting the operator
00102 
00103 // A trait class reverting the operator
00104 template <class xop>
00105 class rev_op
00106 {
00107 public:
00108   typedef xop op;
00109 };
00110 
00111 #define DEFINE_REVERTED_OPERATOR(OP1,OP2) \
00112   template <> \
00113   class rev_op<octave_int_cmp_op::OP1> \
00114   { \
00115   public: \
00116     typedef octave_int_cmp_op::OP2 op; \
00117   }
00118 
00119 DEFINE_REVERTED_OPERATOR(lt,gt);
00120 DEFINE_REVERTED_OPERATOR(gt,lt);
00121 DEFINE_REVERTED_OPERATOR(le,ge);
00122 DEFINE_REVERTED_OPERATOR(ge,le);
00123 
00124 template <class xop>
00125 bool
00126 octave_int_cmp_op::emulate_mop (double x, uint64_t y)
00127 {
00128   typedef typename rev_op<xop>::op rop;
00129   return mop<rop> (y, x);
00130 }
00131 
00132 template <class xop>
00133 bool
00134 octave_int_cmp_op::emulate_mop (double x, int64_t y)
00135 {
00136   typedef typename rev_op<xop>::op rop;
00137   return mop<rop> (y, x);
00138 }
00139 
00140 
00141 // Define handlers for int64 multiplication
00142 
00143 template <>
00144 uint64_t
00145 octave_int_arith_base<uint64_t, false>::mul (uint64_t x, uint64_t y)
00146 {
00147   // Get upper words
00148   uint64_t ux = x >> 32, uy = y >> 32;
00149   uint64_t res;
00150   if (ux)
00151     {
00152       if (uy)
00153         goto overflow;
00154       else
00155         {
00156           uint64_t ly = static_cast<uint32_t> (y), uxly = ux*ly;
00157           if (uxly >> 32)
00158             goto overflow;
00159           uxly <<= 32; // never overflows
00160           uint64_t lx = static_cast<uint32_t> (x), lxly = lx*ly;
00161           res = add (uxly, lxly);
00162         }
00163     }
00164   else if (uy)
00165     {
00166       uint64_t lx = static_cast<uint32_t> (x), uylx = uy*lx;
00167       if (uylx >> 32)
00168         goto overflow;
00169       uylx <<= 32; // never overflows
00170       uint64_t ly = static_cast<uint32_t> (y), lylx = ly*lx;
00171       res = add (uylx, lylx);
00172     }
00173   else
00174     {
00175       uint64_t lx = static_cast<uint32_t> (x);
00176       uint64_t ly = static_cast<uint32_t> (y);
00177       res = lx*ly;
00178     }
00179 
00180   return res;
00181 
00182 overflow:
00183   return max_val ();
00184 }
00185 
00186 template <>
00187 int64_t
00188 octave_int_arith_base<int64_t, true>::mul (int64_t x, int64_t y)
00189 {
00190   // The signed case is far worse. The problem is that
00191   // even if neither integer fits into signed 32-bit range, the result may
00192   // still be OK. Uh oh.
00193 
00194   // Essentially, what we do is compute sign, multiply absolute values
00195   // (as above) and impose the sign.
00196   // FIXME -- can we do something faster if we HAVE_FAST_INT_OPS?
00197 
00198   uint64_t usx = octave_int_abs (x), usy = octave_int_abs (y);
00199   bool positive = (x < 0) == (y < 0);
00200 
00201   // Get upper words
00202   uint64_t ux = usx >> 32, uy = usy >> 32;
00203   uint64_t res;
00204   if (ux)
00205     {
00206       if (uy)
00207         goto overflow;
00208       else
00209         {
00210           uint64_t ly = static_cast<uint32_t> (usy), uxly = ux*ly;
00211           if (uxly >> 32)
00212             goto overflow;
00213           uxly <<= 32; // never overflows
00214           uint64_t lx = static_cast<uint32_t> (usx), lxly = lx*ly;
00215           res = uxly + lxly;
00216           if (res < uxly)
00217             goto overflow;
00218         }
00219     }
00220   else if (uy)
00221     {
00222       uint64_t lx = static_cast<uint32_t> (usx), uylx = uy*lx;
00223       if (uylx >> 32)
00224         goto overflow;
00225       uylx <<= 32; // never overflows
00226       uint64_t ly = static_cast<uint32_t> (usy), lylx = ly*lx;
00227       res = uylx + lylx;
00228       if (res < uylx)
00229         goto overflow;
00230     }
00231   else
00232     {
00233       uint64_t lx = static_cast<uint32_t> (usx);
00234       uint64_t ly = static_cast<uint32_t> (usy);
00235       res = lx*ly;
00236     }
00237 
00238   if (positive)
00239     {
00240       if (res > static_cast<uint64_t> (max_val ()))
00241         {
00242           return max_val ();
00243         }
00244       else
00245         return static_cast<int64_t> (res);
00246     }
00247   else
00248     {
00249       if (res > static_cast<uint64_t> (-min_val ()))
00250         {
00251           return min_val ();
00252         }
00253       else
00254         return -static_cast<int64_t> (res);
00255     }
00256 
00257 
00258 overflow:
00259   return positive ? max_val () : min_val ();
00260 
00261 }
00262 
00263 #define INT_DOUBLE_BINOP_DECL(OP,SUFFIX) \
00264   template <> \
00265   OCTAVE_API octave_ ## SUFFIX \
00266   operator OP (const octave_ ## SUFFIX & x, const double& y)
00267 
00268 #define DOUBLE_INT_BINOP_DECL(OP,SUFFIX) \
00269   template <> \
00270   OCTAVE_API octave_ ## SUFFIX \
00271   operator OP (const double& x, const octave_ ## SUFFIX & y)
00272 
00273 INT_DOUBLE_BINOP_DECL (+, uint64)
00274 {
00275   return (y < 0) ? x - octave_uint64(-y) : x + octave_uint64(y);
00276 }
00277 
00278 DOUBLE_INT_BINOP_DECL (+, uint64)
00279 { return y + x; }
00280 
00281 INT_DOUBLE_BINOP_DECL (+, int64)
00282 {
00283   if (fabs (y) < static_cast<double> (octave_int64::max ()))
00284     return x + octave_int64 (y);
00285   else
00286     {
00287       // If the number is within the int64 range (the most common case,
00288       // probably), the above will work as expected. If not, it's more
00289       // complicated - as long as y is within _twice_ the signed range, the
00290       // result may still be an integer. An instance of such an operation is
00291       // 3*2**62 + (1+intmin('int64')) that should yield int64(2**62) + 1.  So
00292       // what we do is to try to convert y/2 and add it twice. Note that if y/2
00293       // overflows, the result must overflow as well, and that y/2 cannot be a
00294       // fractional number.
00295       octave_int64 y2 (y / 2);
00296       return (x + y2) + y2;
00297     }
00298 }
00299 
00300 DOUBLE_INT_BINOP_DECL (+, int64)
00301 {
00302   return y + x;
00303 }
00304 
00305 INT_DOUBLE_BINOP_DECL (-, uint64)
00306 {
00307   return x + (-y);
00308 }
00309 
00310 DOUBLE_INT_BINOP_DECL (-, uint64)
00311 {
00312   if (x <= static_cast<double> (octave_uint64::max ()))
00313     return octave_uint64(x) - y;
00314   else
00315     {
00316       // Again a trick to get the corner cases right. Things like
00317       // 3**2**63 - intmax('uint64') should produce the correct result, i.e.
00318       // int64(2**63) + 1.
00319       const double p2_64 = std::pow (2.0, 64);
00320       if (y.bool_value ())
00321         {
00322           const uint64_t p2_64my = (~y.value ()) + 1; // Equals 2**64 - y
00323           return octave_uint64 (x - p2_64) + octave_uint64 (p2_64my);
00324         }
00325       else
00326         return octave_uint64 (p2_64);
00327     }
00328 }
00329 
00330 INT_DOUBLE_BINOP_DECL (-, int64)
00331 {
00332   return x + (-y);
00333 }
00334 
00335 DOUBLE_INT_BINOP_DECL (-, int64)
00336 {
00337   static const bool twosc = (std::numeric_limits<int64_t>::min ()
00338                              < -std::numeric_limits<int64_t>::max ());
00339   // In case of symmetric integers (not two's complement), this will probably
00340   // be eliminated at compile time.
00341   if (twosc && y.value () == std::numeric_limits<int64_t>::min ())
00342     {
00343       return octave_int64 (x + std::pow(2.0, 63));
00344     }
00345   else
00346     return x + (-y);
00347 }
00348 
00349 // NOTE:
00350 // Emulated mixed multiplications are tricky due to possible precision loss.
00351 // Here, after sorting out common cases for speed, we follow the strategy
00352 // of converting the double number into the form sign * 64-bit integer* 2**exponent,
00353 // multiply the 64-bit integers to get a 128-bit number, split that number into 32-bit words
00354 // and form 4 double-valued summands (none of which loases precision), then convert these
00355 // into integers and sum them. Though it is not immediately obvious, this should work
00356 // even w.r.t. rounding (none of the summands lose precision).
00357 
00358 // Multiplies two unsigned 64-bit ints to get a 128-bit number represented
00359 // as four 32-bit words.
00360 static void
00361 umul128 (uint64_t x, uint64_t y, uint32_t w[4])
00362 {
00363   uint64_t lx = static_cast<uint32_t> (x), ux = x >> 32;
00364   uint64_t ly = static_cast<uint32_t> (y), uy = y >> 32;
00365   uint64_t a = lx * ly;
00366   w[0] = a; a >>= 32;
00367   uint64_t uxly = ux*ly, uylx = uy*lx;
00368   a += static_cast<uint32_t> (uxly); uxly >>= 32;
00369   a += static_cast<uint32_t> (uylx); uylx >>= 32;
00370   w[1] = a; a >>= 32;
00371   uint64_t uxuy = ux * uy;
00372   a += uxly; a += uylx; a += uxuy;
00373   w[2] = a; a >>= 32;
00374   w[3] = a;
00375 }
00376 
00377 // Splits a double into bool sign, unsigned 64-bit mantissa and int exponent
00378 static void
00379 dblesplit (double x, bool& sign, uint64_t& mtis, int& exp)
00380 {
00381   sign = x < 0; x = fabs (x);
00382   x = frexp (x, &exp);
00383   exp -= 52;
00384   mtis = static_cast<uint64_t> (ldexp (x, 52));
00385 }
00386 
00387 // Gets a double number from a 32-bit unsigned integer mantissa, exponent and sign.
00388 static double
00389 dbleget (bool sign, uint32_t mtis, int exp)
00390 {
00391   double x = ldexp (static_cast<double> (mtis), exp);
00392   return sign ? -x : x;
00393 }
00394 
00395 
00396 INT_DOUBLE_BINOP_DECL (*, uint64)
00397 {
00398   if (y >= 0 && y < octave_uint64::max () && y == xround (y))
00399     {
00400       return x * octave_uint64 (static_cast<uint64_t> (y));
00401     }
00402   else if (y == 0.5)
00403     {
00404       return x / octave_uint64 (static_cast<uint64_t> (2));
00405     }
00406   else if (y < 0 || xisnan (y) || xisinf (y))
00407     {
00408       return octave_uint64 (x.value () * y);
00409     }
00410   else
00411     {
00412       bool sign;
00413       uint64_t my;
00414       int e;
00415       dblesplit (y, sign, my, e);
00416       uint32_t w[4];
00417       umul128 (x.value (), my, w);
00418       octave_uint64 res = octave_uint64::zero;
00419       for (short i = 0; i < 4; i++)
00420         {
00421           res += octave_uint64 (dbleget (sign, w[i], e));
00422           e += 32;
00423         }
00424       return res;
00425     }
00426 }
00427 
00428 DOUBLE_INT_BINOP_DECL (*, uint64)
00429 { return y * x; }
00430 
00431 INT_DOUBLE_BINOP_DECL (*, int64)
00432 {
00433   if (fabs (y) < octave_int64::max () && y == xround (y))
00434     {
00435       return x * octave_int64 (static_cast<int64_t> (y));
00436     }
00437   else if (fabs (y) == 0.5)
00438     {
00439       return x / octave_int64 (static_cast<uint64_t> (4*y));
00440     }
00441   else if (xisnan (y) || xisinf (y))
00442     {
00443       return octave_int64 (x.value () * y);
00444     }
00445   else
00446     {
00447       bool sign;
00448       uint64_t my;
00449       int e;
00450       dblesplit (y, sign, my, e);
00451       uint32_t w[4];
00452       sign = (sign != (x.value () < 0));
00453       umul128 (octave_int_abs (x.value ()), my, w);
00454       octave_int64 res = octave_int64::zero;
00455       for (short i = 0; i < 4; i++)
00456         {
00457           res += octave_int64 (dbleget (sign, w[i], e));
00458           e += 32;
00459         }
00460       return res;
00461     }
00462 }
00463 
00464 DOUBLE_INT_BINOP_DECL (*, int64)
00465 { return y * x; }
00466 
00467 DOUBLE_INT_BINOP_DECL (/, uint64)
00468 {
00469   return octave_uint64 (x / static_cast<double> (y));
00470 }
00471 
00472 DOUBLE_INT_BINOP_DECL (/, int64)
00473 {
00474   return octave_int64 (x / static_cast<double> (y));
00475 }
00476 
00477 INT_DOUBLE_BINOP_DECL (/, uint64)
00478 {
00479   if (y >= 0 && y < octave_uint64::max () && y == xround (y))
00480     {
00481       return x / octave_uint64 (y);
00482     }
00483   else
00484     return x * (1.0/y);
00485 }
00486 
00487 INT_DOUBLE_BINOP_DECL (/, int64)
00488 {
00489   if (fabs (y) < octave_int64::max () && y == xround (y))
00490     {
00491       return x / octave_int64 (y);
00492     }
00493   else
00494     return x * (1.0/y);
00495 }
00496 
00497 #define INSTANTIATE_INT64_DOUBLE_CMP_OP0(OP,T1,T2) \
00498   template OCTAVE_API bool \
00499   octave_int_cmp_op::emulate_mop<octave_int_cmp_op::OP> (T1 x, T2 y)
00500 
00501 #define INSTANTIATE_INT64_DOUBLE_CMP_OP(OP) \
00502   INSTANTIATE_INT64_DOUBLE_CMP_OP0(OP, double, int64_t); \
00503   INSTANTIATE_INT64_DOUBLE_CMP_OP0(OP, double, uint64_t); \
00504   INSTANTIATE_INT64_DOUBLE_CMP_OP0(OP, int64_t, double); \
00505   INSTANTIATE_INT64_DOUBLE_CMP_OP0(OP, uint64_t, double)
00506 
00507 INSTANTIATE_INT64_DOUBLE_CMP_OP(lt);
00508 INSTANTIATE_INT64_DOUBLE_CMP_OP(le);
00509 INSTANTIATE_INT64_DOUBLE_CMP_OP(gt);
00510 INSTANTIATE_INT64_DOUBLE_CMP_OP(ge);
00511 INSTANTIATE_INT64_DOUBLE_CMP_OP(eq);
00512 INSTANTIATE_INT64_DOUBLE_CMP_OP(ne);
00513 
00514 #endif
00515 
00516 //template <class T>
00517 //bool
00518 //xisnan (const octave_int<T>&)
00519 //{
00520 //  return false;
00521 //}
00522 
00523 template <class T>
00524 octave_int<T>
00525 pow (const octave_int<T>& a, const octave_int<T>& b)
00526 {
00527   octave_int<T> retval;
00528 
00529   octave_int<T> zero = static_cast<T> (0);
00530   octave_int<T> one = static_cast<T> (1);
00531 
00532   if (b == zero || a == one)
00533     retval = one;
00534   else if (b < zero)
00535     {
00536       if (a == -one)
00537         retval = (b.value () % 2) ? a : one;
00538       else
00539         retval = zero;
00540     }
00541   else
00542     {
00543       octave_int<T> a_val = a;
00544       T b_val = b; // no need to do saturation on b
00545 
00546       retval = a;
00547 
00548       b_val -= 1;
00549 
00550       while (b_val != 0)
00551         {
00552           if (b_val & 1)
00553             retval = retval * a_val;
00554 
00555           b_val = b_val >> 1;
00556 
00557           if (b_val)
00558             a_val = a_val * a_val;
00559         }
00560     }
00561 
00562   return retval;
00563 }
00564 
00565 template <class T>
00566 octave_int<T>
00567 pow (const double& a, const octave_int<T>& b)
00568 { return octave_int<T> (pow (a, b.double_value ())); }
00569 
00570 template <class T>
00571 octave_int<T>
00572 pow (const octave_int<T>& a, const double& b)
00573 {
00574   return ((b >= 0 && b < std::numeric_limits<T>::digits && b == xround (b))
00575           ? pow (a, octave_int<T> (static_cast<T> (b)))
00576           : octave_int<T> (pow (a.double_value (), b)));
00577 }
00578 
00579 template <class T>
00580 octave_int<T>
00581 pow (const float& a, const octave_int<T>& b)
00582 { return octave_int<T> (pow (a, b.float_value ())); }
00583 
00584 template <class T>
00585 octave_int<T>
00586 pow (const octave_int<T>& a, const float& b)
00587 {
00588   return ((b >= 0 && b < std::numeric_limits<T>::digits && b == xround (b))
00589           ? pow (a, octave_int<T> (static_cast<T> (b)))
00590           : octave_int<T> (pow (a.double_value (), static_cast<double> (b))));
00591 }
00592 
00593 // FIXME: Do we really need a differently named single-precision
00594 //        function integer power function here instead of an overloaded
00595 //        one?
00596 template <class T>
00597 octave_int<T>
00598 powf (const float& a, const octave_int<T>& b)
00599 { return octave_int<T> (pow (a, b.float_value ())); }
00600 
00601 template <class T>
00602 octave_int<T>
00603 powf (const octave_int<T>& a, const float& b)
00604 {
00605   return ((b >= 0 && b < std::numeric_limits<T>::digits && b == xround (b))
00606           ? pow (a, octave_int<T> (static_cast<T> (b)))
00607           : octave_int<T> (pow (a.double_value (), static_cast<double> (b))));
00608 }
00609 
00610 #define INSTANTIATE_INTTYPE(T) \
00611   template class OCTAVE_API octave_int<T>; \
00612   template OCTAVE_API octave_int<T> pow (const octave_int<T>&, const octave_int<T>&); \
00613   template OCTAVE_API octave_int<T> pow (const double&, const octave_int<T>&); \
00614   template OCTAVE_API octave_int<T> pow (const octave_int<T>&, const double&); \
00615   template OCTAVE_API octave_int<T> pow (const float&, const octave_int<T>&);  \
00616   template OCTAVE_API octave_int<T> pow (const octave_int<T>&, const float&);  \
00617   template OCTAVE_API octave_int<T> powf (const float&, const octave_int<T>&); \
00618   template OCTAVE_API octave_int<T> powf (const octave_int<T>&, const float&); \
00619   template OCTAVE_API octave_int<T> \
00620   bitshift (const octave_int<T>&, int, const octave_int<T>&); \
00621 
00622 INSTANTIATE_INTTYPE (int8_t);
00623 INSTANTIATE_INTTYPE (int16_t);
00624 INSTANTIATE_INTTYPE (int32_t);
00625 INSTANTIATE_INTTYPE (int64_t);
00626 
00627 INSTANTIATE_INTTYPE (uint8_t);
00628 INSTANTIATE_INTTYPE (uint16_t);
00629 INSTANTIATE_INTTYPE (uint32_t);
00630 INSTANTIATE_INTTYPE (uint64_t);
00631 
00632 
00633 // Tests follow.
00634 
00635 /*
00636 
00637 %!assert(intmax("int64")/intmin("int64"),int64(-1))
00638 %!assert(intmin("int64")/int64(-1),intmax("int64"))
00639 %!assert(int64(2**63),intmax("int64"))
00640 %!assert(uint64(2**64),intmax("uint64"))
00641 %!test
00642 %! a = 1.9*2^61; b = uint64(a); b++; assert(b > a)
00643 %!test
00644 %! a = -1.9*2^61; b = int64(a); b++; assert(b > a)
00645 %!test
00646 %! a = int64(-2**60) + 2; assert(1.25*a == (5*a)/4)
00647 %!test
00648 %! a = uint64(2**61) + 2; assert(1.25*a == (5*a)/4)
00649 %!assert(int32(2**31+0.5),intmax('int32'))
00650 %!assert(int32(-2**31-0.5),intmin('int32'))
00651 %!assert((int64(2**62)+1)**1, int64(2**62)+1)
00652 %!assert((int64(2**30)+1)**2, int64(2**60+2**31) + 1)
00653 */
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines