fork download 
  1. //Annotated Version
  2. #include<bits/stdc++.h>
  3. using namespace std;
  4.  
  5. #define PI (atan(1)*4)
  6. #define double long double
  7. #define PRECISION 0.000000001
  8.  
  9. complex<double>xx,yy;
  10. vector<int>seq=vector<int>{0};//butterfly sequence transform container
  11.  
  12. //generate the sequence(butterfly transform?) required to do iterative fft
  13. //butterfly transform essentially takes care of all the odd/even ordering in the fft
  14. //the ordering needs the resulting indices differences to be n/2k for the difference of two adjacent blocks of lengths k=2^a
  15. //so instead of thinking of it as sorting as reverse(~x) or whatever the "correct" formula is
  16. //we can recursively generate it, where each time we multiply the first half of the sequence by 2(since length of sequence is doubled, n/2k doubles -> differences also double)
  17. //then the second half is just the first half shifted by one, following the (2n)/(2*(n/2))=1 formula
  18. //note that given a sequence of length n and we need it for a polynomial of length n/2^k,
  19. //we can just divide the first n/2^m terms by 2^m, dividing all the differences to correctly mesh together(n/2k/2^m=(n/2^m)/2k)
  20. //we also know that the first n/2^m terms are divisble by 2^m, due to having sorted in the first half at least m times, gaining a factor of 2 each time it got sorted into the first half
  21. //this ends up being O(2*n) time, which i think is just as fast as the other methods but remains much more intuitive and cleaner
  22. //additionally you only really have to generate it "once"
  23. void seqgen(int n){
  24. 	if(seq.size()<n){seq.reserve(n);}
  25.     while(seq.size()<n){
  26.         //we double the sequence length in each iteration
  27.         int cur=seq.size();
  28.         for(register int a=0;a<cur;++a){
  29.             seq[a]*=2;//add the factor of two to the first number as described earlier
  30.             seq.push_back(seq[a]+1);//the second half of the sequence is just 1 + each corresponding number of the first half,
  31.             //as described earlier the two blocks are of length 2^(n-1) in wihch case n/(2*2^(n-1))=1 and is the valid difference
  32.             //no overlap because first half is even, second half is odd
  33.         }
  34.     }
  35. }
  36.  
  37. vector<complex<double>>cpy;//temp vector to hold contents of poly
  38.  
  39. //returns results of plugging in roots of unity without inverse
  40. //with inverse returns coefficients given results from roots of unity
  41. vector<complex<double>> fft(vector<complex<double>>poly,bool flip=false){//polynomial + whether or not is inverse fft
  42.     //perform preprocessing necessary for butterfly transform
  43.     //i think that preprocessing before ends up being faster than the standard butterfly transform by a slight amount for repeated uses of fft
  44.     seqgen(poly.size());//generate sequence
  45.     for(int a=0;a<poly.size();++a){
  46.         cpy.push_back(poly[((long long)seq[a]*(long long)poly.size())/seq.size()]);//references from earlier(line 16), so we only have to generate at least one sequence greater than or equal to poly length
  47.         //butterfly transform in linear time
  48.     }
  49.     cpy.swap(poly);cpy.clear();
  50.     for(register int l=2;l<=poly.size();l*=2){//block size/root of unity
  51.         complex<double>root=polar(1.0,2*PI/l);//lth root of unity = w
  52.         complex<double>running=1;//running root of unity
  53.         //merge two loops into one, hopefully more efficient?
  54.         for(register int a=0;a<poly.size();++a){
  55.             if(a%l==l/2){a+=l/2;}//we only really want to iterate from [kl,kl+l/2) for all k
  56.             if(a==poly.size()){break;}
  57.             if(!(a%l)){running=1;}//reset running every cycle of length l
  58.             //we use the thing here now
  59.             //A(a)=A^even(a^2)+a*A^odd(a^2) because we simply shift the terms together
  60.             //plug in w^x=a
  61.             //A(w^x)=A^even(w^2x)+(w^x)A^odd(w^2x)
  62.             //A_even(w^2x)=A_old(w_0^x) where w_0 is the root of unity from the previous iteration
  63.             // because w_0=w^2 b/c w_0 is e^iπ/(l/2) but w is e^iπ/l
  64.             //A(w^x)=A_even(w_0^x)+(w^x)A_odd(w_0^x)
  65.             // for numbers in [kl,kl+l/2) we can just iterate over
  66.             //for numbers in [kl+l/2,(k+1)l) we can do it at the same time as [kl,kl+l/2) with the formula A(w^x)=A_even(w_0^x)-(w^x)A_odd(w_0^x)
  67.             // because w^(kl+a)=-w^((k+1)l-a) due to them being conjugates(product is magnitude squared)
  68.             //so here x=A_even(w_0^x) and y=w^x A_odd(w_0^2x) where w^x is held as the value running
  69.             //when l=1 we simply multiply it by the 1st root of unity=1 so we can skip l=1
  70.              xx=poly[a];yy=running*poly[a+l/2];
  71.             poly[a]=xx+yy;poly[a+l/2]=xx-yy;
  72.             running*=root;
  73.         }
  74.     }
  75.     if(flip){
  76.         //inverse fft
  77.         //same as normal fft, but you add this at the end
  78.         //if you think of terms of matrices you need to find inverse of vandermonde matrix
  79.         //which functions the exact same way but you must reverse the nonconstant terms
  80.         //and you must divide by n because of things that happen when simplifying
  81.         reverse(poly.begin()+1,poly.end());
  82.         for(int a=0;a<poly.size();++a){
  83.             poly[a]/=poly.size();
  84.         }
  85.     }
  86.     return poly;
  87. }
  88.  
  89. vector<int>pmul(vector<int>p1,vector<int>p2){
  90.     int cap=p1.size()+p2.size()-1;//here because loj requires you to include leading 0's, is bullshit
  91.     int rn=p1.size()+p2.size()-1;//first find the total size
  92.     rn=(1<<(int)ceil(log2(rn)));//round to the first power of 2 above
  93.     //add cushioning zeroes
  94.     for(int a=p1.size();a<rn;++a){
  95.         p1.push_back(0);
  96.     }
  97.     for(int a=p2.size();a<rn;++a){
  98.         p2.push_back(0);
  99.     }
  100.     vector<complex<double>>f1,f2;//copy over p1 and p2 into a more appropriate datastructure
  101.     //i could use pointers in fft but fft has other uses i think and i dont want to ruin those
  102.     //so this is a bit more general?
  103.     for(int a=0;a<rn;++a){
  104.         f1.push_back(complex<double>(p1[a],0));
  105.         f2.push_back(complex<double>(p2[a],0));
  106.     }
  107.     f1=fft(f1),f2=fft(f2);//transformed point sets
  108.     //merge the two results together
  109.     //C(x)=A(x)B(x)
  110.     for(int a=0;a<rn;++a){
  111.         f1[a]*=f2[a];
  112.     }
  113.     //f1 now holds the point form of the product of p1 and p2
  114.     //so now we can just apply inverse fft
  115.     vector<complex<double>>ret=fft(f1,true);
  116.     //add the contents into an integer vector and prune the ending 0's off
  117.     vector<int>ans;
  118.     for(int a=0;a<ret.size();++a){
  119.         ans.push_back((int)round(ret[a].real()));//results should be integers for product of integer vector
  120.     }
  121.     //used to be prune off leading 0's, but loj needs those so uh
  122.     while(ans.size()>cap){
  123.         ans.pop_back();
  124.     }
  125.     return ans;
  126.  
  127. }
  128.  
Compilation error #stdin compilation error #stdout 0s 0KB
comments (?)
stdin 
Standard input is empty
compilation info 
/usr/bin/ld: /usr/lib/gcc/x86_64-linux-gnu/8/../../../x86_64-linux-gnu/Scrt1.o: in function `_start':
(.text+0x20): undefined reference to `main'
collect2: error: ld returned 1 exit status
stdout 
Standard output is empty
https://ideone.com/MLcdwn
language:
C++14 (gcc 8.3)
created:
2 years ago
visibility:
secret

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!