Rychle násobenie

Komentář K. Vandas : Nechápu, proč to dělat takhle. Pokud to byla jen 2 čísla, myslím, že se dá rozdělit 2. činitel na součty, se kterými vynásobíš vzájemně všechny součty prvního činitele. Je to sice hrozná věc, ale je správná a bez použití něčeho, na co musel pan F. přicházet déle než 3 hodiny :) (Např. 13 * 19 = ( 6 + 7 ) * ( 9 + 10 ) )... Je to kvůli zmenšení počtu cifer tak, aby se to do proměnných Pascalu už dostalo, btw. nevim, jak může (16 000)*(16 000) násobení a pak nějaké sčítání běžet 25 minut :) Sám jsem před chvilkou udělal násobení 50 čísel s 1000 ciframa, který mi běží dost rychle i na 333 MHz / 32 RAM :)
/Michal Demin :Asi som sa nevyjadril spravne,... ale nebolo to 16000 * 16000, ale cislo1*cislo2, kde strlen(cislo1)=strlen(cislo2)=16000
este taky malinky detail, rozdelit cisla na soucet, potom vynasobit, a znova spojit, tak pochybujem ze by si sa dostal na cas mensi ako je O(n^2), tymto algorimtomm som dosiahol rychlost (n log^2 n) co je celkom slusne. BTW, doporucolal by som si precitat knizku (aspon tak nahmatkovo) The Art of Computer Programming Vol. 2, je v nej celkom dopodrobna popisanych dost vela algoritmov. mne osobne sa algoritmus s fourierovou transofrmaciou pacil najviac, je to mozno len vec mojho nazoru. Jedina vec co ma zaraza, je ze pri nejakych vstupoch moja implementacia nefunguje

Piste mi komentare a nejasnosti toto som pisal velmi na rychlo a upravim ak bude treba
nechodte mi ale nadavat ze je to pomale pri 2*2, treba si precitat nejaku tu teoriu, tam sa pise ze tento algoritmus je rychly pri vecsich vstupoch

Pri písaní zápočtového programu som narazil na problém pri násobení. Násobenie 2 čísel s dĺžkou 16000 číslic trvalo 25minut. Toto bolo použitím tzv «násobenia na papieri».
Tak som začal pátrať po algoritme, ktorý by dokázal násobiť číslo v čase rýchlejšom ako O(n^2)
Narazil som na jeden algoritmus a to konkretne Násobenie pomocou fourierovej transformacie.

Nejdem sa moc zaoberať teóriou. Možete si ju nastudovat v nespočetne vela dokumentoch ako je napriklad file:fft.ps alebo knihe The Art of Computer Programming Vol. 2

Konkretna implementacia v C

Stiahni tu file:main.c, doporucujem ale zdrojak co je na stranke, sem tam este mozem upravit nejaky ten detail (bugy, atd)

K implementacii algoritmu som použil FFTW – the Fastest Fourier Transform in the West. Funkcia fft_long_mul má dva argumenty, dva char retazce zakončené '\0' (binárnou nulou). Tieto retazce možu byt teoreticky neobmedzenej velkosti.
Program si potom alokuje dostatok miesta v pamati. (O(n)) Spočíta sa fourierova transformícia oboch vstupov, vynísobia sa vístupy, a vypočíta sa opačná transformacia.
Jediný problem, ktorý môže nastať (a aj nastane) je, že FFTW si pri vytvarani planu robi nejake testy rýchlosti ( Tieto testy trvaju radovo niekolko násobne dlhsie ako celkový beh programu, pri malých vstupných argumentoch) Koli tomuto sa uklada wisdom file a ak existuje sa znova nacita. Dostane sa to zhruba na rychlost bc
Pri extremne velkych vstupoch to pre nejaku pricinu zblbne, tomu sa ale budem venovat neskor :)

Pri kompilovaní treba prilinkovat fftw3 kniznicu, a math knižnicu.


#include <stdlib.h>
#include <stdio.h>

#include <math.h>

#include <string.h>
#include <complex.h>
#include <fftw3.h>

#include <unistd.h>

#define MAX(x,y)        (x > y ? x : y)
#define uint            unsigned int
#define DEFAULT_WISDOM_FILE ".wisdom"


char * fft_long_mul(char *num1,char *num2);

int main(int argc, char *argv[]) {
    char * tmp;
    tmp = fft_long_mul("123456789123456789123456789123456789123456789123456789","123456789123456789123456789123456789123456789123456789987654321");
    printf("\n\n Out = %s \n",tmp);
    free(tmp);
    return 1;
} 

