Форум программистов, компьютерный форум CyberForum.ru

Умножение Карацубы - C++

Восстановить пароль Регистрация
 
Рейтинг: Рейтинг темы: голосов - 12, средняя оценка - 4.67
vlad_light
4 / 4 / 0
Регистрация: 24.09.2012
Сообщений: 178
17.05.2013, 20:26     Умножение Карацубы #1
Помогите, пожалуйста, довести до ума код:
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
#include <iostream>
 
typedef unsigned int        int32;
typedef unsigned long long  int64;
 
const int KARATSUBA_CUTOFF  = 2;
const int MAX_LENGTH        = 1024;
 
const int fit_with_zeros (int32* multiplier1, const int length1, int32* multiplier2, const int length2);
 
void KA (const int32* multiplier1, const int32* multiplier2, int32* product, const int length);
 
void GSA (const int32* multiplier1, const int32* multiplier2, int32* product, const int length);
 
int main ()
{
    int32 a[MAX_LENGTH] = {0xAAAAA1AA, 0xBBB3B5BB,0xABBAAAAA, 0xBBBB1BBB,0xAAAAA1AA, 0xBBB3B5BB,0xABBAAAAA, 0xBBBB1BBB};
    int32 b[MAX_LENGTH] = {0x11159111, 0x33333333,0x11111FF1, 0x33933333,0xAAAAA1AA, 0xBBB3B5BB,0xABBAAAAA, 0xBBBB1BBB};
    int32 c[6 * MAX_LENGTH], d[6 * MAX_LENGTH];
    memset (c, 0, 6 * MAX_LENGTH * sizeof (int32));
    const int a_length = 8;
    const int b_length = 8;
    const int max_length = fit_with_zeros (a, a_length, b, b_length);
 
    GSA (a, b, d, max_length);
 
    std::cout << std::endl << "Result S: ";
    for (int i = 0; i < 2 * max_length; ++i)
        std::cout << d[i] << " ";
 
    KA (a, b, c, max_length);
 
    std::cout << std::endl << "Result K: ";
    for (int i = 0; i < 2 * max_length; ++i)
        std::cout << c[i] << " ";
 
    std::cin.get ();
 
    return 0;
}
 
const int fit_with_zeros (int32* multiplier1, const int length1, int32* multiplier2, const int length2)
{
    int max_length = 1;
 
    while (max_length < std::max (length1, length2))
        max_length <<= 1;
 
    for (int i = length1; i < max_length; ++i)
        multiplier1[i] = 0;
 
    for (int i = length2; i < max_length; ++i)
        multiplier2[i] = 0;
 
    return max_length;
}
 
void KA (const int32* multiplier1, const int32* multiplier2, int32* product, const int length)
{
    if (length < KARATSUBA_CUTOFF)
    {
        GSA (multiplier1, multiplier2, product, length);
        return;
    }
 
    const   int32* multiplier1_low  = &multiplier1[0];
    const   int32* multiplier1_high = &multiplier1[length >> 1];
    const   int32* multiplier2_low  = &multiplier2[0];
    const   int32* multiplier2_high = &multiplier2[length >> 1];
 
            int32* multiplier1_high_low_sub = &product[length * 5];
            int32* multiplier2_low_high_sub = &product[length * 5 + length >> 1];   
 
    for (int i = 0; i < (length >> 1); ++i)
    {
        multiplier1_high_low_sub[i] = multiplier1_high[i] - multiplier1_low[i];
        multiplier2_low_high_sub[i] = multiplier2_low[i] - multiplier2_high[i];
    }
 
    int32* product_low  = &product[0];
    int32* product_high = &product[length];
    int32* product_mid  = &product[length << 1];
 
    KA (multiplier1_low, multiplier2_low, product_low, length >> 1);
    KA (multiplier1_high, multiplier2_high, product_high, length >> 1);
    KA (multiplier1_high_low_sub, multiplier2_low_high_sub, product_mid, length >> 1);
 
    int carry = 0;
    int64 current_sum;
 
    for (int i = 0; i < length; ++i)
    {
        current_sum = (int64) product_mid[i] + product_high[i] + product_low[i] + carry;
        product_mid[i] = current_sum;
        carry = current_sum >> 32;
    }
 
    product_high[length >> 1] += carry;
 
    carry = 0;
 
    for (int i = 0; i < length; ++i)
    {
        current_sum = (int64) product[i + (length >> 1)] + product_mid[i] + carry;
        product[i + (length >> 1)] = current_sum;
        carry = current_sum >> 32;
    }
 
    product[length + (length >> 1)] += carry;
}
 
