Deprecated: Function set_magic_quotes_runtime() is deprecated in /DISK2/WWW/lokiware.info/mff/wakka.php on line 35
#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;
}