src/real.h

   1 #ifndef _REAL_H
   2 #define _REAL_H
   3
   4 #include "common.h"
   5 #include <cmath>
   6 #include <cfloat>
   7
   8
   9 #define REAL_SINGLE 0
  10 #define REAL_DOUBLE 1
  11 #define REAL_LONG_DOUBLE 2
  12 #define REAL_VFPU 3
  13 #define REAL_RATIONAL 4
  14 #define REAL_RATIONAL_ARBINT 5
  15 #define REAL_MPFRCPP 6
  16 #define REAL_IRRAM 7
  17 #define REAL_PARANOIDNUMBER 8
  18 #define REAL_GMPRAT 9
  19
  20 #ifndef REALTYPE
  21         #error "REALTYPE was not defined!"
  22 #endif
  23
  24 #define XSTR(x) STR(x)
  25 #define STR(x) #x
  26 //#pragma message "REALTYPE = " XSTR(REALTYPE)
  27
  28 #if REALTYPE == REAL_VFPU
  29         #include "vfpu.h"
  30 #endif
  31
  32 #if REALTYPE == REAL_RATIONAL
  33         #include "rational.h"
  34 #endif //REALTYPE
  35
  36 #if REALTYPE == REAL_RATIONAL_ARBINT
  37         #include "rational.h"
  38         #include "arbint.h"
  39         #include "gmpint.h"
  40 #endif //REALTYPE
  41
  42 #if REALTYPE == REAL_MPFRCPP
  43         #include <mpreal.h>
  44 #endif //REALTYPE
  45
  46 #if REALTYPE == REAL_IRRAM
  47         #include "../contrib/iRRAM/include/iRRAM/lib.h"
  48 #endif
  49
  50 #if REALTYPE == REAL_PARANOIDNUMBER
  51         #include "paranoidnumber.h"
  52 #endif
  53
  54 #if REALTYPE == REAL_GMPRAT
  55         #include "gmprat.h"
  56 #endif
  57
  58 namespace IPDF
  59 {
  60         extern const char * g_real_name[];
  61
  62 #if REALTYPE == REAL_SINGLE
  63         typedef float Real;
  64         inline Real RealFromStr(const char * str) {return strtof(str, NULL);}
  65         inline std::string Str(const Real & a) {std::stringstream s; s << a; return s.str();}
  66 #elif REALTYPE == REAL_DOUBLE
  67         typedef double Real;
  68         inline Real RealFromStr(const char * str) {return strtod(str, NULL);}
  69         inline std::string Str(const Real & a) {std::stringstream s; s << a; return s.str();}
  70 #elif REALTYPE == REAL_LONG_DOUBLE
  71         typedef long double Real;
  72         inline Real RealFromStr(const char * str) {return strtold(str, NULL);}
  73         inline std::string Str(const Real & a) {std::stringstream s; s << a; return s.str();}
  74 #elif REALTYPE == REAL_VFPU
  75         typedef VFPU::VFloat Real;
  76         inline float Float(const Real & r) {return r.m_value;}
  77         inline double Double(const Real & r) {return r.m_value;}
  78         inline Real RealFromStr(const char * str) {return Real(strtod(str, NULL));}
  79 #elif REALTYPE == REAL_RATIONAL
  80         typedef Rational<int64_t> Real;
  81         inline float Float(const Real & r) {return (float)r.ToDouble();}
  82         inline double Double(const Real & r) {return r.ToDouble();}
  83         inline Real RealFromStr(const char * str) {return Real(strtod(str, NULL));}
  84 #elif REALTYPE == REAL_RATIONAL_ARBINT
  85         #define ARBINT Arbint // Set to Gmpint or Arbint here
  86
  87         typedef Rational<ARBINT> Real;
  88         inline float Float(const Real & r) {return (float)r.ToDouble();}
  89         inline double Double(const Real & r) {return r.ToDouble();}
  90         inline int64_t Int64(const Real & r) {return r.ToInt64();}
  91         inline Rational<ARBINT> Sqrt(const Rational<ARBINT> & r) {return r.Sqrt();}
  92         inline Real RealFromStr(const char * str) {return Real(strtod(str, NULL));}
  93 #elif REALTYPE == REAL_MPFRCPP
  94         typedef mpfr::mpreal Real;
  95         inline double Double(const Real & r) {return r.toDouble();}
  96         inline float Float(const Real & r) {return r.toDouble();}
  97         inline int64_t Int64(const Real & r) {return r.toLong();}
  98         inline Real Sqrt(const Real & r) {return mpfr::sqrt(r, mpfr::mpreal::get_default_rnd());}
  99         inline Real Abs(const Real & r) {return mpfr::abs(r, mpfr::mpreal::get_default_rnd());}
 100         inline Real RealFromStr(const char * str) {return Real(strtod(str, NULL));}
 101 #elif REALTYPE == REAL_IRRAM
 102         typedef iRRAM::REAL Real;
 103         inline double Double(const Real & r) {return r.as_double(53);}
 104         inline float Float(const Real & r) {return r.as_double(53);}
 105         inline int64_t Int64(const Real & r) {return (int64_t)r.as_double(53);}
 106         inline Real Sqrt(const Real & r) {return iRRAM::sqrt(r);}
 107         inline Real RealFromStr(const char * str) {return Real(strtod(str, NULL));}
 108 #elif REALTYPE == REAL_PARANOIDNUMBER
 109         typedef ParanoidNumber Real;
 110         inline double Double(const Real & r) {return r.Digit();}
 111         inline float Float(const Real & r) {return r.Digit();}
 112         inline int64_t Int64(const Real & r) {return (int64_t)r.Digit();}
 113         inline Real Sqrt(const Real & r) {return Real(sqrt(r.Digit()));}
 114         inline Real RealFromStr(const char * str) {return Real(str);}
 115         inline Real Abs(const Real & a) {return Real(fabs(a.Digit()));}
 116 #elif REALTYPE == REAL_GMPRAT
 117         typedef Gmprat Real;
 118         inline double Double(const Real & r) {return r.ToDouble();}
 119         inline float Float(const Real & r) {return (float)(r.ToDouble());}
 120         inline int64_t Int64(const Real & r) {return (int64_t)r.ToDouble();}
 121         inline Real Sqrt(const Real & r) {return Real(sqrt(r.ToDouble()));}
 122         inline Real RealFromStr(const char * str) {return Real(strtod(str, NULL));}
 123         inline Real Abs(const Real & a) {return (a > Real(0)) ? a : Real(0)-a;}
 124
 125
 126 #else
 127         #error "Type of Real unspecified."
 128 #endif //REALTYPE
 129
 130         // Allow us to call Float on the primative types
 131         // Useful so I can template some things that could be either (a more complicated) Real or a primitive type
 132         // Mostly in the testers.
 133         inline float Float(double f) {return (float)f;}
 134         inline float Float(long double f) {return (float)(f);}
 135         inline double Double(float f) {return (double)f;}
 136         inline double Double(double f) {return (double)f;}
 137         inline double Double(long double f) {return (double)(f);}
 138         inline double Sqrt(double f) {return sqrt(f);}
 139         inline double Abs(double a) {return fabs(a);}
 140         inline double Log10(double a) {return log(a)/log(10.0);}
 141         inline size_t Size(double a) {return sizeof(a);}
 142         inline size_t Size(float a) {return sizeof(a);}
 143
 144         // Don't cause an exception
 145         inline float ClampFloat(double d)
 146         {
 147                 float f = (fabs(d) < FLT_MAX) ? ((fabs(d) > FLT_MIN) ? (float)d : FLT_MIN) : FLT_MAX;
 148                 return copysign(f, d);
 149         }
 150
 151         inline int64_t Int64(double a)
 152         {
 153                 if (a < INT64_MIN)
 154                         return INT64_MIN;
 155                 if (a > INT64_MAX)
 156                         return INT64_MAX;
 157                 return (int64_t)(a);
 158         }
 159
 160         inline Real Power(const Real & a, int n)
 161         {
 162                 if (n < 0)
 163                 {
 164                         return Power(Real(1)/a, -n);
 165                 }
 166                 Real r(1);
 167                 for (int i = 0; i < n; ++i)
 168                         r *= a;
 169                 return r;
 170         }
 171         struct Vec2
 172         {
 173                 Real x;
 174                 Real y;
 175                 Vec2() : x(0.0), y(0.0) {}
 176                 Vec2(Real _x, Real _y) : x(_x), y(_y) {}
 177                 Vec2(const std::pair<Real, Real> & p) : x(p.first), y(p.second) {}
 178                 #if REALTYPE != REAL_IRRAM
 179                 Vec2(const std::pair<int64_t, int64_t> & p) : x(p.first), y(p.second) {}
 180                 #else
 181                 Vec2(const std::pair<int64_t, int64_t> & p) : x((int)p.first), y((int)p.second) {}
 182                 // Apparently iRRAM didn't implement -= for a constant argument
 183                 #endif
 184                 bool operator==(const Vec2& other) const { return (x == other.x) && (y == other.y); }
 185                 bool operator!=(const Vec2& other) const { return !(*this == other); }
 186
 187                 Vec2& operator=(const Vec2& other) { x = other.x; y = other.y; return *this; }
 188                 Vec2& operator+=(const Vec2& other) { x += other.x; y += other.y; return *this; }
 189                 Vec2& operator-=(const Vec2& other) { x -= other.x; y -= other.y; return *this; }
 190                 Vec2& operator*=(const Real& lambda) { x *= lambda; y *= lambda; return *this; }
 191                 Vec2& operator/=(const Real& lambda) { x /= lambda; y /= lambda; return *this; }
 192
 193                 Vec2 operator+(const Vec2& other) const { return Vec2(x + other.x, y + other.y); }
 194                 Vec2 operator-(const Vec2& other) const { return Vec2(x - other.x, y - other.y); }
 195                 Vec2 operator*(const Real& lambda) const { return Vec2(x * lambda, y * lambda); }
 196                 Vec2 operator/(const Real& lambda) const { return Vec2(x / lambda, y / lambda); }
 197
 198                 const Real SquareLength() const { return (x*x + y*y); }
 199
 200         };
 201
 202         //TODO: Make sure there is actually a RealFromStr(const char * str) function
 203         //              Or this will recurse infinitely
 204         //              (If you remove this it will also break).
 205         inline Real RealFromStr(const std::string & str) {return RealFromStr(str.c_str());}
 206
 207
 208         // things stolen from wikipedia and googling
 209         inline const char * HumanScale(double f)
 210         {
 211                 if (f < 1e-36)
 212                         return "RATHER SMALL";
 213                 if (f < 1e-35)
 214                         return "Plank Length";
 215                 if (f < 1e-25)
 216                         return "Turtles all the way";
 217                 if (f < 1e-24)
 218                         return "More turtles";
 219                 if (f < 1e-23)
 220                         return "Turtles";
 221                 if (f < 1e-22)
 222                         return "This small";
 223                 if (f < 1e-21)
 224                         return "To find things";
 225                 if (f < 1e-20)
 226                         return "It is pretty difficult";
 227                 if (f < 1e-19)
 228                         return "Not much";
 229                 if (f < 1e-17)
 230                         return "Weak Force";
 231                 if (f < 1e-16)
 232                         return "Proton";
 233                 if (f < 1e-15)
 234                         return "(Classical) Electron";
 235                 if (f < 1e-11)
 236                         return "Inter atomic (still)";
 237                 if (f < 1e-10)
 238                         return "Inter atomic";
 239                 if (f < 1e-9)
 240                         return "Atom";
 241                 if (f < 1e-8)
 242                         return "DNA";
 243                 if (f < 1e-7)
 244                         return "Virus";
 245                 if (f < 1e-6)
 246                         return "Light";
 247                 if (f < 1e-5)
 248                         return "Bacteria";
 249                 if (f < 1e-4)
 250                         return "4004 Transistor";
 251                 if (f < 1e-3)
 252                         return "Ant";
 253                 if (f < 1e-2)
 254                         return "Coin";
 255                 if (f < 1e-1)
 256                         return "iPhone";
 257                 if (f < 1e0)
 258                         return "Person";
 259                 if (f < 1e1)
 260                         return "Building";
 261                 if (f < 1e2)
 262                         return "Football Field";
 263                 if (f < 1e3)
 264                         return "Mountain";
 265                 if (f < 1e4)
 266                         return "Clouds";
 267                 if (f < 1e5)
 268                         return "Countries";
 269                 if (f < 1e6)
 270                         return "Earth";
 271                 if (f < 1e8)
 272                         return "Between Earth and Moon";
 273                 if (f < 1e9)
 274                         return "Solar System";
 275                 if (f < 1e13)
 276                         return "Distance to nearest Star";
 277                 if (f < 1e21)
 278                         return "Milky Way";
 279                 if (f < 1e26)
 280                         return "Universe";
 281                 if (f < 1e27)
 282                         return "A bigger Universe";
 283                 if (f < 1e28)
 284                         return "Really big things";
 285                 if (f < 1e29)
 286                         return "Almost as big as...";
 287                 if (f < 1e30)
 288                         return "Wolfram's Magestic Ego";
 289                 return "QUITE BIG";
 290
 291         }
 292
 293         inline void DebugRealInfo()
 294         {
 295                 Debug("Compiled with REAL = %d => \"%s\" sizeof(Real) == %d bytes", REALTYPE, g_real_name[REALTYPE], sizeof(Real));
 296                 #if REALTYPE == REAL_PARANOIDNUMBER
 297                         #ifdef PARANOID_SIZE_LIMIT
 298                                 Debug("Size limit of %d is being enforced", PARANOID_SIZE_LIMIT);
 299                         #endif
 300                 #endif
 301         }
 302 }
 303
 304 #endif //_REAL_H