First lets prove the left bound of $$$\lfloor \frac{n}{i} \rfloor = q$$$

We can use the classic division algorithm to write $$$n = ql + r$$$ with $$$0 \leq r \lt l$$$ and $$$r \in \mathbb{Z}$$$ which implies $$$0 \lt r+1 \leq l$$$

$$$r = n - ql \iff r + 1 = n - ql + 1 \leq l \iff ql + l \geq n + 1 \iff l(q + 1) \geq n+1 \iff l \geq \frac{n+1}{q+1}$$$

So, the minimum possible $$$l$$$ would be $$$\lceil \frac{n+1}{q+1} \rceil$$$ $$$\square$$$

For the proof of the right bound check Proof of Fact 3 in the blog.

First lets prove the left bound of $$$\lceil \frac{n}{i} \rceil = q$$$

We can use the division algorithm to write $$$n = (q-1)l + r$$$ where $$$0 \lt r \leq l$$$

$$$r = n - (q-1)l \leq l \iff (q-1)l + l \geq n \iff l(q-1+1) \geq n \iff l \geq \frac{n}{q}$$$

So, the minimum possible $$$l$$$ would be $$$\lceil \frac{n}{q} \rceil$$$ $$$\square$$$

The proof of the right bound is kinda similar to the Proof of Fact 3 in the blog.

Suppose $$$l$$$ is the minimum value for which $$$\lceil \frac{n}{l} \rceil = q$$$. Now, suppose that for some other $$$l^\prime$$$ ($$$l^\prime \geq l$$$), we have that $$$\lceil \frac{n}{l^\prime} \rceil = q$$$, and so we write $$$n = (q-1)l^\prime + r^\prime$$$ with $$$0 \lt r^\prime \leq l$$$ and $$$r^\prime \in \mathbb{Z}$$$ which implies $$$0 \leq r^\prime-1 \lt l^\prime$$$.

$$$r^\prime = n - (q-1)l^\prime \iff r^\prime - 1 = n - (q-1)l^\prime - 1 \geq 0 \iff (q-1)l^\prime \leq n - 1 \iff l^\prime \leq \frac{n-1}{q-1}$$$

So, the maximum possible $$$l^\prime$$$ would be $$$\lfloor \frac{n-1}{q-1} \rfloor$$$ $$$\square$$$

#include <iostream>
#include <cmath>
using namespace std;
#define endl '\n'
 
typedef long long ll;
const int MOD = 1e9 + 7;
 
ll sqr(ll x){
    ll ret = max(0LL, (ll)sqrt(x) - 3);
    while(ret*ret <= x)
        ++ret;
    return ret-1;
}
 
ll calc(ll l, ll r){
    if (l > r) return 0;
    l %= MOD, r %= MOD;
    ll a = (l-1) * l / 2 % MOD;
    ll b = r * (r+1) / 2 % MOD;
    return (b-a+MOD)%MOD;
}
 
int main(){
    cin.tie(0)->sync_with_stdio(0);
    cin.exceptions(cin.failbit);
 
    ll n, s, i, ans, l, r;
    cin>>n;
    s = sqr(n);
    ans = 0;
    for (i = 1; i <= s; ++i)
        (ans += n/i * i) %= MOD;
    for (i = 1; i <= s; ++i){
        l = max(s+1, (n+1 + i)/(i+1));
        r = n/i;
        (ans += i * calc(l, r)) %= MOD;
    }
    cout<<ans<<endl;
 
    return 0;
}

I used the max function in the second loop to prevent double counting in the range [1, s] since they are already accounted for in the first loop.

// Author: Sahil Yasar

#include <iostream>
#include <vector>
#include <cstring>
using namespace std;
#define endl '\n'

namespace persistentSegTree{
    const int N = 1e5 + 10;
    const int M = 5e6 + 10;
    typedef int T;

    T arr[M];
    int L[M], R[M], root[N];
    size_t sz;
    int nodes = 0, cnt = 0;
    T val;
    T (*f)(T, T);

    void build(T a[], int v, int tl, int tr){
        if (tl == tr) arr[v] = a[tl];
        else{
            int tm = (tl + tr) / 2;
            build(a, L[v] = nodes++, tl, tm);
            build(a, R[v] = nodes++, tm+1, tr);
            arr[v] = f(arr[L[v]], arr[R[v]]);
        }
    }
    void PSEGtree(T a[], size_t n, T (*func)(T, T), T x){
        sz = n, f = func, val = x;
        build(a, root[cnt++] = nodes++, 0, sz-1);
    }
    void reset(){ nodes = cnt = 0; }

    T query(int v, int l, int r, int tl, int tr){
        if (l > tr || r < tl) return val;
        if (l <= tl && r >= tr) return arr[v];
        int tm = (tl + tr) / 2;
        return f(query(L[v], l, r, tl, tm),
                query(R[v], l, r, tm+1, tr));
    }
    T query(int l, int r, int rt){
        return query(root[rt], l, r, 0, sz-1);
    }

    void update(int cur, int prev, int pos, T x, int tl, int tr){
        if (tl == tr) return void(arr[cur] = x);
        int tm = (tl + tr) / 2;
        if (pos <= tm){
            R[cur] = R[prev], L[cur] = nodes++;
            update(L[cur], L[prev], pos, x, tl, tm);
        }
        else{
            L[cur] = L[prev], R[cur] = nodes++;
            update(R[cur], R[prev], pos, x, tm+1, tr);
        }
        arr[cur] = f(arr[L[cur]], arr[R[cur]]);
    }
    int update(int pos, T x, int rt = -1){
        if (rt < 0) rt = cnt-1;
        update(root[cnt] = nodes++, root[rt], pos, x, 0, sz-1);
        return cnt++;
    }
};
using namespace persistentSegTree;

const int MAX = 1e5 + 10;
vector<int> g[MAX];
int in[MAX], out[MAX], timer;
vector<int> node;
void dfs(int s, vector<int>& nums, int p = -1){
    node.push_back(nums[s]-1);
    in[s] = timer++;
    for (int& c: g[s])
        if (c != p)
            dfs(c, nums, s);
    out[s] = timer-1;
}

T rangeMEX(int v, int x, int tl, int tr){
    if (tl == tr) return tl;
    int tm = (tl + tr) / 2;
    if (arr[L[v]] < x)
        return rangeMEX(L[v], x, tl, tm);
    return rangeMEX(R[v], x, tm+1, tr);
}
int rangeMEX(int rt[], int l, int r){
    return rangeMEX(root[rt[r]], l, 0, sz-1);
}

int main(){
    cin.tie(0)->sync_with_stdio(0);
    cin.exceptions(cin.failbit);

    int n, i;
    cin>>n;
    vector<int> parents(n);
    vector<int> nums(n);
    for (i = 0; i < n; ++i)
        cin>>parents[i];
    for (i = 0; i < n; ++i)
        cin>>nums[i];

    // 0 is root
    for (int i = 1; i < parents.size(); ++i){
        g[i].push_back(parents[i]);
        g[parents[i]].push_back(i);
    }

    int temp[MAX];
    memset(temp, -1, sizeof(temp));
    PSEGtree(temp, MAX, [](int a, int b){ return min(a, b); }, 1e9);
    timer = 0;
    dfs(0, nums);

    int rt[n];
    for (i = 0; i < n; ++i)
        rt[i] = update(node[i], i);
    vector<int> ans;
    for (i = 0; i < n; ++i)
        ans.push_back(rangeMEX(rt, in[i], out[i]) + 1);

    for (i = 0; i < n; ++i)
        cout<<ans[i]<<" ";
    cout<<endl;

    return 0;
}