Быстрое умножение длинных чисел. Опции

[klem4]

Всем привет. Может кто-то подсказать описание хорошего алгоритма умножения длиннных чисел ? Желательно с примерами работы алгоритма. Необходимо за время ~ 1c. перемножить несколько сотен длинных чисел (100-6000 знаков) на 3-4-х значные числа. Простой столбик в этом время не укладывается

[-TarasBer-]

Надо умножить длинное на короткое?
Алгоритм оптимизировать тут некуда.
Если столбик не очень далеко не укладывается, то оптимизируй код: храни все числа в бинарном виде (ну то есть не цифры храни, а блоки по 4 байта от каждого), после каждого умножения запоминай старое значение edx и прибавляй к следующему результату.

[Lapp]

Думаю, кроме столбика все равно ничего нет. Другой вопрос - КАК ты хранишь числа и как организовыван счет. Тут много ресурсов для оптимизации..

[klem4]

Буду рад любым подсказкам, числа храню в массиве (vector), вот пример умножения факториала 1110 на 1111 (получение 1111!)

Время одного умножения слишком большое(для решения моей задачи):
$ g++ -o forum forum.cpp && time ./forum
real 0m0.021s
user 0m0.016s
sys 0m0.004s

Если нужны пояснения к коду, приведу. Комменты общие сейчас поставлю.

#include <vector>
#include <string>
#include <stdio.h>

using namespace std;
typedef vector<int> Int;
typedef vector<Int> Ints;

// создание длинного числа из строки или целого
Int create( string s );
Int create( int n );

// печать длинного числа с переводом строки или без
void dump( Int a, bool cr );

// сложить два длинных с учетом сдвига второго относительно первого на shift влево
Int add(Int, Int, int shift);

// результат сложения массива длинных Ints (каждое последующее суммируется со сдвигом = сдвиг предыдущего + 1)
Int add_multi(Ints adds);

// умножить длинное на целое < 10
Int mult(Int, int);

// умножить длинное на длинное
Int mult(Int, Int);

