agrad_thread_safe.hpp

Go to the documentation of this file.

#ifndef __STAN__AGRAD__AGRAD_HPP__
#define __STAN__AGRAD__AGRAD_HPP__

#include <stdio.h>
#include <stdlib.h>
#include <vector>
#include <cmath>
#include <cstddef>
#include "stan/memory/stack_alloc.hpp"

// FIXME:  Should include common defs, not this huge cut-and-paste!

namespace stan {

  namespace agrad {

    class vari;

    namespace {
      struct var_allocator {
        std::vector<vari*> var_stack_; 
        memory::stack_alloc memalloc_;
        inline void* alloc(size_t nbytes) {
          return memalloc_.alloc(nbytes);
        }
        inline void recover() {
          var_stack_.resize(0);
          memalloc_.recover_all();
        }
        inline void free() {
          memalloc_.free_all();
        }
      };
#ifdef AGRAD_THREAD_SAFE      
      __thread
#endif
      var_allocator* allocator_;
    }

    class vari {
    private:
      friend class var;

    public: 

      const double val_;

      double adj_;
      
      vari(const double x): 
        val_(x),
        adj_(0.0) {
        allocator_->var_stack_.push_back(this);
      }

      virtual void chain() {
      }

      static inline void* operator new(size_t nbytes) {
        if (allocator_ == 0)
          allocator_ = new var_allocator();
        return allocator_->alloc(nbytes);
      }

      static void recover_memory() {
        return allocator_->recover();
        // allocator_.var_stack_.resize(0);
        // allocator_.memalloc_.recover_all();
      }

      static void free_memory() {
        allocator_->free();
        // allocator_.memalloc_.free_all();
      }

    private:
      static void grad(vari* vi) {
        std::vector<vari*>::iterator it = allocator_->var_stack_.end();
        std::vector<vari*>::iterator begin = allocator_->var_stack_.begin();
        // skip to root variable
        for (; (it >= begin) && (*it != vi); --it)
          ;
        vi->adj_ = 1.0; // droot/droot = 1
        // propagate derivates for remaining vars
        for (; it >= begin; --it)
          (*it)->chain();
      }

    };

    namespace {

      class op_v_vari : public vari {
      protected:
        vari* avi_;
      public:
        op_v_vari(double f, vari* avi) :
          vari(f),
          avi_(avi) {
        }
      };

      class op_vv_vari : public vari {
      protected:
        vari* avi_;
        vari* bvi_;
      public:
        op_vv_vari(double f, vari* avi, vari* bvi):
          vari(f),
          avi_(avi),
          bvi_(bvi) {
        }
      };

      class op_vd_vari : public vari {
      protected:
        vari* avi_;
        double bd_;
      public:
        op_vd_vari(double f, vari* avi, double b) :
          vari(f),
          avi_(avi),
          bd_(b) {
        }
      };

      class op_dv_vari : public vari {
      protected:
        double ad_;
        vari* bvi_;
      public:
        op_dv_vari(double f, double a, vari* bvi) :
          vari(f),
          ad_(a),
          bvi_(bvi) {
        }
      };

      class op_vvv_vari : public vari {
      protected:
        vari* avi_;
        vari* bvi_;
        vari* cvi_;
      public:
        op_vvv_vari(double f, vari* avi, vari* bvi, vari* cvi) :
          vari(f),
          avi_(avi),
          bvi_(bvi),
          cvi_(cvi) {
        }
      };

      class op_vvd_vari : public vari {
      protected:
        vari* avi_;
        vari* bvi_;
        double cd_;
      public:
        op_vvd_vari(double f, vari* avi, vari* bvi, double c) :
          vari(f),
          avi_(avi),
          bvi_(bvi),
          cd_(c) {
        }
      };

      class op_vdv_vari : public vari {
      protected:
        vari* avi_;
        double bd_;
        vari* cvi_;
      public:
        op_vdv_vari(double f, vari* avi, double b, vari* cvi) :
          vari(f),
          avi_(avi),
          bd_(b), 
          cvi_(cvi) {
        }
      };

      class op_vdd_vari : public vari {
      protected:
        vari* avi_;
        double bd_;
        double cd_;
      public:
        op_vdd_vari(double f, vari* avi, double b, double c) :
          vari(f),
          avi_(avi),
          bd_(b), 
          cd_(c) {
        }
      };

