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

#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;
    }
}

#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]);
}