https://github.com/ITensor/ITensor
Tip revision: 6928c721813cde1204e8bf29dfb51a2f44dd9c9e authored by Miles Stoudenmire on 13 November 2015, 17:24:13 UTC
Removed outdated IQTensor constructor declarations
Removed outdated IQTensor constructor declarations
Tip revision: 6928c72
real.h
//
// Distributed under the ITensor Library License, Version 1.2
// (See accompanying LICENSE file.)
//
#ifndef __ITENSOR_REAL_H
#define __ITENSOR_REAL_H
#include <cmath>
#include "math.h"
#include "itensor/types.h"
#include "itensor/util/print.h"
namespace itensor {
static const Real Pi = M_PI;
static const Real Sqrt2 = std::sqrt(2.);
static const Real ISqrt2 = 1./std::sqrt(2.);
template <typename T>
T
sqr(T x) { return x*x; }
static const Real maxlogdouble = log(std::numeric_limits<double>::max());
static const Real LogNum_Accuracy = 1E-12;
class TooBigForReal : public ITError
{
public:
using Parent = ITError;
TooBigForReal(const std::string& message)
: Parent(message)
{ }
};
class TooSmallForReal : public ITError
{
public:
using Parent = ITError;
TooSmallForReal(const std::string& message)
: Parent(message)
{ }
};
//
// LogNum
//
//Stores a real number r as lognum_ = log(|r|) and sign_ = sgn(r)
class LogNum
{
//Log of the magnitude of
//the number represented.
Real lognum_;
int sign_;
public:
Real
logNum() const { return lognum_; }
int
sign() const { return sign_; }
bool
isZero() const { return sign_ == 0; }
bool
isRealZero() const
{ return (sign_ == 0 || lognum_ < -maxlogdouble); }
bool
isFiniteReal() const
{ return (lognum_ < maxlogdouble && lognum_ > -maxlogdouble); }
bool
isTooBigForReal() const
{ return (lognum_ > maxlogdouble); }
bool
isTooSmallForReal() const
{ return (lognum_ < -maxlogdouble); }
bool friend inline
isnan(const LogNum& L) { return std::isnan(L.lognum_); }
//Default constructed to NAN
//to catch initialization errors
LogNum()
:
lognum_(NAN),
sign_(1)
{ }
explicit
LogNum(Real r)
{
if(r == 0)
{
sign_ = 0;
lognum_ = 0;
}
else
if(r < 0)
{
sign_ = -1;
lognum_ = log(-r);
}
else
{
sign_ = 1;
lognum_ = log(r);
}
}
LogNum(Real lognum, int sign)
:
lognum_(lognum),
sign_(sign)
{
#ifdef DEBUG
if(!(sign == -1 || sign == 0 || sign == 1))
Error("sign should be -1, 0, or 1");
#endif
}
void
read(std::istream& s)
{
s.read((char*) &lognum_,sizeof(lognum_));
s.read((char*) &sign_,sizeof(sign_));
}
void
write(std::ostream& s) const
{
s.write((char*) &lognum_,sizeof(lognum_));
s.write((char*) &sign_,sizeof(sign_));
}
//operator Real() const too easy to misuse accidentally
Real
real() const
{
if(sign_ == 0) return 0;
#ifdef DEBUG
if(lognum_ > maxlogdouble)
{
println("lognum_ = ",lognum_);
throw TooBigForReal("LogNum too big to convert to Real");
}
if(lognum_ < -maxlogdouble)
{
println("lognum_ = ",lognum_);
throw TooSmallForReal("LogNum too small to convert to Real");
}
#endif
return sign_ * exp(lognum_);
}
Real
real0() const
{
if(isRealZero()) return 0.0;
return real();
}
LogNum&
operator+=(const LogNum& other)
{
try {
*this =
LogNum(this->real()+other.real());
}
catch(const ITError& e)
{
Error("Could not convert to real in LogNum::operator+=");
}
return *this;
}
bool
operator==(const LogNum& other) const
{ return (sign_ == other.sign_) && (lognum_ == other.lognum_); }
bool
operator!=(const LogNum& other) const
{ return (sign_ != other.sign_) || (lognum_ != other.lognum_); }
bool
approxEquals(const LogNum& other) const
{ return (sign_ == other.sign_) && (std::fabs(lognum_-other.lognum_) < LogNum_Accuracy); }
void
negate() { sign_ *= -1; }
LogNum&
operator*=(const LogNum& other)
{
sign_ *= other.sign_;
lognum_ += other.lognum_;
return *this;
}
LogNum&
operator*=(Real other)
{
if(other == -1) // special case for efficiency
{
sign_ *= -1;
return *this;
}
return *this *= LogNum(other);
}
LogNum&
operator/=(const LogNum& other)
{
#ifdef DEBUG
if(other.sign_ == 0) Error("divide by zero in LogNum");
#endif
sign_ *= other.sign_;
lognum_ -= other.lognum_;
return *this;
}
LogNum&
operator/=(Real other) { return *this /= LogNum(other); }
LogNum
operator-() const
{ LogNum res(*this); res.sign_ *= -1; return res; }
LogNum
operator/(const LogNum& other) const
{ LogNum res(*this); res /= other; return res; }
LogNum
operator*(const LogNum& other) const
{ LogNum res(*this); res *= other; return res; }
bool
operator<(const LogNum& other) const
{
if(sign_ != other.sign_)
{ return sign_ < other.sign_; }
else if(sign_ > 0) //this and other are positive
{ return lognum_ < other.lognum_; }
else if(sign_ == 0) //this and other are zero
{ return false; }
//this and other are negative
return lognum_ > other.lognum_;
}
bool
operator<=(const LogNum& other) const
{
return (operator<(other) || operator==(other));
}
bool
operator>=(const LogNum& other) const
{
return !(operator<(other));
}
bool
operator>(const LogNum& other) const
{
return (!operator<(other)) && (!operator==(other));
}
void
swap(LogNum& other)
{
Real sl = lognum_;
lognum_ = other.lognum_;
other.lognum_ = sl;
int si = sign_;
sign_ = other.sign_;
other.sign_ = si;
}
bool
magnitudeLessThan(const LogNum& other) const
{
if(sign_ == 0) return other.sign_ != 0;
if(other.sign_ == 0) return false;
return lognum_ < other.lognum_;
}
const LogNum&
pow(Real p)
{
lognum_ *= p;
return *this;
}
};
//For backwards compatibility:
using LogNumber = LogNum;
LogNum inline
operator*(LogNum L, Real r) { L *= r; return L; }
LogNum inline
operator*(Real r, LogNum L) { L *= r; return L; }
LogNum inline
sqrt(LogNum L)
{
if(L.sign() < 0)
Error("Negative LogNum in sqrt");
return L.pow(0.5);
}
inline std::ostream&
operator<<(std::ostream& s, LogNum const& N)
{
s << "LogNum(" << N.logNum() << ",";
if(N.sign() == 0) s << "0)";
else s << (N.sign() > 0 ? "+)" : "-)");
return s;
}
} //namespace itensor
#endif