      class op_dvv_vari : public vari {
      protected:
        double ad_;
        vari* bvi_;
        vari* cvi_;
      public:
        op_dvv_vari(double f, double a, vari* bvi, vari* cvi) :
          vari(f),
          ad_(a),
          bvi_(bvi),
          cvi_(cvi) {
        }
      };

      class op_dvd_vari : public vari {
      protected:
        double ad_;
        vari* bvi_;
        double cd_;
      public:
        op_dvd_vari(double f, double a, vari* bvi, double c) :
          vari(f),
          ad_(a),
          bvi_(bvi),
          cd_(c) {
        }
      };

      class op_ddv_vari : public vari {
      protected:
        double ad_;
        double bd_;
        vari* cvi_;
      public:
        op_ddv_vari(double f, double a, double b, vari* cvi) :
          vari(f),
          ad_(a),
          bd_(b),
          cvi_(cvi) {
        }
      };

      class neg_vari : public op_v_vari {
      public: 
        neg_vari(vari* avi) :
          op_v_vari(-(avi->val_), avi) {
        }
        void chain() {
          avi_->adj_ -= adj_;
        }
      };


      class add_vv_vari : public op_vv_vari {
      public:
        add_vv_vari(vari* avi, vari* bvi) :
          op_vv_vari(avi->val_ + bvi->val_, avi, bvi) {
        }
        void chain() {
          avi_->adj_ += adj_;
          bvi_->adj_ += adj_;
        }
      };

      class add_vd_vari : public op_vd_vari {
      public:
        add_vd_vari(vari* avi, double b) :
          op_vd_vari(avi->val_ + b, avi, b) {
        }
        void chain() {
          avi_->adj_ += adj_;
        }
      };

      class increment_vari : public op_v_vari {
      public:
        increment_vari(vari* avi) :
          op_v_vari(avi->val_ + 1.0, avi) {
        }
        void chain() {
          avi_->adj_ += adj_;
        }
      };

      class decrement_vari : public op_v_vari {
      public:
        decrement_vari(vari* avi) :
          op_v_vari(avi->val_ - 1.0, avi) {
        }
        void chain() {
          avi_->adj_ += adj_;
        }
      };

      class subtract_vv_vari : public op_vv_vari {
      public:
        subtract_vv_vari(vari* avi, vari* bvi) :
          op_vv_vari(avi->val_ - bvi->val_, avi, bvi) {
        }
        void chain() {
          avi_->adj_ += adj_;
          bvi_->adj_ -= adj_;
        }
      };
    
      class subtract_vd_vari : public op_vd_vari {
      public:
        subtract_vd_vari(vari* avi, double b) :
          op_vd_vari(avi->val_ - b, avi, b) {
        }
        void chain() {
          avi_->adj_ += adj_;
        }
      };

      class subtract_dv_vari : public op_dv_vari {
      public:
        subtract_dv_vari(double a, vari* bvi) :
          op_dv_vari(a - bvi->val_, a, bvi) {
        }
        void chain() {
          bvi_->adj_ -= adj_;
        }
      };

      class multiply_vv_vari : public op_vv_vari {
      public:
        multiply_vv_vari(vari* avi, vari* bvi) :
          op_vv_vari(avi->val_ * bvi->val_, avi, bvi) {
        }
        void chain() {
          avi_->adj_ += bvi_->val_ * adj_;
          bvi_->adj_ += avi_->val_ * adj_;
        }
      };

      class multiply_vd_vari : public op_vd_vari {
      public:
        multiply_vd_vari(vari* avi, double b) :
          op_vd_vari(avi->val_ * b, avi, b) {
        }
        void chain() {
          avi_->adj_ += adj_ * bd_;
        }
      };

      // (a/b)' = a' * (1 / b) - b' * (a / [b * b])
      class divide_vv_vari : public op_vv_vari {
      public:
        divide_vv_vari(vari* avi, vari* bvi) :
          op_vv_vari(avi->val_ / bvi->val_, avi, bvi) {
        }
        void chain() {
          avi_->adj_ += adj_ / bvi_->val_;
          bvi_->adj_ -= adj_ * avi_->val_ / (bvi_->val_ * bvi_->val_);
        }
      };

      class divide_vd_vari : public op_vd_vari {
      public:
        divide_vd_vari(vari* avi, double b) :
          op_vd_vari(avi->val_ / b, avi, b) {
        }
        void chain() {
          avi_->adj_ += adj_ / bd_;
        }
      };

