Большие целые числа и C

Поскольку современные владельцы очень эффективны при работе с числами фиксированной длины в битах, почему бы вам не иметь их массив?

Предположим, вы используете unsigned long long. Они должны иметь ширину 64 бита, поэтому максимально возможное unsigned long long должно быть 2^64 - 1. Давайте представим любое число в виде набора чисел следующим образом:

-big_num = ( n_s, n_0, n_1, ...)

-n_s будет принимать только 0 и 1 для представления знака + и -

-n_0 будет представлять число от 0 до 10 ^ a -1 (показатель степени a будет определяться)

-n_1 будет представлять число от 10^a до 10^(a+1) -1 и так далее, и так далее...

ОПРЕДЕЛЕНИЕ:

Все n_ ДОЛЖНЫ быть ограничены 10^a-1. Таким образом, при добавлении двух big_num это означает, что нам нужно добавить n_ следующим образом:

// A + B = ( (to be determent later),
//          bound(n_A_1 + n_B_1) and carry to next,
//          bound(n_A_2 + n_B_2 + carry) and carry to next,
//        ...)

Ограничение может быть выполнено как:

bound(n_A_i + n_B_i + carry) = (n_A_i + n_B_i + carry)%(10^a)

Следовательно, перенос на i+1 определяется как:

// carry (to be used in i+1) = (n_A_i + n_B_i + carry)/(10^a) 
// (division of unsigned in c++ will floor the result by construction)

Это говорит нам о том, что наихудший случай — carry = 10^a -1, и, таким образом, наихудшее сложение (n_A_i + n_B_i + carry) — это: (worst case) (10^a-1) + (10^a-1) + (10^a-1) = 3*(10^a-1) Поскольку тип unsigned long long, если мы не хотим иметь переполнение при этом сложении, мы должны ограничить нашу экспоненту a так, чтобы:

//    3*(10^a-1) <= 2^64 - 1, and a an positive integer
// => a <= floor( Log10((2^64 - 1)/3 + 1) )
// => a <= 18

Таким образом, теперь зафиксировано максимально возможное значение a=18, и, таким образом, максимально возможное число n_, представленное с помощью unsigned long long, равно 10^18 -1 = 999,999,999,999,999,999. С этой базовой настройкой давайте теперь перейдем к реальному коду. Сейчас я буду использовать std::vector для хранения big_num, о котором мы говорили, но это может измениться:

// Example code with unsigned long long
#include <cstdlib>
#include <vector>
//
// FOR NOW BigNum WILL BE REPRESENTED
// BY std::vector. YOU CAN CHANGE THIS LATTER
// DEPENDING ON WHAT OPTIMIZATIONS YOU WANT
//
using BigNum = std::vector<unsigned long long>;

// suffix ULL garanties number be interpeted as unsigned long long
#define MAX_BASE_10 999999999999999999ULL

// random generate big number
void randomize_BigNum(BigNum &a){
  // assuming MyRandom() returns a random number
  // of type unsigned long long
  for(size_t i=1; i<a.size(); i++)
    a[i] = MyRandom()%(MAX_NUM_BASE_10+1); // cap the numbers
}

// wrapper functions
void add(const BigNum &a, const BigNum &b, BigNum &c); // c = a + b
void add(const BigNum &a, BigNum &b);         // b = a + b

// actual work done here
void add_equal_size(const BigNum &a, const BigNum &b, BigNum &c, size_t &N);
void add_equal_size(const BigNum &a, const BigNum &b, size_t &N);
void blindly_add_one(BigNum &c);
// Missing cases
// void add_equal_size(BigNum &a, BigNum &b, BigNum &c, size_t &Na, size_t &Nb);
// void add_equal_size(BigNum &a, BigNum &b, size_t &Na, size_t &Nb);

int main(){
  size_t n=10;
  BigNum a(n), b(n), c(n);
  randomize_BigNum(a);
  randomize_BigNum(b);
  add(a,b,c);
  return;
}

Функции-оболочки должны выглядеть следующим образом. Они защитят от неправильного размера вызовов массива:

// To do: add support for when size of a,b,c not equal

// c = a + b
void add(const BigNum &a, const BigNum &b, BigNum &c){

  c.resize(std::max(a.size(),b.size()));

  if(a.size()==b.size())
    add_equal_size(a,b,c,a.size());
  else
    // To do: add_unequal_size(a,b,c,a.size(),b.size());

  return;
};
// b = a + b
void add(const BigNum &a, const BigNum &b){

  if(a.size()==b.size())
    add_equal_size(a,b,a.size());
  else{
    b.resize(a.size());
    // To do: add_unequal_size(a,b,a.size());
  }

  return;
};

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

