Is there anything in Sage which solves #SAT ?

edit

Since the solutions of a SAT formula corresponds to integer points of some polytope, you can solve #SAT by asking Sage to construct that polytope and count its integer points. See the code below.

The default integral_points_count for rational polytopes basically enumerates all integer points, so the choice of the backend could have an importance on the speed. I bet normaliz will be faster, it requires the pynormaliz optional package to be installed. Please report about your real-life experiments.

Also, the latte_int optional package could speed things in some cases, please report your experiments, since that could depend on the actual problem (on the code below, you should use latte=True to use it).

That said for further improvements, we really need to know your real problem.

Typically, the polytopes we consider are included in the unit cube of large dimension, hence I do not see how Barvinok theory (which is what latte_int relies on) could speed things up since no set of >=3 points is aligned.

So, depending on your problem, one could imagine to describe it with a "fatter" polytope of smaller dimension, whithout reducing it to SAT first, so that Barvinok theory could be much more efficient.

Other approaches could be thought about too. For example, given a SAT formula, it is easy to enumerate all solution with a SAT solver: for each solution, add its "opposit" as a new constraint and solve it again. I have plans to add such an __iter__ method to the current SAT solvers in Sage.

def sharp_sat(clauses, backend='cdd', latte=False):
    r"""
    Count the number of solutions of the given CNF formula.

    INPUT:

    - clauses is a list of lists of integers representing a CNF formula in the
      DIMACS notation, for example:

      [[1,2,-3], [-1,2]] means

          (x1 or x2 or (not x3)) and ((not 1) or 2)

    - backend : one of `'cdd'`, `'normaliz'`, `'ppl'`, `'polymake'`.

    - latte : if `True`, use `latte_int` (optional package to be installed) for
      the computation.
    """
    p = MixedIntegerLinearProgram(solver='ppl')
    x = p.new_variable(binary=True)
    for clause in clauses:
        p.add_constraint(sum([x[i] if i>0 else (1-x[-i]) for i in clause]) >= 1)

    P = p.polyhedron(backend=backend)
    if latte:
        return P.integral_points_count(explicit_enumeration_threshold=0)
    else:
        return P.integral_points_count()