      class divide_dv_vari : public op_dv_vari {
      public:
        divide_dv_vari(double a, vari* bvi) :
          op_dv_vari(a / bvi->val_, a, bvi) {
        }
        void chain() {
          bvi_->adj_ -= adj_ * ad_ / (bvi_->val_ * bvi_->val_);
        }
      };

      class exp_vari : public op_v_vari {
      public:
        exp_vari(vari* avi) :
          op_v_vari(std::exp(avi->val_),avi) {
        }
        void chain() {
          avi_->adj_ += adj_ * val_;
        }
      };

      class log_vari : public op_v_vari {
      public:
        log_vari(vari* avi) :
          op_v_vari(std::log(avi->val_),avi) {
        }
        void chain() {
          avi_->adj_ += adj_ / avi_->val_;
        }
      };

      double LOG_10 = std::log(10.0);
    
      class log10_vari : public op_v_vari {
      public:
        const double exp_val_;
        log10_vari(vari* avi) :
          op_v_vari(std::log10(avi->val_),avi),
          exp_val_(avi->val_) {
        }
        void chain() {
          avi_->adj_ += adj_ / (LOG_10 * exp_val_);
        }
      };

      class sqrt_vari : public op_v_vari {
      public:
        sqrt_vari(vari* avi) :
          op_v_vari(std::sqrt(avi->val_),avi) {
        }
        void chain() {
          avi_->adj_ += adj_ / (2.0 * val_);
        }
      };

      class pow_vv_vari : public op_vv_vari {
      public:
        pow_vv_vari(vari* avi, vari* bvi) :
          op_vv_vari(std::pow(avi->val_,bvi->val_),avi,bvi) {
        }
        void chain() {
          if (avi_->val_ == 0.0) return; // partials zero, avoids /0 & log(0)
          avi_->adj_ += adj_ * bvi_->val_ * val_ / avi_->val_;
          bvi_->adj_ += adj_ * std::log(avi_->val_) * val_;
        }
      };

      class pow_vd_vari : public op_vd_vari {
      public:
        pow_vd_vari(vari* avi, double b) :
          op_vd_vari(std::pow(avi->val_,b),avi,b) {
        }
        void chain() {
          if (avi_->val_ == 0.0) return; // partials zero, avoids /0 & log(0)
          avi_->adj_ += adj_ * bd_ * val_ / avi_->val_;
        }
      };

      class pow_dv_vari : public op_dv_vari {
      public:
        pow_dv_vari(double a, vari* bvi) :
          op_dv_vari(std::pow(a,bvi->val_),a,bvi) {
        }
        void chain() {
          if (ad_ == 0.0) return; // partials zero, avoids /0 & log(0)
          bvi_->adj_ += adj_ * std::log(ad_) * val_;
        }
      };

      class cos_vari : public op_v_vari {
      public:
        cos_vari(vari* avi) :
          op_v_vari(std::cos(avi->val_),avi) {
        }
        void chain() {
          avi_->adj_ -= adj_ * std::sin(avi_->val_);
        }
      };

      class sin_vari : public op_v_vari {
      public:
        sin_vari(vari* avi) :
          op_v_vari(std::sin(avi->val_),avi) {
        }
        void chain() {
          avi_->adj_ += adj_ * std::cos(avi_->val_);
        }
      };

      class tan_vari : public op_v_vari {
      public:
        tan_vari(vari* avi) :
          op_v_vari(std::tan(avi->val_),avi) {
        }
        void chain() {
          avi_->adj_ += adj_ * (1.0 + val_ * val_); 
        }
      };

      class acos_vari : public op_v_vari {
      public:
        acos_vari(vari* avi) :
          op_v_vari(std::acos(avi->val_),avi) {
        }
        void chain() {
          avi_->adj_ -= adj_ / std::sqrt(1.0 - (avi_->val_ * avi_->val_));
        }
      };

      class asin_vari : public op_v_vari {
      public:
        asin_vari(vari* avi) :
          op_v_vari(std::asin(avi->val_),avi) {
        }
        void chain() {
          avi_->adj_ += adj_ / std::sqrt(1.0 - (avi_->val_ * avi_->val_));
        }
      };

      class atan_vari : public op_v_vari {
      public:
        atan_vari(vari* avi) :
          op_v_vari(std::atan(avi->val_),avi) {
        }
        void chain() {
          avi_->adj_ += adj_ / (1.0 + (avi_->val_ * avi_->val_));
        }
      };