// If a,b,c same size array
// c = a + b
void add_equal_size(const BigNum &a, const BigNum &b, BigNum &c, const size_t &N){
  // start with sign of c is sign of a
  // Specific details follow on whether I need to flip the
  // sign or not
  c[0] = a[0];
  unsigned long long carry=0;

  // DISTINGUISH TWO GRAND CASES:
  //
  // a and b have the same sign
  // a and b have oposite sign
  // no need to check which has which sign (details follow)
  //
  if(a[0]==b[0]){// if a and b have the same sign
    //
    // This means that either +a+b or -a-b=-(a+b)
    // In both cases I just need to add two numbers a and b
    // and I already got the sign of the result c correct form the
    // start
    //
    for(size_t i=1; i<N;i++){
      c[i]  = (a[i] + b[i] + carry)%(MAX_BASE_10+1);
      carry = c[i]/(MAX_BASE_10+1);      
    }
    if(carry){// if carry>0 then I need to extend my array to fit the final carry
      c.resize(N+1);
      c[N]=carry;
    }
  }
  else{// if a and b have opposite sign
    //
    // If I have opposite sign then I am subtracting the
    // numbers. The following is inspired by how
    // you can subtract two numbers with bitwise operations.
    for(size_t i=1; i<N;i++){
      c[i]  = (a[i] + (MAX_BASE_10 - b[i]) + carry)%(MAX_BASE_10+1);
      carry = c[i]/(MAX_BASE_10+1);      
    }
    if(carry){ // I carried? then I got the sign right from the start
      // just add 1 and I am done
      blindly_add_one(c);
    }
    else{ // I didn't carry? then I got the sign wrong from the start
      // flip the sign
      c[0] ^= 1ULL;
      // and take the compliment
      for(size_t i=1; i;<N;i++)
    c[i] = MAX_BASE_10 - c[i];
    }
  }
  return;
};

Ниже приведены некоторые подробности о случае // if a and b have opposite sign: Давайте поработаем с основанием 10. Допустим, мы вычитаем a - b. Давайте преобразуем это в сложение. Определите следующую операцию:

Назовем основание 10 цифр числа di. Тогда любое число будет n = d1 + 10*d2 + 10*10*d3... Дополнение цифры теперь будет определяться как:

     `compliment(d1) = 9-d1`

Тогда комплимент числа n равен:

   compliment(n) =         compliment(d1)
                   +    10*compliment(d2)
                   + 10*10*compliment(d3)
                   ...

Рассмотрим два случая, a>b и a<b:

ПРИМЕР a>b: если не сказать a=830 и b=126. Сделайте следующее 830 - 126 -> 830 + compliment(126) = 830 + 873 = 1703 хорошо, так что если a>b, я опускаю 1 и добавляю 1, результат будет 704!

ПРИМЕР a<b: если не говорить a=126 и b=830. Выполните следующие действия 126 - 830 -> 126 + compliment(830) = 126 + 169 = 295 ...? Ну а если я сделаю комплимент? compliment(295) = 704 !!! так что если a<b у меня уже есть результат... с обратным знаком.

Переходя к нашему случаю, поскольку каждое число в массиве ограничено MAX_BASE_10, дополнение наших чисел равно

compliment(n) = MAX_BASE_10 - n

Таким образом, используя этот комплимент для преобразования вычитания в сложение, мне нужно обратить внимание только на то, есть ли у меня лишняя 1 в конце сложения (случай a>b). Алгоритм теперь такой

• ДЛЯ КАЖДОГО вычитания МАССИВА (i-я итерация):
• na_i - nb_i + перенос(i-1)
• конвертировать -> na_i + комплимент (nb_i) + перенос (i-1)
• связал результат -> (na_i + комплимент (nb_i) + перенос (i-1))% MAX_BASE_10

• найти перенос -> (na_i + комплимент (nb_i) + перенос (i-1))/MAX_BASE_10

• продолжайте добавлять номера массива...

• В конце массива, если я перенес, забудьте перенос и добавьте 1. В противном случае возьмите комплимент результата

Это «и добавить один» выполняется еще одной функцией:

// Just add 1, no matter the sign of c
void blindly_add_one(BigNum &c){
  unsigned long long carry=1;
  for(size_t i=1; i<N;i++){
    c[i]  = carry%(MAX_BASE_10+1);
    carry = c[i]/(MAX_BASE_10+1);
  }
  if(carry){ // if carry>0 then I need to extend my basis to fit the number
    c.resize(N+1);
    c[N]=carry;
  }
};

Хорошо до сюда. Конкретно в этом коде не забываем, что в начале функции мы устанавливаем знак c в знак a. Так что, если я перенесу в конце, это означает, что у меня было |a|>|b|, и я сделал либо +a-b>0, либо -a+b=-(a-b)<0. В любом случае установка знака результатов c на знак a была правильной. Если я не ношу, у меня было |a|<|b| либо с +a-b<0, либо с -a+b=-(a-b)>0. В любом случае установка знака результатов c на знак a была НЕПРАВИЛЬНОЙ, поэтому мне нужно перевернуть знак, если я не ношу.

Следующие функции работают так же, как и предыдущая, только вместо того, чтобы c = a + b делать это, нужно b = a + b

// same logic as above, only b = a + b
void add_equal_size(BigNum &a, BigNum &b, size_t &N){

  unsigned long long carry=0;
  if(a[0]==b[0]){// if a and b have the same sign
    for(size_t i=1; i<N;i++){
      b[i]  = (a[i] + b[i] + carry)%(MAX_BASE_10+1);
      carry = b[i]/(MAX_BASE_10+1);      
    }
    if(carry){// if carry>0 then I need to extend my basis to fit the number
      b.resize(N+1);
      b[N]=carry;
    }
  }
  else{ // if a and b have oposite sign
    b[0] = a[0];
    for(size_t i=1; i<N;i++){
      b[i]  = (a[i] + (MAX_BASE_10 - b[i]) + carry)%(MAX_BASE_10+1);
      carry = b[i]/(MAX_BASE_10+1);      
    }
    if(carry){
      add_one(b);
    }
    else{
      b[0] ^= 1ULL;
      for(size_t i=1; i;<N;i++)
    b[i] = MAX_BASE_10 - b[i];
    }
  }
  return;
};

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

ОТКУДА ИДТИ

С этого момента нужно сделать много вещей для оптимизации кода, я упомяну несколько, которые я мог придумать:

- Попробуйте заменить добавление массивов возможными вызовами BLAS.

-Убедитесь, что вы используете преимущества векторизации. В зависимости от того, как вы пишете свои циклы, вы можете генерировать или не генерировать векторизованный код. Если ваши массивы станут большими, вы можете извлечь из этого выгоду.

- В духе вышеизложенного убедитесь, что вы правильно выровняли массивы в памяти, чтобы действительно воспользоваться преимуществами векторизации. Насколько я понимаю, std::vector доза не гарантирует выравнивания. Ни дозировать вслепую malloc. Я думаю, что библиотеки boost имеют векторную версию, в которой вы можете объявить фиксированное выравнивание, и в этом случае вы можете запросить 64-битный выровненный массив для вашего массива unsigned long long. Другой вариант — иметь свой собственный класс, который управляет необработанным указателем и выравниванием по дозе с помощью специального распределителя. Заимствование aligned_malloc и aligned_free с сайта https://embeddedartistry.com/blog/2017/02/22/generating-aligned-memory/ у вас может быть такой класс для замены std::vector:

// aligned_malloc and aligned_free from:
// https://embeddedartistry.com/blog/2017/02/22/generating-aligned-memory/

// wrapping in absolutly minimal class to handle
// memory allocation and freeing 
class BigNum{
private:
  unsigned long long *ptr;
  size_t size;
public:  
  BigNum() : ptr(nullptr)
       , size(0)
  {};  

  BigNum(const size_t &N) : ptr(nullptr)
              , size(N)
  {
    resize(N);
  }  
  // Defining destructor this will now delete copy and move constructor and assignment. Make your own if you need them
  ~BigNum(){
    aligned_free(ptr);
  }  

  // Access an object in aligned storage
  const unsigned long long& operator[](std::size_t pos) const{
    return *reinterpret_cast<const unsigned long long*>(&ptr[pos]);
  }
  // return my size
  void size(){    
    return size;
  }
  // resize memory allocation
  void resize(const size_t &N){
    size = N;
    if(N){
      void* temp = aligned_malloc(ptr,N+1); // N+1, always keep first entry for the sign of BigNum 
      if(temp!=nullptr)
    ptr = static_cast<unsigned long long>(temp);
      else
    throw std::bad_alloc();
    }
    else{
      aligned_free(ptr);
    }
  }
};