Оптимизируйте меня! (C, производительность) - повторение вопроса о бит-twiddling

Ответ 1

Вы, кажется, тратите много времени на деление по факторам. Гораздо быстрее заменить деление на умножение на обратный делитель (деление: ~ 15-80 (!) циклов, в зависимости от делителя, умножения: ~ 4 цикла), IF, конечно, вы можете прекомпилировать обратные.

Хотя это вряд ли возможно с q, r, s - из-за диапазона этих vars, это очень легко сделать с p, который всегда исходит из небольшого статического массива pr []. Предкоммутируйте обратные числа этих простых чисел и сохраните их в другом массиве. Затем, вместо деления на p, умножьте на обратное, взятое из второго массива. (Или создать один массив структур.)

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

EDIT:

После консультации с Hacker Delight (отличная книга, BTW) по этому вопросу кажется, что вы можете сделать это еще быстрее, используя тот факт, что все деления в вашем коде точны (т.е. остаток равен нулю).

Кажется, что для каждого делителя d, который является нечетным, а база B = 2 ^word_size существует единственный мультипликативный обратный d⃰, который удовлетворяет условиям: d⃰ < B и d·d⃰ ≡ 1 (mod B). Для каждого х, являющегося точным кратным d, это означает x/d ≡ x·d⃰ (mod B). Это означает, что вы можете просто заменить деление на умножение, не добавлять исправления, проверки, проблемы округления, что угодно. (Доказательства этих теорем можно найти в книге.) Примечание, что этот мультипликативный обратный не должен быть равен обратному, как определено предыдущим методом!

Как проверить, является ли данный x точным кратным d - i.e. x mod d = 0? Легко! x mod d = 0 iff x·d⃰ mod B ≤ ⌊(B-1)/d⌋. Обратите внимание, что этот верхний предел можно предварительно вычислить.

Итак, в коде:

unsigned x, d;
unsigned inv_d = mulinv(d);          //precompute this!
unsigned limit = (unsigned)-1 / d;   //precompute this!

unsigned q = x*inv_d;
if(q <= limit)
{
   //x % d == 0
   //q == x/d
} else {
   //x % d != 0
   //q is garbage
}

Предполагая, что массив pr[] станет массивом struct prime:

struct prime {
   ulong p;
   ulong inv_p;  //equal to mulinv(p)
   ulong limit;  //equal to (ulong)-1 / p
}

цикл while(smallprimes) в вашем коде становится:

while (smallprimes) {
    ulong bit = smallprimes & (-smallprimes);
    int bit_ix = __builtin_ffsll(bit);
    ulong p = pr[bit_ix].p;
    ulong inv_p = pr[bit_ix].inv_p;
    ulong limit = pr[bit_ix].limit;
    nu = minuu(nu, vals(p - 1));
    prod *= ugcd(nn, star(p));
    n *= inv_p;
    for(;;) {
        ulong q = n * inv_p;
        if (q > limit)
            break;
        n = q;
    }
    smallprimes ^= bit;
}

И для функции mulinv():

ulong mulinv(ulong d) //d needs to be odd
{
   ulong x = d;
   for(;;)
   {
      ulong tmp = d * x;
      if(tmp == 1)
         return x;
      x *= 2 - tmp;
   }
}

Обратите внимание, что вы можете заменить ulong на любой другой неподписанный тип - используйте один и тот же тип последовательно.

Доказательства, почему s и как s доступны в книге. Сердечно рекомендуется прочитать: -).

Ответ 2

Если ваш компилятор поддерживает атрибуты функции GCC, вы можете пометить свои чистые функции этим атрибутом:

ulong star(ulong n) __attribute__ ((const));

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

Есть ли причина, по которой вы использовали opencoded vals() вместо использования __builtin_ctz()?

Ответ 3

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

Если вы действительно ищете целые числа, которые максимизируют количество не свидетелей для теста MR (т.е. oeis.org/classic/A141768, о котором вы упоминаете), тогда можно было бы использовать, что количество не свидетелей не может быть больше, чем phi (n)/4, и что целые числа, у которых есть это много не свидетелей, либо являются произведением двух простых чисел вида

(к + 1) * (2k + 1)

или они являются числами Кармайчеля с тремя основными факторами. Я думаю, что выше некоторого предела все целые числа в последовательности имеют эту форму и что это можно проверить, доказав верхнюю границу для свидетелей всех других целых чисел. Например. целые числа с 4 или более факторами всегда имеют не более phi (n)/8 свидетелей. Аналогичные результаты могут быть получены из формулы для числа оснований для других целых чисел.

