Самостоятельные числа в С++

Ответ 1

Я построил альтернативное решение C, которое не требует каких-либо операций с модулем или делением:

#include <stdio.h>
#include <string.h>

int main(int argc, char *argv[]) {
   int v[1100000];
   int j1, j2, j3, j4, j5, j6, s, n=0;
   memset(v, 0, sizeof(v));
   for (j6=0; j6<10; j6++) {
      for (j5=0; j5<10; j5++) {
         for (j4=0; j4<10; j4++) {
            for (j3=0; j3<10; j3++) {
               for (j2=0; j2<10; j2++) {
                  for (j1=0; j1<10; j1++) {
                     s = j6 + j5 + j4 + j3 + j2 + j1;
                     v[n + s] = 1;
                     n++;
                  }
               }
            }
         }
      }
   }
   for (n=1; n<=1000000; n++) {
      if (!v[n]) printf("%6d\n", n);
   }
}

Он генерирует 97786 собственных номеров, включая 1 и 1000000.
С выходом требуется

real        0m1.419s
user        0m0.060s
sys         0m0.152s

Когда я перенаправляю вывод в /dev/null, он принимает

real     0m0.030s
user     0m0.024s
sys      0m0.004s

на моей четырехъядерной платформе 3 Ghz.

Для сравнения, ваша версия производит одинаковое количество чисел, поэтому я предполагаю, что мы либо корректны, либо одинаково ошибочны; но ваша версия жует

real    0m0.064s
user    0m0.060s
sys     0m0.000s

при тех же условиях или примерно в 2 раза.

Это или тот факт, что вы используете long s, что не нужно на моей машине. Здесь int достигает 2 миллиардов. Может быть, вы должны проверить INT_MAX на своем?

Обновление

У меня была догадка, что лучше вычислить сумму кусочно. Вот мой новый код:

#include <stdio.h>
#include <string.h>

int main(int argc, char *argv[]) {
   char v[1100000];
   int j1, j2, j3, j4, j5, j6, s, n=0;
   int s1, s2, s3, s4, s5;
   memset(v, 0, sizeof(v));
   for (j6=0; j6<10; j6++) {
      for (j5=0; j5<10; j5++) {
         s5 = j6 + j5;
         for (j4=0; j4<10; j4++) {
            s4 = s5 + j4;
            for (j3=0; j3<10; j3++) {
               s3 = s4 + j3;
               for (j2=0; j2<10; j2++) {
                  s2 = s3 + j2;
                  for (j1=0; j1<10; j1++) {
                     v[s2 + j1 + n++] = 1;
                  }
               }
            }
         }
      }
   }
   for (n=1; n<=1000000; n++) {
      if (!v[n]) printf("%d\n", n);
   }
}

... и что вы знаете, что сократило время для верхнего цикла от 12 мс до 4 мс. Или, может быть, 8, мои часы, кажется, немного дрожат там.

Положение дел, резюме

Фактическое обнаружение номеров самодиагностики до 1 М теперь занимает примерно 4 мс, и у меня возникают проблемы с дальнейшими улучшениями. С другой стороны, до тех пор, пока вывод будет на консоли, он будет продолжать занимать около 1,4 секунды, и я стараюсь использовать буферизацию. Время ввода-вывода так сильно затмевает время вычисления, что любая дальнейшая оптимизация будет бесполезной. Таким образом, хотя и вдохновленный дальнейшими комментариями, я решил оставить достаточно одного.

Все цитируемые времена находятся на моей (довольно быстрой) машине и предназначены для сравнения только друг с другом. Ваш пробег может отличаться.

Ответ 2

Создайте числа один раз, скопируйте вывод в свой код как гигантскую строку. Распечатайте строку.

Ответ 3

Те моды (%) выглядят дорогими. Если вам разрешено перейти на базовую 16 (или даже базовую 2), то вы, вероятно, можете скопировать это намного быстрее. Если вам нужно оставаться в десятичной форме, попробуйте создать массив цифр для каждого места (единицы, десятки, сотни) и создать код прокрутки. Это упростит суммирование чисел.