      class atan2_vv_vari : public op_vv_vari {
      public:
        atan2_vv_vari(vari* avi, vari* bvi) :
          op_vv_vari(std::atan2(avi->val_,bvi->val_),avi,bvi) {
        }
        void chain() {
          double a_sq_plus_b_sq = (avi_->val_ * avi_->val_) + (bvi_->val_ * bvi_->val_);
          avi_->adj_ += bvi_->val_ / a_sq_plus_b_sq;
          bvi_->adj_ -= avi_->val_ / a_sq_plus_b_sq;
        }
      };

      class atan2_vd_vari : public op_vd_vari {
      public:
        atan2_vd_vari(vari* avi, double b) :
          op_vd_vari(std::atan2(avi->val_,b),avi,b) {
        }
        void chain() {
          double a_sq_plus_b_sq = (avi_->val_ * avi_->val_) + (bd_ * bd_);
          avi_->adj_ += bd_ / a_sq_plus_b_sq;
        }
      };

      class atan2_dv_vari : public op_dv_vari {
      public:
        atan2_dv_vari(double a, vari* bvi) :
          op_dv_vari(std::atan2(a,bvi->val_),a,bvi) {
        }
        void chain() {
          double a_sq_plus_b_sq = (ad_ * ad_) + (bvi_->val_ * bvi_->val_);
          bvi_->adj_ -= ad_ / a_sq_plus_b_sq;
        }
      };

      class cosh_vari : public op_v_vari {
      public:
        cosh_vari(vari* avi) :
          op_v_vari(std::cosh(avi->val_),avi) {
        }
        void chain() {
          avi_->adj_ += adj_ * std::sinh(avi_->val_);
        }
      };

      class sinh_vari : public op_v_vari {
      public:
        sinh_vari(vari* avi) :
          op_v_vari(std::sinh(avi->val_),avi) {
        }
        void chain() {
          avi_->adj_ += adj_ * std::cosh(avi_->val_);
        }
      };

      class tanh_vari : public op_v_vari {
      public:
        tanh_vari(vari* avi) :
          op_v_vari(std::tanh(avi->val_),avi) {
        }
        void chain() {
          double cosh = std::cosh(avi_->val_);
          avi_->adj_ += adj_ / (cosh * cosh);
        }
      };


      class floor_vari : public vari {
      public:
        floor_vari(vari* avi) :
          vari(std::floor(avi->val_)) {
        }
      };

      class ceil_vari : public vari {
      public:
        ceil_vari(vari* avi) :
          vari(std::ceil(avi->val_)) {
        }
      };

      class fmod_vv_vari : public op_vv_vari {
      public:
        fmod_vv_vari(vari* avi, vari* bvi) :
          op_vv_vari(std::fmod(avi->val_,bvi->val_),avi,bvi) {
        }
        void chain() {
          avi_->adj_ += adj_;
          bvi_->adj_ -= adj_ * static_cast<int>(avi_->val_ / bvi_->val_);
        }
      };

      class fmod_vd_vari : public op_v_vari {
      public:
        fmod_vd_vari(vari* avi, double b) :
          op_v_vari(std::fmod(avi->val_,b),avi) {
        }
        void chain() {
          avi_->adj_ += adj_;
        }
      };
      
      class fmod_dv_vari : public op_dv_vari {
      public:
        fmod_dv_vari(double a, vari* bvi) :
          op_dv_vari(std::fmod(a,bvi->val_),a,bvi) {
        }
        void chain() {
          int d = static_cast<int>(ad_ / bvi_->val_);
          bvi_->adj_ -= adj_ * d;
        }
      };





    }

    // ********************* vari UP, var DOWN ***********************************************************


    class var {
    public:

      typedef double Scalar;

      vari * vi_;

      explicit var(vari* vi) :
        vi_(vi) {
      }

      var() :
        vi_(0) {
      }

      var(bool b) :
        vi_(new vari(static_cast<double>(b))) {
      }

      var(char c) :
        vi_(new vari(static_cast<double>(c))) {
      }

      var(short n) :
        vi_(new vari(static_cast<double>(n))) {
      }

      var(unsigned short n) :
        vi_(new vari(static_cast<double>(n))) {
      }

      var(int n) :
        vi_(new vari(static_cast<double>(n))) {
      }

      var(unsigned int n) :
        vi_(new vari(static_cast<double>(n))) {
      }

      var(long int n) :
        vi_(new vari(static_cast<double>(n))) {
      }

      var(unsigned long int n) :
        vi_(new vari(static_cast<double>(n))) {
      }

