Pull to refresh

Comments 116

Разве длинная арифметика не будет работать быстрее двоичного умножения?
Практика показывает, что нет. Разделить одно длинное число на другое — не так уж просто
А если делить с помощью вещественных чисел, а потом корректировать? Умножать длинные числа (в пределах int64*int64) не так уж сложно.
Точно не помню, проверял давно. Но вроде тоже медленнее
Да, очень хорошая идея.

long long mul( long long a, long long b, long long m ) {
  long long q = (long long)((long double)a * (long double)b / (long double)m);
  long long r = a * b - q * m;
  return (r + 5 * m) % m;
}


Данная функция умеет перемножать long long по модулю за O(1). Более подробно на Codeforces.
Если long double у вас 80-битный, то это может сработать. Но Visual Studio долгое время интерпретировало его, как double. И тогда вы потеряете старшие 8-12 бит числа r, поскольку q будет найдено с точностью 52 бита. В общем, это сильно зависит от компилятора.
Да и 5*m вызывает подозрения. Они учитывают случай m>2^62?
(1*1%0x7FFFFFFFFFFFFFFFLL) = 9223372036854775804
так и должно быть?

#include <iostream>

using namespace std;

long long mul( long long a, long long b, long long m ) {
  long long q = (long long)((long double)a * (long double)b / (long double)m);
  long long r = a * b - q * m;
  return (r + 5 * m) % m;
}

int main()
{
   long long a = 1;
   long long b = 1;
   long long m = 0x7FFFFFFFFFFFFFFFLL;
    
   cout << mul(a, b, m);
   
   return 0;
}

На тех компиляторах, на которых sizeof(long double)>8, есть хорошие шансы на успех у такого кода:

long long rmul(long long a,long long b,long long m){
	int s=0;
	if(a>=m/2){ a=m-a; s=1; }
	if(b>=m/2){ b=m-b; s^=1; }

	long long q=(long long)((long double)a*(long double)b/(long double)m+0.5);
	long long c=a*b-q*m;
	if(c>=m || c<=(-1LL)<<62) c-=m;
	else if(c<0) c+=m;
	if(s && c!=0) c=m-c;
	return c;
}

К сожалению, VS к таким не относится. И про «старый x87» он уже не знает, а без него 80-битную арифметику, по-видимому, не получить.
Пытаюсь добиться 80-битных чисел с помощью ассемблерных вставок — всё равно получается ерунда:
bool TestPrecision(long long a,long long b){
// test (a+b-b==a)
	long long c=0;
	__asm {
		        fild a
			fild b
			fadd
			fild b
			fsub
			fistp c
	}
	if(c!=a) printf("Error: a=%I64x, b=%I64x, c=%I64x\n",a,b,c);
}

Работает так, словно регистры х87 реализованы, как 8-байтовые double: уже на a=1, b=1LL<<53 выдаётся ошибка.
Длинную арифметику писать нужно, а двоичное умножение очень и очень простое. Но вообще да, можно и длинную арифметику организовать. Она, скорее всего, будет работать несколько быстрее.
И да, еще замечание: что вы будете делать с числами Карамайкла? На серьезных соревнованиях такие числа всегда включают в тестовые данные.
Вас интересует ответ, или то, знает ли его автор? ;)
Если проанализировать числа Карамайкла, то можно заметить, что они состоят из нескольких простых множителей, например 3, 5, 7, 11… Если выбрать 100 случайных чисел, то, с большой вероятностью, эти множители будут входить в некоторые из выбранных чисел.
Тест Ферма в данном случае пользы не даст, но НОД будет отличен от 1. Значит программа все равно выдаст, что число непростое.
У таких чисел, как 410041 и 252601 наименьший простой делитель будет 41. Вероятность того, что он встретится у какого-нибудь из 100 случайных чисел уже всего 92% — маловато для вероятностного теста. А у бОльших чисел наименьшие делители могут быть и ещё больше.
Количество итераций в тесте Ферма может быть и не 100, а например, 200. Но устранять эту проблему можно и другим способом:
Перед проверкой числа N тестом Ферма, можно перебрать все простые числа от 2 до min(N, 300), а их там будет очень мало, и найти НОД каждого из этих чисел и N. На время работы это сильно не скажется, так как массив простых чисел от 2 до 300 можно посчитать в самом начале.
Откуда взялось число 300?
В олимпиадном программировании многие константы берутся интуицией
В олимпиадном программировании числа, взятые интуицией, приводят к WA Test 35
И хорошо если не WA126, когда это последний тест :'(
Хорошо — это когда после успешного отлова ошибки, приводящей к WA126 не прилетает TL127!
UFO just landed and posted this here
Т.Е. в олимпиадном программировании всего 128 попыток?
нет, просто речь идёт о том, что много тестов прошло, а на последних обломинго. очень обидно.
Нет) Зависит от самой олимпиады. А WA126 — значит, что программа дала неверный ответ на 126 тесте, а на 1-125 — верный.
а плохо — когда после отлова WA126 прилетает TL125/
Это за успешный отлов не считается!
После отлова возможной причины WA126 :)
Ну-ну.
1207252621 = 613 * 919 * 2143
54519328481 = 503 * 5021 * 21587
20618724001 = 701 * 2801 * 10501
8976678481 = 1009 * 2521 * 3529
432210655801 = 3011 * 5419 * 26489
949631589089 = 6917 * 10193 * 13469

