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

Вы передаете полную факторизацию 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 бит).

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

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

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

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

Вы пытались передать массив простых чисел вместо того, чтобы разбивать их на 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 s и 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, которые вы можете найти в других ответах, но я думаю, что это того стоит попытка. Я считаю, что это будет выигрыш, особенно для чисел с большим количеством малых простых множителей.

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

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

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

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

5) Вы выполняете множество операций по модулю и деления сразу друг за другом. В C это 2 отдельные команды (сначала «/», а затем «%»).Однако в сборке это одна команда: '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.

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

Пробовали ли вы использовать оптимизацию на основе профиля?

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

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

Попробуйте заменить этот шаблон (также для 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 или если некоторые из них встречаются чаще, вы можете получить что-то, специализируя функцию для этих значений.

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

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

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

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

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

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);

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

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

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

2. Вы пробовали предварительно вычислить star (pr [__ builtin_ffsll (bit)]) и vals (pr [__ builtin_ffsll (bit)] - 1) ? Это потребовало бы некоторой простой работы для поиска по массиву, но, возможно, оно того стоит, если таблицы достаточно малы.

3. Не вычисляйте 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;

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

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

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

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

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

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

n * inverse_q

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

Сначала немного придирок ;-) вам следует быть более внимательным к типам, которые вы используете. В некоторых местах вы, похоже, предполагаете, что ulong имеет ширину 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 и тому подобное...