1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
#include <bits/stdc++.h>

#define x first
#define y second

using namespace std;

typedef unsigned int uint;


vector <pair<int,int>> getVecs(int k) {
    set <int> squares;
    for (int i = 1; i <= k; i++) {
        squares.insert(i * i);
    }

    vector <pair<int,int>> vecs;
    for (int x = 1; x <= k; x++) for (int y = 0; y <= k; y++) if (__gcd(x, y) == 1) {
        int r = x * x + y * y;
        if (squares.count(r)) {
            vecs.push_back({x, y});
        }
    }

    sort(vecs.begin(), vecs.end(), [&](const pair <int,int> &a, const pair <int,int> &b) {
        return a.y * b.x < a.x * b.y;
    });

    return vecs;
}

uint solve(int xMax, int yMax, int k) {
    auto vecs = getVecs(k);

    int xBound = 0, yBound = 0;
    for (auto &e : vecs) {
        xBound += e.x;
        yBound += e.y;
    }

    int bound = min(max(xBound, yBound) + 1, max(xMax, yMax));

    vector <vector<uint>> dp(bound, vector <uint> (bound, 0));
    dp[0][0] = 1;

    for (auto &e : vecs) {
        for (int x = bound - e.x - 1; x >= 0; x--) for (int y = bound - e.y - 1; y >= 0; y--) {
            dp[x + e.x][y + e.y] += dp[x][y];
        }
    }

    vector <vector<vector<uint>>> cnt(2 * bound - 1,
        vector <vector<uint>> (2 * bound - 1, vector <uint> (bound, 0))
    );

    for (int x1 = 0; x1 < bound; x1++) for (int y1 = 0; y1 < bound; y1++) if (dp[x1][y1]) {
        for (int x2 = 0; x2 < bound; x2++) for (int y2 = 0; y2 < bound; y2++) if (dp[x2][y2]) {
            int x = x1 - y2;
            int y = y1 + x2;

            assert(x > -bound && x < bound);
            assert(y >= 0 && y < 2 * bound - 1);

            uint here = dp[x1][y1] * dp[x2][y2];
            cnt[x + bound - 1][y][x1] += here;
        }
    }

    uint ans = 0;
    for (int x = -(bound - 1); x < bound; x++) {
        for (int y = 0; y < min(2 * bound - 1, yMax); y++) {
            for (int rx = 0; rx < bound; rx++) if (cnt[x + bound - 1][y][rx]) {
                for (int lx = 0; lx < bound; lx++) if (cnt[x + bound - 1][y][lx]) {
                    uint here = cnt[x + bound - 1][y][rx] * cnt[x + bound - 1][y][lx];

                    int height = y;
                    int width = rx + (lx - x);

                    if (width < xMax && height < yMax) {
                        uint waysToPlace = (uint) (yMax - height) * (xMax - width);
                        ans += here * waysToPlace;
                    }
                }
            }
        }
    }

    ans -= (uint) xMax * yMax;
    for (auto &e : vecs) {
        int x = e.x;
        int y = e.y;

        for (int rot = 0; rot < 2; rot++) {
            tie(x, y) = make_pair(-y, x);

            int height = abs(y);
            int width = abs(x);

            if (width < xMax && height < yMax) {
                ans -= (uint) (yMax - height) * (xMax - width);
            }
        }
    }

    return ans;
}

int main() {
    ios_base::sync_with_stdio(0);

    int x, y, k;
    cin >> x >> y >> k;

    cout << solve(x, y, k);

    return 0;
}