线段树套线段树套线段树！

设每个原子的翻倍系数为 p，初始值为 p_{i,j,k}=1 即不翻倍。那么所有施加在该原子上的力就是最终的翻倍系数 p=\prod f。且最终质量为 mp。

考虑对 p 维护三维差分。因为要维护的不是加减是乘除，而乘除的单位元素不是 0 而是 1，故差分数组初始值应设为 p=1。每次维护时，在边界乘上 f 或除以 f 即乘以 f 的逆元。并且对差分数组还原为原数组时，进行乘除操作而不是加减操作。

t 次差分操作后对差分数组叠前缀和即得翻倍系数 p。之后便可得到质量数组 mp。再对质量数组求前缀和，即可解决之后的 q 次询问。

考虑到逆元，故复杂度为 O(abc\log p+t\log p+q)。注意细节。

#include <bits/stdc++.h>
using namespace std;
using ll = long long;
const ll mn = 1e2 + 5, mod = 1e9 + 7;
ll qpow(ll x, ll y = mod - 2)
{
    ll r = 1;
    for (; y; y >>= 1)
    {
        if (y & 1)
        {
            r = r * x % mod;
        }
        x = x * x % mod;
    }
    return r;
}
ll a, b, c, t, q, m[mn][mn][mn], p[mn][mn][mn];
signed main()
{
    ios::sync_with_stdio(false), cin.tie(0);
    cin >> a >> b >> c >> t >> q;
    for (ll i = 1; i <= a; ++i)
    {
        for (ll j = 1; j <= b; ++j)
        {
            for (ll k = 1; k <= c; ++k)
            {
                cin >> m[i][j][k];
            }
        }
    }
    for (ll i = 0; i <= a; ++i)
    {
        for (ll j = 0; j <= b; ++j)
        {
            for (ll k = 0; k <= c; ++k)
            {
                p[i][j][k] = 1;
            }
        }
    }
    for (ll xl, xr, yl, yr, zl, zr, f, i = 1; i <= t; ++i)
    {
        cin >> xl >> xr >> yl >> yr >> zl >> zr >> f;
        ll fi = qpow(f);
        (p[xl][yl][zl] *= f) %= mod;
        (p[xr + 1][yl][zl] *= fi) %= mod;
        (p[xl][yr + 1][zl] *= fi) %= mod;
        (p[xl][yl][zr + 1] *= fi) %= mod;
        (p[xl][yr + 1][zr + 1] *= f) %= mod;
        (p[xr + 1][yl][zr + 1] *= f) %= mod;
        (p[xr + 1][yr + 1][zl] *= f) %= mod;
        (p[xr + 1][yr + 1][zr + 1] *= fi) %= mod;
    }
    for (ll i = 1; i <= a; ++i)
    {
        for (ll j = 1; j <= b; ++j)
        {
            for (ll k = 1; k <= c; ++k)
            {
                (p[i][j][k] *= p[i - 1][j][k]) %= mod;
                (p[i][j][k] *= p[i][j - 1][k]) %= mod;
                (p[i][j][k] *= p[i][j][k - 1]) %= mod;
                (p[i][j][k] *= qpow(p[i][j - 1][k - 1])) %= mod;
                (p[i][j][k] *= qpow(p[i - 1][j][k - 1])) %= mod;
                (p[i][j][k] *= qpow(p[i - 1][j - 1][k])) %= mod;
                (p[i][j][k] *= p[i - 1][j - 1][k - 1]) %= mod;
                (m[i][j][k] *= p[i][j][k]) %= mod;
                m[i][j][k] += m[i - 1][j][k];
                m[i][j][k] += m[i][j - 1][k];
                m[i][j][k] += m[i][j][k - 1];
                m[i][j][k] -= m[i][j - 1][k - 1];
                m[i][j][k] -= m[i - 1][j][k - 1];
                m[i][j][k] -= m[i - 1][j - 1][k];
                m[i][j][k] += m[i - 1][j - 1][k - 1];
                m[i][j][k] = (m[i][j][k] + 3 * mod) % mod;
            }
        }
    }
    for (ll xl, xr, yl, yr, zl, zr, i = 1; i <= q; ++i)
    {
        cin >> xl >> xr >> yl >> yr >> zl >> zr;
        ll ans = m[xr][yr][zr];
        ans -= m[xl - 1][yr][zr];
        ans -= m[xr][yl - 1][zr];
        ans -= m[xr][yr][zl - 1];
        ans += m[xr][yl - 1][zl - 1];
        ans += m[xl - 1][yr][zl - 1];
        ans += m[xl - 1][yl - 1][zr];
        ans -= m[xl - 1][yl - 1][zl - 1];
        ans = (ans + 4 * mod) % mod;
        cout << ans << '\n';
    }
    return 0;
}