Длинная арифметика. С (С++)

[vlad_light]

Есть класс BigInt со скрытыми переменными uint32* m_integer и uint32 m_length, которые отвечают за само число и его длину соответственно. Нужно реализовать сложение, умножение и квадрат в этом классе.
Сложение: обычное сложение в столбик, выделяю длину на 1 больше на случай переполнения. Разбил на 3 случая: когда числа одинаковой длины, когда первое больше и когда первое меньше. Из-за этого получился достаточно громоздкий код, хотя, наверное, работает быстро. Можно его как-то улучшить?

(Оффтоп)

BigInt BigInt::operator+(const BigInt& rhs) const
{
    if (m_length == 0 || rhs.m_length == 0)
        throw "Empty number in BigInt::operator+";
    
    BigInt result;
 
    if (m_length == rhs.m_length)
    {       
        result.m_integer = new uint32[m_length + 1];
        result.m_length = m_length + 1;
        memset(result.m_integer, 0x0, (m_length + 1) * sizeof(uint32));
 
        for (uint32 i = 0; i < m_length; ++i)
        {
            result.m_integer[i] += m_integer[i] + rhs.m_integer[i];
            if (result.m_integer[i] < m_integer[i])
                result.m_integer[i + 1] = 1;
        }
 
        return result;
    }
 
    if (m_length > rhs.m_length)
    {
        result.m_integer = new uint32[m_length + 1];
        result.m_length = m_length + 1;
        memset(result.m_integer, 0x0, (m_length + 1) * sizeof(uint32));
 
        for (uint32 i = 0; i < rhs.m_length; ++i)
        {
            result.m_integer[i] += m_integer[i] + rhs.m_integer[i];
            if (result.m_integer[i] < m_integer[i])
                result.m_integer[i + 1] = 1;
        }
 
        for (uint32 i = rhs.m_length; i < m_length; ++i)
        {
            if (result.m_integer[i] == 0)
            {
                memcpy(&result.m_integer[i], &m_integer[i], (m_length - i) * sizeof(uint32));
                break;
            }
            else
            {
                result.m_integer[i] = m_integer[i] + 1;
                if (result.m_integer[i] == 0)
                    result.m_integer[i + 1] = 1;
            }
        }
    }
    else
    {
        result.m_integer = new uint32[rhs.m_length + 1];
        result.m_length = rhs.m_length + 1;
        memset(result.m_integer, 0x0, (rhs.m_length + 1) * sizeof(uint32));
 
        for (uint32 i = 0; i < m_length; ++i)
        {
            result.m_integer[i] += m_integer[i] + rhs.m_integer[i];
            if (result.m_integer[i] < m_integer[i])
                result.m_integer[i + 1] = 1;
        }
 
        for (uint32 i = m_length; i < rhs.m_length; ++i)
        {
            if (result.m_integer[i] == 0)
            {
                memcpy(&result.m_integer[i], &rhs.m_integer[i], (rhs.m_length - i) * sizeof(uint32));
                break;
            }               
            else
            {
                result.m_integer[i] = rhs.m_integer[i] + 1;
                if (result.m_integer[i] == 0)
                    result.m_integer[i + 1] = 1;
            }
        }
    }
    
    return result;
}

Умножение: тоже в столбик. Тут мне не нравится, что я использую uint64. Как этого можно избежать, не используя умножение дважды? (считаем, записываем в uint32, считаем тоже самое ещё раз, сдвигаем, записываем в uint32)

(Оффтоп)

BigInt BigInt::operator*(const BigInt& rhs) const
{
    if (m_length == 0 || rhs.m_length == 0)
        throw "Empty number in BigInt::operator*";
 
    BigInt result;
 
    result.m_integer = new uint32[m_length + rhs.m_length];
    result.m_length = m_length + rhs.m_length;
    memset(result.m_integer, 0x0, (m_length + rhs.m_length) * sizeof(uint32));
 
    for (uint32 i = 0; i < m_length; ++i)
    {
        uint32 carry = 0;
 
        for (uint32 j = 0; j < rhs.m_length; ++j)
        {
            uint64 temp = ((uint64)m_integer[i] * rhs.m_integer[j]) + result.m_integer[i + j] + carry;
            result.m_integer[i + j] = (uint32)temp;
            carry = (uint32)(temp >> 32);
        }
 
        result.m_integer[i + rhs.m_length] = carry;
    }
 
    return result;
}

Квадрат: тут я пользовался формулой:

