A log-density class that counts the number of zeros.

LogDensity.h

#include <limits>
#include <cmath>
#include <cassert>
#include <iostream>

class LogDensity
{
    double zeros_ = 0;

    // Cannot be -Inf.  Can be +Inf.
    double ones_ = 0;
public:

    double zeros() const {return zeros_;}
    double& zeros() {return zeros_;}
    
    double ones() const {return ones_;}
    double& ones() {return ones_;}
    
    void check() const {
	assert(not std::isinf(ones_) or ones_ > 0); // ones_ != -Inf
	assert(not std::isinf(zeros_));             // only a finite number of zeros.
    }

    LogDensity& operator *=(double y)
    {
	// No infinite powers
	assert(not std::isinf(y));

	// 0^0 == 1
	zeros_ *= y; // fractional zeros
	ones_ *= y;

	return *this;
    }

    LogDensity& operator /=(double y)
    {
	// No zeroth roots
	assert(y != 0);

	zeros_ /= y;  // fractional zeros
	ones_ /= y;

	return *this;
    }

    LogDensity operator +=(const LogDensity y)
    {
	zeros_ += y.zeros_;
	ones_ += y.ones_;

	return *this;
    }

    LogDensity operator -=(const LogDensity y)
    {
	zeros_ -= y.zeros_;
	ones_ -= y.ones_;

	return *this;
    }

    bool operator<(const LogDensity& y) const
    {
	if (zeros_ > y.zeros_) return true;
	if (zeros_ < y.zeros_) return false;

	return (ones_ < y.ones_);
    }

    bool operator>(const LogDensity& y) const
    {
	if (zeros_ < y.zeros_) return true;
	if (zeros_ > y.zeros_) return false;

	return (ones_ > y.ones_);
    }

    bool operator==(const LogDensity& y) const
    {
	return zeros_ == y.zeros_ and ones_ == y.ones_;
    }

    bool operator!=(const LogDensity& y) const
    {
	return not (operator==(y));
    }

    bool operator<=(const LogDensity& y) const
    {
	return operator<(y) or operator==(y);
    }
					 
    bool operator>=(const LogDensity& y) const
    {
	return operator>(y) or operator==(y);
    }
    
    explicit operator double() const
    {
	check();

	if (std::isinf(ones_) and zeros_ > 0)
	    return std::nan("1"); // NAN
	else if (zeros_ > 0)
	    return 0;
	else if (zeros_ < 0)
	    return std::numeric_limits<double>::infinity();
	else
	    return ones_;
    }

    LogDensity() = default;

    LogDensity(double y)
    {
	if (std::isinf(y) and y<0)
	    zeros_ = 1;
	else
	    ones_ = y;
    }
};

LogDensity operator+(LogDensity x, const LogDensity& y)
{
    x += y;
    return x;
}

LogDensity operator-(LogDensity x, const LogDensity& y)
{
    x -= y;
    return x;
}

LogDensity operator*(double p, LogDensity x)
{
    x *= p;
    return x;
}

LogDensity operator*(LogDensity x, double p)
{
    x *= p;
    return x;
}

LogDensity operator/(LogDensity x, double p)
{
    x /= p;
    return x;
}

std::ostream& operator<<(std::ostream& o, const LogDensity& x)
{
    o<<x.ones()<<" + "<<x.zeros()<<"*-Inf";
    return o;
}

testDensity.cc

#include <iostream>
#include "LogDensity.h"

int main()
{
    double x1 = log(0.5);

    double x2 = log(0);
    
    LogDensity y1 = x1;
    LogDensity y2 = y1 + x2;
    LogDensity y3 = x1 + x2 + x2; // Addition is of double!!!
    LogDensity y4 = y2 + x2;

    std::cerr<<std::boolalpha;
    std::cerr<<"y1 = "<<double(y1)<<"\n";
    std::cerr<<"y2 = "<<double(y2)<<"\n";
    std::cerr<<"y3 = "<<double(y3)<<"\n";
    std::cerr<<"y4 = "<<double(y4)<<"\n";
    std::cerr<<"y1 = "<<y1<<"\n";
    std::cerr<<"y2 = "<<y2<<"\n";
    std::cerr<<"y3 = "<<y3<<"\n";
    std::cerr<<"y4 = "<<y4<<"\n";
    std::cerr<<"y1 > y2 = "<<(y1 > y2)<<"\n";
    std::cerr<<"y2 > y3 = "<<(y2 > y3)<<"\n";
    std::cerr<<"y2 > y4 = "<<(y2 > y4)<<"\n";
    std::cerr<<"y2 - y2 = "<<(y2-y2)<<"\n";
    
    return 0;
}