Deprecated: Function set_magic_quotes_runtime() is deprecated in /DISK2/WWW/lokiware.info/mff/wakka.php on line 35 Matfiz : Rychle Nasobenie
Přihlášení:  Heslo:  
Matfiz: RychleNasobenie ...
Hlavní Stránka | Seznam Stránek | Poslední Změny | Poslední Komentované | Uživatelé | Registrace |
Toto je stará verze stránky RychleNasobenie z 2006-12-18 20:30:23..

Rychle násobenie

Piste mi komentare a nejasnosti toto som pisal velmi na rychlo a upravim ak bude treba
!!(blue) 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


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)


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);

    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);
    
// vratime co je spravne
    return Out_num;
}


 
Na stránce nejsou žádné soubory. [Zobrazit soubory (formulář)]
Na stránce nejsou žádné komentáře. [Zobrazit komentáře (formulář)]