      var(float x) :
        vi_(new vari(static_cast<double>(x))) {
      }

      var(double x) :
        vi_(new vari(x)) {
      }

      var(long double x) :
        vi_(new vari(static_cast<double>(x))) {
      }

      inline double val() const {
        return vi_->val_;
      }

      void grad(std::vector<var>& x,
                std::vector<double>& g) {
        vari::grad(vi_);
        g.resize(x.size());
        for (size_t i = 0U; i < x.size(); ++i) 
          g[i] = x[i].vi_->adj_;
        vari::recover_memory();
      }

      void grad() {
        vari::grad(vi_);
        vari::recover_memory();
      }

      // COMPOUND ASSIGNMENT OPERATORS
    
      inline var& operator+=(const var& b) {
        vi_ = new add_vv_vari(vi_,b.vi_);
        return *this;
      }

      inline var& operator+=(const double& b) {
        vi_ = new add_vd_vari(vi_,b);
        return *this;
      }

      inline var& operator-=(const var& b) {
        vi_ = new subtract_vv_vari(vi_,b.vi_);
        return *this;
      }

      inline var& operator-=(const double& b) {
        vi_ = new subtract_vd_vari(vi_,b);
        return *this;
      }

      inline var& operator*=(const var& b) {
        vi_ = new multiply_vv_vari(vi_,b.vi_);
        return *this;
      }

      inline var& operator*=(const double& b) {
        vi_ = new multiply_vd_vari(vi_,b);
        return *this;
      }

      inline var& operator/=(const var& b) {
        vi_ = new divide_vv_vari(vi_,b.vi_);
        return *this;
      }

      inline var& operator/=(const double& b) {
        vi_ = new divide_vd_vari(vi_,b);
        return *this;
      };


    };

    // COMPARISON OPERATORS

    inline bool operator==(const var& a, const var& b) {
      return a.val() == b.val();
    }

    inline bool operator==(const var& a, const double& b) {
      return a.val() == b;
    }
  
    inline bool operator==(const double& a, const var& b) {
      return a == b.val();
    }

    inline bool operator!=(const var& a, const var& b) {
      return a.val() != b.val();
    }

    inline bool operator!=(const var& a, const double& b) {
      return a.val() != b;
    }

    inline bool operator!=(const double& a, const var& b) {
      return a != b.val();
    }

    inline bool operator<(const var& a, const var& b) {
      return a.val() < b.val();
    }

    inline bool operator<(const var& a, const double& b) {
      return a.val() < b;
    }

    inline bool operator<(const double& a, const var& b) {
      return a < b.val();
    }

    inline bool operator>(const var& a, const var& b) {
      return a.val() > b.val();
    }

    inline bool operator>(const var& a, const double& b) {
      return a.val() > b;
    }

    inline bool operator>(const double& a, const var& b) {
      return a > b.val();
    }

    inline bool operator<=(const var& a, const var& b) {
      return a.val() <= b.val();
    }

    inline bool operator<=(const var& a, const double& b) {
      return a.val() <= b;
    }

    inline bool operator<=(const double& a, const var& b) {
      return a <= b.val();
    }

    inline bool operator>=(const var& a, const var& b) {
      return a.val() >= b.val();
    }

    inline bool operator>=(const var& a, const double& b) {
      return a.val() >= b;
    }

    inline bool operator>=(const double& a, const var& b) {
      return a >= b.val();
    }

    // LOGICAL OPERATORS

    inline bool operator!(const var& a) {
      return !a.val();
    }

    // ARITHMETIC OPERATORS

    inline var operator+(const var& a) {
      return a;
    }

    inline var operator-(const var& a) {
      return var(new neg_vari(a.vi_));
    }

    inline var operator+(const var& a, const var& b) {    
      return var(new add_vv_vari(a.vi_,b.vi_));
    }


    inline var operator+(const var& a, const double& b) {
      return var(new add_vd_vari(a.vi_,b));
    }

    inline var operator+(const double& a, const var& b) {
      return var(new add_vd_vari(b.vi_,a)); // by symmetry
    }

    inline var operator-(const var& a, const var& b) {
      return var(new subtract_vv_vari(a.vi_,b.vi_));
    }

    inline var operator-(const var& a, const double& b) {
      return var(new subtract_vd_vari(a.vi_,b));
    }

    inline var operator-(const double& a, const var& b) {
      return var(new subtract_dv_vari(a,b.vi_));
    }