Это только среди чисел до триллиона. Дальше наверняка будет ещё хуже.

И где вы планируете остановиться в проверке делимости?
Ну, статья была про числа до 232, так что автор может остановиться, хе-хе, на 613…
Тогда зачем там long long? И двоичное умножение?
long long и двоичное умножение не дают асимптотики, указанной в заголовке статьи
Я думаю, что не все поняли, что данный алгоритм является вероятностным. Чисто математически любой вероятностный алгоритм не верен уже изначально. Тест Ферма может только доказать, что число не является простым.
На реальных олимпиадах этот код работает, я сам его писал.
Значит, это была олимпиада не того уровня. Если Mrrl нашел число 1207252621 за каких-то 15 минут — значит, Станкевич его давно наизусть знает, и вставляет в каждый 35й тест :)
Можете простым перебором убедиться, что данный алгоритм безошибочно работает на всех тестах в пределах типа Int. Для больших чисел есть вероятность ошибки, однако, на олимпиаде может встретиться задача именно на такой алгоритм. Знать его все же нужно.
Какой алгоритм? Первый, который на 410041 падал? Или тот, который перебирает делители до 300 и падает на 1207252621 (это число в int входит, между прочим)? Но командных соревнованиях вы бы сейчас 40 штрафных минут уже получили. Еще по одной задачке такой же фейл — и прощай надежда на третье место.
Мне случалось придумывать тесты, на которых люди и по 60 минусов получали :)
Вероятностный алгоритм хорош, если при любых внешних данных (в данном случае, x) и случайном выборе переменных в алгоритме (a) вероятность ошибки будет не больше чего-то. В данном случае, кармайлово число x повышает вероятность ошибки на отдельно взятом тесте вплоть до 1-3*(x^(-1/3)), что выглядит недопустимо близко к 1 (да ещё и растёт с ростом x).
Я вас удивлю, но любое непростое число состоит из простых множителей
Можно минусы обосновать?
Придираться к придирке — верный путь скатиться во флейм и взаимотроллинг. Отсюда и минус — чтобы к этому «1?» не придрался кто-то еще.
for(long long i=2;i<=sqrt(n);i++)

Даже в первом варианте лучше так
for(long long i=3;i<=sqrt(n); i+=2 )
а проверку на 2 выполнить в начале

В таком случае уже не имеет смысла использовать тест Ферма, если проверять до корня. Иногда на олимпиадах нужна более быстрая проверка, пусть и вероятностная.
UFO just landed and posted this here
да почему, 3-4 раза засабмитил, прошло — ура.
А потом жюри делает, ха-ха, rejudge задачи из-за технических проблем. Или случайно смотрит на решение — и делает rejudge ему индивидуально.

Во всех регламентах соревнований всегда есть фраза — если программа показывает разные результаты от запуска к запуску — жюри имеет право выбрать худший.
Если ответ неверен в одном случае на 2^1000 — можно и рискнуть, Вы не находите? =)
Можно, так и нахожу. Разумеется, если нельзя подобрать тест, который валит решение гарантированно. А вот если вероятность составляет "3-4 раза засабмитил, прошло — ура" — то нет никакой гарантии, что задача случайно не пропадет из числа решенных уже после окончания тура…
А разве нет принципа «плюсы не перепроверяются»?
Нет. Есть принцип — «все посылки после плюса не учитываются», но он затрагивает другие ситуации. А еще есть принцип «жюри всегда право». Не говоря уже о негласном «лучшие тесты всегда пишутся во время соревнований» :)

