cpp_src/number_theory/GaussSum.hpp

View this file on GitHub
Last update: 2023-12-31 13:48:59+09:00
Include: #include "cpp_src/number_theory/GaussSum.hpp"
Link: http://mathforum.org/library/drmath/view/73120.html
Link: https://min-25.hatenablog.com/entry/2018/04/27/225535

Code

// x_i=floor((a*i+b)/c), i=0,1,..n-1
// a,c>0, b>=0
// http://mathforum.org/library/drmath/view/73120.html
// https://min-25.hatenablog.com/entry/2018/04/27/225535

ll gauss_sum(ll n, ll a, ll b, ll c) {
    if (n == 0) return 0;
    ll res = 0;
    {
        ll p = a / c;
        res += n * (n - 1) / 2 * p;
        a %= c;
    }
    {
        ll p = b / c;
        res += n * p;
        b %= c;
    }
    if (a == 0) return res;
    ll top = (a * (n - 1) + b) / c;
    res += top * n;
    ll h = 1;
    if (h <= top) {
        res -= gauss_sum(top - h + 1, c, c * h - (b + 1), a) + top - h + 1;
    }
    return res;
}

#line 1 "cpp_src/number_theory/GaussSum.hpp"
// x_i=floor((a*i+b)/c), i=0,1,..n-1
// a,c>0, b>=0
// http://mathforum.org/library/drmath/view/73120.html
// https://min-25.hatenablog.com/entry/2018/04/27/225535

ll gauss_sum(ll n, ll a, ll b, ll c) {
    if (n == 0) return 0;
    ll res = 0;
    {
        ll p = a / c;
        res += n * (n - 1) / 2 * p;
        a %= c;
    }
    {
        ll p = b / c;
        res += n * p;
        b %= c;
    }
    if (a == 0) return res;
    ll top = (a * (n - 1) + b) / c;
    res += top * n;
    ll h = 1;
    if (h <= top) {
        res -= gauss_sum(top - h + 1, c, c * h - (b + 1), a) + top - h + 1;
    }
    return res;
}