Generic mean function

It seems useful to have a generic function to calculate the mean of all elements in a container, so I wrote one. By default it calculates the arithmetic mean, but should be able to accommodate other means (like geometric) if passed the appropriate functors.

Assumptions:

• value_type can be default-constructed to yield a "zero-value" of its type

I am particularly interested in feedback on:

• The validity of my assumptions
• The mean() interface
• How I validate my iterators
• General code style.

mean.hpp

#ifndef MEAN_HPP
#define MEAN_HPP

#include <algorithm>
#include <functional>
#include <iterator>
#include <numeric>
#include <type_traits>

template <class ForwardIt, class DivOp = std::divides<>, class AddOp = std::plus<>,
          class = std::enable_if_t<
                  std::is_base_of<
                  std::forward_iterator_tag,
                  typename std::iterator_traits<ForwardIt>::iterator_category>::value>>
typename std::iterator_traits<ForwardIt>::value_type mean(ForwardIt begin,
                                                          ForwardIt end,
                                                          DivOp div = DivOp(),
                                                          AddOp add = AddOp())
{
    using value_type = typename std::iterator_traits<ForwardIt>::value_type;
    const auto sum = std::accumulate(begin, end, value_type(), add);
    const auto count = std::distance(begin, end);
    return div(sum, count);
}

#endif

Below is included some example code to demonstrate the usage of mean(). I'm less interested in feedback on the example code, but it is of course still welcome.

main.cpp

#include "mean.hpp"
#include <iostream>
#include <vector>

struct Point
{
    Point() : x(0), y(0) {}
    Point(int x, int y) : x(x), y(y) {}

    int x, y;
};

Point operator+ (const Point& lhs, const Point& rhs)
{
    return Point(lhs.x + rhs.x, lhs.y + rhs.y);
}

int main()
{
    // Works on primitives
    const std::vector<float> scalars{1, 2, 3, 4, 5, 6, 7, 8, 9};
    const auto scalarResult = mean(scalars.begin(), scalars.end());
    std::cout << scalarResult << std::endl;

    // Works on class types with appropriate operators
    const std::vector<Point> points{{1,1}, {2,2}, {3,3}};
    const auto pointDivision = [] (const Point& sum, const float count)
    {
        return Point(sum.x / count, sum.y / count);
    };
    const auto pointResult = mean(points.begin(), points.end(), pointDivision);
    std::cout << pointResult.x << " " << pointResult.y << std::endl;
}