UPD: забавно, в системе Челябинской системе Polypody была стандартная операция «перепроверить все плюсы», а вот «перепроверить только минусы» — не было
На четвертьфиналах в НГТУ регулярно бывало, что обнаруживали неверный тест. Например, 5.
Тогда перепроверяли все отправленные задачи, в том числе пройденные.
А помнишь, когда greedy алгоритм проходил потому что теста на него не было? Было обидно :(
Ага, помню, я отказался решать какую то задачу, т.к. отсек жадный алгоритм, как непроходящий краевые условия, а аналитический алгоритм был слишком сложен за отведенное время. А какой то парень из НГУ от балды заслал и он прошел, т.к. оказалось, что не было достаточно серьезного теста.
а я тогда помню хотел жадный всё равно сделать. но Лёха нас обоих удебил приведя контрпример сразу же :(
Простейшая модификация функции ferma уменьшит вероятность ошибки примерно до 10^(-60) при любом x. Это уже вполне приемлемо, хотя алгоритм и остаётся вероятностным.
Кроме того, разве на олимпиадах не используют quicksort или nthelement? Они тоже вероятностные :)
Быстрая сортировка, хоть и вероятностная, но, в недетерминированом исполнении, не обладает гарантированно валящими ее исходными данными.
Worst-case expected-time? :D
Мне кажется, mayorovp имеет в виду, что здесь не «работающий вероятностный алгоритм», а «алгоритм, работающий почти на всех тестах». Какой тест мы бы ни придумали, мат ожидание времени работы QS будет хорошим, здесь же не так
это смотря как написать qsort =) можно написать так, что будет такой тест
UFO just landed and posted this here
Превратить недерминированное время в вероятностный ответ довольно просто, только нужна подходящая задача.
Допустим, в задаче ищется существование решения или наилучшее решение. Перебором. qsort (или другой метод с похожими свойствами) вызывается на каждом шаге перебора. Если он всегда успевает за n*log(n), то перебор успеет закончиться, если не всегда — можно вылететь по TL. И следим по таймеру, не кончилось ли время. Если кончается — выдаём ответ, не дождавшись конца перебора.
Как-то так. Не знаю, насколько сейчас практикуются такие приёмы, но при подготовке тура их приходится учитывать.
Такие приемы практикуются в «конкурсных» задачах, где нет правильных ответов, и за каждый тест начисляются баллы согласно близости ответа тестируемой программы к лучшему ответу.
UFO just landed and posted this here
тогда уж
for(long long i=3,j=2;i<=sqrt(n); i+=j,j=6-j)
Можно еще i <= sqrt(n) заменить на i*i <= n.
Умножить два целых быстрее, чем взять корень в вещественных числах.

ну я вообще предполагал что любой дурак сделает
long long sqrtn = sqrt(n)
взять один раз корень быстрее, чем N раз умножать
Сейчас многие считают, что этим дураком должен быть компилятор.
ну если не знают что такое «volatile»…
А, кстати, есть ли в современном C/C++ обозначение для «чистой» функции, результат которой зависит только от аргументов?
Нет, #pragma intrinsic — это обозначение функции, которая реализована компилятором специальным образом.
В этом смысле, вызов sqrt может быть соптимизирован: раз компилятор знает, как это устроено, то он может догадаться, что его значение в данном случае не изменится.
Стоп, а модификатор const у функции разве не оно?
const
Many functions do not examine any values except their arguments, and have no effects except the return value. Basically this is just slightly more strict class than the pure attribute below, since function is not allowed to read global memory.

есть еще pure:
pure
Many functions have no effects except the return value and their return value depends only on the parameters and/or global variables. Such a function can be subject to common subexpression elimination and loop optimization just as an arithmetic operator would be. These functions should be declared with the attribute pure.
Да, const похоже. Осталось проверить, есть ли оно у математических функций в стандартных math.h.
А смысла pure я не понял. Если оно может зависеть от глобальных переменных — как его можно вынести из цикла? Вдруг другой поток (или неизвестная функция, вызванная в том же цикле) их поменяет?
Я так понимаю, разница pure и const в том, что pure никуда не пишет кроме результата, а const еще и не читает ничего
Но если в цикле нет ни одной неизвестной не-pure функции — то выносить вычисление из цикла все-таки можно. Что же до других потоков — то их действия не обязаны быть видимы сразу же, если только глобальная переменная не volatile. А такие переменные читать из pure — функции либо нельзя, либо UB.
VS2013 выносит из цикла не только sqrt, но и sin. Допустим, sqrt описан, как intrinsic. Но у sin никаких пометок нет! Откуда он знает, что это const функция?
Во-первых, sqrt тоже не как intrinsic описан. А во-вторых, не пофиг ли компилятору, особенно от M$, забыли или нет указать, что хорошо известная ему функция должна быть реализована им же?..