$$(a_0+a_1x+a_2x^2+\ldots +a_nx^n)^2 = \sum _{i = 0}^n a_ix^{2i} + 2\sum _{k = 1}^n\sum _{i = 0}^{n - k}a_ia_{i+k}x^{2i+k}$$

(формулу придумал сам, если что -- подправьте). Не нравится "кривое" сложение во втором слагаемом. Можно его как-то улучшить?

(Оффтоп)

BigInt BigInt::sqr() const
{
    if (m_length == 0)
        throw "Empty number in BigInt::sqr";
 
    BigInt result;
 
    result.m_integer = new uint32[m_length << 1];
    result.m_length = m_length << 1;
    memset(result.m_integer, 0x0, (m_length << 1) * sizeof(uint32));
 
    // cross-products
    for (uint32 shift = 1; shift < m_length; ++shift)
    {
        uint32 carry = 0;
 
        for (uint32 i = 0; i < (m_length - shift); ++i)
        {
            uint64 temp = ((uint64)m_integer[i] * m_integer[i + shift]) + result.m_integer[shift + (i << 1)] + carry;
            result.m_integer[shift + (i << 1)] = (uint32)temp;
 
            temp >>= 32;
            result.m_integer[shift + (i << 1) + 1] += temp;
            if (result.m_integer[shift + (i << 1) + 1] < temp)
                carry = 1;
            else
                carry = 0;
        }
    }
 
    // shift
    for (uint32 i = result.m_length - 1; i > 0; --i)
    {
        result.m_integer[i] <<= 1;
        result.m_integer[i] |= result.m_integer[i - 1] >> 31;
    }
 
    // squares
    uint32 carry = 0;
    for (uint32 i = 0; i < m_length; ++i)
    {
        uint64 temp = ((uint64)m_integer[i] * m_integer[i]) + result.m_integer[i << 1] + carry;
        result.m_integer[i << 1] = (uint32)temp;
 
        temp >>= 32;
        result.m_integer[(i << 1) + 1] += temp;
        if (result.m_integer[(i << 1) + 1] < temp)
            carry = 1;
        else
            carry = 0;
    }
 
    return result;
}

[arseniiv]

vlad_light в сообщении #834577 писал(а):

Квадрат

А почему не использовать уже определённое умножение?

vlad_light в сообщении #834577 писал(а):

if (m_length == 0 || rhs.m_length == 0)
  throw "Empty number in BigInt::operator+";

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

vlad_light в сообщении #834577 писал(а):

Разбил на 3 случая: когда числа одинаковой длины, когда первое больше и когда первое меньше. Из-за этого получился достаточно громоздкий код, хотя, наверное, работает быстро. Можно его как-то улучшить?

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

Правильность кода не проверял.

[kp9r4d]

Мне очень нравится реализация длинки, которая приведена на emaxx, хотя и константы за счёт использования векторов получаются несколько больше, зато пишется в случае чего такая штука за 3-5 минут, что иногда бывает полезно (АСМ соревнования всякие там).

[Xaositect]

vlad_light в сообщении #834577 писал(а):

Сложение: обычное сложение в столбик, выделяю длину на 1 больше на случай переполнения. Разбил на 3 случая: когда числа одинаковой длины, когда первое больше и когда первое меньше. Из-за этого получился достаточно громоздкий код, хотя, наверное, работает быстро. Можно его как-то улучшить?

Вместо того, чтобы обнулять результат, скопируйте в него большее число и прибавьте к нему меньшее. Можете объявить вспомогательную функцию или оператор +=.

vlad_light в сообщении #834577 писал(а):

Умножение: тоже в столбик. Тут мне не нравится, что я использую uint64. Как этого можно избежать, не используя умножение дважды? (считаем, записываем в uint32, считаем тоже самое ещё раз, сдвигаем, записываем в uint32)

По-моему, uint64 - это лучший вариант.

[Pavia]

vlad_light
Пишите короткие функции. Не бойтесь делать много функций. Используйте, заимствуйте уже написанный протестированный и отработанный код.
А также я бы вам по рекомендовал подумать над тем, как можно уменьшить ваш код. К примеру указатели обменять.

А когда сделаете можете сравнить с
http://cybern.ru/csharp-long-math-add-substract.html

[vlad_light]

arseniiv в сообщении #834583 писал(а):

А почему не использовать уже определённое умножение?

Можно, но моя формула требует $\frac {n(n+1)}{2}$ умножений и $1$ сдвиг (умножение на $2$), а обычное умножение требует $n^2$ умножений, что, примерно, в 2 раза медленнее.

