Грубая сила, однопоточная простая факторизация

Ответ 1

Сравните такой подход к (предварительно сгенерированному) ситу. Modulo стоит дорого, поэтому оба подхода существенно выполняют две вещи: генерируют потенциальные факторы и выполняют операции по модулю. Любая программа должна разумно генерировать новый фактор-кандидат в меньших циклах, чем по модулю, поэтому любая программа связана по модулю.

Данный подход отфильтровывает постоянную долю всех целых чисел, а именно кратные 2 и 3, или 75%. Каждый четвертый (заданный) номер используется в качестве аргумента для оператора modulo. Я назову его фильтром пропуска.

С другой стороны, сито использует только простые числа в качестве аргументов для оператора modulo, а средняя разница между последовательными штрихами регулируется теорема о простом числе - 1/ln (N). Например, e ^ 20 составляет чуть менее 500 миллионов, поэтому число более 500 миллионов человек имеет 5% -ную вероятность быть простым. Если рассматривать все числа до 2 ^ 32, 5% - хорошее эмпирическое правило.

Следовательно, сито будет тратить в 5 раз меньше времени на операции div в качестве вашего пропускающего фильтра. Следующим фактором, который следует учитывать, является скорость, с которой сито создает простые числа, то есть считывает их из памяти или диска. Если выборка одного шага быстрее, чем 4 div s, то сито быстрее. По моим таблицам div пропускная способность моего Core2 не более одного на 12 циклов. Это будут проблемы с жестким разделением, поэтому пусть консервативно бюджет составляет 50 циклов на рассылку. Для процессора 2,5 ГГц это 20 наносекунд.

Через 20 нс жесткий диск емкостью 50 МБ/с может читать около одного байта. Простое решение состоит в том, чтобы использовать 4 байта на рассылку, поэтому диск будет медленнее. Но мы можем быть более умными. Если мы хотим закодировать все простые числа в порядке, мы можем просто кодировать их различия. Опять же, ожидаемая разница равна 1/ln (N). Кроме того, они все четные, что экономит дополнительный бит. И они никогда не равны нулю, что делает расширение для многобайтовой кодировки бесплатным. Таким образом, используя один байт на рассылку, различия до 512 могут храниться в одном байте, что дает нам до 303371455241 в соответствии с

Производительность на MacBook Pro 2,2 ГГц:

dkrauss$ time ./multibyte_sieve g
4294967291

real    2m8.845s
user    1m15.177s
sys    0m2.446s
dkrauss$ time ./multibyte_sieve 18446743721522234449
4294967231 4294967279 

real    0m5.405s
user    0m4.773s
sys 0m0.458s
dkrauss$ time ./mike 18446743721522234449
4294967231 4294967279
real    0m25.147s
user    0m24.170s
sys 0m0.096s

Ответ 2

Мой другой ответ довольно длинный и сильно отличается от этого, так что здесь что-то еще.

Вместо того, чтобы просто отфильтровать кратность первых двух простых чисел или закодировать все соответствующие простые числа по одному байту, эта программа отфильтровывает кратность всех простых чисел, которые вписываются в восемь бит, в частности от 2 до 211. Поэтому вместо передачи 33% от числа, это переходит на 10% к оператору деления.

Простые числа хранятся в трех регистрах SSE, а их модули с работающим счетчиком хранятся в трех других. Если модуль любого простого с счетчиком равен нулю, счетчик не может быть простым. Кроме того, если какой-либо модуль равен единице, то (счетчик + 2) не может быть простым и т.д. До (счетчик + 30). Даже числа игнорируются, а смещения, такие как +3, +6 и +5, пропускаются. Векторная обработка позволяет одновременно обновлять шестнадцатибайтные переменные.

После того, как вы выбрали кухонную раковину с полной оптимизацией (но не более специфичной для платформы, чем встроенная директива), я получил повышение производительности 1,78x (24 с против 13,4 с на моем ноутбуке). Если вы используете библиотеку bignum (даже очень быструю), преимущество больше. Ниже приведена более читаемая версия предварительной оптимизации.

/*
 *  factorize_sse.cpp
 *  Filter out multiples of the first 47 primes while factorizing a number.
 *
 *  Created by David Krauss on 10/14/10.
 *
 */

#include <cmath>
#include <sstream>
#include <iostream>
#include <xmmintrin.h>
using namespace std;

inline void remove_factor( unsigned long &n, unsigned long factor ) __attribute__((always_inline));
void remove_factor( unsigned long &n, unsigned long factor ) {
    while ( n % factor == 0 ) {
        n /= factor;
        cout << factor << " ";
    }
}