До кучи — если выключить стандартную библиотеку, то компилятор жалуется на отсутствие не то memset, не то memcpy, но при этом запрещает реализовывать в своем коде.
Во-первых, sqrt тоже не как intrinsic описан.


Он описан так:
_CRT_JIT_INTRINSIC  double  __cdecl sqrt(_In_ double _X);

Возможно, правда, что «JIT_INTRINSIC» означает, что это действует только в managed коде — я не разбирался с её расшифровкой.
www.cdsan.com/Src_VsInc.php?fid=122&ln=35
/* jit64 instrinsic stuff */
#ifndef _CRT_JIT_INTRINSIC
#if defined(_M_CEE) && (defined(_M_AMD64) || defined(_M_IA64))
/* This is only needed when managed code is calling the native APIs, targeting the 64-bit runtime */
#define _CRT_JIT_INTRINSIC __declspec(jitintrinsic)
#else
#define _CRT_JIT_INTRINSIC 
#endif
#endif
</code>
По крайней мере, Станкевич, Лопатин и иже с ними никогда не были против
Он на всех числах работает. Менее 25% ошибки на 1 тест для любого x.
Я, конечно, понимаю, что здесь собрался кружок по спортивному программированию, но, блин, укажите нормальную асимптотику в заголовке.
AKS проверяет числа за O(L^{6 + eps}), Рабин-Миллера и Ферма можно написать за O(k * L^2 \log L), где L — длина числа.

