Алгоритм целочисленного разделения

Я думал об алгоритме в делении больших чисел: делении с остатком bigint C на bigint D, где мы знаем представление C в базе b, а D имеет вид b ^ k-1. Вероятно, это проще всего показать на примере. Попробуем делить C = 21979182173 на D = 999.

Мы записываем число как наборы из трех цифр: 21 979 182 173
Мы берем суммы (по модулю 999) последовательных множеств, начиная слева: 21 001 183 356
Добавим 1 к тем наборам, которые предшествуют тем, где мы "перешли 999": 22 001 183 356

Действительно, 21979182173/999 = 22001183 и остаток 356.

Я вычислил сложность и, если я не ошибаюсь, алгоритм должен работать в O (n), n - количество цифр C в представлении base b. Я также сделал очень грубую и неоптимизированную версию алгоритма (только для b = 10) в С++, протестировал его против универсального алгоритма общего GMP, и он действительно выглядит лучше, чем GMP. Я не мог найти ничего подобного в любом месте, где бы я ни выглядел, поэтому мне пришлось прибегнуть к тестированию его против общего деления.

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

Таким образом, мой вопрос будет: есть ли статья или книга или что-то там, где описан вышеописанный алгоритм, возможно, обсуждая реализации? Если нет, было бы ли для меня разумным попытаться реализовать и протестировать такой алгоритм, скажем, в C/С++, или этот алгоритм каким-то образом изначально плох?

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

Спасибо большое!

Дальнейшее разъяснение вопросов, поднятых в комментариях/ответах:

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

Я прекрасно понимаю, что работа на основаниях 2 ^ n, вообще говоря, является, несомненно, наиболее эффективным способом ведения дел. Практически все библиотеки bigint используют 2 ^ 32 или что-то еще. Однако, что, если (и, я подчеркиваю, это было бы полезно только для этого конкретного алгоритма!), Мы реализуем bigints как массив цифр в базе b? Разумеется, здесь мы должны быть "разумными": b = 10, наиболее естественный случай, представляется достаточно разумным. Я знаю, что он более или менее неэффективен, учитывая память и время, учитывая, как числа хранятся внутри, но я был в состоянии, если мои (основные и, возможно, как-то ошибочные) тесты верны, дают результаты быстрее, чем общее разделение GMP, что имело бы смысл в реализации такого алгоритма.

Замечания о Ninefingers Мне пришлось бы использовать в этом случае дорогостоящую операцию modulo. Надеюсь, что нет: я вижу, если старый + новый скрещенный, скажем, 999, просто посмотрев на количество цифр старого + нового + 1. Если у нас 4 цифры, мы закончили. Более того, начиная с старых < 999 и новых <= 999, мы знаем, что если старый + новый + 1 имеет 4 цифры (он не может иметь больше), то (старый + новый)% 999 равен удалению самой левой цифры ( старый + новый + 1), который, я полагаю, мы можем сделать дешево.

Конечно, я не оспариваю очевидных ограничений этого алгоритма и не утверждаю, что он не может быть улучшен - он может делиться только с определенным классом чисел, и мы должны априори знать представление дивиденда в базе b, Однако, например, для b = 10 последнее кажется естественным.

Теперь, скажем, мы реализовали бонусы, как я изложил выше. Скажем C = (a_1a_2... a_n) в базе b и D = b ^ k-1. Алгоритм (который может быть, вероятно, гораздо более оптимизирован) будет идти следующим образом. Я надеюсь, что не так много опечаток.

если k > n, мы, очевидно, выполняем
добавьте нуль (т.е. a_0 = 0) в начале C (на всякий случай мы попытаемся разделить, скажем, 9999 с 99)
l = n% k (mod для "правильных" целых чисел - не должно быть слишком дорого)
old = (a_0... a_l) (первый набор цифр, возможно с меньшим числом символов)
для (i = l + 1; я < n; я = я + k) (Мы будем иметь пол (n/k) или так итерации)
- Новый = (a_i... а_ (я + к-1))
- new = new + old (это добавление bigint, таким образом, O (k))
- aux = new + 1 (опять же, добавление bigint - O (k) - что мне не нравится)
- если aux имеет более чем k цифр
  - удалить первую цифру aux
  - old = old + 1 (добавление bigint еще раз)
  - заполнить старый нулями в начале, чтобы он имел столько же цифр, сколько и должно
  - (a_ (i-k)... a_ (i-1)) = old (если я = l + 1, (a _ 0... a _ l) = old)
  - Новый = Окс
