Policzenie wszystkich kombinacji liczby danej wzorem

#include <iostream>
#include <iomanip>
using namespace std;
#include <conio.h>
#include <string>
#include <vector>
#include <cmath>
#include <limits>

//====================================================================================================================

double computeW( int pX, int pY, int pP, const vector < vector < int >>& pXn, const vector < int >& pXIndxs );
bool isNaturalNum( double value );
double newton( unsigned int n, unsigned int k );

//====================================================================================================================

int main()
{
    int x = 73;
    int y = 5;
    int p = 49667;
   
    cout << "Podaj x: ";
    cin >> x;
    cout << "Podaj y: ";
    cin >> y;
    cout << "Podaj p: ";
    cin >> p;
   
    vector < vector < int >> xn( y - 1, vector < int >( x + 1 ) );
   
    for( auto & vec: xn )
    {
        for( int i = 0; i <= x; ++i )
        {
            vec[ i ] = i;
        }
    }
   
    vector < int > xIndxs( xn.size(), 0 );
    double w;
    bool exit = false;
    double counter = 0;
    unsigned long long howManyFound = 0;
   
    cout << "Start " << endl;
   
    while( !exit )
    {
        counter++;
       
        w = computeW( x, y, p, xn, xIndxs );
        if( isNaturalNum( w ) )
        {
            ++howManyFound;
           
            cout << "w= " << fixed << w << "; ";
            for( int i = 0; i < xn.size(); ++i )
            {
                cout << "x" << i + 1 << "=" << xn[ i ][ xIndxs[ i ] ] <<( i == xn.size() - 1 ? ""
                    : ", " );
            }
            cout << endl;
        }
       
        ++xIndxs[ 0 ];
       
        for( int j = xIndxs.size() - 1; j > 0; --j )
        {
            if( xIndxs[ j ] == x )
            {
                if( j == xIndxs.size() - 1 ) { exit = true; break; }
               
                ++xIndxs[ j + 1 ];
               
                for( int pom = j; pom >= 0; --pom )
                {
                    xIndxs[ pom ] = xIndxs[ j + 1 ];
                }
               
                break;
            }
        }
       
        if( xIndxs[ 0 ] == x + 1 )
        {
            if( xIndxs.size() >= 2 )
            {
                ++xIndxs[ 1 ];
                xIndxs[ 0 ] = xIndxs[ 1 ];
            }
            else
            {
                exit = true;
            }
        }
    }
   
    cout << "Koniec" << endl;
    cout << "wszystkich iteracji: " << fixed << setprecision( 0 ) << counter << endl;
    cout << "Policzonych, calkowitych liczb w: " << howManyFound << endl;
   
    getch();
    return 0;
}
//====================================================================================================================
double computeW( int pX, int pY, int pP, const vector < vector < int >>& pXn, const vector < int >& pXIndxs )
{
    double retVal = 0;
   
    int lastPow;
    for( int i = 0; i < pXn.size(); ++i )
    {
        retVal += pow( 2, pXn[ i ][ pXIndxs[ i ] ] ) * pow( pP, i ) / pow( 2, i );
       
        lastPow = i + 1;
    }
   
    retVal += pow( pP, pY - 1 ) / pow( 2, pY - 1 );
    retVal *= pow( 2, pY - 1 );
    retVal /= pow( 2, pX + pY ) - pow( pP, pY );
   
    return retVal;
}

//*****************************************************************************
bool isNaturalNum( double value )
{
    unsigned long long intPart = value;
    double rest = value - intPart;
   
    return rest == 0;
}
//*****************************************************************************
double newton( unsigned int n, unsigned int k ) // not working
{
    double Wynik = 1;
   
    for( unsigned int i = 1; i <= k; i++ )
    {
        Wynik = Wynik *( n - i + 1 ) / i; // Obliczanie ze wzoru iteracyjnego
    }
   
    return Wynik;
}