Функция Эйлера

    for (i = 3; i <= 4000000; i++)
    {
        s = 0; 
        m = static_cast <long long> (sqrt(double(i)));
 
        for (j = 1; j <= m; j++)
            if (i % j == 0)
            {
                t = i / j;
                s += j * ph[t] + (t) * ph[j];
            }
 
        if (m * m == i)
            s -= m * ph[m];
 
        d[i] = d[i-1] + s;
    }

#include <iostream>
#include <vector>
 
using namespace std;
 
// Генерация простых чисел меньших, чем n
// Используется решето Эратосфена
vector<char> gen_primes(size_t n)
{
    vector<char> is_prime(n);
    for (size_t i = 2; i < n; ++i)
        is_prime[i] = 1;
    for (size_t i = 1; n - i*i > 2*i + 1; )
    {
        ++i;
        if (is_prime[i])
            for (size_t j = i*i; j < n; j += i)
                is_prime[j] = 0;
    }
    return is_prime;
}
 
int main()
{
    unsigned const N = 10000;
    vector<int> euler(N + 1);
 
    vector<char> is_prime = gen_primes(euler.size());
 
    // Вычисление функций Эйлера
    // Инициализация
    for (size_t i = 0; i < euler.size(); ++i)
        euler[i] = i;
    // Расчёт по формуле произведения Эйлера
    for (unsigned i = 2; i < is_prime.size(); ++i)
    {
        if (is_prime[i])
            for (unsigned j = i; j < euler.size(); j += i)
                euler[j] = euler[j] / i * (i - 1);
    }
 
    for (unsigned i = 1; i < euler.size(); ++i)
        cout << "y(" << i << ") = " << euler[i] << endl;
 
    return 0;
}

#include <iostream>
 
#pragma comment (linker, "/stack:65000000")
 
#define MAXN 4000000
 
int main()
 {
    unsigned long long f[MAXN + 1] = {};
    int n, i, j, minp[MAXN + 1] = {}, mult[MAXN + 1] = {};
 
    for (i = 2; (j = i * i) < MAXN + 1; i++) 
        if (!minp[i])
            for (; j < MAXN + 1; j += i) 
                if (!minp[j]) minp[j] = i;
 
    for (i = 0; i < MAXN + 1; i++)
        if (minp[i] == 0)
            minp[i] = i;
 
    f[1] = 1;
        
    for (i = 2; i <= MAXN; i++)
        if ((i / minp[i]) % minp[i])
        {
            mult[i] = 2 * minp[i] - 1;
            j = i / minp[i];
            f[i] = f[j] * mult[i];
        }
        else
        {
            j = i / minp[i];
            mult[i] = mult[j] + minp[i] - 1;
            f[i] = minp[i] * (f[j] / mult[j]) * mult[i];
        }
 
    f[1] = 0;
 
    for (i = 2; i <= MAXN; i++)
        f[i] = f[i] + f[i-1] - i;
 
    while (scanf("%I64d", &n) == 1 && n != 0)
        printf("%I64d\n", f[n]);
 
    return 0;
 }