int main()
{
    Int a =  create (& quot;145685245062266955405470957253929395052590595389527229775646520190327903725
 54007898400184159901236772947210799341232610833194173685864708721316793931909617
 22026542039712903994026600771128190953332907351098095093325712235742366518965010
 36571646308576256679272529325332526486576675624833429362316874934756651946649134
 21410633095260060537202807613431275295936480324590769021182190080881874327665650
 82583218005587067768716119803433798349035964310403606789737274766911641119521737
 46180656971020016031834607226074608594179294098389656489920957213783176615357614
 14359115669876891073985209131632267336569535231402438829077076942241254266873239
 33851466476850382651242921683126512786170363342152348041274427388172196550649455
 04401632115036773889846308182465643456178540901589017145380088890788421309831606
 43164226461351388108484835111945923176511690540408121566911147935224253920676551
 18647308475078508028834253397857726143901906354042126876485278456654893656062499
 25314050447252977272210708850546544885414450253552644077861793760271638336853880
 90906042155615771122351797746664895677687827680045944606819654993445463775105030
 62444412141832970108764127956861734717364630676445026481879245943442426165844096
 07611262493119341070227211745826733087220725629450421909884069575797597455107505
 82124621146336278229630516458056839008369980416547287938512590094701669234821173
 23384201492314647650739275485131342519309934362226353036296318568807152391082874
 98012708422343976153016889113943927273447360011554333301306939682348088452729905
 75843876283904126047988479220409814640069306760537593290117911314545746452949340
 53382041885300571225804204342431303730575560849376831834711525483668892375430232
 65026366216328132203543393012960988307772319594027441112056698283444243265160316
 44427279643772694590856982768531258528854972002877502225592362340489050009670717
 21577995711126231664571344342816789492287780579060212551757865802215913547054447
 11155740006547334931711056824872216863537059642129801584847265791171095033358865
 64558956250708640776789567388972216276588794189464973671190022877952494827564789
 75282440342710765029244416033234680346588871788663330292877711951125270788410677
 93183006671733289181877681317209401908042137605998301921907449818635479782500922
 19571874012818071327198208465830561262476520827812910313314823857879770562447910
 90777264932140189384617585751625700155612698447315032474397695537090973801981148
 67832246304778402509061034345682207861345822729165491122453539190230833073670813
 17568990801499084824509871236338438005351423039831796266968985125938052504797785
 95802226159971365933989904645726991898398910965318798719888460374605824000000000
 00000000000000000000000000000000000000000000000000000000000000000000000000000000
 00000000000000000000000000000000000000000000000000000000000000000000000000000000
 00000000000000000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000");
    Int b = create("1111");

    Int r = mult(a, b);
    //dump(r, true);

    return 0;
}

Int create( string s )
{
    Int result;
    for (long long i = 0; i < s.length(); i++ )
	result.push_back((long long)(s[i]) - (long long)'0');
    return result;
}

void dump( Int a, bool cr )
{
    for (long long i = 0; i < a.size(); i++ )
	printf("%d", a[i]);
	
    if ( cr )
	printf ("\n");
}

Int mult(Int a, int n)
{
    Int result;
    int a_size = a.size();

    int bonus = 0;
    for (int i = a_size - 1; i >= 0; i-- )
    {
	int value = bonus + a[i] * n;
	if ( value >= 10 )
	{
	    bonus = (int)(value / 10);
	    value %= 10;
	}
	else
	{
	    bonus = 0;
	}
	result.insert(result.begin(), value);
    }

    if ( bonus > 0 )
    {
	result.insert(result.begin(), bonus);
    }

    return result;
}

Int mult(Int a, Int b)
{
    Int result;
    Ints adds;

    int power = 0;

    while ( a[a.size() - 1] == 0 )
    {
	a.pop_back();
	++power;
    }

    while ( b[b.size() - 1] == 0 )
    {
	b.pop_back();
	++power;
    }

    int a_size = a.size();
    int b_size = b.size();


    int t = 0;
    for (int i = b_size - 1; i >= 0; i-- )
	adds.push_back( mult(a, b[i]) );

    result = add_multi( adds );

    while ( power-- > 0 )
	result.push_back(0);

    return result;
}

Int add_multi(Ints adds )
{
    Int result;

    int ads_size = adds.size();

    int idxs[ ads_size ];

    for (int i = 0; i < ads_size; i++ )
	idxs[i] = -i;

    int digits_count = adds[0].size();
    int bonus = 0;

    bool add;
    do
    {
	int digit_value = bonus;
	add = false;

	for (int j = 0; j < ads_size; j++ )
	{
	    if ( idxs[j] >= 0 && idxs[j] < adds[j].size() )
	    {
		digit_value += adds[j][adds[j].size() - idxs[j] - 1];
		add = true;
	    }

	    ++idxs[j];
	}
	if ( digit_value >= 0 )
	{
	    bonus = digit_value / 10;
	    digit_value = digit_value % 10;
	}
	else
	{
	    bonus = 0;
	}

	if ( digit_value > 0 || add )
	result.insert(result.begin(), digit_value);
    } while ( add );

    if ( bonus > 0 )
    {
	result.push_back(bonus);
    }

    return result;
}

Int add(Int a, Int b, int shift)
{
    Int result;

    int a_size = a.size() - 1;
    int b_size = b.size() - 1;

    while ( shift-- > 0 )
    {
	result.insert(result.begin(), a[a_size--]);
    }
    
    int bonus = 0;
    while (( a_size >= 0 ) && ( b_size >= 0 ))
    {
	int value = bonus + a[a_size--] + b[b_size--];
	if ( value >= 10 )
	{
	    value -= 10;
	    bonus = 1;
	}
	else
	{
	    bonus = 0;
	}

	result.insert(result.begin(), value);
    }

    while ( a_size >= 0 )
    {
	int value = a[a_size];

	if ( bonus == 1 )
	{
	    ++value;
	    bonus = 0;
	}

	if ( value >= 10 )
	{
	    value -= 10;
	    bonus = 1;
	}

	result.insert(result.begin(), value);
	--a_size;
    }

    while ( b_size >= 0 )
    {
	int value = b[b_size];

	if ( bonus == 1 )
	{
	    ++value;
	    bonus = 0;
	}

	if ( value >= 10 )
	{
	    value -= 10;
	    bonus = 1;
	}

	result.insert(result.begin(), value);
	--b_size;
    }

    if ( bonus > 0 )
	result.insert(result.begin(), 1);

    return result;
}

Int create( int n )
{
    Int result;
    do
    {
	result.insert(result.begin(), n % 10);
	n /= 10;
    } while ( n > 0 );
    return result;
}

ps редактор подсветки не подружился с первой кавычкой(заменил на ")

[Гость]

> for (int i = a_size - 1; i >= 0; i-- )
{
int value = bonus + a[i] * n;
if ( value >= 10 )
{
bonus = (int)(value / 10);
value %= 10;
}
else


> 10

То есть отнование системы счисления - 10? Это плохо.

Самое быстрое - основанием брать 2^32.
И в процессоре есть готовая операция, которая сразу получает (a+b) mod 2^32 и (a+b) div 2^32.
Правда, компилятор Борланда, например, никак не заставить это применить, насчёт компиляторов от Микрософта и Интела не знаю.

Если тебе надо, чтобы ещё и быстро выводилось, то пойди на компромисс - в качестве основния возьми 1000000000. Считаться будет медленнее, чем с 2^32, но в 9 раз быстрее, чем с 10. И, опять же, перемножить два 32-разрядных числа и поделить 64-разрядный результат на 1000000000 - это очень просто записывается на асме (тупо mul, div подряд, промежуточный 64-разрядный результат хранится в 2х регистрах, edx и eax), но невозможно заставить компилятор Борланда сделать это по-простому.

[volvo]

Кстати, Андрей, в твоем конкретном случае (при умножении на 1111), ты делаешь 4 раза одну и ту же работу, вот тут:

    for (int i = b_size - 1; i >= 0; i-- )
	 adds.push_back( mult(a, b[i]) );

Сделай чуть более хитро: опиши вектор из Int-ов, от 0 до 9, и когда приходишь в mult, проверяй, не делал ли ты уже этого умножения. Если делал - то сразу возвращай соотв. вектор. Не делал - делай, и запоминай, в следующий раз вернешь, когда понадобится. Так сэкономишь кучу времени. А если умножаешь на 1 - то можно вообще сразу возвращать исходный a, без холостых действий.

И в процессоре есть готовая операция, которая сразу получает (a+b) mod 2^32 и (a+b) div 2^32.

Как называется?

[-TarasBer-]

> Как называется?

Там опечатка, не +, а * надо.

[klem4]

Спасибо за советы, volvo - хорошее замечание, воспользуюсь. Но видимо без перевода в СС с большим основанием никак сильно не ускорить. Хотя тут придется реализовывать ф-ю деления большого на большое, для перевода большого в новую СС, попробую на днях.

Сообщение отредактировано: klem4 - 9.03.2011 22:53

[TarasBer]

> Хотя тут придется реализовывать ф-ю деления большого на большое

Надо только функцию, которая говорит частное и остаток деления a*b на 1000000000
Функция из 3 строчек на асме.

Для Д7 это выглядит так:

type TDivMod = record
  d, m: cardinal;
end;

function GetDivMod(a, b: cardinal): TDivMod;
asm
  mul edx
  mov ecx, 1000000000
  div ecx
  // эти две строчки из-за того, что результат почему-то возвращается не в edx:eax, а на стеке
  mov ebp - $18, eax
  mov ebp - $14, edx
end;

(Можно и через int64 написать, но это тормозить будет).

Ещё можно её модифицировать, чтобы она сразу для a*b+c говорила, только не забудь про случай переноса при сложении.

Сообщение отредактировано: TarasBer - 10.03.2011 11:21

[klem4]

Дык не получится так просто, в cardinal ты число из 5000 цифр не запишешь я думаю.

[-TarasBer-]

Это просто для замены

Код

int value = bonus + a[i] * n;
    if ( value >= 10 )
    {
        bonus = (int)(value / 10);
        value %= 10;
    }

На

Код

(value, bonus) = MulDiv(a[i], n);

Не надо основание системы счисления менять на 10^500, это для 1000000000.

А функция деления длинного на длинное тебе не нужна, ты ж не собрался, надеюсь, основание длинным делать.
Число храни как массив cardinal (или u32 или uint32, или что там), каждый элемент от 0 до 999999999. Ну просто чтобы было проще выводить результат. Иначе бы можно было бы тупо взять в качестве основания 2^32 и не париться.

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

[klem4]

Большое спасибо за подсказки, все получилось, использовал хранение в СС с основанием 10^8. Задача собственно заключалась в следующем: дано число, являющееся факториалом N, где 1 <= N <= 2000. Определить надо N.