Вычисление функции Эйлера

#include <omp.h>
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <vector>
 
using namespace std;
long long n, a, b, d;
long double ed;
vector<long long> muls;
vector<long long> emuls;
 
double phi(long double a, long double b, long double c, long double d)
{
    ed = (c / d);
    return a * b * ed;
}
 
int gcd(long long a, long long b)
{
    long nod = 1L;
 
    if (a > b) {
        long tmp = a;
        a = b;
        b = tmp;
    }
 
    while (a > 1L && b > 1L) {
        for (long i = 2; i <= a; i++) {
            if (a % i == 0 && b % i == 0) {
                nod *= i;
                a /= i;
                b /= i;
                break;
            }
            if (a % i == 0) {
                a /= i;
                break;
            }
            if (b % i == 0) {
                b /= i;
                break;
            }
        }
    }
    return nod;
}
 
int main()
{
    muls = vector<long long>(2);
    emuls = vector<long long>(2);
    n = (2LL * 2 * 2 * 2) * (3 * 3 * 3) * (11 * 11 * 11 * 11 * 11) * (13 * 13) * (29 * 29) * 37;
    a = (2LL * 2) * (3 * 3) * (11 * 11) * 13 * 29;
    b = n / a;
    d = gcd(a, b);
    cout << "n = " << n << endl;
    cout << "a = " << a << endl;
    cout << "b = " << b << endl;
    cout << "d = " << d << endl;
 
    if (a * b == n)
    {
        muls.push_back(a);
        muls.push_back(b);
        muls.push_back(d);
        emuls.push_back(a);
        emuls.push_back(b);
        emuls.push_back(d);
    }
    else
        cout << "Error";
    
#pragma omp parallel
    {
        int i;
#pragma omp for private(i)
        for (i = 2; i <= 4; i++)
        {
            for (int j = 2; j * j <= muls[i]; j++)
            {
                if (muls[i] % j == 0)
                {
                    while (muls[i] % j == 0)
                    {
#pragma omp critical
                        {
                        muls[i] /= j;
                        emuls[i] -= emuls[i] / i;
                        }
                    }
                }
            }
            if (n > 1)
            {
                emuls[i] -= emuls[i] / muls[i];
            }
        
            cout << emuls[i] << '\n';
        }
 
        cout << "phi = " << phi(emuls[2], emuls[3], d, emuls[4]) << endl;
 
    }
 
}