воскресенье, 29 августа 2010 г.

Длинное DIV длинное & Длинное MOD длинное

[Вся длинная арифметика] 

Пришло время познакомиться с самой “страшной” операцией на длинную арифметику. Первый раз я ее реализовал на Accepted 31 июля 2008. До этого момента мне было даже страшно подумать над ее реализацией. Но как оказалось не так уж это деление и страшно, как его малюют. Мы не будем изобретать велосипед, а просто грамотно, максимально понятно и быстро реализуем деление столбиком. Рассмотрим нижепредставленный рисунок.

Как видно osn = 100.

Заметим факт, что количество разрядов у частного от деления не превосходит количества разрядов у делимого. Т.е. начинаем формировать ответ со старшего разряда. На каждой итерации имеем текущее значение, которое пытаемся уменьшить на максимально большое количество раз делимым. Итак: как будем искать это “максимально большое количество”? За ответом отсылаю в нашу любимую первую главу “Программирование в алгоритмах”. Ответ прост. Вместо того, чтобы линейно перебирать цифры в интервале от 0 до osn в поисках нужного ответа, мы будем использовать бинарный поиск.

Корректность этого утверждения вытекает из того, что рассматриваемая функция, которую можно представить в виде b * x(где b – делимое, x – подбираемое значение) является монотонной. В нашем случае возрастающей.

Оценим выигрыш в скорости относительно “лобовой реализации”. Для osn = 10000 предложенная реализация будет работать на 10000/log(10000) ~ 769 раз быстрее.

Рассмотрим конкретную реализацию целочисленного деления двух длинных чисел:

  1. void LevelUp()
  2. {
  3.   for (int i = amount;i>=1;i--)
  4.     digits[i] = digits[i-1];
  5.   if (digits[amount])
  6.     amount++;
  7. }
  8. ...
  9. BigInt operator / (const BigInt &a, const BigInt &b)
  10. {
  11.   BigInt res;
  12.   BigInt curValue;
  13.   for (int i = a.amount-1; i>=0; i--)
  14.   {
  15.     curValue.LevelUp(); // * osn
  16.     curValue.digits[0] = a.digits[i];
  17.     // подбираем максимальное число x, такое что b * x <= curValue
  18.     int x = 0;
  19.     int l = 0, r = osn;
  20.     while (l <= r)
  21.     {
  22.       int m = (l + r) >> 1;
  23.       BigInt cur = b * m;
  24.       if (cur <= curValue)
  25.       {
  26.         x = m;
  27.         l = m+1;
  28.       }
  29.       else
  30.         r = m-1;
  31.     }
  32.     res.digits[i] = x;
  33.     curValue = curValue - b * x;
  34.   }
  35.   // избавляемся от лидирующих нулей
  36.   int pos = a.amount;
  37.   while (pos>=0 && !res.digits[pos])
  38.     pos--;
  39.   res.amount = pos+1;
  40.  
  41.   return res;
  42. }
* This source code was highlighted with Source Code Highlighter.

Обратите внимание на строчку 20. Строгая проверка дает неправильное решение. Советую разобраться почему.

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

  1. BigInt operator % (const BigInt &a, const BigInt &b)
  2. {
  3.   BigInt res;
  4.   BigInt curValue;
  5.   for (int i = a.amount-1; i>=0; i--)
  6.   {
  7.     curValue.LevelUp(); // * osn
  8.     curValue.digits[0] = a.digits[i];
  9.     // подбираем максимальное число x, такое что b * x <= curValue
  10.     int x = 0;
  11.     int l = 0, r = osn;
  12.     while (l <= r)
  13.     {
  14.       int m = (l + r) >> 1;
  15.       BigInt cur = b * m;
  16.       if (cur <= curValue)
  17.       {
  18.         x = m;
  19.         l = m+1;
  20.       }
  21.       else
  22.         r = m-1;
  23.     }
  24.     res.digits[i] = x;
  25.     curValue = curValue - b * x;
  26.   }
  27.   return curValue;
  28. }
* This source code was highlighted with Source Code Highlighter.