int main( int argc, char *argv[] ) {
    unsigned long n;

    if ( argc != 2
        || ( istringstream( argv[1] ) >> n >> ws ).rdstate() != ios::eofbit ) {
        cerr << "Usage: " << argv[0] << " <number>\n";
        return 1;
    }

    int primes[] = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47,
                     53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127,
                     131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211 };
    for ( int *p = primes; p < primes + sizeof primes/sizeof *primes; ++ p ) {
        remove_factor( n, * p );
    }

    //int histo[8] = {}, total = 0;

    enum { bias = 15 - 128 };
    __m128i const prime1 =       _mm_set_epi8( 21, 21, 21, 22, 22, 26, 26, 17, 19, 23, 29, 31, 37, 41, 43, 47 ),
            prime2 =             _mm_set_epi8( 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127 ),
            prime3 =             _mm_set_epi8( 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211 ),
            vbias = _mm_set1_epi8( bias ),
            v3 = _mm_set1_epi8( 3+bias ), v5 = _mm_set1_epi8( 5+bias ), v6 = _mm_set1_epi8( 6+bias ), v8 = _mm_set1_epi8( 8+bias ),
            v9 = _mm_set1_epi8( 9+bias ), v11 = _mm_set1_epi8( 11+bias ), v14 = _mm_set1_epi8( 14+bias ), v15 = _mm_set1_epi8( 15+bias );
    __m128i mod1 = _mm_add_epi8( _mm_set_epi8(  3, 10, 17,  5, 16,  6, 19,  8,  9, 11, 14, 15, 18, 20, 21, 23 ), vbias ),
            mod2 = _mm_add_epi8( _mm_set_epi8( 26, 29, 30, 33, 35, 36, 39, 41, 44, 48,  50,  51,  53,  54,  56,  63 ), vbias ),
            mod3 = _mm_add_epi8( _mm_set_epi8(  65,  68,  69,  74,  75,  78,  81,  83,  86,  89,  90,  95,  96,  98,  99, 105 ), vbias );

    for ( unsigned long factor = 1, limit = sqrtl( n ); factor <= limit + 30; factor += 30 ) {
        if ( n == 1 ) goto done;

        // up to 2^32, distribution of number candidates produced (0 up to 7) is
        // 0.010841     0.0785208   0.222928    0.31905     0.246109    0.101023    0.0200728   0.00145546 
        unsigned candidates[8], *cand_pen = candidates;
        * cand_pen = 6;
        cand_pen += !( _mm_movemask_epi8( _mm_cmpeq_epi8( mod1,  v3 ) ) | _mm_movemask_epi8( _mm_or_si128( _mm_cmpeq_epi8( mod2,  v3 ), _mm_cmpeq_epi8( mod3,  v3 ) ) ) );
        * cand_pen = 10;                                                                                                                                            
        cand_pen += !( _mm_movemask_epi8( _mm_cmpeq_epi8( mod1,  v5 ) ) | _mm_movemask_epi8( _mm_or_si128( _mm_cmpeq_epi8( mod2,  v5 ), _mm_cmpeq_epi8( mod3,  v5 ) ) ) );
        * cand_pen = 12;                                                                                                                            
        cand_pen += !( _mm_movemask_epi8( _mm_cmpeq_epi8( mod1,  v6 ) ) | _mm_movemask_epi8( _mm_or_si128( _mm_cmpeq_epi8( mod2,  v6 ), _mm_cmpeq_epi8( mod3,  v6 ) ) ) );
        * cand_pen = 16;                                                                                                                            
        cand_pen += !( _mm_movemask_epi8( _mm_cmpeq_epi8( mod1,  v8 ) ) | _mm_movemask_epi8( _mm_or_si128( _mm_cmpeq_epi8( mod2,  v8 ), _mm_cmpeq_epi8( mod3,  v8 ) ) ) );
        * cand_pen = 18;                                                                                                                            
        cand_pen += !( _mm_movemask_epi8( _mm_cmpeq_epi8( mod1,  v9 ) ) | _mm_movemask_epi8( _mm_or_si128( _mm_cmpeq_epi8( mod2,  v9 ), _mm_cmpeq_epi8( mod3,  v9 ) ) ) );
        * cand_pen = 22;                                                                                                                            
        cand_pen += !( _mm_movemask_epi8( _mm_cmpeq_epi8( mod1, v11 ) ) | _mm_movemask_epi8( _mm_or_si128( _mm_cmpeq_epi8( mod2, v11 ), _mm_cmpeq_epi8( mod3, v11 ) ) ) );
        * cand_pen = 28;                                                                                                                            
        cand_pen += !( _mm_movemask_epi8( _mm_cmpeq_epi8( mod1, v14 ) ) | _mm_movemask_epi8( _mm_or_si128( _mm_cmpeq_epi8( mod2, v14 ), _mm_cmpeq_epi8( mod3, v14 ) ) ) );
        * cand_pen = 30;                                                                                                                            
        cand_pen += !( _mm_movemask_epi8( _mm_cmpeq_epi8( mod1, v15 ) ) | _mm_movemask_epi8( _mm_or_si128( _mm_cmpeq_epi8( mod2, v15 ), _mm_cmpeq_epi8( mod3, v15 ) ) ) );

        /*++ total;
        ++ histo[ cand_pen - candidates ];

        cout << "( ";
        while ( cand_pen != candidates ) cout << factor + * -- cand_pen << " ";
        cout << ")" << endl; */

        mod1 = _mm_sub_epi8( mod1, _mm_set1_epi8( 15 ) ); // update residuals
        __m128i mask1 = _mm_cmplt_epi8( mod1, _mm_set1_epi8( 1+bias ) );
        mask1 = _mm_and_si128( mask1, prime1 ); // residual goes to zero or negative?
        mod1 = _mm_add_epi8( mask1, mod1 ); // combine reset into zero or negative

        mod2 = _mm_sub_epi8( mod2, _mm_set1_epi8( 15 ) );
        __m128i mask2 = _mm_cmplt_epi8( mod2, _mm_set1_epi8( 1+bias ) );
        mask2 = _mm_and_si128( mask2, prime2 );
        mod2 = _mm_add_epi8( mask2, mod2 );

        mod3 = _mm_sub_epi8( mod3, _mm_set1_epi8( 15 ) );
        __m128i mask3 = _mm_cmplt_epi8( mod3, _mm_set1_epi8( 1+bias ) );
        mask3 = _mm_and_si128( mask3, prime3 );
        mod3 = _mm_add_epi8( mask3, mod3 );

        if ( cand_pen - candidates == 0 ) continue;
        remove_factor( n, factor + candidates[ 0 ] );
        if ( cand_pen - candidates == 1 ) continue;
        remove_factor( n, factor + candidates[ 1 ] );
        if ( cand_pen - candidates == 2 ) continue;
        remove_factor( n, factor + candidates[ 2 ] );
        if ( cand_pen - candidates == 3 ) continue;
        remove_factor( n, factor + candidates[ 3 ] );
        if ( cand_pen - candidates == 4 ) continue;
        remove_factor( n, factor + candidates[ 4 ] );
        if ( cand_pen - candidates == 5 ) continue;
        remove_factor( n, factor + candidates[ 5 ] );
        if ( cand_pen - candidates == 6 ) continue;
        remove_factor( n, factor + candidates[ 6 ] );
    }

    cout << n;