Вы же не будете составлять RSA-ключ из чисел, влезающих в long long?
Я не буду, но за остальных я бы не стал ручаться… А реальная асимптотика тут O(k L3)
А почему не O(k L2 log(L) log(log(L))?

(у алгоритма, реализованного через двоичное умножение, вообще получается O(k L4) — он же вычисляет остаток от деления на P после каждого сдвига. К счастью, деление там можно заменить вычитанием, и вернуть потерянный порядок).
Как — почему? Возведение в степень — это log N умножений, каждое умножение — это L2 операций. L = log N. Отсюда и L3. Откуда вообще в подобном алгоритме, как описанный в этой статье, может появиться log L?
Если L будет большим, то и умножение будет браться из длинной арифметики. А там оно уже быстрое (если подобрать подходящую библиотеку или язык, в котором длинная арифметика входит в стандартные библиотеки).
Быстрое длинное умножение — это уже оптимизация, не подразумеваемая явно… Хотя, возможно, я просто слишком много раз писал длинное умножение без сторонних библиотек за ограниченное время.

Да, с быстрым длинным умножением асимптотика выходит другая, согласен.

PS УПС! А ведь там в алгоритме еще требуется длинное деление по модулю… Надеюсь, у него тоже существуют быстрые версии? Иначе асимптотика падает до O(k L4)…
Gcd можно искать через двоичного Евклида. Это O(L^2).
А так, конечно, есть. Например, метод касательных (Newton-Raphson method). Мне кажется, он работает за O(\log L M(L)), где M(L) — время умножения двух чисел.
Если через длинную арифметику (а не «двоичным умножением»), то оно делается один раз на умножение, то есть, остаётся то же O(k L3) (при классической реализации). Поскольку делим мы всегда на одно и то же число, то можем заранее посчитать 2^(2L)/x, и искать остаток за два умножения (умножили на обратный — взяли старшее «слово» — умножили его на x — вычли из a*b — при необходимости прибавили или вычли x). Существует ли настоящее быстрое деление, я, к сожалению, не знаю. За O(L*log(L)^2*log(log(L)) — скорее всего, существует (для вычисления a/b считаем (a/M)/(b/M), где M=2^(L/2), потом корректируем… или что-то в этом роде).
А зачем проверять чётные?

bool SimplicityTest(int number) {
  if (number%2 == 0) { return false; }
  int temp = 3;
  do { if (number%temp == 0) { return false; }
  temp=temp+2; }
  while (temp<number);
  return true; }


В моём случае числа заведомо меньше int, если что.
использовать надо +2 +4 — это еще на треть сокращает поле
Ваш код работает за O(N), а код из первого примера — за O(N1/2), не говоря уже о правильном варианте, упомянутом выше в комментариях.
То есть при простом N около 106 вашему коду понадобится порядка 5 105 операций, а тому коду — всего 1000. Не надо заниматься микрооптимизациями прежде, чем достигните нормальной асимптотики.
Не надо заниматься микрооптимизациями прежде, чем достигните нормальной асимптотики.

Не совсем уверен. Очень аккуратный O(n*log(n)), учитывающий особенности кэша памяти, скорее всего, обгонит надёжного, но громоздкого O(n). В том числе, за счёт микрооптимизаций.
Конечно, к сравнению O(n) и O(sqrt(n)) это не относится, если только у нас нет уверенности, что n достаточно мало.
писать рекурсивные функции для решения такого рода задач не есть правильно. Даёшь только нерекурсивные версии. Самое позднее при переходе на длинную арифметику станет больно вызывать много раз одну и ту же функцию с тремя параметрами.
Тест Ферма говорит «n – простое с вероятностью 1– е^t», где e<= Fi(n)/n, Fi(n) — функция Эйлера, t — количество итераций. В случае составного числа имеющего только большые делители e приблизительно равно единице. То есть качество проверки практически равно 0. И тут уже не особо важно 100 или 1000 итераций вы сделаете, вероятность-то все равно близка к 0. Я это студентам каждый год рассказываю. Данный тест надо знать, но использовать его в реальных приложениях не нужно.
Не хотел никого задеть. Ну вот вам пример:
A = 18446743979220271189 = 4294967291*4294967279
Число А — составное, равное произведению двух простых чисел. Число A чуть меньше 2^64.
Тогда вероятность того, что тест ферма выдаст вам правильный ответ равна p = 1 — (0.9999999995)^t. При 100 итерациях p = 5 x 10^(-8), при 1000 итерациях p = 5 x 10^(-7).
На соревнованиях легко могут подсунуть такое число. В итоге ваша программа не выиграет. И дело даже не дойдет до чисел Кармайкла.
Если очень хочется использовать вероятностный алгоритма присмотритесь к тесту Миллера-Рабина. Для него вероятность e меньше в 4 раза. Как следствие, число необходимых итераций для получения заданной точности в сотни раз меньше. К сожалению оценку сложности на память не помню.
Тогда вероятность того, что тест ферма выдаст вам правильный ответ равна p = 1 — (0.9999999995)^t.

Проверил на Maple. Для первых 50 простых чисел все значения p^(A-1) mod A отличаются от 1… Какова вероятность такого события?
Провел собственное небольшое расследование. Похоже на то, что точное значение e в оценке e<= Fi(n)/n сильно отличается от верхней границы. Получается она реально годна только для сравнения эффективности тестов. Например аналогичная вероятность для теста Соловея-Штрассена в 2 раза ниже этой e<= Fi(n)/(2n), а для теста Миллера-Рабина в 4 раза e<= Fi(n)/(4n). Которые, похоже, так же сильно завышены.
Ноги у казанной оценки вроде бы растут отсюда «Молдовян Н.А., Молдовян А.А. Введение в криптосистемы с открытым ключом. – СПб.: БХВ-Петербург, 2005. 288 с.», но сам не читал… Авторы этой книги в некоторых кругах имеют дурную славу.
Мне кажется, вы оба сейчас перепутали условную вероятность того, что число простое при пройденном тесте, и условную вероятность того, что тест будет пройден, для составного числа.
Никто не перепутал. Изначально предполагается, что n составное.
В случае простого числа тест ошибиться не может, будут всегда выполняться все t =100 итераций.
Про «условную вероятность того, что число простое при пройденном тесте» можно говорить только когда есть какое-то априорное распределение тестируемых чисел. А его взять неоткуда, мы не знаем, откуда берутся тесты…
Оценка Fi(n)/n в тесте Ферма вполне может достигаться. А именно, на кармайкловых числах. Для остальных вероятность ложного положительного результата заметно меньше.
Нету ли у Вас ссылок на литературу (касаемо теста Ферма), лучше бумажную? Я о точности оценок.
К сожалению, нет — никогда не интересовался этим вопросом в такой степени, чтобы была нужна литература. Но в простейшем случае, когда n=p*q (p,q — простые) всё легко посчитать и так: число решений уравнения a^(n-1)=1 (mod n) равно gcd(p-1,q-1)^2. Например, для n=91 тест Ферма пройдёт в 36 случаях. А для больших p=2*q-1 вероятность ошибки будет стремиться к 1/2. Но для двух «случайно выбранных» простых чисел gcd(p-1,q-1) будет довольно маленьким, и вероятность ошибки тоже.

...if(b==1)

return a;

...

А почему не "return a%m;"? По-моему a*1 mod m = a mod m = a%m. Это ошибка автора или я что-то не так понял?

Выбран крайне неудачный ГПСЧ. Дело в том, что rand() возвращает число в диапазоне от 0 до RAND_MAX = 32767.

Sign up to leave a comment.

Articles