29 комментариев:

  1. ааа можете объяснить строчки 24-30?

    ОтветитьУдалить
  2. в строчках 20-31 происходит поиск числа x, на которое накладываются следующие ограничения:
    1) b * x <= curValue
    2) Оно должно быть максимальным и удовлетворять условию 1).

    Метод, с помощью которого ищем само число называется бинарным. Его реализацию в чистом виде можно посмотреть здесь: http://cppalgo.blogspot.com/2010/07/binarysearch-lowerbound-upperbound.html

    Итак, чтобы понять строчки 24-30 предлагаю разобрать пошагово алгоритм на конкретном примере.

    На второй итерации алгоритма имеем:
    curValue = 1-23, где 1 - старшый разряд, 23 - младший
    b = 53.

    Задача найти такое x, которое:
    53 * x <= 123 - неравенство (1)

    На каждой итерации вместо числа x будет подставлять переменную m в неравенство (1)

    Итерации бинарного поиска:
    1) l = 0, r = 100, m = 50 - понятно, что много. При этом нужно сдвинуть правую границу, чтобы сузить границы поиска. Можно просто присвоить r = 50, но мы же уже точно знаем, что 50 - не ответ на наш вопрос, поэтому можно смело присвоить r = m - 1, не боясь потерять правильный корень.
    2) l = 0, r = 49, m = 24
    3) l = 0, r = 23, m = 11
    4) l = 0, r = 10, m = 5
    5) l = 0, r = 4, m = 2 - первый раз выполнится условие cur <= curValue. Это говорит о том, что m - число, которое удовлетворяет неравенству (1), поэтому мы запоминаем это число в итоговую переменную x. Теперь наша задача грамотно подвинуть левую границу. Как и в случае с правой границей мы можем сделать просто l = m. Но мы уже обработали случай с числом m, поэтому можно также смело присвоить l = m + 1. После этой итерации не факт, что найденное число x будет максимальным. Для этого нужно выполнить остальные итерации бинарного поиска.
    6) l = 3, r = 4, m = 3.
    7) l = 3, r = 2 - конец итераций.

    Как можно было заменить "приемчик" с r = m - 1 и l = m + 1 позволяет избежать зацикливания на итерации 6).

    Найденное число x = 2. Оно полностью удовлетворяет поставленной задачи.

    ОтветитьУдалить
  3. ага, понятно, спасибо огромное

    ОтветитьУдалить
  4. Объясние пожалуйста, как работает функция LevelUp и куда ее нужно вставить в программе

    ОтветитьУдалить
  5. Наверное название функции LevelUp было не очень удачное, т.к. оно не в полной мере отображает ее смысл.

    Во-первых, эта функция относится к типу(структуре) BigInt и применяется для объектов этого типа. Это видно в строчке 15 реализации операции "/". Отсюда делаем вывод, что эту функцию нужно поместить внутрь структуры BigInt.

    Во-вторых, вполне ризонный вопрос - что делает эта функция? Название функции ответа на этот вопрос не дает. Глядим на ее реализацию. Если все еще остались вопросы. Поясню: эта функция умножает длинное число на основание системы счисления. Допустим мы имеем число 12345 и osn = 100. В памяти мы храним это число в виде массива с первыми треями элементами: 45, 23, 1. После выполнения функции получаем число 1234500, и массив с первыми 4 элементами: 0,45,23,1.
    Пример, когда osn = 10 не приводил. Надеюсь читатель это сможет сделать самостоятельно.

    Думаю это внесло определенную долю понимания.
    Если все же остались вопросы, этот исходник должен их прояснить: http://www.everfall.com/paste/id.php?3jq29e52lbwo

    ОтветитьУдалить
  6. Здравствуйте.. смысл деления вроде понятен. Такой вопросик только : почему при вычислениях есть неточность.. например 4200/42 = 99 а при числах побольше результат к примеру 34...04 а правильный 34...4 , т.е. дополнительный 0 в разряде десятков вылезает. Я по-другому немного храню число (с использованием вектора), но код функции не менял. Не подскажете в чем может быть дело? У Вас такого не было? Считает как надо?

    ОтветитьУдалить
  7. 1. По поводу 4200/42. У меня получается правильный ответ.
    2. Если вылазеет лишний ноль, то думаю что здесь скорее всего проблема с выбором основании системы счисления. Так например, в своем текущем решении у меня выставлены следующие настройки:
    int len = 4;
    const char * digitFormat = "%.4d";
    int osn = 1e4;

    функция output выглядит так:
    void output()
    {
    printf("%d",digits[amount-1]);
    for (int i= amount-2;i>=0;i--)
    printf(digitFormat,digits[i]);
    }

    Следовательно, если будет выводиться цифра числа отличная от старшего разряда, то необходимо предусмотреть возможность вывода лидирующих нулей. Лишний ноль может вылезти, если например osn = 10, а digitFormat = "%.2d";

    3. Данную реализацию тестировал здесь:
    http://informatics.mccme.ru/moodle/mod/statements/view3.php?id=251&chapterid=141
    http://informatics.mccme.ru/moodle/mod/statements/view3.php?id=251&chapterid=142
    там проблем не было.

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

    ОтветитьУдалить
  9. http://www.everfall.com/paste/id.php?3fn39vbuk8aj - вот мое решение

    ОтветитьУдалить
  10. Кажется понял в чем дело.. вместо условия (cur <= curValue) я писал (cur < curValue).. ща посмотрим что получится после замены.

    ОтветитьУдалить
  11. Почини рисунок, будь котиком.

    ОтветитьУдалить
  12. Использовал Ваш % для реализации RSA. Спасибо! Не знаю на сколько существенно, но мне пришлось добавить следующие проверки перед вычисление:

    if(a < b){
    return a;
    }

    if(a == b){
    return zero;
    }

    Может кому пригодится

    ОтветитьУдалить
    Ответы
    1. RSA прекрасно реализуется без длинной арифметики

      Удалить
    2. ссылка из википедии:
      "Выбираются два различных случайных простых числа p и q заданного размера (например, 1024 бита каждое)."

      Удалить
  13. а почему не работает рисунок(?

    ОтветитьУдалить
  14. Пожалуйста почините рисунок.
    И еще, могли бы вы показать адаптацию алгоритма для osn = 10

    ОтветитьУдалить
    Ответы
    1. рисунок за 2 года пока не починили(. Надеюсь исправимся.
      Что касается деления на 10. Не совсем понятно что вас интересует. вы можете смело исправлять osn на 10 и получите валидную реализацию деления длинного на длинное, но по osn = 10

      Удалить
  15. деление 678 / 4 дает не то что надо. В чем проблема?

    ОтветитьУдалить
    Ответы
    1. Так что там? Решение то не верное.

      Удалить
    2. Я бы не был так уверен, что решение неверное - скорее всего, проблема на вашей стороне.

      Удалить
    3. Dan Sagunov как на моей? Я делю 678 на 4 и получаю 1010. Если делитель состоит из одной цифры то ответы не верные в большинстве случаев...

      Удалить
    4. Ну во-первых, чтобы делать утверждения об ошибочной реализации, наверное, нужно показать исходный код, с которого вы запускаетесь.
      Во-вторых, данный цикл статей не преследовал задачи предоставления целостной реализации под копипаст. Все-таки в некоторых вещах нужно поразбираться. Скорее всего у вас ошибка при вводе/инициализации/выводе.
      В-третьих, ваше поведение считаю непрофессиональным. Если вы новичок и сидите на диване в ожидании готового решения, то это не есть путь война. Вам никто, ничего не должен! Если же вы бывалый и действительно нашли ошибку, то могли уже на нее указать, получив от нас благодарный фидбек.
      У авторов блога нет задачи искать ошибки в чужих реализациях, даже если они основаны на фрагментах наших исходников. Но чтобы закрыть тему, для вас сделаем исключение

      Удалить
    5. Ок)
      http://paste.ubuntu.com/24362452/ .сpp
      http://paste.ubuntu.com/24362465/ .h

      Удалить
    6. Ошибка в конструкторе по умолчанию: http://ideone.com/bpDQII

      Удалить
  16. По коду:
    1) curValue не инициализирован, но домножается на osn в методе LevelUp - как выглядит конструктор по умолчанию для BigInt?
    2) Не очень понятна логика работы алгоритма в месте, где LevelUp делается только один раз. В таком случае лидирующие нули получаются, но мы лишний раз выполняем бинарный поиск для их получения. Почему бы сначала не определить минимальное число curValue, которое >= делителя?

    ОтветитьУдалить
    Ответы
    1. 1) Конструктор по умолчанию BigInt инициализирует его так, чтобы он соответствовал значению 0.
      2) Уверен, что как раз логика вам понятна, поскольку вы написали, что мы лишь осуществим лишний проход бинпоиском. Вы правы - в случае, если curValue слишком мал, бинпоиск можно не выполнять.

      Удалить