Что касается микрооптимизаций: всякий раз, когда вы знаете, что целое число делится на какое-то частное, тогда можно заменить деление на умножение на обратное к фактору по модулю 2 ^ 64. И тесты n% q == 0 можно заменить тестом

n * inverse_q < max_q,

где inverse_q = q ^ (- 1) mod 2 ^ 64 и max_q = 2 ^ 64/q. Очевидно, что inverse_q и max_q должны быть предварительно вычислены, чтобы быть эффективными, но поскольку вы используете сито, я предполагаю, что это не должно быть препятствием.

Ответ 4

Небольшая оптимизация, но:

ulong f;
ulong nn;
ulong nu = 0xFFFF;  // "Infinity" for the purpose of minimum
ulong prod = 1;

if (!smallprimes & !q)
    return 0;

// no need to do this operations before because of the previous return
f = __builtin_popcountll(smallprimes) + (q > 1) + (r > 1) + (s > 1);
nn = star(n);

BTW: вы должны отредактировать сообщение, чтобы добавить star() и другие функции, которые вы используете для определения

Ответ 5

Попробуйте заменить этот шаблон (для r и q тоже):

n /= p; 
while (n % p == 0) 
  n /= p; 

При этом:

ulong m;
  ...
m = n / p; 
do { 
  n = m; 
  m = n / p; 
} while ( m * p == n); 

В моих ограниченных тестах я получил небольшое ускорение (10%) от устранения модуля.

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

Ответ 6

1) Я бы заставил компилятор выплюнуть собранную им сборку и попытаться вывести, если то, что она делает, это лучшее, что она может сделать... и если вы заметите проблемы, измените код, чтобы сборка выглядела лучше. Таким образом, вы также можете убедиться, что функции, на которые вы надеетесь, будут встроены (например, звезда и vals). (Возможно, вам придется добавить прагму или даже превратить их в макросы)

2) Замечательно, что вы пытаетесь это сделать на многоядерной машине, но этот цикл является однопроходным. Я предполагаю, что есть зонтичные функции, которые разбивают нагрузку на несколько потоков, чтобы использовать больше ядер?

3) Невозможно предложить ускорение, если то, что пытается вычислить фактическая функция, неясно. Как правило, наиболее впечатляющие ускорения не достигаются с помощью бит-скручивания, но с изменением алгоритма. Так может быть немного комментариев; ^)

4) Если вы действительно хотите увеличить скорость до 10 * или больше, проверьте CUDA или openCL, который позволяет запускать программы на вашем графическом оборудовании. Он сияет такими функциями!

5) Вы выполняете множество по модулю и делите сразу после друг друга. В C это две отдельные команды (сначала "/", а затем "%" ). Однако в сборке это 1 команда: "DIV" или "IDIV", которая возвращает как остаток, так и коэффициент в один ход:

B.4.75 IDIV: Signed Integer Divide

IDIV r/m8                     ; F6 /7                [8086]
IDIV r/m16                    ; o16 F7 /7            [8086]
IDIV r/m32                    ; o32 F7 /7            [386]

IDIV performs signed integer division. The explicit operand provided is the divisor; the dividend and destination operands are implicit, in the following way:

For IDIV r/m8, AX is divided by the given operand; the quotient is stored in AL and the remainder in AH.

For IDIV r/m16, DX:AX is divided by the given operand; the quotient is stored in AX and the remainder in DX.

For IDIV r/m32, EDX:EAX is divided by the given operand; the quotient is stored in EAX and the remainder in EDX.

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

Ответ 7

• Убедитесь, что ваши функции включены. Если они вне очереди, накладные расходы могут скомпенсироваться, особенно в первом цикле while. Лучший способ быть уверенным в том, чтобы проверить сборку.

• Пробовали ли вы предварительные вычисления star( pr[__builtin_ffsll(bit)] ) и vals( pr[__builtin_ffsll(bit)] - 1)? Это должно было бы проделать некоторую простую работу для поиска массива, но это может стоить того, если таблицы достаточно малы.

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

BASES_END:
ulong addToF = 0;
if (n > 1) {
    nu = minuu(nu, vals(n - 1));
    prod *= ugcd(nn, star(n));
    addToF = 1;
}
// ... early out if nu == 1...
// ... compute f ...
f += addToF;

