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