void GSA (const int32* multiplier1, const int32* multiplier2, int32* product, const int length)
{
    memset (product, 0, length * sizeof (int32));
 
    for (int i = 0; i < length; ++i)
    {
        int32 carry = 0;
 
        for (int j = 0; j < length; ++j)
        {
            int64 current_product = (int64) multiplier1[i] * multiplier2[j] + product[i + j] + carry;
            product[i + j] = current_product;
            carry = current_product >> 32; 
        }
 
        product[i + length] = carry;
    }
}


Оргигинал брал отсюда:
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
#include <stdlib.h>
 
#include <stdio.h>
 
#include <time.h>
 
#include <iostream>
 
 
#define MAX_DIGITS 1024
 
 
#define KARAT_CUTOFF 4
 
void           
 karatsuba(const int *a, int *b, int *ret, const int d);
void   
 gradeSchool(const int *a, int *b, int *ret, int d);
void   
 doCarry(int *a, int d);
void      
 getNum(int *a, int *d_a);
void 
 printNum(int *a, int d);
 
int
main() {
    int             a[MAX_DIGITS]; // first multiplicand
    int             b[MAX_DIGITS]; // second multiplicand
    int             r[6 * MAX_DIGITS]; // result goes here
    int             d_a; // length of a
    int             d_b; // length of b
    int             d; // maximum length
    int             i; // counter
    clock_t         start; // for timing
    clock_t         stop; // for timing
 
    getNum(a, &d_a);
    getNum(b, &d_b);
 
    if(d_a < 0 || d_b < 0) {
        printf("0\n");
        exit(0);
        return 0;
    }
 
    i = (d_a > d_b) ? d_a : d_b;
    for(d = 1; d < i; d *= 2);
    for(i = d_a; i < d; i++) a[i] = 0;   // fill both a and b  with zeros
    for(i = d_b; i < d; i++) b[i] = 0;   // to nearest power of 2
 
    // do the trials, first for Karatsuba, then for grade-school.
    // For each trial, we print the result, followed by the time
    // taken per multiplication, followed by the number of
    // multiplications done. We do as many multiplications as we
    // can until we pass away an entire second.
    start = clock();
    stop = start + CLOCKS_PER_SEC;
    for(i = 0; clock() < stop; i++) {
        karatsuba(a, b, r, d); // compute product w/o regard to carry
        doCarry(r, 2 * d); // now do any carrying
    }
    start = clock() - start;
    printNum(r, 2 * d);
    printf(" %f ms (%d trials)\n", 1000 * (double) start / CLOCKS_PER_SEC / i, i);
 
    start = clock();
    stop = start + CLOCKS_PER_SEC;
    for(i = 0; clock() < stop; i++) {
        gradeSchool(a, b, r, d); // compute product in old way
        doCarry(r, 2 * d); // now do any carrying
    }
    start = clock() - start;
    printNum(r, 2 * d);
    printf(" %f ms (%d trials)\n", 1000 * (double) start / CLOCKS_PER_SEC / i, i);
 
    std::cin.get ();
}
 
// ret must have space for 6d digits.
// the result will be in only the first 2d digits
// my use of the space in ret is pretty creative.
// | ar*br  | al*bl  | asum*bsum | lower-recursion space | asum | bsum |
//  d digits d digits  d digits     3d digits              d/2    d/2
void
karatsuba(const int *a, int *b, int *ret, const int d) {
    int             i;
    int             *ar = &a[0]; // low-order half of a
    int             *al = &a[d/2]; // high-order half of a
    int             *br = &b[0]; // low-order half of b
    int             *bl = &b[d/2]; // high-order half of b
    int             *asum = &ret[d * 5]; // sum of a's halves
    int             *bsum = &ret[d * 5 + d/2]; // sum of b's halves
    int             *x1 = &ret[d * 0]; // ar*br's location
    int             *x2 = &ret[d * 1]; // al*bl's location
    int             *x3 = &ret[d * 2]; // asum*bsum's location
 
    // when d is small, we're better off just reverting to
    // grade-school multiplication, since it's faster at this point.
    if(d <= KARAT_CUTOFF) {
        gradeSchool(a, b, ret, d);
        return;
    }
 
    // compute asum and bsum
    for(i = 0; i < d / 2; i++) {
        asum[i] = al[i] + ar[i];
        bsum[i] = bl[i] + br[i];
    }
 
    // do recursive calls (I have to be careful about the order,
    // since the scratch space for the recursion on x1 includes
    // the space used for x2 and x3)
    karatsuba(ar, br, x1, d/2);
    karatsuba(al, bl, x2, d/2);
    karatsuba(asum, bsum, x3, d/2);
 
    // combine recursive steps
    for(i = 0; i < d; i++) x3[i] = x3[i] - x1[i] - x2[i];
    for(i = 0; i < d; i++) ret[i + d/2] += x3[i];
}
 
void
gradeSchool(const int *a, int *b, int *ret, int d) {
    int             i, j;
 
    for(i = 0; i < 2 * d; i++) ret[i] = 0;
    for(i = 0; i < d; i++) {
        for(j = 0; j < d; j++) ret[i + j] += a[i] * b[j];
    }
}
 