done:
    /*cout << endl;
    for ( int hx = 0; hx < 8; ++ hx ) cout << (float) histo[hx] / total << " ";*/
    cout << endl;
}

.

dkrauss$ /usr/local/bin/g++ main.cpp -o factorize_sse -O3 --profile-use
dkrauss$ time ./factorize_sse 18446743721522234449
4294967231 4294967279 

real    0m13.437s
user    0m13.393s
sys 0m0.011s

Ниже приведен первый черновик. Оптимизация включена

• Сделать циклический счетчик безусловным слиянием (исключить ветку).
• Получите ILP, разворачивая цикл в 15 раз, увеличивая шаг до 30.
  • Вдохновленный вашей оптимизацией.
  • 30 кажется сладким пятном, так как он бесплатно удаляет множественные числа 2, 3 и 5.
  • Прерывания между 5 и 15 могут иметь несколько кратных в один шаг, поэтому поместите несколько копий в различной фазе в вектор.
• Фактор out remove_factor.
• Измените условные, непредсказуемые вызовы remove_factor на нерасширяющиеся записи массивов.
• Полностью разверните окончательный цикл с вызовом remove_factor и убедитесь, что функция всегда включена.
  • Устранить окончательную развернутую итерацию, так как среди кандидатов всегда много.
• Добавьте еще один вектор, содержащий все оставшиеся простые числа, которые достаточно малы.
• Сделайте больше места, добавив смещение к счетчикам и добавьте еще один вектор. Теперь есть только шесть простых чисел, которые могут быть отфильтрованы без столкновения до 16 бит, и у меня также закончились регистры: цикл требует 3 простых вектора, 3 вектора модуля, 8 констант для поиска и одну константу каждый для увеличения и для проверки диапазона. Это составляет 16.
  • В этом приложении коэффициент усиления минимален (но положителен), но первоначальной целью этого метода было отфильтровать простые числа для сита в другом ответе. Так что следите за обновлениями...