Надеюсь, что это поможет.

Ответ 8

Сначала немного придирки;-) вам следует более осторожно относиться к типам, которые вы используете. В некоторых местах вам кажется, что улонг имеет ширину 64 бит, используйте uint64_t там. А также для всех других типов тщательно переосмысливайте, что вы ожидаете от них, и используйте соответствующий тип.

Оптимизация, которую я мог видеть, - это целочисленное деление. Ваш код делает это много, это, наверное, самая дорогая вещь, которую вы делаете. Разделение небольших целых чисел (uint32_t) может быть намного более эффективным, чем большие. В частности, для uint32_t существует инструкция ассемблера, которая делит и по модулю за один раз, называется divl.

Если вы используете соответствующие типы, ваш компилятор может сделать это для вас. Но лучше проверить ассемблер (опция -S на gcc), как уже сказал кто-то. В противном случае легко добавить некоторые небольшие фрагменты ассемблера здесь и там. Я нашел что-то вроде этого в моем коде:

register uint32_t a asm("eax") = 0;
register uint32_t ret asm("edx") = 0;
asm("divl %4"
    : "=a" (a), "=d" (ret)
    : "0" (a), "1" (ret), "rm" (divisor));

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

Ответ 9

Вы пытались использовать оптимизацию с помощью профиля?

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

Затем перекомпилируйте и свяжите его с опцией -fprofile-use.

Ответ 10

Пробовал ли вы просмотр таблицы для первого цикла while? Вы можете разделить smallprimes на 4 16-битные значения, посмотреть их вклад и объединить их. Но, возможно, вам нужны побочные эффекты.

Ответ 11

Попробовали ли вы передать массив простых чисел, а не разбивать их на smallprimes, q, r и s? Поскольку я не знаю, что делает внешний код, я, вероятно, ошибаюсь, но есть вероятность, что у вас также есть функция для преобразования некоторых простых чисел в растровое изображение smallprimes, и внутри этой функции вы преобразуете растровое изображение обратно в массив простых чисел, эффективно. Кроме того, вы, похоже, выполняете идентичную обработку для элементов smallprimes, q, r и s. Он должен сэкономить вам крошечный объем обработки за звонок.

Кроме того, вы, похоже, знаете, что переданные в простых числах разделяют n. У вас есть достаточное знание вне пределов власти каждого шага, которое делит n? Вы могли бы сэкономить много времени, если сможете устранить операцию по модулю, передав эту информацию этой функции. Другими словами, если n - pow(p_0,e_0)*pow(p_1,e_1)*...*pow(p_k,e_k)*n_leftover, и если вы знаете больше об этих e_i и n_leftover, то их передача будет означать много вещей, которые вам не нужно делать в этой функции.

Возможно, существует способ обнаружить n_leftover (неактивную часть n) с меньшим количеством операций по модулю, но это только догадка, поэтому вам может понадобиться немного поэкспериментировать с ней. Идея состоит в том, чтобы использовать gcd для удаления известных факторов из n повторно, пока вы не избавитесь от всех известных простых факторов. Позвольте мне дать почти-c-код:

factors=p_0*p_1*...*p_k*q*r*s;
n_leftover=n/factors;
do {
    factors=gcd(n_leftover, factors);
    n_leftover = n_leftover/factors;
} while (factors != 1);

Я вовсе не уверен, что это будет лучше, чем код, который у вас есть, не говоря уже о комбинированных предложениях mod/div, которые вы можете найти в других ответах, но я думаю, что стоит попробовать. Я чувствую, что это будет победа, особенно для чисел с большим количеством небольших простых факторов.

Ответ 12

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

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

edit должен был восстановить код здесь:

#include 
#define SIZE (1024*1024) //must be 2^n
#define MASK (SIZE-1)
typedef struct {
    int p;
    int next;
} p_type;

p_type primes[SIZE];
int sieve[SIZE];

void init_sieve()
{
    int i,n;
    int count = 1;
    primes[1].p = 3;
    sieve[1] = 1;
    for (n=5;SIZE>n;n+=2)
    {
        int flag = 0;
        for (i=1;count>=i;i++)
        {
            if ((n%primes[i].p) == 0)
            {
                flag = 1;
                break;
            }
        }
        if (flag==0)
        {
            count++;
            primes[count].p = n;
            sieve[n>>1] = count;
        }
    }
}