arseniiv в сообщении #834583 писал(а):

Разве плохо рассматривать числа без цифр как ноль?

Не подумал, спасибо!

arseniiv в сообщении #834583 писал(а):

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

Всё верно, спасибо! Про коммутативность забыл

kp9r4d в сообщении #834588 писал(а):

Мне очень нравится реализация длинки, которая приведена на emaxx

Ну, там по-детски

Xaositect в сообщении #834593 писал(а):

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

Спасибо!

Xaositect в сообщении #834593 писал(а):

Можете объявить вспомогательную функцию или оператор +=.

Я сделал наоборот: += через +

BigInt& BigInt::operator+=(const BigInt& rhs)
{
   *this = *this + rhs;
   return *this;
}

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

Xaositect в сообщении #834593 писал(а):

По-моему, uint64 - это лучший вариант.

Я имел ввиду, что не везде он есть.

Pavia в сообщении #834631 писал(а):

Пишите короткие функции. Не бойтесь делать много функций. Используйте, заимствуйте уже написанный протестированный и отработанный код.
А также я бы вам по рекомендовал подумать над тем, как можно уменьшить ваш код.

Стараюсь

Pavia в сообщении #834631 писал(а):

К примеру указатели обменять.

Можно по-подробнее, пожалуйста?

Pavia в сообщении #834631 писал(а):

А когда сделаете можете сравнить с http://cybern.ru/csharp-long-math-add-substract.html