Считываемая версия:

/*
 *  factorize_sse.cpp
 *  Filter out multiples of the first 17 primes while factorizing a number.
 *
 *  Created by David Krauss on 10/14/10.
 *
 */

#include <cmath>
#include <sstream>
#include <iostream>
#include <xmmintrin.h>
using namespace std;

int main( int argc, char *argv[] ) {
    unsigned long n;

    if ( argc != 2
        || ( istringstream( argv[1] ) >> n >> ws ).rdstate() != ios::eofbit ) {
        cerr << "Usage: " << argv[0] << " <number>\n";
        return 1;
    }

    int primes[] = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59 };
    for ( int *p = primes; p < primes + sizeof primes/sizeof *primes; ++ p ) {
        while ( n % * p == 0 ) {
            n /= * p;
            cout << * p << " ";
        }
    }

    if ( n != 1 ) {
        __m128i       mod   = _mm_set_epi8( 1, 2, 3,  5,  6,  8,  9, 11, 14, 15, 18, 20, 21, 23, 26, 29 );
        __m128i const prime = _mm_set_epi8( 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59 ),
                      one = _mm_set1_epi8( 1 );

        for ( unsigned long factor = 1, limit = sqrtl( n ); factor < limit; ) {
            factor += 2;
            __m128i mask = _mm_cmpeq_epi8( mod, one ); // residual going to zero?
            mod = _mm_sub_epi8( mod, one ); // update other residuals
            if ( _mm_movemask_epi8( mask ) ) {
                mask = _mm_and_si128( mask, prime ); // reset cycle if going to zero
                mod = _mm_or_si128( mask, mod ); // combine reset into zeroed position

            } else while ( n % factor == 0 ) {
                n /= factor;
                cout << factor << " ";
                if ( n == 1 ) goto done;
            }
        }
        cout << n;
    }
done:
    cout << endl;
}

Ответ 3

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

... без использования вероятностного подхода (например, Миллера-Рабина) для определения первичности

При равномерном распределении 75% ваших входов понадобится миллиард циклов итераций, поэтому сначала стоит потратить миллион операций на менее детерминированные методы, даже если вы получите неубедительный ответ и вернетесь к пробному делению.

Я обнаружил, что метод Брентта по методу Полларда Ро очень хорош, хотя более сложный код и понимание. Лучший пример, который я видел, - это обсуждение форума. Метод полагается на удачу, но помогает достаточно часто, чтобы быть полезной.

Тест на примитивность Миллера-Рабина на самом деле детерминирован до 10 ^ 15, что может избавить вас от бесполезного поиска.

Я попробовал несколько десятков вариаций и остановился на следующем для факторинга значений int64:

• Судебное деление на небольшие факторы. (Я использую первые 8000 предварительно вычисленных простых чисел.)
• 10 попыток с Поллардом Ро, каждый из которых использует 16 итераций
• Пробное деление на sqrt (n).

Обратите внимание, что Поллард Ро обнаруживает факторы, которые не обязательно являются первичными, поэтому рекурсия может использоваться для их определения.

Ответ 4

Этот код довольно медленный, и я уверен, что понимаю, почему. Это не невероятно медленнее, но определенно медленнее в диапазоне 10-20%. Разделение не должно выполняться один раз для каждого цикла, но единственный способ сделать это - вызвать sqrt или что-то подобное.

// system specific typedef for ulong should go here (or use boost::uint64_t)
typedef std::vector<ulong> ULongVector;

void GetFactors(ULongVector &factors, ulong num)
{
  if (num<2)
    return;

  ulong workingNum=num;
  ulong nextOffset=2;

  while (workingNum%2==0)
  {
    factors.push_back(2);
    workingNum/=2;
  }

  while (workingNum%3==0)
  {
    factors.push_back(3);
    workingNum/=3;
  }

  ulong n = 5;
  while ((workingNum != 1) && ((workingNum / n) >= n)) {
    // Is workingNum divisible by n?
    if (workingNum%n==0)
    {
      // n is a factor!
      // so is the number multiplied by n to get workingNum

      // Insert n into the list of factors
      factors.push_back(n);

      // Divide working number by n
      workingNum/=n;
    } else {
      n+=nextOffset;
      nextOffset=(nextOffset==2UL ? 4UL : 2UL);
    }
  }

  if (workingNum != 1) {
    // workingNum is prime, add it to the list of factors        
    factors.push_back(workingNum);
  }
}

Ответ 5

Включение некоторых идей из Omnifarious и других улучшений:

// system specific typedef for ulong should go here (or use boost::uint64_t)
typedef unsigned __int64 ulong;
typedef std::vector<ulong> ULongVector;

// Caller needs to pass in an empty factors vector
void GetFactors(ULongVector &factors, ulong num)
{
  if (num<2)
    return;

  ulong workingNum=num;

  // Factor out factors of 2
  while (workingNum%2==0)
  {
    factors.push_back(2);
    workingNum/=2;
  }

  // Factor out factors of 3
  while (workingNum%3==0)
  {
    factors.push_back(3);
    workingNum/=3;
  }

  if (workingNum==1)
    return;

  // Factor out factors >=5
  ulong nextOffset=2;
  char nextShift = 1;
  ulong n = 5;
  ulong nn = 25;
  do {
    // Is workingNum divisible by n?
    if (workingNum%n==0)
    {
      // n is a factor!
      // so is the number multiplied by n to get workingNum

      // Insert n into the list of factors
      factors.push_back(n);

      // Divide working number by n
      workingNum/=n;

      // Test for done...
      if (workingNum==1)
        return;

      // Try n again
    }  
    else {
      nn += (n << (nextShift+1)) + (1<<(nextShift*2)); // (n+b)^2 = n^2 + 2*n*b + b*2
      n += nextOffset;
      nextOffset ^= 6;
      nextShift ^= 3;
      // invariant: nn == n*n
      if (n & 0x100000000LL) break; // careful of integer wraparound in n^2
    }
  } while (nn <= workingNum);

  // workingNum is prime, add it to the list of factors        
  factors.push_back(workingNum);
}

Ответ 6

Естественное обобщение состоит в том, чтобы предкомпрометировать пропуски с использованием более известных простых чисел, чем только 2 и 3. Как 2, 3, 5, 7, 11, для периода шаблона 2310 (да, славное число). И, возможно, больше, но он имеет уменьшенные прибыли - график времени выполнения может точно определить, где прекомпция начинает оказывать отрицательное воздействие, но, конечно, она зависит от количества чисел, которые необходимо учитывать...

Да, я оставляю детали кодирования для вас.: -)

Приветствия и hth.,

- Alf

Ответ 7

Я не уверен, насколько они эффективны, но вместо

while (workingNum%2==0)

вы могли бы сделать

while (workingNum & 1 == 0)

Я не уверен, что gcc или msvc (или любой другой компилятор, который вы используете) достаточно умны, чтобы изменить выражение workingNum% 2, но вероятность того, что он выполняет разделение и смотрит на модуль в edx...

Мое следующее предложение может быть совершенно ненужным в зависимости от вашего компилятора, но вы можете попробовать поставить workNum/= 3 перед вызовом метода. g++ может быть достаточно умным, чтобы увидеть ненужное разделение и просто использовать частное в eax (вы можете сделать это и внутри своей большой петли). Или более тщательный (но болезненный) подход заключался бы в сборке следующего кода.

while (workingNum%3==0)
{
  factors.push_back(3);
  workingNum/=3;
}

Компилятор, вероятно, переводит модульную операцию в divison, а затем смотрит на модуль в edx. Проблема в том, что вы снова выполняете разделение (и я сомневаюсь, что компилятор видит, что вы просто неявно выполняли деление в состоянии цикла). Итак, вы можете встроить это. Это представляет две проблемы:

1) Метод вызывает push_back (3). Это может испортить регистры, что делает это совершенно ненужным.

2) Получение регистра для работыNum, но это можно определить, выполнив начальную модульную проверку (чтобы заставить ее в% eax), или в текущий момент она будет/должна быть в eax.

Вы можете написать цикл как (предполагая, что workNum находится в eax, и это 32-битный AT & T-синтаксис, только потому, что я не знаю 64-битную сборку или синтаксис Intel)

asm( "
     movl     $3, %ebx
  WorkNumMod3Loop: 
     movl     %eax, %ecx    # just to be safe, backup workingNUm
     movl     $0, %edx      # zero out edx
     idivl    $3            # divide by 3. quotient in eax, remainder in edx
     cmpl     $0, %edx      # compare if it 0
     jne      AfterNumMod3Loop    # if 0 is the remainder, jump out

     # no need to perform division because new workingNum is already in eax
     #factors.push_back(3) call

     je       WorkNumMod3Loop
  AfterNumMod3Loop: 
     movl     %ecx, %eax"
);

Вы должны посмотреть на вывод сборки для этих циклов. Возможно, ваш компилятор уже делает эти оптимизации, но я в этом сомневаюсь. Я бы не удивился, если бы поместить рабочийNum/= n прежде, чем вызов метода немного улучшит производительность в некоторых случаях.