- заполнить ноль нулями в начале, чтобы он имел столько же цифр, сколько и должно
- (a_i... а_ (я + к-1) = новый
Quot = (A_0... а_ (п-к + 1))
бэр = новый

Там, спасибо за обсуждение этого со мной - как я уже сказал, мне кажется, что это интересный алгоритм "частного случая", который пытается реализовать, протестировать и обсудить, если никто не видит в нем каких-либо фатальных изъянов. Если это не так широко обсуждается до сих пор, еще лучше. Пожалуйста, дай мне знать, что ты думаешь. Извините за длинный пост.

Кроме того, еще несколько личных комментариев:

@Ninefingers: у меня есть некоторые (очень простые!) знания о том, как работает GMP, что он делает и общих алгоритмов разделения bigint, поэтому я смог понять большую часть ваших аргументов. Я также знаю, что GMP очень оптимизирован и в какой-то мере настраивает себя для разных платформ, поэтому я, конечно же, не пытаюсь "победить" в целом - это кажется столь же плодотворным, как атака танком с заостренной палкой. Однако это не идея этого алгоритма - он работает в очень особых случаях (что GMP, похоже, не охватывает). На несвязанной ноте, уверены ли вы, что общие деления выполняются в O (n)? Самое большее, что я видел, это M (n). (И это может, если я правильно понимаю, на практике (Schönhage-Strassen и т.д.) Не достигает O (n). Алгоритм Фюрера, который до сих пор не достигает O (n), является, если я прав, почти чисто теоретическое.)

@Avi Berger: На самом деле это не похоже на то, что "выкидывать девятки", хотя идея схожа. Однако вышеупомянутый алгоритм должен работать все время, если я не ошибаюсь.

Ответ 1

Ваш алгоритм представляет собой вариацию базового 10 алгоритма, известного как "изгнание девяток". В вашем примере используется база 1000 и "выкапывание" 999 (на один меньше, чем у базы). Раньше это учили в начальной школе, чтобы быстро провести расчеты. У меня был учитель математики средней школы, который был в ужасе, узнав, что его больше не учат, и наполнили его.

Отбрасывание 999 в базе 1000 не будет работать как алгоритм общего деления. Он будет генерировать значения, которые сравниваются по модулю 999 с фактическим фактором и остатком, а не с фактическими значениями. Ваш алгоритм немного отличается, и я не проверял, работает ли он, но он основан на эффективном использовании базы 1000, а делитель на 1 меньше базы. Если вы хотите попробовать его для разделения на 47, вам придется сначала преобразовать в базовую систему счисления 48.

Google "вычеркивает девятки" для получения дополнительной информации.

Изменить: я изначально читал ваш пост слишком быстро, и вы знаете об этом как о рабочем алгоритме. Поскольку @Ninefingers и @Karl Bielefeldt заявили более четко, чем я, в своих комментариях, то, что вы не включаете в свою оценку эффективности, является преобразование в базу, подходящую для конкретного делителя.

Ответ 2

Я чувствую необходимость добавить к этому на основе моего комментария. Это не ответ, а объяснение фона.

Библиотека bignum использует так называемые конечности - поиск mp_limb_t в источнике gmp, обычно это целочисленное поле фиксированного размера.

Когда вы делаете что-то вроде добавления, один из способов (хотя и неэффективный) подходит к нему:

doublelimb r = limb_a + limb_b + carryfrompreviousiteration

Эта двойная конечность ловит переполнение limb_a + limb_b в случае, если сумма больше размера конечности. Так что если сумма больше 2 ^ 32, если мы используем uint32_t в качестве нашего размера конечности, перехват может быть пойман.

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

Это может показаться неэффективным, но это всего лишь способ сделать что-то. Чтобы действительно разбить большие пушки, x86 имеет adc в качестве инструкции - добавить с переносом. То, что это делает, является арифметикой и в ваших полях и устанавливает бит переноса, если арифметика переполняет размер регистра. В следующий раз, когда вы делаете add или adc, факторы процессора также переносятся. В вычитании он называется флагом заимствования.

Это также относится к операциям сдвига. Таким образом, эта особенность процессора имеет решающее значение для того, что делает bignums быстрыми. Так что дело в том, что в чипе есть электронные схемы для этого - делать это в ПО всегда будет медленнее.

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

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

Я также хотел бы добавить, что, хотя у вас есть, вероятно, быстро для этого случая, имейте в виду, что в качестве библиотеки bignum gmp для вас работает справедливая работа, например управление памятью. Если вы используете mpz_, вы используете абстракцию выше того, что я здесь описал, для стартеров. Наконец, gmp использует ручную оптимизированную сборку с развернутыми контурами почти для каждой платформы, о которой вы когда-либо слышали, и еще больше. Там очень хорошая причина, по которой он поставляется с Mathematica, Maple и др.

Теперь, только для справки, некоторые материалы для чтения.

Современная компьютерная арифметика - это работа, подобная Кнуту, для библиотек произвольной точности.
Дональд Кнут, Семинумерные алгоритмы ( "Искусство программирования", том II).
блог Уильяма Харта об алгоритме реализации bsdnt в котором он обсуждает различные алгоритмы деления. Если вас интересуют библиотеки bignum, это отличный ресурс. Я считал себя хорошим программистом, пока не начал следить за такими вещами...

Подводя итог вам: инструкции сборки сборок сосать, поэтому люди обычно вычисляют инверсии и умножаются вместо этого, как и при определении деления в модульной арифметике. Различные методы, которые существуют (см. MCA), в основном O (n).

Изменить: Хорошо, не все методы O (n). Большинство методов, называемых div1 (разделение на что-то не больше, чем конечность O (n). Когда вы идете больше, вы заканчиваете сложность O (n ^ 2), чего трудно избежать.

Теперь вы можете реализовать bigints как массив цифр? Ну да, конечно, можно. Однако рассмотрим идею только при добавлении

/* you wouldn't do this just before add, it just to 
   show you the declaration.
 */
uint32_t* x = malloc(num_limbs*sizeof(uint32_t));
uint32_t* y = malloc(num_limbs*sizeof(uint32_t));
uint32_t* a = malloc(num_limbs*sizeof(uint32_t));
uint32_t m;

for ( i = 0; i < num_limbs; i++ )
{
    m = 0;
    uint64_t t = x[i] + y[i] + m;
    /* now we need to work out if that overflowed at all */
    if ( (t/somebase) >= 1 ) /* expensive division */
    {
        m = t % somebase; /* get the overflow */
    }
}

/* frees somewhere */

Это приблизительный эскиз того, что вы ищете для добавления по вашей схеме. Поэтому вам нужно запустить преобразование между базами. Таким образом, вам понадобится преобразование в ваше представление для базы, а затем назад, когда вы закончите, потому что эта форма просто очень медленная везде. Мы не говорим о различии между O (n) и O (n ^ 2) здесь, но мы говорим о дорогостоящей инструкции деления на конечность или дорогостоящем преобразовании каждый раз, когда вы хотите разделить. См. это.

Далее, как вы расширяете свое разделение для разделения общего дела? Под этим я имею в виду, когда вы хотите разделить эти два числа x и y из приведенного выше кода. Вы не можете, это ответ, не прибегая к средствам на основе бигума, которые дороги. См. Кнут. Принимая по модулю число, большее вашего размера, не работает.

Позвольте мне объяснить. Попробуйте 21979182173 mod 1099. Предположим здесь для простоты, что наибольшее поле размера может иметь три цифры. Это надуманный пример, но самый большой размер поля, который я знаю, использует 128 бит с использованием расширений gcc. Во всяком случае, дело в том, что вы:

21 979 182 173

Разделите свой номер на конечности. Затем вы принимаете по модулю и сумме:

21 1000 1182 1355

Это не работает. Именно здесь Avi является правильным, потому что это форма извлечения девяток или их адаптация, но здесь это не работает, потому что наши поля переполнены для начала - вы используете модуль, чтобы обеспечить, чтобы каждое поле находилось внутри его размер конечности/поля.

Итак, какое решение? Разделите свой номер на ряд подходящих размеров? И начните использовать функции bignum, чтобы вычислить все, что вам нужно? Это будет намного медленнее, чем любой существующий способ непосредственного манипулирования полями.

Теперь, возможно, вы предлагаете только этот случай для деления на конечность, а не на бонусы, и в этом случае он может работать, но деление на генсе и предварительно сконвертированные обратные и т.д. делают без требования преобразования. Я понятия не имею, будет ли этот алгоритм быстрее, чем говорить о подразделении hensel; это было бы интересным сравнением; проблема связана с общим представлением через библиотеку bignum. Представление, выбранное в существующих библиотеках bignum, по причинам, которые я расширил, - это имеет смысл на уровне сборки, где это было сделано впервые.

В качестве побочного примечания; вам не нужно использовать uint32_t для представления ваших конечностей. Размер в идеале определяется размером регистров системы (например, uint64_t), чтобы вы могли использовать версии, оптимизированные для сборки. Таким образом, в 64-битной системе adc rax, rbx устанавливается только переполнение (CF), если результат перегружает 2 ^ 64 бит.

tl; dr версия: проблема не в вашем алгоритме или идее; это проблема преобразования баз, поскольку представление, которое вам нужно для вашего алгоритма, не самый эффективный способ сделать это в add/sub/mul и т.д. Чтобы перефразировать knuth: Это показывает вам разницу между математической элегантностью и вычислительной эффективностью.

Ответ 3

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

Вы можете использовать базу 999, если хотите; нет ничего особенного в использовании базы с мощностью-10, за исключением того, что она делает преобразование в десятичное целое очень дешевым. (Вы можете работать на одном конечности за раз, вместо того, чтобы выполнять полное деление на целое целое число. Это похоже на разницу между преобразованием двоичного целого в десятичное число и умножением на каждые 4 бита на шестую цифру. с наиболее значимыми битами, но преобразование в базы без питания 2 должно быть LSB-первым с использованием деления.)

Например, чтобы вычислить первые 1000 десятичных цифр Фибоначчи (10 ⁹) для вопроса с кодом-гольфа с требованием производительности, мои 105 байтов ответа машинного кода x86 использовал тот же алгоритм, что и этот ответ на Python: обычная итерация a+=b; b+=a Fibonacci, но делятся на (мощность) 10 каждый время a становится слишком большим.

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

Деление на мощность 2 не работает, если вы не отслеживаете, сколько полномочий 2 вы отбросили, потому что возможное двоичное → десятичное преобразование в конце будет зависеть от этого.

Итак, для этого алгоритма вам необходимо выполнить добавление с расширенной точностью и деление на 10 (или любую мощность, которую вы хотите).

Я сохранил базы-10 ⁹ конечности в 32-битных целочисленных элементах. Деление на 10 ⁹ тривиально дешево: просто приращение указателя, чтобы пропустить нижнюю конечность. Вместо фактического выполнения memmove я просто компенсирую указатель, используемый следующей итерацией.

Я думаю, что деление на 10, кроме 10 ^ 9, будет несколько дешевым, но потребует фактического деления на каждую конечность и распространения остатка на следующую конечность.

Добавление с добавленной точностью несколько более дорогое, чем с бинарными конечностями, потому что я должен выполнить выполнение вручную с помощью сравнения: sum[i] = a[i] + b[i]; carry = sum < a; (сравнение без знака). А также вручную обернуть 10 ^ 9 на основе этого сравнения с инструкцией условного перемещения. Но я смог использовать это выполнение как вход в adc (инструкция x86 add-with-carry).

Вам не нужен полный modulo для обработки упаковки при добавлении, потому что вы знаете, что вы завернули не более одного раза.

Это отбрасывает чуть более 2 бит каждой 32-разрядной конечности: 10 ^ 9 вместо 2^32 = 4.29... * 10^9. Хранение базы-10 цифр по одному на байт было бы значительно менее экономичным по площади и намного хуже для производительности, потому что 8-битное двоичное добавление стоит так же, как и 64-битное двоичное дополнение на современном 64-битном процессоре.

Я стремился к размеру кода: для чистой производительности я бы использовал 64-битные конечности, содержащие base-10 ^ 19 "цифр. (2^64 = 1.84... * 10^19, поэтому это составляет менее 1 бит на 64.) Это позволяет получить в два раза больше работы с каждой аппаратной командой add. Хм, на самом деле это может быть проблемой: сумма двух конечностей может обернуть 64-битное целое число, поэтому просто проверка на > 10^19 уже недостаточна. Вы можете работать в базе 5*10^18 или в базе 10^18 или выполнять более сложное обнаружение выполнения, которое проверяет двоичную перенос, а также ручную перенос.

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

В целом моя версия работала примерно на 10 раз быстрее, чем версия расширенной точности Python на одном и том же оборудовании (но у нее была возможность для значительной оптимизации скорости, делясь реже). (70 секунд или 80 секунд против 12 минут)

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