void
doCarry(int *a, int d) {
    int             c;
    int             i;
 
    c = 0;
    for(i = 0; i < d; i++) {
        a[i] += c;
        if(a[i] < 0) {
            c = -(-(a[i] + 1) / 10 + 1);
        } else {
            c = a[i] / 10;
        }
        a[i] -= c * 10;
    }
    if(c != 0) fprintf(stderr, "Overflow %d\n", c);
}
 
void
getNum(int *a, int *d_a) {
    int             c;
    int             i;
 
    *d_a = 0;
    while(true) {
        c = getchar();
        if(c == '\n' || c == EOF) break;
        if(*d_a >= MAX_DIGITS) {
            fprintf(stderr, "using only first %d digits\n", MAX_DIGITS);
            while(c != '\n' && c != EOF) c = getchar();
        }
        a[*d_a] = c - '0';
        ++(*d_a);
    }
    // reverse the number so that the 1's digit is first
    for(i = 0; i * 2 < *d_a - 1; i++) {
        c = a[i], a[i] = a[*d_a - i - 1], a[*d_a - i - 1] = c;
    }
}
 
void
printNum(int *a, int d) {
    int i;
    for(i = d - 1; i > 0; i--) if(a[i] != 0) break;
    for(; i >= 0; i--) printf("%d", a[i]);
}
Similar
Эксперт
41792 / 34177 / 6122
Регистрация: 12.04.2006
Сообщений: 57,940
17.05.2013, 20:26     Умножение Карацубы
Посмотрите здесь:

Умножение матриц C++
Умножение матриц C++
Умножение матриц C++
C++ умножение
Метод Карацубы умножения длинных чисел C++
После регистрации реклама в сообщениях будет скрыта и будут доступны все возможности форума.
ViktorKozlov
133 / 125 / 2
Регистрация: 13.12.2012
Сообщений: 293
17.05.2013, 21:25     Умножение Карацубы #2
А что не так?
taras atavin
Ушёл с форума.
 Аватар для taras atavin
3569 / 1752 / 91
Регистрация: 24.11.2009
Сообщений: 27,619
17.05.2013, 21:28     Умножение Карацубы #3
Задачу в студию.
vlad_light
4 / 4 / 0
Регистрация: 24.09.2012
Сообщений: 178
19.05.2013, 13:47  [ТС]     Умножение Карацубы #4
Задачу в студию.
Нужно реализовать алгоритм умножения Карацубы.
Я написал умножение "в столбик" (функция GSA (const int*,const int*, int*, const int)) и нашёл в интернете реализацию Алгоритма Карацубы (2-ой cut). Отличия в том, что я использую 32-ух разрядные числа (вместо 0...9) в одной ячейке, т.е. мне нужно подогнать алгоритм Карацубы под 32-ух разрядные числа.
А что не так?
Приведенный мною алгоритм выдаёт неправильный результат.
Основной упор делается на скорость.

Добавлено через 17 часов 39 минут
Помогите, пожалуйста, найти ошибку
Вродь алгоритм всё правильно должен считать, но, к сожалению, не считает
У меня уже мозг

Добавлено через 17 часов 41 минуту
Урааа!!!
Я нашел ошибку. В строках 76-77 и 93 несогласовано работают "+" и "-". Теперь подскажите, как их согласовать?
Если их поменять местами то в строках 76-77 может возникнуть переполнение на 1 бит, следовательно число не поместится в соответствующую ячейку в строке 86. Если их перемножить отдельно, то у нас добавится ещё одно умножение и столбик станет работать быстрее (поскольку использует меньше сложений и столько же умножений).

Не по теме:

На всякий случай напишу, чем отличаются "+" и "-":
http://www.cyberforum.ru/cgi-bin/latex.cgi?A\times B=(A_{low}+A_{high}x)\times (B_{low}+B_{high}x)=
http://www.cyberforum.ru/cgi-bin/latex.cgi?A_{low}\times B_{low}+A_{high}\times B_{high}x^2+((A_{low}+A_{high})\times (B_{low}+B_{high})-A_{low}\times B_{low}-A_{high}\times B_{high})\times x=
http://www.cyberforum.ru/cgi-bin/latex.cgi?A_{low}\times B_{low}+A_{high}\times B_{high}x^2+((A_{low}-A_{high})\times (B_{high}-B_{low})+A_{low}\times B_{low}+A_{high}\times B_{high})\times x

Yandex
Объявления
19.05.2013, 13:47     Умножение Карацубы
Ответ Создать тему
Опции темы

Текущее время: 14:55. Часовой пояс GMT +3.
КиберФорум - форум программистов, компьютерный форум, программирование
Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2016, vBulletin Solutions, Inc.
Рейтинг@Mail.ru