int main()
{
    int ptr,n;
    init_sieve();
    printf("init_done\n");

    // factor odd numbers starting with 3
    for (n=1;1000000000>n;n++)
    {
        ptr = sieve[n&MASK];
        if (ptr == 0) //prime
        {
//            printf("%d is prime",n*2+1);
        }
        else //composite
        {
//            printf ("%d has divisors:",n*2+1);
            while(ptr!=0)
            {
//                printf ("%d ",primes[ptr].p);
                sieve[n&MASK]=primes[ptr].next;
                //move the prime to the next number it divides
                primes[ptr].next = sieve[(n+primes[ptr].p)&MASK];
                sieve[(n+primes[ptr].p)&MASK] = ptr;
                ptr = sieve[n&MASK];
            }
        }
//        printf("\n");
    }
    return 0;
}

Функция init создает базу факторов и инициализирует сито. Это занимает около 13 секунд на моем ноутбуке. Затем все цифры до 1 миллиарда учитываются или определяются как простые в течение еще 25 секунд. Числа меньше SIZE никогда не сообщаются как простые, потому что они имеют 1 фактор в базе факторов, но это может быть изменено.

Идея состоит в том, чтобы поддерживать связанный список для каждой записи в сите. Числа факторизуются, просто вытаскивая их факторы из связанного списка. Когда они вытаскиваются, они вставляются в список для следующего числа, которое будет делиться на это число. Это тоже очень кэш. Размер сита должен быть больше, чем наибольшее простое в факторной базе. Точно так же это сито может достигать 2 ** 40 примерно за 7 часов, что кажется вашей целью (за исключением того, что n должно быть 64 бита).

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

Надеюсь, что это поможет.

Кстати, это первый раз, когда я публиковал этот алгоритм публично.

Ответ 13

просто мысль, но, возможно, использование ваших вариантов оптимизации компиляторов поможет, если вы еще этого не сделали. другая мысль заключалась бы в том, что если деньги не будут проблемой, вы можете использовать компилятор Intel C/С++, предполагая, что вы используете процессор Intel. Я также предполагаю, что другие производители процессоров (AMD и т.д.) Будут иметь похожие компиляторы

Ответ 14

Если вы собираетесь немедленно выйти на (!smallprimes&!q), почему бы не выполнить этот тест перед тем, как вызвать функцию, и сохранить служебные данные вызова функции?

Кроме того, кажется, что у вас фактически есть 3 различные функции, которые являются линейными, за исключением цикла smallprimes.   bases1(s,n,q), bases2(s,n,q,r) и bases3(s,n,q,r,s).

Это может быть победа на самом деле создать их как 3 отдельные функции без ветвей и gotos и вызвать соответствующий:

if (!(smallprimes|q)) { r = 0;}
else if (s) { r = bases3(s,n,q,r,s);}
else if (r) { r = bases2(s,n,q,r); }
else        { r = bases1(s,n,q);

Это было бы наиболее эффективно, если предыдущая обработка уже дала вызывающему коду "знание", какую функцию выполнить, и вам не нужно ее проверять.

Ответ 15

Если деления, которые вы используете, имеют числа, которые известны во время компиляции, но часто используются во время выполнения (делятся на одно и то же число много раз), тогда я бы предложил использовать libdivide библиотека, которая в основном реализует во время выполнения оптимизацию, которую компиляторы выполняют для констант времени компиляции (используя маски сдвигов и т.д.). Это может принести огромную пользу. Кроме того, избегая использования x% y == 0 для чего-то вроде z = x/y, z * y == x как ergosys, предложенного выше, также должно иметь измеримое улучшение.

Ответ 16

Является ли код на вашем верхнем посту оптимизированной версией? Если да, есть еще слишком много операций разделения, которые в значительной степени ожидают циклы процессора.

Этот код чрезмерно перегружает бит

if (!smallprimes & !q)
    return 0;

изменить на логическое и &&

if (!smallprimes && !q)
    return 0;

сделает его коротким, закороченным быстрее, без выделения q

И следующий код

ulong bit = smallprimes & (-smallprimes);
ulong p = pr[__builtin_ffsll(bit)];

который используется для поиска последнего набора бит мелких символов. Почему бы вам не использовать более простой способ

ulong p = pr[__builtin_ctz(smallprimes)];

Другим виновником снижения производительности может быть слишком много разветвлений программ. Вы можете рассмотреть возможность перехода на другие эквивалентные менее ветвящиеся или неконкурентные эквиваленты