char * fft_long_mul(char *num1,char *num2) {
    
    uint Num1_size, Num2_size;
    fftw_complex *FFT_In1, *FFT_In2;
    fftw_complex *FFT_Out1, *FFT_Out2;
    fftw_plan FFT_Plan,FFT_Plan_Out;
    FILE *FFT_Wisdom;    
    char *Out_num;
    uint Num;
    uint i,cr,tmp;
// najskor alokujeme dostatocne velike polia ... polia musia byt velke aspon ako MAX(len(num1),len(num2)) * 2 ( a dokonca to cislo musi byt pow of 2) , 
    Num1_size = strlen(num1);
    Num2_size = strlen(num2);

    Num = 2;
    while (Num<Num1_size+Num2_size)
        Num *= 2;

//printf("%d \n",Num);

    // alokujem si potrebne miesto v pamati ...
    FFT_In1 = fftw_malloc(Num * sizeof(fftw_complex));
    FFT_In2 = fftw_malloc(Num * sizeof(fftw_complex));
    FFT_Out1 = fftw_malloc(Num * sizeof(fftw_complex));
    FFT_Out2 = fftw_malloc(Num * sizeof(fftw_complex));

    Out_num = malloc(Num+1);

    if (!(FFT_In1 && FFT_In2 && FFT_Out1 && FFT_Out2 && Out_num)) {
                fftw_free(FFT_In1);
                fftw_free(FFT_In2);
                fftw_free(FFT_Out1);
                fftw_free(FFT_Out2);
                free(Out_num);
                fprintf(stderr,"Error allocating memory\n");
                exit(-1);
    }

// if wisdom file exists, load it... 
    if ((FFT_Wisdom=fopen(DEFAULT_WISDOM_FILE,"r"))<=0) { 
//        fprintf(stderr,"FFT_Procesor: No wisdom file ... \n");  
    } else {    
        fftw_import_wisdom_from_file(FFT_Wisdom);
        fclose(FFT_Wisdom);
    }
// vygenerujem si plan k FFTW,      
    FFT_Plan = fftw_plan_dft_1d(Num, FFT_In1, FFT_Out1, FFTW_FORWARD, FFTW_ESTIMATE);//ESTIMATE MEASURE PATIENT 
    FFT_Plan_Out = fftw_plan_dft_1d(Num, FFT_Out1, FFT_In1, FFTW_BACKWARD, FFTW_ESTIMATE);//ESTIMATE MEASURE PATIENT 
//    fprintf(stderr,"FFT_Procesor: Done building plan ...\n");

// Save new wisdom file
    if ((FFT_Wisdom=fopen(DEFAULT_WISDOM_FILE,"w"))<0) { 
        fprintf(stderr,"Error saving wisdom file ... \n");  
    } else {
        fftw_export_wisdom_to_file(FFT_Wisdom);
        fclose(FFT_Wisdom);
    }
    
// naplnim si vstupne polia, mojimi cislami, ktore som dostal
    for (i=0; i<strlen(num1); i++) {
        FFT_In1[i] = (fftw_complex)(num1[Num1_size-i-1]-'0');
    }
    for (i=0; i<strlen(num2); i++) {
        FFT_In2[i] = (fftw_complex)(num2[Num2_size-i-1]-'0');
    }
            
// spracujem vstupy fourierovou transformaciou ...
    fftw_execute_dft(FFT_Plan, FFT_In1, FFT_Out1);
    fftw_execute_dft(FFT_Plan, FFT_In2, FFT_Out2);

//vynasobim jednotlive koeficienty z fourierovej transformacie
    for (i=0; i<Num; i++) {
        FFT_Out1[Num-i-1] *= FFT_Out2[Num-i-1];
    }
    

// spracujem FFT_Out1 fourierovou transformaciou naspet , a vysledok by mal byt nasobok num1 * num2 (po preusporiadani)
    fftw_execute_dft(FFT_Plan_Out, FFT_Out1, FFT_Out2);

// FFT_Out2 obsahuje moj vystup, ktory este treba podelit num

    for (i=0; i<Num; i++) {
        FFT_Out2[Num-i-1] /= (fftw_complex)Num;
    }


    cr = 0;
    for (i=0; i<Num; i++) {
        tmp = (int)nearbyint(creal(FFT_Out2[i])); // nearbyint je potrebne lebo pri double 16.0000000 != 16.00000000 a toto mi robilo celkom velike problemy
        tmp +=cr;
        Out_num[Num-i-1] = tmp%10 + '0';
        cr = tmp/10;

    }
/// zbavim sa nepotrebneho bordelu 
    fftw_free(FFT_In1);
    fftw_free(FFT_Out1);    
    fftw_free(FFT_In2);
    fftw_free(FFT_Out2);    
    fftw_destroy_plan(FFT_Plan);
    fftw_destroy_plan(FFT_Plan_Out);
    Out_num[Num] = '\0';
// vratime co je spravne
    return Out_num;
}