    inline var operator*(const var& a, const var& b) {
      return var(new multiply_vv_vari(a.vi_,b.vi_));
    }

    inline var operator*(const var& a, const double& b) {
      return var(new multiply_vd_vari(a.vi_,b));
    }

    inline var operator*(const double& a, const var& b) {
      return var(new multiply_vd_vari(b.vi_,a)); // by symmetry
    }

    inline var operator/(const var& a, const var& b) {
      return var(new divide_vv_vari(a.vi_,b.vi_));
    }

    inline var operator/(const var& a, const double& b) {
      return var(new divide_vd_vari(a.vi_,b));
    }

    inline var operator/(const double& a, const var& b) {
      return var(new divide_dv_vari(a,b.vi_));
    }

    inline var& operator++(var& a) {
      a.vi_ = new increment_vari(a.vi_);
      return a;
    }

    inline var operator++(var& a, int /*dummy*/) {
      var temp(a);
      a.vi_ = new increment_vari(a.vi_);
      return temp;
    }

    inline var& operator--(var& a) {
      a.vi_ = new decrement_vari(a.vi_);
      return a;
    }

    inline var operator--(var& a, int /*dummy*/) {
      var temp(a);
      a.vi_ = new decrement_vari(a.vi_);
      return temp;
    }

    // CMATH EXP AND LOG
  
    inline var exp(const var& a) {
      return var(new exp_vari(a.vi_));
    }

    inline var log(const var& a) {
      return var(new log_vari(a.vi_));
    }

    inline var log10(const var& a) {
      return var(new log10_vari(a.vi_));
    }


    // POWER FUNCTIONS

    inline var sqrt(const var& a) {
      return var(new sqrt_vari(a.vi_));
    }

    inline var pow(const var& base, const var& exponent) {
      return var(new pow_vv_vari(base.vi_,exponent.vi_));
    }
  
    inline var pow(const var& base, const double& exponent) {
      return var(new pow_vd_vari(base.vi_,exponent));
    }

    inline var pow(const double& base, const var& exponent) {
      return var(new pow_dv_vari(base,exponent.vi_));
    }


    // TRIG FUNCTIONS

    inline var cos(const var& a) {
      return var(new cos_vari(a.vi_));
    }

    inline var sin(const var& a) {
      return var(new sin_vari(a.vi_));
    }

    inline var tan(const var& a) {
      return var(new tan_vari(a.vi_));
    }

    inline var acos(const var& a) {
      return var(new acos_vari(a.vi_));
    }

    inline var asin(const var& a) {
      return var(new asin_vari(a.vi_));
    }

    inline var atan(const var& a) {
      return var(new atan_vari(a.vi_));
    }

    inline var atan2(const var& a, const var& b) {
      return var(new atan2_vv_vari(a.vi_,b.vi_));
    }

    inline var atan2(const var& a, const double& b) {
      return var(new atan2_vd_vari(a.vi_,b));
    }

    inline var atan2(const double& a, const var& b) {
      return var(new atan2_dv_vari(a,b.vi_));
    }

    // HYPERBOLIC FUNCTIONS

    inline var cosh(const var& a) {
      return var(new cosh_vari(a.vi_));
    }

    inline var sinh(const var& a) {
      return var(new sinh_vari(a.vi_));
    }

    inline var tanh(const var& a) {
      return var(new tanh_vari(a.vi_));
    }
  

    // ROUNDING FUNCTIONS

    inline var fabs(const var& a) {
      // cut-and-paste from abs()
      if (a.val() > 0.0)
        return a;
      if (a.val() < 0.0)
        return var(new neg_vari(a.vi_));
      return var(new vari(0.0));
    }

    inline var floor(const var& a) {
      return var(new floor_vari(a.vi_));
    }

    inline var ceil(const var& a) {
      return var(new ceil_vari(a.vi_));
    }

    inline var fmod(const var& a, const var& b) {
      return var(new fmod_vv_vari(a.vi_,b.vi_));
    }
  
    inline var fmod(const var& a, const double& b) {
      return var(new fmod_vd_vari(a.vi_,b));
    }

    inline var fmod(const double& a, const var& b) {
      return var(new fmod_dv_vari(a,b.vi_));
    }


    // STD LIB FUNCTIONS

    inline var abs(const var& a) {   
      // cut-and-paste from fabs()
      if (a.val() > 0.0)
        return a;
      if (a.val() < 0.0)
        return var(new neg_vari(a.vi_));
      return var(new vari(0.0));
    }

  }

}

#endif