From 89b8c368ea62218ae429633aede551c84972abce Mon Sep 17 00:00:00 2001 From: Adam Welc <adam.welc@oracle.com> Date: Thu, 19 Jun 2014 16:57:31 -0700 Subject: [PATCH] Implemented support for "fft" function via C code extracted from GNU R. --- com.oracle.truffle.r.native/Makefile | 12 +- .../lib/darwin/libRDerived.dylib | Bin 0 -> 25064 bytes com.oracle.truffle.r.native/src/fft.c | 880 ++++++++++++++++++ .../r/nodes/builtin/base/FFTFunctions.java | 61 -- .../nodes/builtin/base/ForeignFunctions.java | 84 ++ .../truffle/r/nodes/builtin/stats/R/fft.R | 47 + .../truffle/r/nodes/builtin/stats/R/init.R | 23 + .../r/nodes/builtin/stats/R/package-info.java | 28 + .../src/com/oracle/truffle/r/runtime/DCF.java | 1 - .../r/runtime/data/RComplexVector.java | 8 + .../truffle/r/runtime/ffi/RDerivedRFFI.java | 2 +- .../r/runtime/ffi/jnr/JNR_RFFIFactory.java | 6 +- 12 files changed, 1084 insertions(+), 68 deletions(-) create mode 100755 com.oracle.truffle.r.native/lib/darwin/libRDerived.dylib create mode 100644 com.oracle.truffle.r.native/src/fft.c delete mode 100644 com.oracle.truffle.r.nodes/src/com/oracle/truffle/r/nodes/builtin/base/FFTFunctions.java create mode 100644 com.oracle.truffle.r.nodes/src/com/oracle/truffle/r/nodes/builtin/stats/R/fft.R create mode 100644 com.oracle.truffle.r.nodes/src/com/oracle/truffle/r/nodes/builtin/stats/R/init.R create mode 100644 com.oracle.truffle.r.nodes/src/com/oracle/truffle/r/nodes/builtin/stats/R/package-info.java diff --git a/com.oracle.truffle.r.native/Makefile b/com.oracle.truffle.r.native/Makefile index 027569e390..fcc82719fd 100644 --- a/com.oracle.truffle.r.native/Makefile +++ b/com.oracle.truffle.r.native/Makefile @@ -20,9 +20,17 @@ # or visit www.oracle.com if you need additional information or have any # questions. # -# A placeholder to keep mx happy -# The only native code at this stage is in the form of binary Lapack/Blas libraries copied from GnuR all: +ifneq ($(shell uname), Darwin) + gcc -fPIC -shared -o ./lib/linux/libRDerived.so ./src/fft.c +else + gcc -fPIC -dynamiclib -o ./lib/darwin/libRDerived.dylib ./src/fft.c +endif clean: +ifneq ($(shell uname), Darwin) + rm -f ./lib/linux/libRDerived.* +else + rm -f ./lib/darwin/libRDerived.* +endif diff --git a/com.oracle.truffle.r.native/lib/darwin/libRDerived.dylib b/com.oracle.truffle.r.native/lib/darwin/libRDerived.dylib new file mode 100755 index 0000000000000000000000000000000000000000..8d58a619d4c489fe40556023a711d7e24dc88d08 GIT binary patch literal 25064 zcmX^A>+L^w1_nlE1_lN;1_lOx1_p)*RtAPv1_1^TkYr$BxWUN45Fa1n8W92#LBj#q z`Roh~46INLQV}0tl3Gy$VnHyvd7nWRGcYjxwS&mR_z)fg0|N^~1DwqOa!+wdX;Lv< z2;IC7eh?`h2*m*6GeX%QiUVpM$j=~td_2g!_>9z?g480g`KbQ=12s<tq?mz$0mNs5 zvcVM8JTM!qKRG|I7~x@b_bK>8R6;_TfdRy4fe2wG<Kt5^;?s%}b5kMG=;j$f%{u{+ zW&rU)aSG*PP@s720LQC~V~8V?G%P-04nyOE;uT~r%&qb9c{%aLmAOgzIq?N0MGW!r z*v*3(#}L540MZY$5fuK&;TRvElZeke0jPUm=7IRg=7HP?a~sHx`1qvaVp9+cg3;Zl z05#7BA`K=%=@!gjU|>K_Gaxarc}5@>1f!c512s<sA`K>y{fBNIC>}h(>BYy>&)e0- z6PCvkz*-rg*^q$&%6|bh0OU_76RePzfq}t-fq_8+&CPoHIhjfNDTzhpnRy^)kV|S& zW?5>AUP@(7W)cH~5(7hn0my^|1_llW1_q7<h#4Ru5LSTlEfN?QG#D5dKz;$)u_1wh z0hG5o5+P3GW?%qCEl3Lk$Za4DQp?GJ8sA_ANJK#*q})Lg0yCwdemVD{H89uP^uoU5 z)7(~yaa(}$A2{wgxEL5XxEQn;5J3de2yzoF96<InFs}rMh<<5t5y<}_99&sklA5dM z1PcS0J_o3sV13}2XJpXNFDTK^$xlwq0qF$c^yFkH6<@4p2s1+s>W&0xIBG)q8PIss zM&g6w8)WAw9u0xf5Eu=C(GVC7fzc2c4S~@R7!3hBgg~fA=TncyHwlak3?7|(PyGA; z-=lNufq(!1cTU{{B71lI`~TmivvmuY>fP|~|Nk8zXL)qGZh$dXz!)7cM#J&ekpKVx zgB3Qf1)0j>(YzL9EQd#Pt^$bA01*Zt!U9A%uy{1*da!i2f^0bWL%#En@l+7AbE?h% z|Nk{xL1Z_W<X_(kc3G>#|NsBtLLm2bwp#rE|KDSWgsw}-Z^f7P9-XcWJi1*Mc<cf> zgn{9SeVPY<-T@GC>Jxw5!88y4x&t1atrZ|$-CznNci_bnTLuQ0qyf6*StLoY&K;|K zB3@K&dkwO0Li1jbs~PHDJs<{t;ul~&YRkaj!N2Q=J=`wW2_W&c2yxd5AcukkXCMhK z=mz_!v(*Rg)$XYvZ+BkQJlJ`nvsDKy4)S&9p@YBVzwZUbR%dI@|NsAww`%<V|6hTD zfuS4B;$III?rc>7OO`;R=Xh%dl-2Dj(H$xPG1+5=N3TdLDAYW<T|GP)FLX0>wtD>k z|G&9bf}zf!+f@Qng+Mo0h2<fB-$M|^uh}QOW|{Dsal$8lEwJSt;Pk=3px^)s=M$g! zBTs;W`2<KF<VBF^g-`r(7eJyHJUUw?z=1y%OhNTN07v2>kW!EVkVpjaE_~vT3|#=$ z3F01jp<;uGMFFrutsn|4`H4T`fJZl+4>Aej6pv1@phq{@Zjf3GbI({K%;f+FX2Add z|3Url=2lR!GIoP0%n)K|{=vlGYR3hN5|DR04;_3Z&v@bD4_E|)1e*^qb+*QU;}yj1 zo(eYnLT76PSPVphO+E^VomNoh@?gB+(aqq|>H47?tPDi~DA&OiT)?6LV#1;BR*+2~ ztvsk^pYZ4g<%YZipZH^3KY#-GfJY}27K0!*VKWFK12gEvC;k{;a2%fi8w9gg1l3+_ z20`2kGw8x6{)kYpi!NX{2x2cbgCH_6gARP+k8lOM2;wzJ$b-U~@j@p!Tp*!i4bE}= zeF^_SMGq+VX`;x2G9G`Q6Id4HD{FAh<?pXkXJCLO{tw^;3<`r@XmYGO;nC^3!=t+u zl$ty`U1xNI696K49O?$=TFXQHeup|+UH<?7|8f!o1H%Ng6zmH#1C(qJeBzIEo$-l3 z;-E*j>kg2}g-`r3U_V{}r4*=jAPKNikdnFs9^F%+j><a(HUK0Ka_fmSehrXGGe82M zaEF^T!=u~v1IPof87Fju{fnMvK_LN3S7`AIO0(TlL4oX%e95Di=PD#cPld$sDUWU; zkX^9iItG%{K~9In1&WD~=)`3rL>6Wu$mJm4Ap+w978Aj)=;dL;B?ZcI0*Fjtf|_<e zK=aHgkQYF~1!92WnGuV5VCP}C3M_@oDu^V^Dp-*i0ShQbgzHd4zL$pyDNRJea!_|D zIC8=94h|zoynFCN-3^W@P{jcAD#ZRn9^I~Bji5r*_XEsiaH6`9#;<Yl6MqEAxH?E4 z#htZ46&rug&cC3n1!6;z5u`$b6n`MIL3Lz_8N5Pu-2w4s_f$}s2~!M8agczv1{a9@ z{fS^>LHUyrCC7p?6Mx^PKOjp%x*+On5+Dr=h^dT+AUsf60I~&^%0X<XjqsBEP&c@y z11SX=2jv~G23yMC^ZLjC|1UrNhS}Eqf)QGAgIxBBUkjAl>p(;vC^UOP6(6X;1=US? zNEH#NkbnjIffuHx@Cpf3j6h=JzzY>5;Tez;q<bpJzaamB-G1=kOL@i<7k?lZES7)% z{|D9CFi(Pf0rCw<156CW28n@E1SFVIG9x(afD$dVN&*>yq5zU`K)KbU)Aa(Rv_nw< zHvuV!;i#1`lD+E>EC#}D#AYBw5>^*n0%ciHg#b&mt~YS{5tQ~p^$#evB7A><qY7_L z<k1bTW<YTn5qblp9+WFQx<h}!s)q;rLH5IZ6$vVH!8U*j5x9#!cyx!}z@r1C5NrpC z0m^vzLz#i(Vi@E`aC)2o^#lG;$7%s2*kPqIIMJW_#1Ad!F)g@(C0Zfr1Re|E77932 zA=Nl2ouZ^Gh?OAqpn$=MHBfp+QIBRmC=L*@1}S(@)T5aXb`LmJK-`17>;p9=_<LHw z#Sn;%+N=QS<L~d3W?*=k3d){%8kj8nt;$0124*V@xX}pG-y6W(dVt^K(7_)}kmlxV z#tE<#^a0d$*#T}~gGv@iHn#+~Ir#m+6tvah!bozf1Gy4$-GJ2U04E=CtK$Nht&R<# zRtKm~03{Z1eFHY`fCq9-2CJH4KxHVXCIdMF(jWqBnUTh?0g6{pf&k}~4v%iv4VWf& zfP}zxBc_QfU|K+hwlAnb1Xc|;e+OC>=-Pm#odim#kWh!#IPms9DC!`&6I9@WI@LG} zqYED0*h)7@8iti_u!05NN^pIE&p@mh7h)nT<AMSXwK??yXHJDU5#$R{w!u-NgPKa< z!VAOz#aQbPaDy5YNFZ&fjTvw~(>)bzCMnIx)(>DqG3q9W=Rifo15k{ElEze!{m^m_ z)N%q#fburjVIT?6C;~_kmY{?b2@sFN$``oDA(F7T0r?lj<E{tr8HhD*ASS}%2CY!_ z=yW}SC2k;kAWj5_;|z~p98Scj99v(2gApTc;AS#liyMe6%$1;=nRn?Ee+(!9z!Clc zY$!(DK(vA42JA6VB@MP88U>)Z0ZV|&oL-Q_KoTcBI>CydBQyLREX)iH9^I}lpwSyM z6=WMoC8%HY044)A942$Xqto?;N4G23r6)e|NA!YH4@eKF>vjMp1L|Q-g~@<SI{_|J zK^7fI<JWVYQ6&KvIFZJ$=R1SHg`1gy0bC+r1R*#jLfY!!6o*n^hJrncy*y&X8ONZY z1DgdesKCZxEvT^d@DWaeC4JDSKqR=l1a}^xfetDUB3;3SASf<C!3K&~Y&{ZqoWe~6 zyRE~c8ysukumkl(H=vcHOx@riM$7^M5@O)k1^a6SD4IZZ0>~8*hc$pqMROP|;loQ8 zxWgbR1?;K@SWH7r#?iY0nFT7YK?xX=gRQ}%mHd57zyALR$2OvaqlwXF0BN%Zck}uC z`hLLr3ed5Ho=c3d>K)QS0_g{hj-p!#>QV6bNr3f(Y7uL2$ArH%`_2FVFF%TbD{_d} zFp@ev6M({~7nTV?NgbXEFp@ev6M#*FWdcyfL1Y4Gn(yNK3mruU#lwX(e!b8c{H;I# zK!XN*MqmV|U2x+XLjfYHAnhzzISY?VaJmJV08V6JLon394FRV+uo_T2VkR(<Ij%3T zD}Wo~3K`l01w5nxM2{X&0^sle^!NY&m+!v+|33kePN>r%R8hy;ArwXu#^?}&T#DQe zT=MPzf6xFEX#5)#O`sIt3gv*aI%w=0BmpXe;A7{Y1Ov)vpxh2>Xu))Vt2I#NQ3q}= zfDMBh02+dK?Eq&AkQ*S08WdWfXePc<1knqMF_2ycaJv^|0+#rL)N~*vAW3*l2T3Hj zE8_=PidBdnP*1jdD%e0$njo&&8cL89gu{stNthEs;Rz~#;lnGS_5y~VA$o8*5h9DT zCI+<^Fmejq(*lSx7)K~#G&jK+5KGSqZYI|D0z?+(N^lZ7gfyB7ZOdUe6yiBZSb|ap zN_0acVPOesdVmuthyhBy;8FmKLm}EA!3j1KCB`ALFo%LGS45xIqq7xMmO^qm$Sr8) z6Ub1H&ej|F8b44q7qF{=gwr99Zm=Ps+zD#GhJvbmc*ctWM**Z#1P$YZY8w>CK>PqQ z1Y`@y5KxQ97vz~bP~HXy(^P0p59%_3+Fc+kK>g{+UQij~(G9Z%ls-Z21LT?#WF)wj z0QnT^5LeK^A;=08L!g%6H3aHNa4~WL)_Mcg*1ezt#-kgYIG|R7+9b#>@#u8@0jtuH zwSb}tw-z!+BEV$>^4KXnrob~QAWuj5f=4bO^#+by45|acsUIAz;E@H;_)&K&hWYSF zgp4<VdMmDA^`IKkquUi8U%sF|IjHXf^Dm^xMh+#YOav@ZfrA<)zQBecM;0DKaK@J} zsJ&Qs!lN73ga$X(;8hE##R2j=D8(US4OEkW8=Bxc2^Md->!cf4vNT#;fGk2D<pJk% zJaGp~%C6wJ1BERX^U>l0oSVRL0f}o2^U>6U-2+w+iEBtr2a0h}O$W-pkP-ztIn}cn z+~mVGisZ)0!0>Y74{(bEG+PPI`sh==3E=7zWKk#55E68#8YBrBLTcU%Vl(r%oO}8I zKZL`|-*W6bsL2YFu?9~r^S8J#ftsu!kxsA(XucQ3gXDea>^x*r7+eH^XXil^;-Hj= zJR$4~8e0L6zZ`gBFNd{}qlY9sqZ@1;*o)9f3s(&f#zV;CUQAfjSa>iV0(aCrI+>2Q z7C>irUuJ?^-XJc9iQTOr|A2&H6(PtXSmgpLT6}k4QDfo32vY-2$e_U}wD}SY>nuFL zQ$;ZQV08${K2YTV3Nol!pydtSQ^5*Bm6Qo|f-Tg-gYgq|b`6wKAcJAZEjPyIwIE%L z2p1#Ua;O^|B-lI)ZWM#Gf+XR3K*5Pc57hIZ37y3<@ChBJX0Yju5DnJgaOLl5{rCSr zd^rImxx;43!9$#Av*e&mf$Rj7S#oeO1~CIg4J6xv!qB6$6;w(=t5MLDqenMHq|@~R zB0)eJr;OmqT6p_a6FlLB;t6D1TooYEhf%;o3L!{a1=f(kF{y~r?t=yvnmyo=F^_KW z^j;pQSp%-PK@3p)2DDZHBW)Sra2z~nkOmmxj)O?zn})<<56p3(5<L&YL#`Jn2ted< zF<6?1PAbEaJ-8WzJmQ2c1PKNPR*N8^0ka5JoIz(bJUSVX6Ah$80hw~-@7WEWqQsGO zK{fveaQ4Db0I?CI08&4I6yVMrh~x`TvXB%EaTp}|LNZ@xYYJr945*!kmJ6W8qw5bC z<Hm7tRs}7Cfv`XuCtgkg7nv9(f&zF!ia<BS0;J(aP$>p05y5>2(4ZH{Yv4IBP|XR7 zI#~L}FxP-&bHV8k#oP<c5RYIf{BS!PT9$zdUleD9s%>!N4({w9Sj@%kY+UAonj5I* z-hjjzBREJQ0Sl>q>v$kFoF-TX#d8n^{5?_N-UKL@c25N@2ZK&(gTn!mMv;paXkiME z6>!@U<RDPHkO_-=qy!2IAW&!{&y~VTJ6L^zZa%oJ3pO8;zA@CJnGbdkDB!{FK`BUJ z)d@2=tAT0`)YSl>MiYPEjGv$e42X?o8m_PQ=l}nZWo?jYIBQn{{+0#IpcW2j$gTMR zBV@%8qLBghbPURRA5gO!vEB#N)WS012AWrbPIArg=ydIX<WXcdV~H3@4F)RTK_vmE z0>}U<c!~>@yTJ8B-37?F3wR|CywHRMEUYX8w-rEbbWp;BPp9Er&;lMNgDH>n{eY-( zVNIk~i21Pm4sNwz8Ro@|Ch*u<-*3=jRp{CiuwIa5pg0HjCBdykaKmc@$`oM-+BmW+ zXuc0LZwTvTf?7?m81n@cJK)l0LpOK?w|gomrXcwPIU20N&1C+*bZ|}rxwv~OC^cYH z30i`|-`~RknuP!_HiX(1NnrgJatjSSwgyT*AidyN1bLy3teMmaFQ<aD11PmX2AOMY zAqfqVCO`wsAisc?ji3%RgGJ$MTSTw@|NoK=ZWMKfG9!f$t$a{Kf}{1o3lAh=f<u|x zKy&7x>F#b<2anzW7U-&4CeUJ5{uU8728M1|3C089oC~T_x?Kez1w}18dPBJR03*`M zSO*VKn1Dh7rJryJyyzBGoO;agK-U3oyhG&>Yiu3BYizrxf+fI>^{J56aG@R^jGwxN zI;XAyt;U8d@&z-I8|Sd@8^k>zFTgD3L6SPr%>k(+Kt+NHOv1y1@k2KQ|N4WSu$DL6 z&mhM@0)oE>R2sr-E_6RZ)(3<8vE8l?-3}a(?jdCH6F68wX$z93z+6xfj%bK<w}O@} zTZ5N%^SAzEgieTo^njd-Rxc1L-at_TG9RQBUeZ8X*PxOH+=hb5qZtjZ(!dHZ?1mIo zP`k0M7lbG{gy{am?S{zUwi_Z3vKvyrf%6!6y%;3K7{NUXh{4w21>yWX_dxw2(0WTy z9D#faU#$&ENswk7O5F!b5@;%6rsFUO(m1V|h8$a<z=4SK_k@58Ji<<{IRPpEVIFV< zk1&9e3}m_b%RVLs2Iz8eNV<U(KA<5EP=bc&2KyA+lZQ4FpizQUy`V=bB;-I^;ZZVF zlLI82L8%uK9<baAN)E741Mda^?+<`gP@n>;dn&RxvI*eIw5Jgqk03T`6$eoX5eJ!r zTA_htVKE9>^^7AK(Jejh2iHF!m%@Axs<2^+9N}K{^yq|^9zmTWXr%yZ;PvJ}Y6%>f z2s%oJw1Wb+Hv?2|VrT;Gs5lPZQ{D|;P~Hh1@^~GID4yY&)B{|_z(N&XoO^gMA{Xc2 z0vT(Oz80hv+Pk%Om8hT7>?*;Cb^itfsKN$2=Qy~h0IGi=xgNxNT|EJkk09L$Nb-ZM z%0=l?fU+5~<H0*eKm`UUK*5C?Y7YezBehGgcpOrffjXd_OxWEEvAx?>0$im)3JDO2 z(+McWA#}8~FAcnd0~AVNZ)U=L1~Lg=;h|U5kQP2JY4Bjw_gkF;5VNeov0J(f!s`St zTX@Yl0mU^aZiW=$5IdUpf&!1B;vz&h#4P^SFmQ4Oxf)X^xDf%d9;5)`(b9g%_8k88 z;7wGJ4LUCanL&jv<Ai49WP)iTL+gQxJCO1PGA00yOi<idg7^IJ`++G~tq;m^uw(*q z0Tu&0p?ynGM(sSrzaBhF3F%LP>Ti(IppqZtcSt3Os0~0PXUD-+B517PWh1z;1xiL> zYoL{bs{%ORfYgG7i7ohD4X|s2RMGG%98?Tzf`|Q@Ar66bdZCNFTXVsY0}2IHSK~@A z0uTwt3(!4UkTrmik`|iC!NCRc2|RB>63O=)u#kj#mGs&VLnp+yL^Nd(;Q(164=LqQ z10x9>7@+upBxs0VVd?lPxZ1%CDvSuma3rJ($Kps()Wc&6LmN)NfubI+4dF3t5!?@o z;O4!cM9RS5`soj7?gNz08Tnh+gIDl^nl<1;2_)JL?iy-Bx{u(@hrEXmtPC^{)4Uf{ z$ujVF=0n!OoIDKPn%P<V!lSeFfk$WP4ba?XH$$iEm*(0p4E4&#S(k#&%K>*IT^}3= zm0qA_Z8spSmn%Rc60To5S)17y7>={nLe+M=e(7}m(Omn3vDD#Z3rG!X8XE(HNAd-a zULHuFqgNEvB|6R`2s%0k<UUpxaOL--H-I0ygqWq%^#`a7aQ)HzfWP%X$*$wBe?SNC zy!P#86^9vT0yFL?D+5Ef@1IT{F|cXhKwQ@!klCc#4{)R40vB2jlx%Ny{lnDl`iFnL z>mSg*HrEf<u0QH-G}r!Msdw*Y-Neek(8;<1WD@H#Rt5%8)__)yy{w@i;m{wQB8(t^ z^0#~h6%nlUtPBib`^#7v7+&55iGiHY+<KtYvYRy#WC%(i@Ut>79A^;(4P%0}*@H}V z{eirW+Vw{_s|MI5Y9L8b1rW(9#md0Ic%e6dsoV8O^9x4EhDcT(ka4bmK*n{$J^!ck zLhFH2({9#BEYLWB2igP{lt4QO;v(sD{n73E2jni;&J#%g)%8!e>mN`^&F*HM#KOSP z8^GAj+RDPfV0oxUyj!#yWD08~3&?>W&8%f00nro?$(jTa{=r|Te4N!26#2N)mJeuD zm(_@cfuZ$4Nnkgt1`7j2r|X~B4&AJhEDQ|IxqleDL;rxslHh6T04S~T2!azD+>@)B z85m&6;yyD2Lv!sPhSFW#qL)F6SucWQ!7bV57mT17h8BLkti>RC7N$<tO&|s<D6O|1 zsMG6aod+@mO|ctDjR;dWYb!G-96;MsS*t)|tOX#|KNvw8Q$frZOenVSf)s%?`XDqq zg2Y(OL8^Z+fi&uXm@k-7G+qJCcY-tuBQ$b=G>S4o)}B7-6zCTH#RM7wa=p<h5DU)h zj-ier!5)ooR)B)Nvv!3?XXyfu&d?bioxT%%I$b+_I(<Rw4;O&SfIFbwKOhFERG0u_ zT=D33?Eo=AYamv5G#}squl8kaVFJZcDHnL0{^cJ~vS%$|Vqmazy-_OETziM1j;Fcy z3PT-RbL|BNm@?2zVYlmzPS+#du1EOSyI$#ZJ<#oXfPcN~ole(1-K^Firzc<V=)B;; zcmOo!e((XaM<=KY+u+ly!_>_x1yV2rBy{i*i$}L>gHNvtQ#UI!M7S4LNx3$7^vW=~ zbhADKr6}7wj0_Ae-LV&3x<e0i`yS~I-Q&?+dcvdIb%jqaZ#^RegHNZAiiJ<NuK<W) z;M19-qT$n-qN3o@8KWZM)9u>;5*KjkbW!2xv{7;Gyy@C`!K3r|3$cU$|L*`b7d#q| zfPy#LF~%`A_V9mIDMki{|Ee+|qkCmp85tNHga5031?{cO02Qb&pZ)*;KLd1})yq53 zZPE-34Dhq(VCQ-<Rs}IIRtPW(@UU}CU}Ru0VPIg8VPIfbeBeJA%(es>20FJ6<OV*0 zHYU&-WCmU~mT8s@3_R?hv$jAgK>9#OezE-mYy1BORHlL?SkU!7XGYcsI<E{gp1A8T zL|++5IaYn36R~=kdsv!TSyCC1jR2{}z&!pCAu%Wo;f<1`Aut*OqaiRF0;3@?8Umvs zFd71*Aut*OqaiRF0;3@?8UjN(1OycYf<W#UbV#YpNr4=KpOl%Gl3K(NSjaGep-@23 zAwD_37{rtibcip`%ww1U(!dy>mR1taz{$XrmY7_UU&J6;o?n#BAkV_kH<^)vh2hsM z@QB1``Ps#9*cljp-2;t=GPppG5C@q68ZH5=`E?I8Mgcxv9K;6=TZ6_-xpOj;U`KU> z)PqMxKoWe63=FXQ0YK-A^KmdRz|L_8T_eE9$-n?Rhh2<;fdRx%fS$`P3+1Ok`HE0} z29&Q3<>x^8x}ZA|7#SG)JV8!iU;qv4@G~+nRQW>q&fwF<!4tb+0npjrj7Z1EgOnlj zLGI;ZU|;|ZFflMNIKcP{450ALtzd}HOG`{<h|kYSiO*w*&y5H3vP&4^a}z7lz+?ef YDM$ho>>zFN;GhRb6F7n}qLP6D02YudsQ>@~ literal 0 HcmV?d00001 diff --git a/com.oracle.truffle.r.native/src/fft.c b/com.oracle.truffle.r.native/src/fft.c new file mode 100644 index 0000000000..ce69cd5cae --- /dev/null +++ b/com.oracle.truffle.r.native/src/fft.c @@ -0,0 +1,880 @@ +/* + * R : A Computer Language for Statistical Data Analysis + * Copyright (C) 1995, 1996, 1997 Robert Gentleman and Ross Ihaka + * Copyright (C) 1998--2000 R Core Team + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, a copy is available at + * http://www.r-project.org/Licenses/ + */ + +// modifications required for standalone compilation + +//#ifdef HAVE_CONFIG_H +//#include <config.h> +//#endif + +#include <stdlib.h> /* for abs */ +#include <math.h> + +// modifications required for standalone compilation + +//#include <Rmath.h> /* for imax2(.),..*/ +//#include <R_ext/Applic.h> +#define imax2(_x,_y) ((_x<_y) ? _y : _x) +#define imin2(_x,_y) ((_x<_y) ? _x : _y) +#define Rboolean int +#define FALSE 0 +#define TRUE 1 +#define M_SQRT_3 1.732050807568877293527446341506 + +/* Fast Fourier Transform + * + * These routines are based on code by Richard Singleton in the + * book "Programs for Digital Signal Processing" put out by IEEE. + * + * I have translated them to C and moved the memory allocation + * so that it takes place under the control of the algorithm + * which calls these; for R, see ../main/fourier.c + * + * void fft_factor(int n, int *maxf, int *maxp) + * + * This factorizes the series length and computes the values of + * maxf and maxp which determine the amount of scratch storage + * required by the algorithm. + * + * If maxf is zero on return, an error occured during factorization. + * The nature of the error can be determined from the value of maxp. + * If maxp is zero, an invalid (zero) parameter was passed and + * if maxp is one, the internal nfac array was too small. This can only + * happen for series lengths which exceed 12,754,584. + * + * PROBLEM (see fftmx below): nfac[] is overwritten by fftmx() in fft_work() + * ------- Consequence: fft_factor() must be called way too often, + * at least from do_mvfft() [ ../main/fourier.c ] + * + * The following arrays need to be allocated following the call to + * fft_factor and preceding the call to fft_work. + * + * work double[4*maxf] + * iwork int[maxp] + * + * int fft_work(double *a, double *b, int nseg, int n, int nspn, + * int isn, double *work, int *iwork) + * + * The routine returns 1 if the transform was completed successfully and + * 0 if invalid values of the parameters were supplied. + * + * Ross Ihaka + * University of Auckland + * February 1997 + * ========================================================================== + * + * Header from the original Singleton algorithm: + * + * -------------------------------------------------------------- + * subroutine: fft + * multivariate complex fourier transform, computed in place + * using mixed-radix fast fourier transform algorithm. + * -------------------------------------------------------------- + * + * arrays a and b originally hold the real and imaginary + * components of the data, and return the real and + * imaginary components of the resulting fourier coefficients. + * multivariate data is indexed according to the fortran + * array element successor function, without limit + * on the number of implied multiple subscripts. + * the subroutine is called once for each variate. + * the calls for a multivariate transform may be in any order. + * + * n is the dimension of the current variable. + * nspn is the spacing of consecutive data values + * while indexing the current variable. + * nseg nseg*n*nspn is the total number of complex data values. + * isn the sign of isn determines the sign of the complex + * exponential, and the magnitude of isn is normally one. + * the magnitude of isn determines the indexing increment for a&b. + * + * if fft is called twice, with opposite signs on isn, an + * identity transformation is done...calls can be in either order. + * the results are scaled by 1/n when the sign of isn is positive. + * + * a tri-variate transform with a(n1,n2,n3), b(n1,n2,n3) + * is computed by + * call fft(a,b,n2*n3,n1, 1, -1) + * call fft(a,b,n3 ,n2,n1, -1) + * call fft(a,b,1, n3,n1*n2,-1) + * + * a single-variate transform of n complex data values is computed by + * call fft(a,b,1,n,1,-1) + * + * the data may alternatively be stored in a single complex + * array a, then the magnitude of isn changed to two to + * give the correct indexing increment and a(2) used to + * pass the initial address for the sequence of imaginary + * values, e.g. + * call fft(a,a(2),nseg,n,nspn,-2) + * + * nfac[15] (array) is working storage for factoring n. the smallest + * number exceeding the 15 locations provided is 12,754,584. + * + */ + +static void fftmx(double *a, double *b, int ntot, int n, int nspan, int isn, + int m, int kt, double *at, double *ck, double *bt, double *sk, + int *np, int *nfac) +{ +/* called from fft_work() */ + +/* Design BUG: One purpose of fft_factor() would be to compute + * ---------- nfac[] once and for all; and fft_work() [i.e. fftmx ] + * could reuse the factorization. + * However: nfac[] is `destroyed' currently in the code below + */ + double aa, aj, ajm, ajp, ak, akm, akp; + double bb, bj, bjm, bjp, bk, bkm, bkp; + double c1, c2=0, c3=0, c72, cd; + double dr, rad; + double s1, s120, s2=0, s3=0, s72, sd; + int i, inc, j, jc, jf, jj; + int k, k1, k2, k3=0, k4, kk, klim, ks, kspan, kspnn; + int lim, maxf, mm, nn, nt; + + a--; b--; at--; ck--; bt--; sk--; + np--; + nfac--;/*the global one!*/ + + inc = abs(isn); + nt = inc*ntot; + ks = inc*nspan; + rad = M_PI_4;/* = pi/4 =^= 45 degrees */ + s72 = rad/0.625;/* 72 = 45 / .625 degrees */ + c72 = cos(s72); + s72 = sin(s72); + s120 = 0.5*M_SQRT_3;/* sin(120) = sqrt(3)/2 */ + if(isn <= 0) { + s72 = -s72; + s120 = -s120; + rad = -rad; + } else { +#ifdef SCALING + /* scale by 1/n for isn > 0 */ + ak = 1.0/n; + for(j=1 ; j<=nt ; j+=inc) { + a[j] *= ak; + b[j] *= ak; + } +#endif + } + + kspan = ks; + nn = nt - inc; + jc = ks/n; + + /* sin, cos values are re-initialized each lim steps */ + + lim = 32; + klim = lim*jc; + i = 0; + jf = 0; + maxf = nfac[m - kt]; + if(kt > 0) maxf = imax2(nfac[kt],maxf); + + /* compute fourier transform */ + +L_start: + dr = (8.0*jc)/kspan; + cd = sin(0.5*dr*rad); + cd = 2.0*cd*cd; + sd = sin(dr*rad); + kk = 1; + i++; + if( nfac[i] != 2) goto L110; + +/* transform for factor of 2 (including rotation factor) */ + + kspan /= 2; + k1 = kspan + 2; + do { + do { + k2 = kk + kspan; + ak = a[k2]; + bk = b[k2]; + a[k2] = a[kk] - ak; + b[k2] = b[kk] - bk; + a[kk] += ak; + b[kk] += bk; + kk = k2 + kspan; + } while(kk <= nn); + kk -= nn; + } while(kk <= jc); + + if(kk > kspan) goto L_fin; +L60: + c1 = 1.0 - cd; + s1 = sd; + mm = imin2(k1/2,klim); + goto L80; + +L70: + ak = c1 - (cd*c1+sd*s1); + s1 = (sd*c1-cd*s1) + s1; + +/* the following three statements compensate for truncation error. */ +/* if rounded arithmetic is used (nowadays always ?!), substitute c1=ak */ +#ifdef TRUNCATED_ARITHMETIC + c1 = 0.5/(ak*ak+s1*s1) + 0.5; + s1 = c1*s1; + c1 = c1*ak; +#else + c1 = ak; +#endif + +L80: + do { + k2 = kk + kspan; + ak = a[kk] - a[k2]; + bk = b[kk] - b[k2]; + a[kk] += a[k2]; + b[kk] += b[k2]; + a[k2] = c1*ak - s1*bk; + b[k2] = s1*ak + c1*bk; + kk = k2 + kspan; + } while(kk < nt); + k2 = kk - nt; + c1 = -c1; + kk = k1 - k2; + if( kk > k2) goto L80; + kk += jc; + if(kk <= mm) goto L70; + if(kk >= k2) { + k1 = k1 + inc + inc; + kk = (k1-kspan)/2 + jc; + if( kk <= jc+jc) goto L60; + goto L_start; + } + + s1 = ((kk-1)/jc)*dr*rad; + c1 = cos(s1); + s1 = sin(s1); + mm = imin2(k1/2,mm+klim); + goto L80; + +/* transform for factor of 3 (optional code) */ + +L100: + k1 = kk + kspan; + k2 = k1 + kspan; + ak = a[kk]; + bk = b[kk]; + aj = a[k1] + a[k2]; + bj = b[k1] + b[k2]; + a[kk] = ak + aj; + b[kk] = bk + bj; + ak = -0.5*aj + ak; + bk = -0.5*bj + bk; + aj = (a[k1]-a[k2])*s120; + bj = (b[k1]-b[k2])*s120; + a[k1] = ak - bj; + b[k1] = bk + aj; + a[k2] = ak + bj; + b[k2] = bk - aj; + kk = k2 + kspan; + if( kk < nn) goto L100; + kk = kk - nn; + if( kk <= kspan) goto L100; + goto L290; + +/* transform for factor of 4 */ + +L110: + if( nfac[i] != 4) goto L_f_odd; + kspnn = kspan; + kspan /= 4; +L120: + c1 = 1.0; + s1 = 0; + mm = imin2(kspan,klim); + goto L150; +L130: + c2 = c1 - (cd*c1+sd*s1); + s1 = (sd*c1-cd*s1) + s1; + +/* the following three statements compensate for truncation error. */ +/* if rounded arithmetic is used (nowadays always ?!), substitute c1=c2 */ +#ifdef TRUNCATED_ARITHMETIC + c1 = 0.5/(c2*c2+s1*s1) + 0.5; + s1 = c1*s1; + c1 = c1*c2; +#else + c1 = c2; +#endif + +L140: + c2 = c1*c1 - s1*s1; + s2 = c1*s1*2.0; + c3 = c2*c1 - s2*s1; + s3 = c2*s1 + s2*c1; + +L150: + k1 = kk + kspan; + k2 = k1 + kspan; + k3 = k2 + kspan; + akp = a[kk] + a[k2]; + akm = a[kk] - a[k2]; + ajp = a[k1] + a[k3]; + ajm = a[k1] - a[k3]; + a[kk] = akp + ajp; + ajp = akp - ajp; + bkp = b[kk] + b[k2]; + bkm = b[kk] - b[k2]; + bjp = b[k1] + b[k3]; + bjm = b[k1] - b[k3]; + b[kk] = bkp + bjp; + bjp = bkp - bjp; + if( isn < 0) goto L180; + akp = akm - bjm; + akm = akm + bjm; + bkp = bkm + ajm; + bkm = bkm - ajm; + if( s1 == 0.0) goto L190; +L160: + a[k1] = akp*c1 - bkp*s1; + b[k1] = akp*s1 + bkp*c1; + a[k2] = ajp*c2 - bjp*s2; + b[k2] = ajp*s2 + bjp*c2; + a[k3] = akm*c3 - bkm*s3; + b[k3] = akm*s3 + bkm*c3; + kk = k3 + kspan; + if( kk <= nt) goto L150; +L170: + kk = kk - nt + jc; + if( kk <= mm) goto L130; + if( kk < kspan) goto L200; + kk = kk - kspan + inc; + if(kk <= jc) goto L120; + if(kspan == jc) goto L_fin; + goto L_start; +L180: + akp = akm + bjm; + akm = akm - bjm; + bkp = bkm - ajm; + bkm = bkm + ajm; + if( s1 != 0.0) goto L160; +L190: + a[k1] = akp; + b[k1] = bkp; + a[k2] = ajp; + b[k2] = bjp; + a[k3] = akm; + b[k3] = bkm; + kk = k3 + kspan; + if( kk <= nt) goto L150; + goto L170; +L200: + s1 = ((kk-1)/jc)*dr*rad; + c1 = cos(s1); + s1 = sin(s1); + mm = imin2(kspan,mm+klim); + goto L140; + +/* transform for factor of 5 (optional code) */ + +L_f5: + c2 = c72*c72 - s72*s72; + s2 = 2.0*c72*s72; +L220: + k1 = kk + kspan; + k2 = k1 + kspan; + k3 = k2 + kspan; + k4 = k3 + kspan; + akp = a[k1] + a[k4]; + akm = a[k1] - a[k4]; + bkp = b[k1] + b[k4]; + bkm = b[k1] - b[k4]; + ajp = a[k2] + a[k3]; + ajm = a[k2] - a[k3]; + bjp = b[k2] + b[k3]; + bjm = b[k2] - b[k3]; + aa = a[kk]; + bb = b[kk]; + a[kk] = aa + akp + ajp; + b[kk] = bb + bkp + bjp; + ak = akp*c72 + ajp*c2 + aa; + bk = bkp*c72 + bjp*c2 + bb; + aj = akm*s72 + ajm*s2; + bj = bkm*s72 + bjm*s2; + a[k1] = ak - bj; + a[k4] = ak + bj; + b[k1] = bk + aj; + b[k4] = bk - aj; + ak = akp*c2 + ajp*c72 + aa; + bk = bkp*c2 + bjp*c72 + bb; + aj = akm*s2 - ajm*s72; + bj = bkm*s2 - bjm*s72; + a[k2] = ak - bj; + a[k3] = ak + bj; + b[k2] = bk + aj; + b[k3] = bk - aj; + kk = k4 + kspan; + if( kk < nn) goto L220; + kk = kk - nn; + if( kk <= kspan) goto L220; + goto L290; + +/* transform for odd factors */ + +L_f_odd: + k = nfac[i]; + kspnn = kspan; + kspan /= k; + if(k == 3) goto L100; + if(k == 5) goto L_f5; + if(k == jf) goto L250; + jf = k; + s1 = rad/(k/8.0); + c1 = cos(s1); + s1 = sin(s1); + ck[jf] = 1.0; + sk[jf] = 0.0; + + for(j = 1; j < k; j++) { /* k is changing as well */ + ck[j] = ck[k]*c1 + sk[k]*s1; + sk[j] = ck[k]*s1 - sk[k]*c1; + k--; + ck[k] = ck[j]; + sk[k] = -sk[j]; + } + +L250: + k1 = kk; + k2 = kk + kspnn; + aa = a[kk]; + bb = b[kk]; + ak = aa; + bk = bb; + j = 1; + k1 = k1 + kspan; +L260: + k2 = k2 - kspan; + j++; + at[j] = a[k1] + a[k2]; + ak = at[j] + ak; + bt[j] = b[k1] + b[k2]; + bk = bt[j] + bk; + j++; + at[j] = a[k1] - a[k2]; + bt[j] = b[k1] - b[k2]; + k1 = k1 + kspan; + if( k1 < k2) goto L260; + a[kk] = ak; + b[kk] = bk; + k1 = kk; + k2 = kk + kspnn; + j = 1; +L270: + k1 += kspan; + k2 -= kspan; + jj = j; + ak = aa; + bk = bb; + aj = 0.0; + bj = 0.0; + k = 1; + for(k=2; k < jf; k++) { + ak += at[k]*ck[jj]; + bk += bt[k]*ck[jj]; + k++; + aj += at[k]*sk[jj]; + bj += bt[k]*sk[jj]; + jj += j; + if(jj > jf) jj -= jf; + } + k = jf - j; + a[k1] = ak - bj; + b[k1] = bk + aj; + a[k2] = ak + bj; + b[k2] = bk - aj; + j++; + if( j < k) goto L270; + kk = kk + kspnn; + if( kk <= nn) goto L250; + kk = kk - nn; + if( kk <= kspan) goto L250; + +/* multiply by rotation factor (except for factors of 2 and 4) */ + +L290: + if(i == m) goto L_fin; + kk = jc + 1; +L300: + c2 = 1.0 - cd; + s1 = sd; + mm = imin2(kspan,klim); + + do { /* L320: */ + c1 = c2; + s2 = s1; + kk += kspan; + do { /* L330: */ + do { + ak = a[kk]; + a[kk] = c2*ak - s2*b[kk]; + b[kk] = s2*ak + c2*b[kk]; + kk += kspnn; + } while(kk <= nt); + ak = s1*s2; + s2 = s1*c2 + c1*s2; + c2 = c1*c2 - ak; + kk += -nt + kspan; + } while(kk <= kspnn); + kk += -kspnn + jc; + if(kk <= mm) { /* L310: */ + c2 = c1 - (cd*c1+sd*s1); + s1 = s1 + (sd*c1-cd*s1); +/* the following three statements compensate for truncation error.*/ +/* if rounded arithmetic is used (nowadays always ?!), they may be deleted. */ +#ifdef TRUNCATED_ARITHMETIC + c1 = 0.5/(c2*c2+s1*s1) + 0.5; + s1 = c1*s1; + c2 = c1*c2; +#endif + continue/* goto L320*/; + } + if(kk >= kspan) { + kk = kk - kspan + jc + inc; + if( kk <= jc+jc) goto L300; + goto L_start; + } + s1 = ((kk-1)/jc)*dr*rad; + c2 = cos(s1); + s1 = sin(s1); + mm = imin2(kspan,mm+klim); + } while(1); + +/*------------------------------------------------------------*/ + + +/* permute the results to normal order---done in two stages */ +/* permutation for square factors of n */ + +L_fin: + np[1] = ks; + if( kt == 0) goto L440; + k = kt + kt + 1; + if( m < k) k--; + np[k+1] = jc; + for(j = 1; j < k; j++, k--) { + np[j+1] = np[j]/nfac[j]; + np[k] = np[k+1]*nfac[j]; + } + k3 = np[k+1]; + kspan = np[2]; + kk = jc + 1; + k2 = kspan + 1; + j = 1; + + if(n == ntot) { + + /* permutation for single-variate transform (optional code) */ + + L370: + do { + ak = a[kk]; a[kk] = a[k2]; a[k2] = ak; + bk = b[kk]; b[kk] = b[k2]; b[k2] = bk; + kk += inc; + k2 += kspan; + } while(k2 < ks); + L380: + do { k2 -= np[j]; j++; k2 += np[j+1]; } while(k2 > np[j]); + j = 1; + do { + if(kk < k2) goto L370; + kk += inc; + k2 += kspan; + } while(k2 < ks); + if( kk < ks) goto L380; + jc = k3; + + } else { + + /* permutation for multivariate transform */ + + L400: + k = kk + jc; + do { + ak = a[kk]; a[kk] = a[k2]; a[k2] = ak; + bk = b[kk]; b[kk] = b[k2]; b[k2] = bk; + kk += inc; + k2 += inc; + } while( kk < k); + kk += ks - jc; + k2 += ks - jc; + if(kk < nt) goto L400; + k2 += - nt + kspan; + kk += - nt + jc; + if( k2 < ks) goto L400; + + do { + do { k2 -= np[j]; j++; k2 += np[j+1]; } while(k2 > np[j]); + j = 1; + do { + if(kk < k2) goto L400; + kk += jc; + k2 += kspan; + } while(k2 < ks); + } while(kk < ks); + jc = k3; + } + +L440: + if( 2*kt+1 >= m) return; + kspnn = np[kt+1]; + +/* permutation for square-free factors of n */ + + /* Here, nfac[] is overwritten... -- now CUMULATIVE ("cumprod") factors */ + nn = m - kt; + nfac[nn+1] = 1; + for(j = nn; j > kt; j--) + nfac[j] *= nfac[j+1]; + kt++; + nn = nfac[kt] - 1; + jj = 0; + j = 0; + goto L480; +L460: + jj -= k2; + k2 = kk; + k++; + kk = nfac[k]; +L470: + jj += kk; + if( jj >= k2) goto L460; + np[j] = jj; +L480: + k2 = nfac[kt]; + k = kt + 1; + kk = nfac[k]; + j++; + if( j <= nn) goto L470; + +/* determine the permutation cycles of length greater than 1 */ + + j = 0; + goto L500; + + do { + do { k = kk; kk = np[k]; np[k] = -kk; } while(kk != j); + k3 = kk; + L500: + do { j++; kk = np[j]; } while(kk < 0); + } while(kk != j); + np[j] = -j; + if( j != nn) goto L500; + maxf *= inc; + goto L570; + +/* reorder a and b, following the permutation cycles */ + +L_ord: + do j--; while(np[j] < 0); + jj = jc; + +L520: + kspan = imin2(jj,maxf); + jj -= kspan; + k = np[j]; + kk = jc*k + i + jj; + + for(k1= kk + kspan, k2= 1; k1 != kk; + k1 -= inc, k2++) { + at[k2] = a[k1]; + bt[k2] = b[k1]; + } + + do { + k1 = kk + kspan; + k2 = k1 - jc*(k+np[k]); + k = -np[k]; + do { + a[k1] = a[k2]; + b[k1] = b[k2]; + k1 -= inc; + k2 -= inc; + } while( k1 != kk); + kk = k2; + } while(k != j); + + for(k1= kk + kspan, k2= 1; k1 > kk; + k1 -= inc, k2++) { + a[k1] = at[k2]; + b[k1] = bt[k2]; + } + + if(jj != 0) goto L520; + if( j != 1) goto L_ord; + +L570: + j = k3 + 1; + nt = nt - kspnn; + i = nt - inc + 1; + if( nt >= 0) goto L_ord; +} /* fftmx */ + +static int old_n = 0; + +static int nfac[15]; +static int m_fac; +static int kt; +static int maxf; +static int maxp; + +/* At the end of factorization, + * nfac[] contains the factors, + * m_fac contains the number of factors and + * kt contains the number of square factors */ + +/* non-API, but used by package RandomFields */ +// signature modification (pointers to scalars) required to enable JNR call +void fft_factor(int *nptr, int *pmaxf, int *pmaxp) +{ +/* fft_factor - factorization check and determination of memory + * requirements for the fft. + * + * On return, *pmaxf will give the maximum factor size + * and *pmaxp will give the amount of integer scratch storage required. + * + * If *pmaxf == 0, there was an error, the error type is indicated by *pmaxp: + * + * If *pmaxp == 0 There was an illegal zero parameter among nseg, n, and nspn. + * If *pmaxp == 1 There we more than 15 factors to ntot. */ + + int n = *nptr; + int j, jj, k; + + /* check series length */ + + if (n <= 0) { + old_n = 0; *pmaxf = 0; *pmaxp = 0; + return; + } + else old_n = n; + + /* determine the factors of n */ + + m_fac = 0; + k = n;/* k := remaining unfactored factor of n */ + if (k == 1) + return; + + /* extract square factors first ------------------ */ + + /* extract 4^2 = 16 separately + * ==> at most one remaining factor 2^2 = 4, done below */ + while(k % 16 == 0) { + nfac[m_fac++] = 4; + k /= 16; + } + + /* extract 3^2, 5^2, ... */ + for(j = 3; (jj= j*j) <= k; j += 2) { + while(k % jj == 0) { + nfac[m_fac++] = j; + k /= jj; + } + } + + if(k <= 4) { + kt = m_fac; + nfac[m_fac] = k; + if(k != 1) m_fac++; + } + else { + if(k % 4 == 0) { + nfac[m_fac++] = 2; + k /= 4; + } + + /* all square factors out now, but k >= 5 still */ + + kt = m_fac; + maxp = imax2(kt+kt+2, k-1); + j = 2; + do { + if (k % j == 0) { + nfac[m_fac++] = j; + k /= j; + } + j = ((j+1)/2)*2 + 1; + } + while(j <= k); + } + + if (m_fac <= kt+1) + maxp = m_fac+kt+1; + if (m_fac+kt > 15) { /* error - too many factors */ + old_n = 0; *pmaxf = 0; *pmaxp = 0; + return; + } + else { + if (kt != 0) { + j = kt; + while(j != 0) + nfac[m_fac++] = nfac[--j]; + } + maxf = nfac[m_fac-kt-1]; +/* The last squared factor is not necessarily the largest PR#1429 */ + if (kt > 0) maxf = imax2(nfac[kt-1], maxf); + if (kt > 1) maxf = imax2(nfac[kt-2], maxf); + if (kt > 2) maxf = imax2(nfac[kt-3], maxf); + } + *pmaxf = maxf; + *pmaxp = maxp; +} + + +// signature modification (pointers to scalars) required to enable JNR call +// we also need to do pointer shift for imaginary parts on the callee side (below) +// rather than on the caller side as in GNU R +Rboolean fft_work(double *a, int *nsegptr, int *nptr, int *nspnptr, int *isnptr, + double *work, int *iwork) +{ + int nseg = *nsegptr; + int n = *nptr; + int nspn = *nspnptr; + int isn = *isnptr; + double *b=&(a[1]); + int nf, nspan, ntot; + + /* check that factorization was successful */ + + if(old_n == 0) return FALSE; + + /* check that the parameters match those of the factorization call */ + + if(n != old_n || nseg <= 0 || nspn <= 0 || isn == 0) + return FALSE; + + /* perform the transform */ + + nf = n; + nspan = nf * nspn; + ntot = nspan * nseg; + + fftmx(a, b, ntot, nf, nspan, isn, m_fac, kt, + &work[0], &work[maxf], &work[2*maxf], &work[3*maxf], + iwork, nfac); + + return TRUE; +} diff --git a/com.oracle.truffle.r.nodes/src/com/oracle/truffle/r/nodes/builtin/base/FFTFunctions.java b/com.oracle.truffle.r.nodes/src/com/oracle/truffle/r/nodes/builtin/base/FFTFunctions.java deleted file mode 100644 index 83b1fc15f1..0000000000 --- a/com.oracle.truffle.r.nodes/src/com/oracle/truffle/r/nodes/builtin/base/FFTFunctions.java +++ /dev/null @@ -1,61 +0,0 @@ -/* - * Copyright (c) 2014, 2014, Oracle and/or its affiliates. All rights reserved. - * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER. - * - * This code is free software; you can redistribute it and/or modify it - * under the terms of the GNU General Public License version 2 only, as - * published by the Free Software Foundation. - * - * This code is distributed in the hope that it will be useful, but WITHOUT - * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License - * version 2 for more details (a copy is included in the LICENSE file that - * accompanied this code). - * - * You should have received a copy of the GNU General Public License version - * 2 along with this work; if not, write to the Free Software Foundation, - * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA. - * - * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA - * or visit www.oracle.com if you need additional information or have any - * questions. - */ -package com.oracle.truffle.r.nodes.builtin.base; - -import com.oracle.truffle.api.dsl.*; -import com.oracle.truffle.r.nodes.*; -import com.oracle.truffle.r.nodes.access.*; -import com.oracle.truffle.r.nodes.builtin.*; -import com.oracle.truffle.r.runtime.*; -import com.oracle.truffle.r.runtime.data.*; -import com.oracle.truffle.r.runtime.data.model.*; -import com.oracle.truffle.r.runtime.ffi.*; - -public abstract class FFTFunctions { - - private abstract static class FFTAdapter extends RBuiltinNode { - private static final String[] PARAMETER_NAMES = new String[]{"z", "inverse"}; - - @Override - public Object[] getParameterNames() { - return PARAMETER_NAMES; - } - - @Override - public RNode[] getParameterValues() { - return new RNode[]{ConstantNode.create(RMissing.instance), ConstantNode.create(RRuntime.LOGICAL_FALSE)}; - } - - } - - @RBuiltin(name = "fft", kind = RBuiltinKind.SUBSTITUTE) - public abstract static class FFT extends FFTAdapter { - @Specialization - public Object doFFT(RAbstractVector zIn, byte inverse) { - @SuppressWarnings("unused") - RDerivedRFFI ffi = RFFIFactory.getRFFI().getRDerivedRFFI(); -// ffi.fft_work(a, b, nseg, n, nspn, isn, work, iwork) - return RNull.instance; - } - } -} diff --git a/com.oracle.truffle.r.nodes/src/com/oracle/truffle/r/nodes/builtin/base/ForeignFunctions.java b/com.oracle.truffle.r.nodes/src/com/oracle/truffle/r/nodes/builtin/base/ForeignFunctions.java index b35e12cdd8..43f4d51c2c 100644 --- a/com.oracle.truffle.r.nodes/src/com/oracle/truffle/r/nodes/builtin/base/ForeignFunctions.java +++ b/com.oracle.truffle.r.nodes/src/com/oracle/truffle/r/nodes/builtin/base/ForeignFunctions.java @@ -24,11 +24,14 @@ package com.oracle.truffle.r.nodes.builtin.base; import com.oracle.truffle.api.*; import com.oracle.truffle.api.dsl.*; +import com.oracle.truffle.api.frame.*; import com.oracle.truffle.r.nodes.*; import com.oracle.truffle.r.nodes.access.*; import com.oracle.truffle.r.nodes.builtin.*; +import com.oracle.truffle.r.nodes.unary.*; import com.oracle.truffle.r.runtime.*; import com.oracle.truffle.r.runtime.data.*; +import com.oracle.truffle.r.runtime.data.model.*; import com.oracle.truffle.r.runtime.ffi.*; import com.oracle.truffle.r.runtime.ffi.DLL.*; @@ -289,4 +292,85 @@ public class ForeignFunctions { } + /** + * For now, just some special case functions that are built in to the implementation. + */ + @RBuiltin(name = ".Call", kind = RBuiltinKind.PRIMITIVE, isCombine = true) + @NodeField(name = "argNames", type = String[].class) + public abstract static class Call extends Adapter { + + @Child private CastComplexNode castComplex; + @Child private CastLogicalNode castLogical; + @Child private CastToVectorNode castVector; + + private Object castComplex(VirtualFrame frame, Object operand) { + if (castComplex == null) { + CompilerDirectives.transferToInterpreterAndInvalidate(); + castComplex = insert(CastComplexNodeFactory.create(null, true, false, false)); + } + return castComplex.executeCast(frame, operand); + } + + private Object castLogical(VirtualFrame frame, Object operand) { + if (castLogical == null) { + CompilerDirectives.transferToInterpreterAndInvalidate(); + castLogical = insert(CastLogicalNodeFactory.create(null, true, false, false)); + } + return castLogical.executeCast(frame, operand); + } + + private RAbstractVector castVector(VirtualFrame frame, Object value) { + if (castVector == null) { + CompilerDirectives.transferToInterpreterAndInvalidate(); + castVector = insert(CastToVectorNodeFactory.create(null, false, false, false, false)); + } + return castVector.executeRAbstractVector(frame, value); + } + + // TODO: handle more argumet types (this is sufficient to run the b25-matfunc1 benchmark + @SuppressWarnings("unused") + @Specialization(order = 1, guards = "fft") + public RComplexVector callFFT(VirtualFrame frame, RList f, Object[] args) { + controlVisibility(); + RComplexVector z = (RComplexVector) castComplex(frame, castVector(frame, args[0])); + RComplexVector res = z; + if (res.isShared()) { + res = (RComplexVector) z.copy(); + } + RLogicalVector inverse = (RLogicalVector) castLogical(frame, castVector(frame, args[1])); + int inv = RRuntime.isNA(inverse.getDataAt(0)) || inverse.getDataAt(0) == RRuntime.LOGICAL_FALSE ? -2 : 2; + int retCode = 7; + if (res.getLength() > 1) { + if (z.getDimensions() == null) { + int n = res.getLength(); + int[] maxf = new int[1]; + int[] maxp = new int[1]; + RFFIFactory.getRFFI().getRDerivedRFFI().fft_factor(n, maxf, maxp); + if (maxf[0] == 0) { + throw RError.getGenericError(getEncapsulatingSourceSection(), "fft factorization error"); + } + double[] work = new double[4 * maxf[0]]; + int[] iwork = new int[maxp[0]]; + retCode = RFFIFactory.getRFFI().getRDerivedRFFI().fft_work(res.getDataWithoutCopying(), 1, n, 1, inv, work, iwork); + } + } + + return res; + } + + public boolean fft(RList f) { + if (f.getNames() == RNull.instance) { + return false; + } + RStringVector names = (RStringVector) f.getNames(); + for (int i = 0; i < names.getLength(); i++) { + if (names.getDataAt(i).equals("name")) { + return f.getDataAt(i).equals("fft") ? true : false; + } + } + return false; + } + + } + } diff --git a/com.oracle.truffle.r.nodes/src/com/oracle/truffle/r/nodes/builtin/stats/R/fft.R b/com.oracle.truffle.r.nodes/src/com/oracle/truffle/r/nodes/builtin/stats/R/fft.R new file mode 100644 index 0000000000..dee7a73499 --- /dev/null +++ b/com.oracle.truffle.r.nodes/src/com/oracle/truffle/r/nodes/builtin/stats/R/fft.R @@ -0,0 +1,47 @@ +# File src/library/stats/R/fft.R +# Part of the R package, http://www.R-project.org +# +# Copyright (C) 1995-2012 The R Core Team +# +# This program is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 2 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# A copy of the GNU General Public License is available at +# http://www.r-project.org/Licenses/ + +fft <- function(z, inverse=FALSE) .Call(C_fft, z, inverse) + +#mvfft <- function(z, inverse=FALSE) .Call(C_mvfft, z, inverse) +# +#nextn <- function(n, factors=c(2,3,5)) .Call(C_nextn, n, factors) +# +#convolve <- function(x, y, conj=TRUE, type=c("circular","open","filter")) +#{ +# type <- match.arg(type) +# n <- length(x) +# ny <- length(y) +# Real <- is.numeric(x) && is.numeric(y) +# ## switch(type, circular = ..., ) +# if(type == "circular") { +# if(ny != n) +# stop("length mismatch in convolution") +# } +# else { ## "open" or "filter": Pad with zeros +# n1 <- ny - 1 +# x <- c(rep.int(0, n1), x) +# n <- length(y <- c(y, rep.int(0, n - 1)))# n = nx+ny-1 +# } +# x <- fft(fft(x)* (if(conj)Conj(fft(y)) else fft(y)), inverse=TRUE) +# if(type == "filter") +# (if(Real) Re(x) else x)[-c(1L:n1, (n-n1+1L):n)]/n +# else +# (if(Real) Re(x) else x)/n +#} + diff --git a/com.oracle.truffle.r.nodes/src/com/oracle/truffle/r/nodes/builtin/stats/R/init.R b/com.oracle.truffle.r.nodes/src/com/oracle/truffle/r/nodes/builtin/stats/R/init.R new file mode 100644 index 0000000000..d9663a01b0 --- /dev/null +++ b/com.oracle.truffle.r.nodes/src/com/oracle/truffle/r/nodes/builtin/stats/R/init.R @@ -0,0 +1,23 @@ +# Copyright (c) 2014, Oracle and/or its affiliates. All rights reserved. +# DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER. +# +# This code is free software; you can redistribute it and/or modify it +# under the terms of the GNU General Public License version 2 only, as +# published by the Free Software Foundation. +# +# This code is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +# version 2 for more details (a copy is included in the LICENSE file that +# accompanied this code). +# +# You should have received a copy of the GNU General Public License version +# 2 along with this work; if not, write to the Free Software Foundation, +# Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA. +# +# Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA +# or visit www.oracle.com if you need additional information or have any +# questions. + +# object defined in the stats package describing native fft function +C_fft <- list(name="fft") \ No newline at end of file diff --git a/com.oracle.truffle.r.nodes/src/com/oracle/truffle/r/nodes/builtin/stats/R/package-info.java b/com.oracle.truffle.r.nodes/src/com/oracle/truffle/r/nodes/builtin/stats/R/package-info.java new file mode 100644 index 0000000000..6227322707 --- /dev/null +++ b/com.oracle.truffle.r.nodes/src/com/oracle/truffle/r/nodes/builtin/stats/R/package-info.java @@ -0,0 +1,28 @@ +/* + * Copyright (c) 2014, Oracle and/or its affiliates. All rights reserved. + * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER. + * + * This code is free software; you can redistribute it and/or modify it + * under the terms of the GNU General Public License version 2 only, as + * published by the Free Software Foundation. + * + * This code is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * version 2 for more details (a copy is included in the LICENSE file that + * accompanied this code). + * + * You should have received a copy of the GNU General Public License version + * 2 along with this work; if not, write to the Free Software Foundation, + * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA. + * + * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA + * or visit www.oracle.com if you need additional information or have any + * questions. + */ +/** + * This "package" contains R sources that correspond to (some of) the R functions + * in the "stats" package. They are loaded using the {@link java.lang.Class#getResource} + * mechanism on system startup. + */ +package com.oracle.truffle.r.nodes.builtin.stats.R; \ No newline at end of file diff --git a/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/DCF.java b/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/DCF.java index 9ce4ec2d8a..d3ada37db2 100644 --- a/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/DCF.java +++ b/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/DCF.java @@ -22,7 +22,6 @@ */ package com.oracle.truffle.r.runtime; -import java.io.*; import java.util.*; /** diff --git a/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/data/RComplexVector.java b/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/data/RComplexVector.java index 0dc306dcd6..396406b989 100644 --- a/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/data/RComplexVector.java +++ b/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/data/RComplexVector.java @@ -84,6 +84,14 @@ public final class RComplexVector extends RVector implements RAbstractComplexVec return copy; } + /** + * Intended for external calls where a copy is not needed. WARNING: think carefully before using + * this method rather than {@link #getDataCopy()}. + */ + public double[] getDataWithoutCopying() { + return data; + } + public RComplexVector copyWithNewDimensions(int[] newDimensions) { return RDataFactory.createComplexVector(data, isComplete(), newDimensions); } diff --git a/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/ffi/RDerivedRFFI.java b/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/ffi/RDerivedRFFI.java index aa698758ee..14ed8df6fc 100644 --- a/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/ffi/RDerivedRFFI.java +++ b/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/ffi/RDerivedRFFI.java @@ -36,5 +36,5 @@ public interface RDerivedRFFI { // Checkstyle: stop method name void fft_factor(int n, int[] pmaxf, int[] pmaxp); - int fft_work(double[] a, double[] b, int nseg, int n, int nspn, int isn, double[] work, int[] iwork); + int fft_work(double[] a, int nseg, int n, int nspn, int isn, double[] work, int[] iwork); } diff --git a/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/ffi/jnr/JNR_RFFIFactory.java b/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/ffi/jnr/JNR_RFFIFactory.java index 8fe543a7d4..ebaf7a689f 100644 --- a/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/ffi/jnr/JNR_RFFIFactory.java +++ b/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/ffi/jnr/JNR_RFFIFactory.java @@ -448,7 +448,7 @@ public class JNR_RFFIFactory extends RFFIFactory implements RFFI, BaseRFFI, RDer // @formatter:off void fft_factor(@In int[] n, int[] pmaxf, int[] pmaxp); - int fft_work(double[] a, double[] b, @In int[] nseg, @In int[] n, @In int[] nspn, @In int[] isn, double[] work, int[] iwork); + int fft_work(double[] a, @In int[] nseg, @In int[] n, @In int[] nspn, @In int[] isn, double[] work, int[] iwork); } private static class FFTProvider { @@ -487,12 +487,12 @@ public class JNR_RFFIFactory extends RFFIFactory implements RFFI, BaseRFFI, RDer static int[] isn = new int[1]; } - public int fft_work(double[] a, double[] b, int nseg, int n, int nspn, int isn, double[] work, int[] iwork) { + public int fft_work(double[] a, int nseg, int n, int nspn, int isn, double[] work, int[] iwork) { RefScalars_fft_work.n[0] = n; RefScalars_fft_work.nseg[0] = nseg; RefScalars_fft_work.nspn[0] = nspn; RefScalars_fft_work.isn[0] = isn; - return fft().fft_work(a, b, RefScalars_fft_work.nseg, RefScalars_fft_work.n, RefScalars_fft_work.nspn, RefScalars_fft_work.isn, work, iwork); + return fft().fft_work(a, RefScalars_fft_work.nseg, RefScalars_fft_work.n, RefScalars_fft_work.nspn, RefScalars_fft_work.isn, work, iwork); } -- GitLab