В качестве альтернативы вы можете распознать поведение активной функции ядра (позвоните на нее s):

s = n + f(b,n)

где f(b,n) - сумма цифр числа n в базе b.

Для основания 10 ясно, что, поскольку числа (также известные как наименее значимые) цифры перемещаются из 0,1,2,..., 9, что n и f(b,n) продолжают блокировку при переходе от n до n+1, это только то, что 10% времени, которое 9 катится до 0, что это не так, так:

f(b,n+1) = f(b,n) + 1  // 90% of the time

Таким образом, основная функция self s продвигается как

n+1 + f(b,n+1) = n + 1 + f(b,n) + 1 = n + f(b,n) + 2

s(n+1) = s(n) + 2 // again, 90% of the time

В оставшейся (и легко идентифицируемой) 10% времени 9 возвращается к нулю и добавляет одну к следующей цифре, которая в простейшем случае вычитает (9-1) из общей суммы, но может каскадировать вверх через серию из 9s, чтобы вычесть 99-1, 999-1 и т.д.

Итак, первая оптимизация может удалить большую часть работы из 90% ваших циклов!

if ((n % 10) != 0) 
{
  n + f(b,n) = n-1 + f(b,n-1) + 2;
}

или

if ((n % 10) != 0)
{
  s = old_s + 2;
}

Этого должно быть достаточно, чтобы существенно увеличить производительность без изменения вашего алгоритма.

Если вам нужно больше, то выработайте простой алгоритм для изменения между итерациями для оставшихся 10%.

Ответ 4

Если вы хотите, чтобы ваш результат был быстрым, возможно, стоит изучить замену вывода iostream на простой старый printf() - зависит от правил выигрыша конкурса, важно ли это.

Ответ 5

Многопоточность (используйте разные массивы/диапазоны для каждого потока). Кроме того, не используйте больше потоков, чем количество ядер процессора =)

Ответ 6

cout или printf в цикле будет медленным. Если вы можете удалить любые отпечатки из цикла, вы увидите значительное увеличение производительности.

Ответ 7

Поскольку диапазон ограничен (от 1 до 1000000), максимальная сумма цифр не превышает 9 * 6 = 54. Это означает, что для реализации сита круговой буфер из 54 элементов должен быть совершенно достаточным (и размер сит растет очень медленно по мере увеличения диапазона).

У вас уже есть решение на основе сита, но оно основано на предварительном создании полноразмерного буфера (сито из 1000000 элементов), которое довольно неэлегантно (если не полностью неприемлемо). Производительность вашего решения также страдает от нелокальности доступа к памяти.

Например, это возможно очень простая реализация

#define N 1000000U

void print_self_numbers(void)
{
  #define NMARKS 64U /* make it 64 just in case (and to make division work faster :) */

  unsigned char marks[NMARKS] = { 0 };
  unsigned i, imark;

  for (i = 1, imark = i; i <= N; ++i, imark = (imark + 1) % NMARKS)
  {
    unsigned digits, sum;

    if (!marks[imark])
      printf("%u ", i);
    else
      marks[imark] = 0;

    sum = i;
    for (digits = i; digits > 0; digits /= 10)
      sum += digits % 10;

    marks[sum % NMARKS] = 1;
  }
}

(Я не собираюсь получать наилучшую производительность с точки зрения часов процессора здесь, просто иллюстрируя ключевую идею с помощью циклического буфера.)

Конечно, диапазон можно легко превратить в параметр функции, а размер буфера-буфера можно легко вычислить во время выполнения из диапазона.

Что касается "оптимизаций"... Нет смысла пытаться оптимизировать код, содержащий операции ввода-вывода. Такая оптимизация не будет достигнута. Если вы хотите проанализировать производительность самого алгоритма, вам придется поместить сгенерированные числа в выходной массив и напечатать их позже.

Ответ 8

Для такой простой задачи лучшим вариантом было бы думать об альтернативных алгоритмах для получения того же результата. % 10 обычно не считается быстрой операцией.

