# HG changeset patch # User Chris Cannam # Date 1447069607 0 # Node ID ebc87a62321d5ed5617818ea83eacf63fafc6037 # Parent 59a1ee198dca4cd6d3f458046cb0cf3df8c15f9a Add Nayuki fft.c compiled to JS diff -r 59a1ee198dca -r ebc87a62321d fft/index.html --- a/fft/index.html Mon Nov 09 11:46:35 2015 +0000 +++ b/fft/index.html Mon Nov 09 11:46:47 2015 +0000 @@ -12,6 +12,8 @@ + + @@ -38,6 +40,8 @@
>>16&4;E=E<>>16&2;i=14-(p|q|i)+(E<>>15)|0;i=o>>>(i+7|0)&1|i<<1}else i=0;b=c[360+(i<<2)>>2]|0;a:do if(!b){a=0;b=0;E=86}else{f=d;a=0;g=o<<((i|0)==31?0:25-(i>>>1)|0);h=b;b=0;while(1){e=c[h+4>>2]&-8;d=e-o|0;if(d>>>0
>>0)if((e|0)==(o|0)){a=h;b=h;E=90;break a}else b=h;else d=f;e=c[h+20>>2]|0;h=c[h+16+(g>>>31<<2)>>2]|0;a=(e|0)==0|(e|0)==(h|0)?a:e;e=(h|0)==0;if(e){E=86;break}else{f=d;g=g<<(e&1^1)}}}while(0);if((E|0)==86){if((a|0)==0&(b|0)==0){a=2<>>12&16;q=q>>>m;l=q>>>5&8;q=q>>>l;n=q>>>2&4;q=q>>>n;p=q>>>1&2;q=q>>>p;a=q>>>1&1;a=c[360+((l|m|n|p|a)+(q>>>a)<<2)>>2]|0}if(!a){i=d;j=b}else E=90}if((E|0)==90)while(1){E=0;q=(c[a+4>>2]&-8)-o|0;e=q>>>0 >>0;d=e?q:d;b=e?a:b;e=c[a+16>>2]|0;if(e){a=e;E=90;continue}a=c[a+20>>2]|0;if(!a){i=d;j=b;break}else E=90}if((j|0)!=0?i>>>0<((c[16]|0)-o|0)>>>0:0){f=c[18]|0;if(j>>>0 >>0)fa();h=j+o|0;if(j>>>0>=h>>>0)fa();g=c[j+24>>2]|0;d=c[j+12>>2]|0;do if((d|0)==(j|0)){b=j+20|0;a=c[b>>2]|0;if(!a){b=j+16|0;a=c[b>>2]|0;if(!a){s=0;break}}while(1){d=a+20|0;e=c[d>>2]|0;if(e){a=e;b=d;continue}d=a+16|0;e=c[d>>2]|0;if(!e)break;else{a=e;b=d}}if(b>>>0 >>0)fa();else{c[b>>2]=0;s=a;break}}else{e=c[j+8>>2]|0;if(e>>>0 >>0)fa();a=e+12|0;if((c[a>>2]|0)!=(j|0))fa();b=d+8|0;if((c[b>>2]|0)==(j|0)){c[a>>2]=d;c[b>>2]=e;s=d;break}else fa()}while(0);do if(g){a=c[j+28>>2]|0;b=360+(a<<2)|0;if((j|0)==(c[b>>2]|0)){c[b>>2]=s;if(!s){c[15]=c[15]&~(1<>>0<(c[18]|0)>>>0)fa();a=g+16|0;if((c[a>>2]|0)==(j|0))c[a>>2]=s;else c[g+20>>2]=s;if(!s)break}b=c[18]|0;if(s>>>0>>0)fa();c[s+24>>2]=g;a=c[j+16>>2]|0;do if(a)if(a>>>0>>0)fa();else{c[s+16>>2]=a;c[a+24>>2]=s;break}while(0);a=c[j+20>>2]|0;if(a)if(a>>>0<(c[18]|0)>>>0)fa();else{c[s+20>>2]=a;c[a+24>>2]=s;break}}while(0);do if(i>>>0>=16){c[j+4>>2]=o|3;c[h+4>>2]=i|1;c[h+i>>2]=i;a=i>>>3;if(i>>>0<256){d=96+(a<<1<<2)|0;b=c[14]|0;a=1<>2]|0;if(b>>>0<(c[18]|0)>>>0)fa();else{u=a;v=b}}else{c[14]=b|a;u=d+8|0;v=d}c[u>>2]=h;c[v+12>>2]=h;c[h+8>>2]=v;c[h+12>>2]=d;break}a=i>>>8;if(a)if(i>>>0>16777215)d=31;else{K=(a+1048320|0)>>>16&8;L=a< >>16&4;L=L< >>16&2;d=14-(J|K|d)+(L< >>15)|0;d=i>>>(d+7|0)&1|d<<1}else d=0;e=360+(d<<2)|0;c[h+28>>2]=d;a=h+16|0;c[a+4>>2]=0;c[a>>2]=0;a=c[15]|0;b=1< >2]=h;c[h+24>>2]=e;c[h+12>>2]=h;c[h+8>>2]=h;break}f=i<<((d|0)==31?0:25-(d>>>1)|0);a=c[e>>2]|0;while(1){if((c[a+4>>2]&-8|0)==(i|0)){d=a;E=148;break}b=a+16+(f>>>31<<2)|0;d=c[b>>2]|0;if(!d){E=145;break}else{f=f<<1;a=d}}if((E|0)==145)if(b>>>0<(c[18]|0)>>>0)fa();else{c[b>>2]=h;c[h+24>>2]=a;c[h+12>>2]=h;c[h+8>>2]=h;break}else if((E|0)==148){a=d+8|0;b=c[a>>2]|0;L=c[18]|0;if(b>>>0>=L>>>0&d>>>0>=L>>>0){c[b+12>>2]=h;c[a>>2]=h;c[h+8>>2]=b;c[h+12>>2]=d;c[h+24>>2]=0;break}else fa()}}else{L=i+o|0;c[j+4>>2]=L|3;L=j+L+4|0;c[L>>2]=c[L>>2]|1}while(0);L=j+8|0;return L|0}}}else o=-1;while(0);d=c[16]|0;if(d>>>0>=o>>>0){a=d-o|0;b=c[19]|0;if(a>>>0>15){L=b+o|0;c[19]=L;c[16]=a;c[L+4>>2]=a|1;c[L+a>>2]=a;c[b+4>>2]=o|3}else{c[16]=0;c[19]=0;c[b+4>>2]=d|3;L=b+d+4|0;c[L>>2]=c[L>>2]|1}L=b+8|0;return L|0}a=c[17]|0;if(a>>>0>o>>>0){J=a-o|0;c[17]=J;L=c[20]|0;K=L+o|0;c[20]=K;c[K+4>>2]=J|1;c[L+4>>2]=o|3;L=L+8|0;return L|0}do if(!(c[132]|0)){a=ka(30)|0;if(!(a+-1&a)){c[134]=a;c[133]=a;c[135]=-1;c[136]=-1;c[137]=0;c[125]=0;c[132]=(ia(0)|0)&-16^1431655768;break}else fa()}while(0);h=o+48|0;g=c[134]|0;i=o+47|0;f=g+i|0;g=0-g|0;j=f&g;if(j>>>0<=o>>>0){L=0;return L|0}a=c[124]|0;if((a|0)!=0?(u=c[122]|0,v=u+j|0,v>>>0<=u>>>0|v>>>0>a>>>0):0){L=0;return L|0}b:do if(!(c[125]&4)){a=c[20]|0;c:do if(a){d=504;while(1){b=c[d>>2]|0;if(b>>>0<=a>>>0?(r=d+4|0,(b+(c[r>>2]|0)|0)>>>0>a>>>0):0){e=d;d=r;break}d=c[d+8>>2]|0;if(!d){E=173;break c}}a=f-(c[17]|0)&g;if(a>>>0<2147483647){b=ha(a|0)|0;if((b|0)==((c[e>>2]|0)+(c[d>>2]|0)|0)){if((b|0)!=(-1|0)){h=b;f=a;E=193;break b}}else E=183}}else E=173;while(0);do if((E|0)==173?(t=ha(0)|0,(t|0)!=(-1|0)):0){a=t;b=c[133]|0;d=b+-1|0;if(!(d&a))a=j;else a=j-a+(d+a&0-b)|0;b=c[122]|0;d=b+a|0;if(a>>>0>o>>>0&a>>>0<2147483647){v=c[124]|0;if((v|0)!=0?d>>>0<=b>>>0|d>>>0>v>>>0:0)break;b=ha(a|0)|0;if((b|0)==(t|0)){h=t;f=a;E=193;break b}else E=183}}while(0);d:do if((E|0)==183){d=0-a|0;do if(h>>>0>a>>>0&(a>>>0<2147483647&(b|0)!=(-1|0))?(w=c[134]|0,w=i-a+w&0-w,w>>>0<2147483647):0)if((ha(w|0)|0)==(-1|0)){ha(d|0)|0;break d}else{a=w+a|0;break}while(0);if((b|0)!=(-1|0)){h=b;f=a;E=193;break b}}while(0);c[125]=c[125]|4;E=190}else E=190;while(0);if((((E|0)==190?j>>>0<2147483647:0)?(x=ha(j|0)|0,y=ha(0)|0,x>>>0 >>0&((x|0)!=(-1|0)&(y|0)!=(-1|0))):0)?(z=y-x|0,z>>>0>(o+40|0)>>>0):0){h=x;f=z;E=193}if((E|0)==193){a=(c[122]|0)+f|0;c[122]=a;if(a>>>0>(c[123]|0)>>>0)c[123]=a;i=c[20]|0;do if(i){e=504;do{a=c[e>>2]|0;b=e+4|0;d=c[b>>2]|0;if((h|0)==(a+d|0)){A=a;B=b;C=d;D=e;E=203;break}e=c[e+8>>2]|0}while((e|0)!=0);if(((E|0)==203?(c[D+12>>2]&8|0)==0:0)?i>>>0 >>0&i>>>0>=A>>>0:0){c[B>>2]=C+f;L=i+8|0;L=(L&7|0)==0?0:0-L&7;K=i+L|0;L=f-L+(c[17]|0)|0;c[20]=K;c[17]=L;c[K+4>>2]=L|1;c[K+L+4>>2]=40;c[21]=c[136];break}a=c[18]|0;if(h>>>0>>0){c[18]=h;j=h}else j=a;d=h+f|0;a=504;while(1){if((c[a>>2]|0)==(d|0)){b=a;E=211;break}a=c[a+8>>2]|0;if(!a){b=504;break}}if((E|0)==211)if(!(c[a+12>>2]&8)){c[b>>2]=h;l=a+4|0;c[l>>2]=(c[l>>2]|0)+f;l=h+8|0;l=h+((l&7|0)==0?0:0-l&7)|0;a=d+8|0;a=d+((a&7|0)==0?0:0-a&7)|0;k=l+o|0;g=a-l-o|0;c[l+4>>2]=o|3;do if((a|0)!=(i|0)){if((a|0)==(c[19]|0)){L=(c[16]|0)+g|0;c[16]=L;c[19]=k;c[k+4>>2]=L|1;c[k+L>>2]=L;break}b=c[a+4>>2]|0;if((b&3|0)==1){i=b&-8;f=b>>>3;e:do if(b>>>0>=256){h=c[a+24>>2]|0;e=c[a+12>>2]|0;do if((e|0)==(a|0)){d=a+16|0;e=d+4|0;b=c[e>>2]|0;if(!b){b=c[d>>2]|0;if(!b){J=0;break}}else d=e;while(1){e=b+20|0;f=c[e>>2]|0;if(f){b=f;d=e;continue}e=b+16|0;f=c[e>>2]|0;if(!f)break;else{b=f;d=e}}if(d>>>0 >>0)fa();else{c[d>>2]=0;J=b;break}}else{f=c[a+8>>2]|0;if(f>>>0 >>0)fa();b=f+12|0;if((c[b>>2]|0)!=(a|0))fa();d=e+8|0;if((c[d>>2]|0)==(a|0)){c[b>>2]=e;c[d>>2]=f;J=e;break}else fa()}while(0);if(!h)break;b=c[a+28>>2]|0;d=360+(b<<2)|0;do if((a|0)!=(c[d>>2]|0)){if(h>>>0<(c[18]|0)>>>0)fa();b=h+16|0;if((c[b>>2]|0)==(a|0))c[b>>2]=J;else c[h+20>>2]=J;if(!J)break e}else{c[d>>2]=J;if(J)break;c[15]=c[15]&~(1<>>0 >>0)fa();c[J+24>>2]=h;b=a+16|0;d=c[b>>2]|0;do if(d)if(d>>>0 >>0)fa();else{c[J+16>>2]=d;c[d+24>>2]=J;break}while(0);b=c[b+4>>2]|0;if(!b)break;if(b>>>0<(c[18]|0)>>>0)fa();else{c[J+20>>2]=b;c[b+24>>2]=J;break}}else{d=c[a+8>>2]|0;e=c[a+12>>2]|0;b=96+(f<<1<<2)|0;do if((d|0)!=(b|0)){if(d>>>0 >>0)fa();if((c[d+12>>2]|0)==(a|0))break;fa()}while(0);if((e|0)==(d|0)){c[14]=c[14]&~(1< >>0 >>0)fa();b=e+8|0;if((c[b>>2]|0)==(a|0)){G=b;break}fa()}while(0);c[d+12>>2]=e;c[G>>2]=d}while(0);a=a+i|0;g=i+g|0}a=a+4|0;c[a>>2]=c[a>>2]&-2;c[k+4>>2]=g|1;c[k+g>>2]=g;a=g>>>3;if(g>>>0<256){d=96+(a<<1<<2)|0;b=c[14]|0;a=1<>2]|0;if(b>>>0>=(c[18]|0)>>>0){K=a;L=b;break}fa()}while(0);c[K>>2]=k;c[L+12>>2]=k;c[k+8>>2]=L;c[k+12>>2]=d;break}a=g>>>8;do if(!a)d=0;else{if(g>>>0>16777215){d=31;break}K=(a+1048320|0)>>>16&8;L=a< >>16&4;L=L< >>16&2;d=14-(J|K|d)+(L< >>15)|0;d=g>>>(d+7|0)&1|d<<1}while(0);e=360+(d<<2)|0;c[k+28>>2]=d;a=k+16|0;c[a+4>>2]=0;c[a>>2]=0;a=c[15]|0;b=1< >2]=k;c[k+24>>2]=e;c[k+12>>2]=k;c[k+8>>2]=k;break}f=g<<((d|0)==31?0:25-(d>>>1)|0);a=c[e>>2]|0;while(1){if((c[a+4>>2]&-8|0)==(g|0)){d=a;E=281;break}b=a+16+(f>>>31<<2)|0;d=c[b>>2]|0;if(!d){E=278;break}else{f=f<<1;a=d}}if((E|0)==278)if(b>>>0<(c[18]|0)>>>0)fa();else{c[b>>2]=k;c[k+24>>2]=a;c[k+12>>2]=k;c[k+8>>2]=k;break}else if((E|0)==281){a=d+8|0;b=c[a>>2]|0;L=c[18]|0;if(b>>>0>=L>>>0&d>>>0>=L>>>0){c[b+12>>2]=k;c[a>>2]=k;c[k+8>>2]=b;c[k+12>>2]=d;c[k+24>>2]=0;break}else fa()}}else{L=(c[17]|0)+g|0;c[17]=L;c[20]=k;c[k+4>>2]=L|1}while(0);L=l+8|0;return L|0}else b=504;while(1){a=c[b>>2]|0;if(a>>>0<=i>>>0?(F=a+(c[b+4>>2]|0)|0,F>>>0>i>>>0):0){b=F;break}b=c[b+8>>2]|0}g=b+-47|0;d=g+8|0;d=g+((d&7|0)==0?0:0-d&7)|0;g=i+16|0;d=d>>>0 >>0?i:d;a=d+8|0;e=h+8|0;e=(e&7|0)==0?0:0-e&7;L=h+e|0;e=f+-40-e|0;c[20]=L;c[17]=e;c[L+4>>2]=e|1;c[L+e+4>>2]=40;c[21]=c[136];e=d+4|0;c[e>>2]=27;c[a>>2]=c[126];c[a+4>>2]=c[127];c[a+8>>2]=c[128];c[a+12>>2]=c[129];c[126]=h;c[127]=f;c[129]=0;c[128]=a;a=d+24|0;do{a=a+4|0;c[a>>2]=7}while((a+4|0)>>>0>>0);if((d|0)!=(i|0)){h=d-i|0;c[e>>2]=c[e>>2]&-2;c[i+4>>2]=h|1;c[d>>2]=h;a=h>>>3;if(h>>>0<256){d=96+(a<<1<<2)|0;b=c[14]|0;a=1<>2]|0;if(b>>>0<(c[18]|0)>>>0)fa();else{H=a;I=b}}else{c[14]=b|a;H=d+8|0;I=d}c[H>>2]=i;c[I+12>>2]=i;c[i+8>>2]=I;c[i+12>>2]=d;break}a=h>>>8;if(a)if(h>>>0>16777215)d=31;else{K=(a+1048320|0)>>>16&8;L=a< >>16&4;L=L< >>16&2;d=14-(J|K|d)+(L< >>15)|0;d=h>>>(d+7|0)&1|d<<1}else d=0;f=360+(d<<2)|0;c[i+28>>2]=d;c[i+20>>2]=0;c[g>>2]=0;a=c[15]|0;b=1< >2]=i;c[i+24>>2]=f;c[i+12>>2]=i;c[i+8>>2]=i;break}e=h<<((d|0)==31?0:25-(d>>>1)|0);a=c[f>>2]|0;while(1){if((c[a+4>>2]&-8|0)==(h|0)){d=a;E=307;break}b=a+16+(e>>>31<<2)|0;d=c[b>>2]|0;if(!d){E=304;break}else{e=e<<1;a=d}}if((E|0)==304)if(b>>>0<(c[18]|0)>>>0)fa();else{c[b>>2]=i;c[i+24>>2]=a;c[i+12>>2]=i;c[i+8>>2]=i;break}else if((E|0)==307){a=d+8|0;b=c[a>>2]|0;L=c[18]|0;if(b>>>0>=L>>>0&d>>>0>=L>>>0){c[b+12>>2]=i;c[a>>2]=i;c[i+8>>2]=b;c[i+12>>2]=d;c[i+24>>2]=0;break}else fa()}}}else{L=c[18]|0;if((L|0)==0|h>>>0 >>0)c[18]=h;c[126]=h;c[127]=f;c[129]=0;c[23]=c[132];c[22]=-1;a=0;do{L=96+(a<<1<<2)|0;c[L+12>>2]=L;c[L+8>>2]=L;a=a+1|0}while((a|0)!=32);L=h+8|0;L=(L&7|0)==0?0:0-L&7;K=h+L|0;L=f+-40-L|0;c[20]=K;c[17]=L;c[K+4>>2]=L|1;c[K+L+4>>2]=40;c[21]=c[136]}while(0);a=c[17]|0;if(a>>>0>o>>>0){J=a-o|0;c[17]=J;L=c[20]|0;K=L+o|0;c[20]=K;c[K+4>>2]=J|1;c[L+4>>2]=o|3;L=L+8|0;return L|0}}c[(ya()|0)>>2]=12;L=0;return L|0}function Aa(a){a=a|0;var b=0,d=0,e=0,f=0,g=0,h=0,i=0,j=0,k=0,l=0,m=0,n=0,o=0,p=0,q=0;if(!a)return;d=a+-8|0;h=c[18]|0;if(d>>>0 >>0)fa();a=c[a+-4>>2]|0;b=a&3;if((b|0)==1)fa();e=a&-8;m=d+e|0;do if(!(a&1)){a=c[d>>2]|0;if(!b)return;k=d+(0-a)|0;j=a+e|0;if(k>>>0 >>0)fa();if((k|0)==(c[19]|0)){a=m+4|0;b=c[a>>2]|0;if((b&3|0)!=3){q=k;g=j;break}c[16]=j;c[a>>2]=b&-2;c[k+4>>2]=j|1;c[k+j>>2]=j;return}e=a>>>3;if(a>>>0<256){b=c[k+8>>2]|0;d=c[k+12>>2]|0;a=96+(e<<1<<2)|0;if((b|0)!=(a|0)){if(b>>>0 >>0)fa();if((c[b+12>>2]|0)!=(k|0))fa()}if((d|0)==(b|0)){c[14]=c[14]&~(1< >>0 >>0)fa();a=d+8|0;if((c[a>>2]|0)==(k|0))f=a;else fa()}else f=d+8|0;c[b+12>>2]=d;c[f>>2]=b;q=k;g=j;break}f=c[k+24>>2]|0;d=c[k+12>>2]|0;do if((d|0)==(k|0)){b=k+16|0;d=b+4|0;a=c[d>>2]|0;if(!a){a=c[b>>2]|0;if(!a){i=0;break}}else b=d;while(1){d=a+20|0;e=c[d>>2]|0;if(e){a=e;b=d;continue}d=a+16|0;e=c[d>>2]|0;if(!e)break;else{a=e;b=d}}if(b>>>0 >>0)fa();else{c[b>>2]=0;i=a;break}}else{e=c[k+8>>2]|0;if(e>>>0 >>0)fa();a=e+12|0;if((c[a>>2]|0)!=(k|0))fa();b=d+8|0;if((c[b>>2]|0)==(k|0)){c[a>>2]=d;c[b>>2]=e;i=d;break}else fa()}while(0);if(f){a=c[k+28>>2]|0;b=360+(a<<2)|0;if((k|0)==(c[b>>2]|0)){c[b>>2]=i;if(!i){c[15]=c[15]&~(1<>>0<(c[18]|0)>>>0)fa();a=f+16|0;if((c[a>>2]|0)==(k|0))c[a>>2]=i;else c[f+20>>2]=i;if(!i){q=k;g=j;break}}d=c[18]|0;if(i>>>0 >>0)fa();c[i+24>>2]=f;a=k+16|0;b=c[a>>2]|0;do if(b)if(b>>>0 >>0)fa();else{c[i+16>>2]=b;c[b+24>>2]=i;break}while(0);a=c[a+4>>2]|0;if(a)if(a>>>0<(c[18]|0)>>>0)fa();else{c[i+20>>2]=a;c[a+24>>2]=i;q=k;g=j;break}else{q=k;g=j}}else{q=k;g=j}}else{q=d;g=e}while(0);if(q>>>0>=m>>>0)fa();a=m+4|0;b=c[a>>2]|0;if(!(b&1))fa();if(!(b&2)){if((m|0)==(c[20]|0)){p=(c[17]|0)+g|0;c[17]=p;c[20]=q;c[q+4>>2]=p|1;if((q|0)!=(c[19]|0))return;c[19]=0;c[16]=0;return}if((m|0)==(c[19]|0)){p=(c[16]|0)+g|0;c[16]=p;c[19]=q;c[q+4>>2]=p|1;c[q+p>>2]=p;return}g=(b&-8)+g|0;e=b>>>3;do if(b>>>0>=256){f=c[m+24>>2]|0;a=c[m+12>>2]|0;do if((a|0)==(m|0)){b=m+16|0;d=b+4|0;a=c[d>>2]|0;if(!a){a=c[b>>2]|0;if(!a){n=0;break}}else b=d;while(1){d=a+20|0;e=c[d>>2]|0;if(e){a=e;b=d;continue}d=a+16|0;e=c[d>>2]|0;if(!e)break;else{a=e;b=d}}if(b>>>0<(c[18]|0)>>>0)fa();else{c[b>>2]=0;n=a;break}}else{b=c[m+8>>2]|0;if(b>>>0<(c[18]|0)>>>0)fa();d=b+12|0;if((c[d>>2]|0)!=(m|0))fa();e=a+8|0;if((c[e>>2]|0)==(m|0)){c[d>>2]=a;c[e>>2]=b;n=a;break}else fa()}while(0);if(f){a=c[m+28>>2]|0;b=360+(a<<2)|0;if((m|0)==(c[b>>2]|0)){c[b>>2]=n;if(!n){c[15]=c[15]&~(1<>>0<(c[18]|0)>>>0)fa();a=f+16|0;if((c[a>>2]|0)==(m|0))c[a>>2]=n;else c[f+20>>2]=n;if(!n)break}d=c[18]|0;if(n>>>0 >>0)fa();c[n+24>>2]=f;a=m+16|0;b=c[a>>2]|0;do if(b)if(b>>>0 >>0)fa();else{c[n+16>>2]=b;c[b+24>>2]=n;break}while(0);a=c[a+4>>2]|0;if(a)if(a>>>0<(c[18]|0)>>>0)fa();else{c[n+20>>2]=a;c[a+24>>2]=n;break}}}else{b=c[m+8>>2]|0;d=c[m+12>>2]|0;a=96+(e<<1<<2)|0;if((b|0)!=(a|0)){if(b>>>0<(c[18]|0)>>>0)fa();if((c[b+12>>2]|0)!=(m|0))fa()}if((d|0)==(b|0)){c[14]=c[14]&~(1< >>0<(c[18]|0)>>>0)fa();a=d+8|0;if((c[a>>2]|0)==(m|0))l=a;else fa()}else l=d+8|0;c[b+12>>2]=d;c[l>>2]=b}while(0);c[q+4>>2]=g|1;c[q+g>>2]=g;if((q|0)==(c[19]|0)){c[16]=g;return}}else{c[a>>2]=b&-2;c[q+4>>2]=g|1;c[q+g>>2]=g}a=g>>>3;if(g>>>0<256){d=96+(a<<1<<2)|0;b=c[14]|0;a=1<>2]|0;if(b>>>0<(c[18]|0)>>>0)fa();else{o=a;p=b}}else{c[14]=b|a;o=d+8|0;p=d}c[o>>2]=q;c[p+12>>2]=q;c[q+8>>2]=p;c[q+12>>2]=d;return}a=g>>>8;if(a)if(g>>>0>16777215)d=31;else{o=(a+1048320|0)>>>16&8;p=a< >>16&4;p=p< >>16&2;d=14-(n|o|d)+(p< >>15)|0;d=g>>>(d+7|0)&1|d<<1}else d=0;e=360+(d<<2)|0;c[q+28>>2]=d;c[q+20>>2]=0;c[q+16>>2]=0;a=c[15]|0;b=1< >>1)|0);a=c[e>>2]|0;while(1){if((c[a+4>>2]&-8|0)==(g|0)){d=a;e=130;break}b=a+16+(f>>>31<<2)|0;d=c[b>>2]|0;if(!d){e=127;break}else{f=f<<1;a=d}}if((e|0)==127)if(b>>>0<(c[18]|0)>>>0)fa();else{c[b>>2]=q;c[q+24>>2]=a;c[q+12>>2]=q;c[q+8>>2]=q;break}else if((e|0)==130){a=d+8|0;b=c[a>>2]|0;p=c[18]|0;if(b>>>0>=p>>>0&d>>>0>=p>>>0){c[b+12>>2]=q;c[a>>2]=q;c[q+8>>2]=b;c[q+12>>2]=d;c[q+24>>2]=0;break}else fa()}}else{c[15]=a|b;c[e>>2]=q;c[q+24>>2]=e;c[q+12>>2]=q;c[q+8>>2]=q}while(0);q=(c[22]|0)+-1|0;c[22]=q;if(!q)a=512;else return;while(1){a=c[a>>2]|0;if(!a)break;else a=a+8|0}c[22]=-1;return}function Ba(){}function Ca(b,d,e){b=b|0;d=d|0;e=e|0;var f=0,g=0,h=0,i=0;f=b+e|0;if((e|0)>=20){d=d&255;h=b&3;i=d|d<<8|d<<16|d<<24;g=f&~3;if(h){h=b+4-h|0;while((b|0)<(h|0)){a[b>>0]=d;b=b+1|0}}while((b|0)<(g|0)){c[b>>2]=i;b=b+4|0}}while((b|0)<(f|0)){a[b>>0]=d;b=b+1|0}return b-e|0}function Da(b,d,e){b=b|0;d=d|0;e=e|0;var f=0;if((e|0)>=4096)return ja(b|0,d|0,e|0)|0;f=b|0;if((b&3)==(d&3)){while(b&3){if(!e)return f|0;a[b>>0]=a[d>>0]|0;b=b+1|0;d=d+1|0;e=e-1|0}while((e|0)>=4){c[b>>2]=c[d>>2];b=b+4|0;d=d+4|0;e=e-4|0}}while((e|0)>0){a[b>>0]=a[d>>0]|0;b=b+1|0;d=d+1|0;e=e-1|0}return f|0} + +// EMSCRIPTEN_END_FUNCS +return{_free:Aa,_memset:Ca,_malloc:za,_precalc:va,_dispose:wa,_memcpy:Da,_transform_radix2_precalc:xa,runPostSets:Ba,stackAlloc:ma,stackSave:na,stackRestore:oa,establishStackSpace:pa,setThrew:qa,setTempRet0:ta,getTempRet0:ua}}) + + +// EMSCRIPTEN_END_ASM +(Module.asmGlobalArg,Module.asmLibraryArg,buffer);var _free=Module["_free"]=asm["_free"];var runPostSets=Module["runPostSets"]=asm["runPostSets"];var _memset=Module["_memset"]=asm["_memset"];var _malloc=Module["_malloc"]=asm["_malloc"];var _precalc=Module["_precalc"]=asm["_precalc"];var _dispose=Module["_dispose"]=asm["_dispose"];var _memcpy=Module["_memcpy"]=asm["_memcpy"];var _transform_radix2_precalc=Module["_transform_radix2_precalc"]=asm["_transform_radix2_precalc"];Runtime.stackAlloc=asm["stackAlloc"];Runtime.stackSave=asm["stackSave"];Runtime.stackRestore=asm["stackRestore"];Runtime.establishStackSpace=asm["establishStackSpace"];Runtime.setTempRet0=asm["setTempRet0"];Runtime.getTempRet0=asm["getTempRet0"];function ExitStatus(status){this.name="ExitStatus";this.message="Program terminated with exit("+status+")";this.status=status}ExitStatus.prototype=new Error;ExitStatus.prototype.constructor=ExitStatus;var initialStackTop;var preloadStartTime=null;var calledMain=false;dependenciesFulfilled=function runCaller(){if(!Module["calledRun"])run();if(!Module["calledRun"])dependenciesFulfilled=runCaller};Module["callMain"]=Module.callMain=function callMain(args){assert(runDependencies==0,"cannot call main when async dependencies remain! (listen on __ATMAIN__)");assert(__ATPRERUN__.length==0,"cannot call main when preRun functions remain to be called");args=args||[];ensureInitRuntime();var argc=args.length+1;function pad(){for(var i=0;i<4-1;i++){argv.push(0)}}var argv=[allocate(intArrayFromString(Module["thisProgram"]),"i8",ALLOC_NORMAL)];pad();for(var i=0;i 0){return}preRun();if(runDependencies>0)return;if(Module["calledRun"])return;function doRun(){if(Module["calledRun"])return;Module["calledRun"]=true;if(ABORT)return;ensureInitRuntime();preMain();if(Module["onRuntimeInitialized"])Module["onRuntimeInitialized"]();if(Module["_main"]&&shouldRunNow)Module["callMain"](args);postRun()}if(Module["setStatus"]){Module["setStatus"]("Running...");setTimeout((function(){setTimeout((function(){Module["setStatus"]("")}),1);doRun()}),1)}else{doRun()}}Module["run"]=Module.run=run;function exit(status,implicit){if(implicit&&Module["noExitRuntime"]){return}if(Module["noExitRuntime"]){}else{ABORT=true;EXITSTATUS=status;STACKTOP=initialStackTop;exitRuntime();if(Module["onExit"])Module["onExit"](status)}if(ENVIRONMENT_IS_NODE){process["stdout"]["once"]("drain",(function(){process["exit"](status)}));console.log(" ");setTimeout((function(){process["exit"](status)}),500)}else if(ENVIRONMENT_IS_SHELL&&typeof quit==="function"){quit(status)}throw new ExitStatus(status)}Module["exit"]=Module.exit=exit;var abortDecorators=[];function abort(what){if(what!==undefined){Module.print(what);Module.printErr(what);what=JSON.stringify(what)}else{what=""}ABORT=true;EXITSTATUS=1;var extra="\nIf this abort() is unexpected, build with -s ASSERTIONS=1 which can give more information.";var output="abort("+what+") at "+stackTrace()+extra;if(abortDecorators){abortDecorators.forEach((function(decorator){output=decorator(output,what)}))}throw output}Module["abort"]=Module.abort=abort;if(Module["preInit"]){if(typeof Module["preInit"]=="function")Module["preInit"]=[Module["preInit"]];while(Module["preInit"].length>0){Module["preInit"].pop()()}}var shouldRunNow=true;if(Module["noInitialRun"]){shouldRunNow=false}run() + + + + + + return Module; +}; diff -r 59a1ee198dca -r ebc87a62321d fft/nayukic/fft.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/fft/nayukic/fft.c Mon Nov 09 11:46:47 2015 +0000 @@ -0,0 +1,374 @@ +/* + * Free FFT and convolution (C) + * + * Copyright (c) 2014 Project Nayuki + * http://www.nayuki.io/page/free-small-fft-in-multiple-languages + * + * (MIT License) + * Permission is hereby granted, free of charge, to any person obtaining a copy of + * this software and associated documentation files (the "Software"), to deal in + * the Software without restriction, including without limitation the rights to + * use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of + * the Software, and to permit persons to whom the Software is furnished to do so, + * subject to the following conditions: + * - The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * - The Software is provided "as is", without warranty of any kind, express or + * implied, including but not limited to the warranties of merchantability, + * fitness for a particular purpose and noninfringement. In no event shall the + * authors or copyright holders be liable for any claim, damages or other + * liability, whether in an action of contract, tort or otherwise, arising from, + * out of or in connection with the Software or the use or other dealings in the + * Software. + */ + +#include +#include +#include +#include +#include "fft.h" + + +// Private function prototypes +static size_t reverse_bits(size_t x, unsigned int n); +static void *memdup(const void *src, size_t n); + +#define SIZE_MAX ((size_t)-1) + + +int transform(double real[], double imag[], size_t n) { + if (n == 0) + return 1; + else if ((n & (n - 1)) == 0) // Is power of 2 + return transform_radix2(real, imag, n); + else // More complicated algorithm for arbitrary sizes + return transform_bluestein(real, imag, n); +} + + +int inverse_transform(double real[], double imag[], size_t n) { + return transform(imag, real, n); +} + +tables *precalc(size_t n) { + unsigned int levels; + // Compute levels = floor(log2(n)) + { + size_t temp = n; + levels = 0; + while (temp > 1) { + levels++; + temp >>= 1; + } + if (1u << levels != n) + return 0; // n is not a power of 2 + } + if (SIZE_MAX / sizeof(double) < n / 2) return 0; + tables *tables = malloc(sizeof(tables)); + if (!tables) return tables; + tables->levels = levels; + size_t size = (n / 2) * sizeof(double); + tables->cos = malloc(size); + if (!tables->cos) { + free(tables); + return 0; + } + tables->sin = malloc(size); + if (!tables->sin) { + free(tables->cos); + free(tables); + return 0; + } + int i; + for (i = 0; i < n / 2; i++) { + tables->cos[i] = cos(2 * M_PI * i / n); + tables->sin[i] = sin(2 * M_PI * i / n); + } + return tables; +} + +void dispose(tables *tables) { + if (!tables) return; + free(tables->cos); + free(tables->sin); + free(tables); +} + +void transform_radix2_precalc(double real[], double imag[], int n, tables *tables) { + double *cos_table, *sin_table; + int size; + int i; + + // Trignometric tables + cos_table = tables->cos; + sin_table = tables->sin; + + // Bit-reversed addressing permutation + for (i = 0; i < n; i++) { + int j = reverse_bits(i, tables->levels); + if (j > i) { + double temp = real[i]; + real[i] = real[j]; + real[j] = temp; + temp = imag[i]; + imag[i] = imag[j]; + imag[j] = temp; + } + } + + // Cooley-Tukey decimation-in-time radix-2 FFT + for (size = 2; size <= n; size *= 2) { + int halfsize = size / 2; + int tablestep = n / size; + for (i = 0; i < n; i += size) { + int j; + int k; + for (j = i, k = 0; j < i + halfsize; j++, k += tablestep) { + double tpre = real[j+halfsize] * cos_table[k] + imag[j+halfsize] * sin_table[k]; + double tpim = -real[j+halfsize] * sin_table[k] + imag[j+halfsize] * cos_table[k]; + real[j + halfsize] = real[j] - tpre; + imag[j + halfsize] = imag[j] - tpim; + real[j] += tpre; + imag[j] += tpim; + } + } + if (size == n) // Prevent overflow in 'size *= 2' + break; + } +} + +int transform_radix2(double real[], double imag[], size_t n) { + // Variables + int status = 0; + unsigned int levels; + double *cos_table, *sin_table; + size_t size; + size_t i; + + // Compute levels = floor(log2(n)) + { + size_t temp = n; + levels = 0; + while (temp > 1) { + levels++; + temp >>= 1; + } + if (1u << levels != n) + return 0; // n is not a power of 2 + } + + // Trignometric tables + if (SIZE_MAX / sizeof(double) < n / 2) + return 0; + size = (n / 2) * sizeof(double); + cos_table = malloc(size); + sin_table = malloc(size); + if (cos_table == NULL || sin_table == NULL) + goto cleanup; + for (i = 0; i < n / 2; i++) { + cos_table[i] = cos(2 * M_PI * i / n); + sin_table[i] = sin(2 * M_PI * i / n); + } + + // Bit-reversed addressing permutation + for (i = 0; i < n; i++) { + size_t j = reverse_bits(i, levels); + if (j > i) { + double temp = real[i]; + real[i] = real[j]; + real[j] = temp; + temp = imag[i]; + imag[i] = imag[j]; + imag[j] = temp; + } + } + + // Cooley-Tukey decimation-in-time radix-2 FFT + for (size = 2; size <= n; size *= 2) { + size_t halfsize = size / 2; + size_t tablestep = n / size; + for (i = 0; i < n; i += size) { + size_t j; + size_t k; + for (j = i, k = 0; j < i + halfsize; j++, k += tablestep) { + double tpre = real[j+halfsize] * cos_table[k] + imag[j+halfsize] * sin_table[k]; + double tpim = -real[j+halfsize] * sin_table[k] + imag[j+halfsize] * cos_table[k]; + real[j + halfsize] = real[j] - tpre; + imag[j + halfsize] = imag[j] - tpim; + real[j] += tpre; + imag[j] += tpim; + } + } + if (size == n) // Prevent overflow in 'size *= 2' + break; + } + status = 1; + +cleanup: + free(sin_table); + free(cos_table); + return status; +} + + +int transform_bluestein(double real[], double imag[], size_t n) { + // Variables + int status = 0; + double *cos_table, *sin_table; + double *areal, *aimag; + double *breal, *bimag; + double *creal, *cimag; + size_t m; + size_t size_n, size_m; + size_t i; + + // Find a power-of-2 convolution length m such that m >= n * 2 + 1 + { + size_t target; + if (n > (SIZE_MAX - 1) / 2) + return 0; + target = n * 2 + 1; + for (m = 1; m < target; m *= 2) { + if (SIZE_MAX / 2 < m) + return 0; + } + } + + // Allocate memory + if (SIZE_MAX / sizeof(double) < n || SIZE_MAX / sizeof(double) < m) + return 0; + size_n = n * sizeof(double); + size_m = m * sizeof(double); + cos_table = malloc(size_n); + sin_table = malloc(size_n); + areal = calloc(m, sizeof(double)); + aimag = calloc(m, sizeof(double)); + breal = calloc(m, sizeof(double)); + bimag = calloc(m, sizeof(double)); + creal = malloc(size_m); + cimag = malloc(size_m); + if (cos_table == NULL || sin_table == NULL + || areal == NULL || aimag == NULL + || breal == NULL || bimag == NULL + || creal == NULL || cimag == NULL) + goto cleanup; + + // Trignometric tables + for (i = 0; i < n; i++) { + double temp = M_PI * (size_t)((unsigned long long)i * i % ((unsigned long long)n * 2)) / n; + // Less accurate version if long long is unavailable: double temp = M_PI * i * i / n; + cos_table[i] = cos(temp); + sin_table[i] = sin(temp); + } + + // Temporary vectors and preprocessing + for (i = 0; i < n; i++) { + areal[i] = real[i] * cos_table[i] + imag[i] * sin_table[i]; + aimag[i] = -real[i] * sin_table[i] + imag[i] * cos_table[i]; + } + breal[0] = cos_table[0]; + bimag[0] = sin_table[0]; + for (i = 1; i < n; i++) { + breal[i] = breal[m - i] = cos_table[i]; + bimag[i] = bimag[m - i] = sin_table[i]; + } + + // Convolution + if (!convolve_complex(areal, aimag, breal, bimag, creal, cimag, m)) + goto cleanup; + + // Postprocessing + for (i = 0; i < n; i++) { + real[i] = creal[i] * cos_table[i] + cimag[i] * sin_table[i]; + imag[i] = -creal[i] * sin_table[i] + cimag[i] * cos_table[i]; + } + status = 1; + + // Deallocation +cleanup: + free(cimag); + free(creal); + free(bimag); + free(breal); + free(aimag); + free(areal); + free(sin_table); + free(cos_table); + return status; +} + + +int convolve_real(const double x[], const double y[], double out[], size_t n) { + double *ximag, *yimag, *zimag; + int status = 0; + ximag = calloc(n, sizeof(double)); + yimag = calloc(n, sizeof(double)); + zimag = calloc(n, sizeof(double)); + if (ximag == NULL || yimag == NULL || zimag == NULL) + goto cleanup; + + status = convolve_complex(x, ximag, y, yimag, out, zimag, n); +cleanup: + free(zimag); + free(yimag); + free(ximag); + return status; +} + + +int convolve_complex(const double xreal[], const double ximag[], const double yreal[], const double yimag[], double outreal[], double outimag[], size_t n) { + int status = 0; + size_t size; + size_t i; + double *xr, *xi, *yr, *yi; + if (SIZE_MAX / sizeof(double) < n) + return 0; + size = n * sizeof(double); + xr = memdup(xreal, size); + xi = memdup(ximag, size); + yr = memdup(yreal, size); + yi = memdup(yimag, size); + if (xr == NULL || xi == NULL || yr == NULL || yi == NULL) + goto cleanup; + + if (!transform(xr, xi, n)) + goto cleanup; + if (!transform(yr, yi, n)) + goto cleanup; + for (i = 0; i < n; i++) { + double temp = xr[i] * yr[i] - xi[i] * yi[i]; + xi[i] = xi[i] * yr[i] + xr[i] * yi[i]; + xr[i] = temp; + } + if (!inverse_transform(xr, xi, n)) + goto cleanup; + for (i = 0; i < n; i++) { // Scaling (because this FFT implementation omits it) + outreal[i] = xr[i] / n; + outimag[i] = xi[i] / n; + } + status = 1; + +cleanup: + free(yi); + free(yr); + free(xi); + free(xr); + return status; +} + + +static size_t reverse_bits(size_t x, unsigned int n) { + size_t result = 0; + unsigned int i; + for (i = 0; i < n; i++, x >>= 1) + result = (result << 1) | (x & 1); + return result; +} + + +static void *memdup(const void *src, size_t n) { + void *dest = malloc(n); + if (dest != NULL) + memcpy(dest, src, n); + return dest; +} diff -r 59a1ee198dca -r ebc87a62321d fft/nayukic/fft.h --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/fft/nayukic/fft.h Mon Nov 09 11:46:47 2015 +0000 @@ -0,0 +1,74 @@ +/* + * Free FFT and convolution (C) + * + * Copyright (c) 2014 Project Nayuki + * http://www.nayuki.io/page/free-small-fft-in-multiple-languages + * + * (MIT License) + * Permission is hereby granted, free of charge, to any person obtaining a copy of + * this software and associated documentation files (the "Software"), to deal in + * the Software without restriction, including without limitation the rights to + * use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of + * the Software, and to permit persons to whom the Software is furnished to do so, + * subject to the following conditions: + * - The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * - The Software is provided "as is", without warranty of any kind, express or + * implied, including but not limited to the warranties of merchantability, + * fitness for a particular purpose and noninfringement. In no event shall the + * authors or copyright holders be liable for any claim, damages or other + * liability, whether in an action of contract, tort or otherwise, arising from, + * out of or in connection with the Software or the use or other dealings in the + * Software. + */ + + +/* + * Computes the discrete Fourier transform (DFT) of the given complex vector, storing the result back into the vector. + * The vector can have any length. This is a wrapper function. Returns 1 (true) if successful, 0 (false) otherwise (out of memory). + */ +int transform(double real[], double imag[], size_t n); + +/* + * Computes the inverse discrete Fourier transform (IDFT) of the given complex vector, storing the result back into the vector. + * The vector can have any length. This is a wrapper function. This transform does not perform scaling, so the inverse is not a true inverse. + * Returns 1 (true) if successful, 0 (false) otherwise (out of memory). + */ +int inverse_transform(double real[], double imag[], size_t n); + +/* + * Computes the discrete Fourier transform (DFT) of the given complex vector, storing the result back into the vector. + * The vector's length must be a power of 2. Uses the Cooley-Tukey decimation-in-time radix-2 algorithm. + * Returns 1 (true) if successful, 0 (false) otherwise (n is not a power of 2, or out of memory). + */ +int transform_radix2(double real[], double imag[], size_t n); + +typedef struct { + double *cos; + double *sin; + int levels; +} tables; + +tables *precalc(size_t n); +void dispose(tables *); +void transform_radix2_precalc(double real[], double imag[], int n, tables *tables); + +/* + * Computes the discrete Fourier transform (DFT) of the given complex vector, storing the result back into the vector. + * The vector can have any length. This requires the convolution function, which in turn requires the radix-2 FFT function. + * Uses Bluestein's chirp z-transform algorithm. Returns 1 (true) if successful, 0 (false) otherwise (out of memory). + */ +int transform_bluestein(double real[], double imag[], size_t n); + +/* + * Computes the circular convolution of the given real vectors. Each vector's length must be the same. + * Returns 1 (true) if successful, 0 (false) otherwise (out of memory). + */ +int convolve_real(const double x[], const double y[], double out[], size_t n); + +/* + * Computes the circular convolution of the given complex vectors. Each vector's length must be the same. + * Returns 1 (true) if successful, 0 (false) otherwise (out of memory). + */ +int convolve_complex(const double xreal[], const double ximag[], const double yreal[], const double yimag[], double outreal[], double outimag[], size_t n); + diff -r 59a1ee198dca -r ebc87a62321d fft/nayukic/ffttest.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/fft/nayukic/ffttest.c Mon Nov 09 11:46:47 2015 +0000 @@ -0,0 +1,218 @@ +/* + * FFT and convolution test (C) + * + * Copyright (c) 2014 Project Nayuki + * http://www.nayuki.io/page/free-small-fft-in-multiple-languages + * + * (MIT License) + * Permission is hereby granted, free of charge, to any person obtaining a copy of + * this software and associated documentation files (the "Software"), to deal in + * the Software without restriction, including without limitation the rights to + * use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of + * the Software, and to permit persons to whom the Software is furnished to do so, + * subject to the following conditions: + * - The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * - The Software is provided "as is", without warranty of any kind, express or + * implied, including but not limited to the warranties of merchantability, + * fitness for a particular purpose and noninfringement. In no event shall the + * authors or copyright holders be liable for any claim, damages or other + * liability, whether in an action of contract, tort or otherwise, arising from, + * out of or in connection with the Software or the use or other dealings in the + * Software. + */ + +#include +#include +#include +#include +#include +#include "fft.h" + + +// Private function prototypes +static void test_fft(int n); +static void test_convolution(int n); +static void naive_dft(const double *inreal, const double *inimag, double *outreal, double *outimag, int inverse, int n); +static void naive_convolve(const double *xreal, const double *ximag, const double *yreal, const double *yimag, double *outreal, double *outimag, int n); +static double log10_rms_err(const double *xreal, const double *ximag, const double *yreal, const double *yimag, int n); +static double *random_reals(int n); +static void *memdup(const void *src, size_t n); + +static double max_log_error = -INFINITY; + + +/* Main and test functions */ + +int main(int argc, char **argv) { + int i; + int prev; + srand(time(NULL)); + + // Test power-of-2 size FFTs + for (i = 0; i <= 12; i++) + test_fft(1 << i); + + // Test small size FFTs + for (i = 0; i < 30; i++) + test_fft(i); + + // Test diverse size FFTs + prev = 0; + for (i = 0; i <= 100; i++) { + int n = (int)lround(pow(1500, i / 100.0)); + if (n > prev) { + test_fft(n); + prev = n; + } + } + + // Test power-of-2 size convolutions + for (i = 0; i <= 12; i++) + test_convolution(1 << i); + + // Test diverse size convolutions + prev = 0; + for (i = 0; i <= 100; i++) { + int n = (int)lround(pow(1500, i / 100.0)); + if (n > prev) { + test_convolution(n); + prev = n; + } + } + + printf("\n"); + printf("Max log err = %.1f\n", max_log_error); + printf("Test %s\n", max_log_error < -10 ? "passed" : "failed"); + return 0; +} + + +static void test_fft(int n) { + double *inputreal, *inputimag; + double *refoutreal, *refoutimag; + double *actualoutreal, *actualoutimag; + + inputreal = random_reals(n); + inputimag = random_reals(n); + + refoutreal = malloc(n * sizeof(double)); + refoutimag = malloc(n * sizeof(double)); + naive_dft(inputreal, inputimag, refoutreal, refoutimag, 0, n); + + actualoutreal = memdup(inputreal, n * sizeof(double)); + actualoutimag = memdup(inputimag, n * sizeof(double)); + transform(actualoutreal, actualoutimag, n); + + printf("fftsize=%4d logerr=%5.1f\n", n, log10_rms_err(refoutreal, refoutimag, actualoutreal, actualoutimag, n)); + + free(inputreal); + free(inputimag); + free(refoutreal); + free(refoutimag); + free(actualoutreal); + free(actualoutimag); +} + + +static void test_convolution(int n) { + double *input0real, *input0imag; + double *input1real, *input1imag; + double *refoutreal, *refoutimag; + double *actualoutreal, *actualoutimag; + + input0real = random_reals(n); + input0imag = random_reals(n); + input1real = random_reals(n); + input1imag = random_reals(n); + + refoutreal = malloc(n * sizeof(double)); + refoutimag = malloc(n * sizeof(double)); + naive_convolve(input0real, input0imag, input1real, input1imag, refoutreal, refoutimag, n); + + actualoutreal = malloc(n * sizeof(double)); + actualoutimag = malloc(n * sizeof(double)); + convolve_complex(input0real, input0imag, input1real, input1imag, actualoutreal, actualoutimag, n); + + printf("convsize=%4d logerr=%5.1f\n", n, log10_rms_err(refoutreal, refoutimag, actualoutreal, actualoutimag, n)); + + free(input0real); + free(input0imag); + free(input1real); + free(input1imag); + free(refoutreal); + free(refoutimag); + free(actualoutreal); + free(actualoutimag); +} + + +/* Naive reference computation functions */ + +static void naive_dft(const double *inreal, const double *inimag, double *outreal, double *outimag, int inverse, int n) { + double coef = (inverse ? 2 : -2) * M_PI; + int k; + for (k = 0; k < n; k++) { // For each output element + double sumreal = 0; + double sumimag = 0; + int t; + for (t = 0; t < n; t++) { // For each input element + double angle = coef * ((long long)t * k % n) / n; + sumreal += inreal[t]*cos(angle) - inimag[t]*sin(angle); + sumimag += inreal[t]*sin(angle) + inimag[t]*cos(angle); + } + outreal[k] = sumreal; + outimag[k] = sumimag; + } +} + + +static void naive_convolve(const double *xreal, const double *ximag, const double *yreal, const double *yimag, double *outreal, double *outimag, int n) { + int i; + for (i = 0; i < n; i++) { + double sumreal = 0; + double sumimag = 0; + int j; + for (j = 0; j < n; j++) { + int k = (i - j + n) % n; + sumreal += xreal[k] * yreal[j] - ximag[k] * yimag[j]; + sumimag += xreal[k] * yimag[j] + ximag[k] * yreal[j]; + } + outreal[i] = sumreal; + outimag[i] = sumimag; + } +} + + +/* Utility functions */ + +static double log10_rms_err(const double *xreal, const double *ximag, const double *yreal, const double *yimag, int n) { + double err = 0; + int i; + for (i = 0; i < n; i++) + err += (xreal[i] - yreal[i]) * (xreal[i] - yreal[i]) + (ximag[i] - yimag[i]) * (ximag[i] - yimag[i]); + + err /= n > 0 ? n : 1; + err = sqrt(err); // Now this is a root mean square (RMS) error + err = err > 0 ? log10(err) : -99.0; + if (err > max_log_error) + max_log_error = err; + return err; +} + + +static double *random_reals(int n) { + double *result = malloc(n * sizeof(double)); + int i; + for (i = 0; i < n; i++) + result[i] = (rand() / (RAND_MAX + 1.0)) * 2 - 1; + return result; +} + + +static void *memdup(const void *src, size_t n) { + void *dest = malloc(n); + if (dest != NULL) + memcpy(dest, src, n); + return dest; +} diff -r 59a1ee198dca -r ebc87a62321d fft/test.js --- a/fft/test.js Mon Nov 09 11:46:35 2015 +0000 +++ b/fft/test.js Mon Nov 09 11:46:47 2015 +0000 @@ -103,6 +103,35 @@ report("nayukiobj", start, middle, end, total); } +function testNayukiC(size) { + + var fft = new FFTNayukiC(size); + + var start = performance.now(); + var middle = start; + var end = start; + + total = 0.0; + + for (var i = 0; i < 2*iterations; ++i) { + if (i == iterations) { + middle = performance.now(); + } + var real = inputReals(size); + var imag = zeroReals(size); + fft.forward(real, imag); + for (var j = 0; j < size; ++j) { + total += Math.sqrt(real[j] * real[j] + imag[j] * imag[j]); + } + } + + var end = performance.now(); + + fft.dispose(); + + report("nayukic", start, middle, end, total); +} + function testNockert(size) { var fft = new FFT.complex(size, false); @@ -252,7 +281,8 @@ } var sizes = [ 512, 2048 ]; -var tests = [ testNayuki, testNayukiObj, testKissFFT, testCross, testFFTW, testNockert, testDntj ]; +var tests = [ testNayuki, testNayukiObj, testNayukiC, testKissFFT, + testCross, testFFTW, testNockert, testDntj ]; var nextTest = 0; var nextSize = 0; var interval;