Deprecated: Function set_magic_quotes_runtime() is deprecated in /DISK2/WWW/lokiware.info/mff/wakka.php on line 35 Matfiz : RychleNasobenie

Matfiz: RychleNasobenie

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;

}