Ответ 9

Почему бы не использовать рекуррентное отношение, указанное на странице wikipedia? Это должно быть невероятно быстро.

РЕДАКТИРОВАТЬ: Игнорируйте это. отношение повторения создает некоторые, но не все числа собственных чисел. На самом деле их очень мало. Это не особенно ясно видно на странице thewikipedia: (

Ответ 10

Это может помочь ускорить вывод С++ iostreams:

cin.tie(0);
ios::sync_with_stdio(false);

Поместите их в основное меню, прежде чем вы начнете писать в cout.

Ответ 11

Я создал решение на основе CUDA, основанное на втором алгоритме Carl Smotricz. Код для самоисчисления Self Numbers очень быстрый - на моей машине он выполняется в ~ 45 наносекундах; это примерно в 150 раз быстрее, чем алгоритм Карла Смотрича, который работал на моей машине в 7 миллисекунд.

Однако есть узкое место, и это похоже на интерфейс PCIe. Мне понадобилось 43 миллисекунды для переноса вычисленных данных с графической карты обратно в ОЗУ. Это может быть оптимизировано, и я буду смотреть на это.

Тем не менее, 45 наносезонов довольно быстро прокручиваются. На самом деле это ужасно быстро, и я добавил код в свою программу, которая запускает алгоритм Карла Смотрича и сравнивает результаты для точности. Результаты являются точными. Вот выход программы (скомпилированный в VS2008 64-разрядный, Windows7):

UPDATE

Я перекомпилировал этот код в режиме выпуска с полной оптимизацией и использовал статические библиотеки времени выполнения, со значительными результатами. Оптимизатор, похоже, очень хорошо справился с алгоритмом Карла, сократив время работы от 7 мс до 1 мс. Реализация CUDA ускорилась, от 35 до 20 человек. Копия памяти с видеокарты в ОЗУ не была затронута.

Выход программы:

Running on device: 'Quadro NVS 295'
Reference Implementation Ran In 15603 ticks (7 ms)
Kernel Executed in 40 ms -- Breakdown:
  [kernel] : 35 us (0.09%)
  [memcpy] : 40 ms (99.91%)
CUDA Implementation Ran In 111889 ticks (51 ms)
Compute Slots: 1000448 (1954 blocks X 512 threads)
Number of Errors: 0

Код выглядит следующим образом:

файл: main.h

#pragma once

#include <cstdlib>
#include <functional>

typedef std::pair<int*, size_t> sized_ptr;
static sized_ptr make_sized_ptr(int* ptr, size_t size)
{
    return make_pair<int*, size_t>(ptr, size);
}

__host__ void ComputeSelfNumbers(sized_ptr hostMem, sized_ptr deviceMemory, unsigned const blocks, unsigned const threads);

inline std::string format_elapsed(double d) 
{
    char buf[256] = {0};

    if( d < 0.00000001 )
    {
        // show in ps with 4 digits
        sprintf(buf, "%0.4f ps", d * 1000000000000.0);
    }
    else if( d < 0.00001 )
    {
        // show in ns
        sprintf(buf, "%0.0f ns", d * 1000000000.0);
    }
    else if( d < 0.001 )
    {
        // show in us
        sprintf(buf, "%0.0f us", d * 1000000.0);
    }
    else if( d < 0.1 )
    {
        // show in ms
        sprintf(buf, "%0.0f ms", d * 1000.0);
    }
    else if( d <= 60.0 )
    {
        // show in seconds
        sprintf(buf, "%0.2f s", d);
    }
    else if( d < 3600.0 )
    {
        // show in min:sec
        sprintf(buf, "%01.0f:%02.2f", floor(d/60.0), fmod(d,60.0));
    }
    // show in h:min:sec
    else 
        sprintf(buf, "%01.0f:%02.0f:%02.2f", floor(d/3600.0), floor(fmod(d,3600.0)/60.0), fmod(d,60.0));

    return buf;
}

inline std::string format_pct(double d)
{
    char buf[256] = {0};
    sprintf(buf, "%.2f", 100.0 * d);
    return buf;
}

файл: main.cpp

#define _CRT_SECURE_NO_WARNINGS 

#include <windows.h>
#include "C:\CUDA\include\cuda_runtime.h"
#include <cstdlib>
#include <iostream>
#include <string>
using namespace std;
#include <cmath>
#include <map>
#include <algorithm>
#include <list>

#include "main.h"

int main()
{
    unsigned numVals = 1000000;
    int* gold = new int[numVals];
    memset(gold, 0, sizeof(int)*numVals);

    LARGE_INTEGER li = {0}, li2 = {0};
    QueryPerformanceFrequency(&li);
    __int64 freq = li.QuadPart;

    // get cuda properties...
    cudaDeviceProp cdp = {0};
    cudaError_t err = cudaGetDeviceProperties(&cdp, 0);
cout << "Running on device: '" << cdp.name << "'" << endl;

    // first run the reference implementation
    QueryPerformanceCounter(&li);
    for( int j6=0, n = 0; j6<10; j6++ ) 
    {
        for( int j5=0; j5<10; j5++ ) 
        {
            for( int j4=0; j4<10; j4++ ) 
            {
                for( int j3=0; j3<10; j3++ ) 
                {
                    for( int j2=0; j2<10; j2++ ) 
                    {
                        for( int j1=0; j1<10; j1++ )  
                        {
                            int s = j6 + j5 + j4 + j3 + j2 + j1;
                            gold[n + s] = 1;
                            n++;
                        }
                    }
                }
            }
        }
    }
    QueryPerformanceCounter(&li2);
    __int64 ticks = li2.QuadPart-li.QuadPart;
    cout << "Reference Implementation Ran In " << ticks << " ticks" << " (" << format_elapsed((double)ticks/(double)freq) << ")" << endl;

    // now run the cuda version...
    unsigned threads = cdp.maxThreadsPerBlock;
    unsigned blocks = numVals/threads;
    if( numVals%threads ) ++blocks;
    unsigned computeSlots = blocks * threads;   // this may be != the number of vals since we want 32-thread warps

    // allocate device memory for test
    int* deviceTest = 0;
    err = cudaMalloc(&deviceTest, sizeof(int)*computeSlots);
    err = cudaMemset(deviceTest, 0, sizeof(int)*computeSlots);

    int* hostTest = new int[numVals];   // the repository for the resulting data on the host
    memset(hostTest, 0, sizeof(int)*numVals);

    // run the CUDA code...
    LARGE_INTEGER li3 = {0}, li4={0};
    QueryPerformanceCounter(&li3);
    ComputeSelfNumbers(make_sized_ptr(hostTest, numVals), make_sized_ptr(deviceTest, computeSlots), blocks, threads);
    QueryPerformanceCounter(&li4);

    __int64 ticksCuda = li4.QuadPart-li3.QuadPart;
    cout << "CUDA Implementation Ran In " << ticksCuda << " ticks" << " (" << format_elapsed((double)ticksCuda/(double)freq) << ")" << endl;
    cout << "Compute Slots: " << computeSlots << " (" << blocks << " blocks X " << threads << " threads)" << endl;


    unsigned errorCount = 0;
    for( size_t i = 0; i < numVals; ++i )
    {
        if( gold[i] != hostTest[i] )
        {
            ++errorCount;
        }
    }

    cout << "Number of Errors: " << errorCount << endl;

    return 0;
}

файл: self.cu

#pragma warning( disable : 4231)
#include <windows.h>
#include <cstdlib>
#include <vector>
#include <iostream>
#include <string>
#include <iomanip>
using namespace std;
#include "main.h"

__global__ void SelfNum(int * slots)
{
    __shared__ int N;
    N = (blockIdx.x * blockDim.x) + threadIdx.x;

    const int numDigits = 10;

    __shared__ int digits[numDigits];
    for( int i = 0, temp = N; i < numDigits; ++i, temp /= 10 )
    {
        digits[numDigits-i-1] = temp - 10 * (temp/10)      /*temp % 10*/;
    }

    __shared__ int s;
    s = 0;
    for( int i = 0; i < numDigits; ++i )
        s += digits[i];

    slots[N+s] = 1;
}

__host__ void ComputeSelfNumbers(sized_ptr hostMem, sized_ptr deviceMem, const unsigned  blocks, const unsigned threads)
{
    LARGE_INTEGER li = {0};
    QueryPerformanceFrequency(&li);
    double freq = (double)li.QuadPart;

    LARGE_INTEGER liStart = {0};
    QueryPerformanceCounter(&liStart);

    // run the kernel
    SelfNum<<<blocks, threads>>>(deviceMem.first);
    LARGE_INTEGER liKernel = {0};
    QueryPerformanceCounter(&liKernel);

    cudaMemcpy(hostMem.first, deviceMem.first, hostMem.second*sizeof(int), cudaMemcpyDeviceToHost); // dont copy the overflow - just throw it away
    LARGE_INTEGER liMemcpy = {0};
    QueryPerformanceCounter(&liMemcpy);

    // display performance stats
    double e = double(liMemcpy.QuadPart - liStart.QuadPart)/freq,
        eKernel = double(liKernel.QuadPart - liStart.QuadPart)/freq,
        eMemcpy = double(liMemcpy.QuadPart - liKernel.QuadPart)/freq;

    double pKernel = eKernel/e,
        pMemcpy = eMemcpy/e;

    cout << "Kernel Executed in " << format_elapsed(e) << " -- Breakdown: " << endl
        << "  [kernel] : " << format_elapsed(eKernel) << " (" << format_pct(pKernel) << "%)" << endl
        << "  [memcpy] : " << format_elapsed(eMemcpy) << " (" << format_pct(pMemcpy) << "%)" << endl;



}

UPDATE2:

Я реорганизовал свою реализацию CUDA, чтобы попытаться немного ускорить ее. Я сделал это, развернув петли вручную, установив некоторое сомнительное использование памяти __shared__, которая могла быть ошибкой, и избавиться от некоторой избыточности.

Выход моего нового ядра:

Reference Implementation Ran In 69610 ticks (5 ms)
Kernel Executed in 2 ms -- Breakdown:
  [kernel] : 39 us (1.57%)
  [memcpy] : 2 ms (98.43%)
CUDA Implementation Ran In 62970 ticks (4 ms)
Compute Slots: 1000448 (1954 blocks X 512 threads)
Number of Errors: 0

Единственный код, который я изменил, - это ядро, так что все, что я опубликую здесь:

__global__ void SelfNum(int * slots)
{
    int N = (blockIdx.x * blockDim.x) + threadIdx.x;

    int s = 0;

    int temp = N;
    s += temp - 10 * (temp/10)      /*temp % 10*/;
    s += temp - 10 * ((temp/=10)/10)      /*temp % 10*/;
    s += temp - 10 * ((temp/=10)/10)      /*temp % 10*/;
    s += temp - 10 * ((temp/=10)/10)      /*temp % 10*/;
    s += temp - 10 * ((temp/=10)/10)      /*temp % 10*/;
    s += temp - 10 * ((temp/=10)/10)      /*temp % 10*/;
    s += temp - 10 * ((temp/=10)/10)      /*temp % 10*/;
    s += temp - 10 * ((temp/=10)/10)      /*temp % 10*/;
    s += temp - 10 * ((temp/=10)/10)      /*temp % 10*/;
    s += temp - 10 * ((temp/=10)/10)      /*temp % 10*/;

    slots[N+s] = 1;
}

Ответ 12

Интересно, поможет ли многопоточность. Этот алгоритм выглядит так, что он хорошо поддается многопоточности. (Плохое испытание: создайте две копии программы и запустите их одновременно. Если она работает менее чем в 200% случаев, многопоточность может помочь).

Ответ 13

Меня действительно удивило, что код ниже был быстрее, чем любой другой, размещенный здесь. Я, вероятно, измерил это неправильно, но, возможно, это помогает; или, по крайней мере, интересно.

#include <iostream>
#include <boost/progress.hpp>

class SelfCalc
{
private:
    bool    array[1000000];
    int     non_self;

public:
    SelfCalc()
    {
        memset(&array, 0, sizeof(array));
    }

    void operator()(const int i)
    {
        if (!(array[i]))
            std::cout << i << '\n';

        non_self = i + (i%10) + (i/10)%10 + (i/100)%10 + (i/1000)%10 + (i/10000)%10 +(i/100000)%10;
        array[non_self] = true;
    }
};

class IntIterator
{
private:
    int value;

public: 
    IntIterator(const int _value):value(_value){}

    int operator*(){ return value; } 
    bool operator!=(const IntIterator &v){ return value != v.value; }
    int operator++(){ return ++value; }
};

int main()
{
    boost::progress_timer t;

    SelfCalc selfCalc;
    IntIterator i(1), end(100000);

    std::for_each(i, end, selfCalc);

    std::cout << 100000 << std::endl;
    return 0;
}

Ответ 14

Забавная проблема. В заявленной задаче не указывается, в какой базе она должна быть. Я поиграл с ней и написал версию base-2. Он генерирует дополнительные несколько тысяч записей, потому что точка завершения 1 000 000 не так естественна с базой-2. Это предварительно подсчитывает количество бит в байте для поиска в таблице. Генерация набора результатов (без ввода-вывода) заняла 2,4 мс.

Одна интересная вещь (при условии, что я написал ее правильно) заключается в том, что версия base-2 имеет около 250 000 "номеров собственных номеров" до 1 000 000, в то время как в этом диапазоне имеется чуть меньше 100 000 базовых чисел базового 10.

#include <windows.h>
#include <stdio.h>
#include <string.h>

void StartTimer( _int64 *pt1 )
{
   QueryPerformanceCounter( (LARGE_INTEGER*)pt1 );
}

double StopTimer( _int64 t1 )
{
   _int64 t2, ldFreq;

   QueryPerformanceCounter( (LARGE_INTEGER*)&t2 );
   QueryPerformanceFrequency( (LARGE_INTEGER*)&ldFreq );
   return ((double)( t2 - t1 ) / (double)ldFreq) * 1000.0;
}

#define RANGE 1000000

char sn[0x100000 + 32];

int bitCount[256];

 // precompute bitcounts for each byte
void PreCountBits()
{
    int i;

    // generate count of bits in each byte
    memset( bitCount, 0, sizeof( bitCount ));
    for ( i = 0; i < 256; i++ )
        {
        int tmp = i;

        while ( tmp )
            {
            if ( tmp & 0x01 )
                bitCount[i]++;
            tmp >>= 1;
            }
        }
}


void GenBase2( )
{
    int i;
    int *b1, *b2, *b3;
    int b1sum, b2sum, b3sum;

    i = 0;
    for ( b1 = bitCount; b1 < bitCount + 256; b1++ )
        {
        b1sum = *b1;
        for ( b2 = bitCount; b2 < bitCount + 256; b2++ )
            {
            b2sum = b1sum + *b2;
            for ( b3 = bitCount; b3 < bitCount + 256; b3++ )
                {
                sn[i++ + *b3 + b2sum] = 1;
                }
            }

        // 1000000 does not provide a great termination number for base 2.  So check
        // here.  Overshoots the target some but avoids repeated checks
        if ( i > RANGE )
            return;
        }
}


int main( int argc, char* argv[] )
{
    int i = 0;
    __int64 t1;


    memset( sn, 0, sizeof( sn ));
    StartTimer( &t1 );
    PreCountBits();
    GenBase2();
    printf( "Generation time = %.3f\n", StopTimer( t1 ));

    #if 1
    for ( i = 1; i <= RANGE; i++ )
        if ( !sn[i] ) printf( "%d\n", i );
    #endif
    return 0;
}

Ответ 15

Может быть, попробуйте просто вычислить рекуррентное отношение, определенное ниже?

http://en.wikipedia.org/wiki/Self_number