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