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