Сравнил. По-моему, там не очень интересная реализация
К сожалению, на квадрат у меня пока времени не было, зато плюс доделал. Мало того, что код уменьшился, так ещё и работает немного быстрее (скорее всего из-за рекомендации Xaositect'а).

(Оффтоп)

BigInt BigInt::operator+(const BigInt& rhs) const
{
   if (m_length >= rhs.m_length)
   {
      if (m_length == 0)
         return rhs;
      if (rhs.m_length == 0)
         return *this;

      BigInt result(m_integer, m_length + 1); // !!!
      result.m_integer[m_length] = 0;

      for (uint32 i = 0; i < rhs.m_length; ++i)
         result.m_integer[i] += rhs.m_integer[i];

      for (uint32 i = 0; i < rhs.m_length; ++i)
         if (result.m_integer[i] < m_integer[i])
            while (++result.m_integer[i + 1] == 0)
               ++i;      

      return result;
   }
   else
      return rhs + *this;
}

Но есть одна проблема в помеченном месте: я считываю память за пределами выделенной -- это нормально или так делать нельзя?

[arseniiv]

vlad_light в сообщении #835731 писал(а):

Можно, но моя формула требует $\frac {n(n+1)}{2}$ умножений и $1$ сдвиг (умножение на $2$), а обычное умножение требует $n^2$ умножений, что, примерно, в 2 раза медленнее.

Асимптотически обе $O(n^2)$, но с «самодельной» реализацией больше способов внести ошибку (и в коде умножения, и в не использующем его коде квадрата vs. в коде умножения и в мааленьком коде квадрата, состоящем из вызова умножения и возврата результата). Это просто заметка.

vlad_light в сообщении #835731 писал(а):

Я имел ввиду, что не везде он есть.

Как так? (Вы же не думаете, что 64-разрядные целые есть только на x64?)

vlad_light в сообщении #835731 писал(а):

Но есть одна проблема в помеченном месте: я считываю память за пределами выделенной -- это нормально или так делать нельзя?

Даже если и получится, что вы надеетесь там (за пределами выделенной области) найти — нули? А они не обязаны там быть.

-- Ср мар 12, 2014 03:27:00 --

Только не понял, где вы в помеченном месте считываете память за пределами выделенной. Там же конструктор BigInt вызывается, а в нём, наверняка, выделяется столько, сколько вы указали.

-- Ср мар 12, 2014 03:29:29 --

И неясно, чему служит result.m_integer[m_length] = 0; — разве в конструкторе не зануляется весь массив m_integer? (Было бы странно, если бы он ставил нули по всем индексам кроме последнего — ведь это легко исправляется.)

[vlad_light]

arseniiv в сообщении #835747 писал(а):

Как так? (Вы же не думаете, что 64-разрядные целые есть только на x64?)

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

arseniiv в сообщении #835747 писал(а):

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

Простите, забыл написать конструктор

BigInt(const uint32* integer, const uint32 length)
{
    m_integer = new uint32[length];
    m_length = length;
    memcpy(m_integer, integer, length * sizeof(uint32));
}

Считывается "плохая" память в функции memcpy, поскольку я выделяю length + 1 памяти, а в массиве this->m_integer содержится только length элементов. Таким образом, с помощью result.m_integer[this->m_length] = 0; я обнуляю "мусор", который скопировался.

[Xaositect]

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

[Joker_vD]

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

Почему вообще нельзя было вместо пары uint32 *m_integer; uint32 m_length; взять обычный std::vector<uint32>? И кстати, кто такой вообще uint32?

[vlad_light]

Joker_vD в сообщении #835826 писал(а):

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

Придется следить за этим "руками".

Joker_vD в сообщении #835826 писал(а):

Почему вообще нельзя было вместо пары uint32 *m_integer; uint32 m_length; взять обычный std::vector<uint32>?

Потому что из-за бОльшей безопасности, std::vector<uint32> работает в разы медленнее.

Joker_vD в сообщении #835826 писал(а):

И кстати, кто такой вообще uint32?

Это typedef unsigned int uint32;

[EtCetera]

vlad_light в сообщении #835976 писал(а):

Joker_vD в сообщении #835826 писал(а):

Почему вообще нельзя было вместо пары uint32 *m_integer; uint32 m_length; взять обычный std::vector<uint32>?

Потому что из-за бОльшей безопасности, std::vector<uint32> работает в разы медленнее.

Если Вы имеете в виду проверку индексов на выход за пределы, то просто не используйте функцию at, которая, действительно, осуществляет такую проверку и бросает исключение std::out_of_range в случае неудачи. Используйте обычный operator[], который таковую проверку не осуществляет.

[vlad_light]

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

[Joker_vD]

(Оффтоп)

vlad_light в сообщении #835995 писал(а):

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

Бенчмарк в студию? http://stackoverflow.com/a/16446991/1440433 показывает обратное.

vlad_light в сообщении #835976 писал(а):

Это typedef unsigned int uint32;

Есть системы, где unsigned int занимает 2 байта, а есть и такие, где все 8.

Насчет исключений в конструкторе в данном случае я поторопился, извините.

[vlad_light]

Joker_vD в сообщении #836072 писал(а):

Бенчмарк в студию?

Как будет время -- обязательно напишу!

Joker_vD в сообщении #836072 писал(а):

unsigned int занимает 2 байта, а есть и такие, где все 8.

Да, знаю. Пусть будет uint32_t из <cstdint>.

Joker_vD в сообщении #836072 писал(а):

Насчет исключений в конструкторе в данном случае я поторопился

Добавил проверку на nullptr.
Написал квадрат, гляньте, пожалуйста, может что-то посоветуете. Формулу писал в стартовом посте.

(Оффтоп)

BigInt BigInt::sqr() const
{
   if (m_length == 0)
      return *this;

   BigInt result;

   result.m_integer = new uint32[m_length << 1];
   result.m_length = m_length << 1;
   memset(result.m_integer, 0x0, (m_length << 1) * sizeof(uint32));

   // cross-products
   for (uint32 shift = 1; shift < m_length; ++shift)
   {
      uint64 temp = 0;

      for (uint32 i = 0; i < (m_length - shift); ++i)
      {
         temp += ((uint64)m_integer[i] * m_integer[i + shift]) + result.m_integer[(i << 1) + shift];
         result.m_integer[(i << 1) + shift] = (uint32)temp;

         temp = (temp >> 32) + result.m_integer[(i << 1) + shift + 1];
         result.m_integer[(i << 1) + shift + 1] = (uint32)temp;

         temp >>= 32;
      }

      if (temp != 0)
         if (++result.m_integer[(m_length << 1) - shift] == 0)
            ++result.m_integer[(m_length << 1) - shift + 1];
   }
   
   // shift
   for (uint32 i = result.m_length - 1; i > 0; --i)
   {
      result.m_integer[i] <<= 1;
      result.m_integer[i] |= result.m_integer[i - 1] >> 31;
   }

   // squares
   uint64 temp = 0;

   for (uint32 i = 0; i < m_length; ++i)
   {
      temp += ((uint64)m_integer[i] * m_integer[i]) + result.m_integer[i << 1];
      result.m_integer[i << 1] = (uint32)temp;

      temp = (temp >> 32) + result.m_integer[(i << 1) + 1];
      result.m_integer[(i << 1) + 1] = (uint32)temp;

      temp >>= 32;
   }

   return result;
}