From 2198ea2ab503c9c647803503e1a4c6bbf0155d9b Mon Sep 17 00:00:00 2001
From: Mick Jordan <mick.jordan@oracle.com>
Date: Fri, 10 Apr 2015 07:18:35 -0700
Subject: [PATCH] move fft.c to stats library, modify source from GnuR with sed
---
.hgignore | 1 +
com.oracle.truffle.r.native/Makefile | 3 +-
.../builtinlibs/src/fft.c | 887 ------------------
com.oracle.truffle.r.native/gnur/Makefile | 1 -
.../library/base/Makefile | 4 +-
com.oracle.truffle.r.native/library/lib.mk | 7 +-
.../library/stats/Makefile | 8 +
.../library/stats/src/ed_fft | 44 +
.../r/runtime/ffi/jnr/JNR_RFFIFactory.java | 5 +-
.../rpackages/rffi/testrffi/src/Makefile | 2 +-
10 files changed, 66 insertions(+), 896 deletions(-)
delete mode 100644 com.oracle.truffle.r.native/builtinlibs/src/fft.c
create mode 100644 com.oracle.truffle.r.native/library/stats/src/ed_fft
diff --git a/.hgignore b/.hgignore
index 8af9c89293..b106ca0ca5 100644
--- a/.hgignore
+++ b/.hgignore
@@ -91,6 +91,7 @@ R.tokens
findbugs.html
com.oracle.truffle.r.native/builtinlibs/lib/*
com.oracle.truffle.r.native/library/*/lib/*
+com.oracle.truffle.r.native/library/stats/src/fft.c
com.oracle.truffle.r.native/platform.mk
com.oracle.truffle.r.native/gnur/Makeconf.done
com.oracle.truffle.r.native/include/jni/include
diff --git a/com.oracle.truffle.r.native/Makefile b/com.oracle.truffle.r.native/Makefile
index c6abbfd958..d78f9fec75 100644
--- a/com.oracle.truffle.r.native/Makefile
+++ b/com.oracle.truffle.r.native/Makefile
@@ -1,5 +1,5 @@
#
-# Copyright (c) 2015, Oracle and/or its affiliates. All rights reserved.
+# Copyright (c) 2014, 2015, 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
@@ -42,4 +42,3 @@ clean:
$(MAKE) -C library clean
$(MAKE) -C run clean
$(MAKE) -C gnur clean
-
\ No newline at end of file
diff --git a/com.oracle.truffle.r.native/builtinlibs/src/fft.c b/com.oracle.truffle.r.native/builtinlibs/src/fft.c
deleted file mode 100644
index 194aab1ca8..0000000000
--- a/com.oracle.truffle.r.native/builtinlibs/src/fft.c
+++ /dev/null
@@ -1,887 +0,0 @@
-/*
- * This material is distributed under the GNU General Public License
- * Version 2. You may review the terms of this license at
- * http://www.gnu.org/licenses/gpl-2.0.html
- *
- * Copyright (c) 1995, 1996, 1997 Robert Gentleman and Ross Ihaka
- * Copyright (c) 1995-2014, The R Core Team
- * Copyright (c) 2002-2008, The R Foundation
- * Copyright (c) 2014, 2015, Oracle and/or its affiliates
- *
- * All rights reserved.
- */
-
-// modifications required for standalone compilation
-
-//#ifdef HAVE_CONFIG_H
-//#include <config.h>
-//#endif
-
-#include <limits.h> /* for INT_MAX */
-#include <stddef.h> /* for size_t */
-#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.
- *
- * Update in R 3.1.0: nfac[20], increased array size. It is now possible to
- * factor any positive int n, up to 2^31 - 1.
- */
-
-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[20];
-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, sqrtk, kchanged;
-
- /* 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, ... */
- kchanged = 0;
- sqrtk = (int)sqrt(k);
- for(j = 3; j <= sqrtk; j += 2) {
- jj = j * j;
- while(k % jj == 0) {
- nfac[m_fac++] = j;
- k /= jj;
- kchanged = 1;
- }
- if (kchanged) {
- kchanged = 0;
- sqrtk = (int)sqrt(k);
- }
- }
-
- 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;
- }
- if (j > INT_MAX - 2)
- break;
- j = ((j+1)/2)*2 + 1;
- }
- while(j <= k);
- }
-
- if (m_fac <= kt+1)
- maxp = m_fac+kt+1;
- if (m_fac+kt > 20) { /* 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*(size_t)maxf], &work[3*(size_t)maxf],
- iwork, nfac);
-
- return TRUE;
-}
diff --git a/com.oracle.truffle.r.native/gnur/Makefile b/com.oracle.truffle.r.native/gnur/Makefile
index a4468ffc2a..9772ef1737 100644
--- a/com.oracle.truffle.r.native/gnur/Makefile
+++ b/com.oracle.truffle.r.native/gnur/Makefile
@@ -32,4 +32,3 @@ clean:
$(MAKE) -f Makefile.gnur clean
$(MAKE) -f Makefile.libs clean
$(MAKE) -f Makefile.platform clean
-
\ No newline at end of file
diff --git a/com.oracle.truffle.r.native/library/base/Makefile b/com.oracle.truffle.r.native/library/base/Makefile
index 7ca188ab76..40d2328561 100644
--- a/com.oracle.truffle.r.native/library/base/Makefile
+++ b/com.oracle.truffle.r.native/library/base/Makefile
@@ -23,14 +23,14 @@
RPROFILE := $(FASTR_LIBDIR)/base/R/Rprofile
RPROFILE_ORIG := $(RPROFILE).orig
-PKG_EXTRAS = $(RPROFILE)
+LIB_PKG_POST = $(RPROFILE)
include ../lib.mk
# edit the Rprofile to add fastr as a default package
# sed's edit in place option with backup is not portable
-$(PKG_EXTRAS): $(RPROFILE_ORIG)
+$(LIB_PKG_POST): $(RPROFILE_ORIG)
$(RPROFILE_ORIG):
cp $(RPROFILE) $(RPROFILE_ORIG)
diff --git a/com.oracle.truffle.r.native/library/lib.mk b/com.oracle.truffle.r.native/library/lib.mk
index cccf72cc7e..4e129f6588 100644
--- a/com.oracle.truffle.r.native/library/lib.mk
+++ b/com.oracle.truffle.r.native/library/lib.mk
@@ -27,6 +27,9 @@
# and overwrites the default. The libraries are stored in the directory denoted
# FASTR_LIBDIR.
+# A package that requires special processing before the library is built should
+# define LIB_PKG_PRE and for post processing define LIB_PKG_POST
+
include $(TOPDIR)/platform.mk
.PHONY: all clean cleanlib cleanobj force libr libcommon
@@ -59,9 +62,9 @@ INCLUDES := $(JNI_INCLUDES) $(FFI_INCLUDES)
PKGDIR := $(FASTR_LIBDIR)/$(PKG)
ifneq ($(C_SOURCES),)
-all: libcommon $(LIB_PKG) $(PKG_EXTRAS)
+all: $(LIB_PKG_PRE) libcommon $(LIB_PKG) $(LIB_PKG_POST)
else
-all: libcommon
+all: $(LIB_PKG_PRE) libcommon $(LIB_PKG_POST)
endif
libcommon: $(PKGDIR)
diff --git a/com.oracle.truffle.r.native/library/stats/Makefile b/com.oracle.truffle.r.native/library/stats/Makefile
index f128fcf8a4..d973b76747 100644
--- a/com.oracle.truffle.r.native/library/stats/Makefile
+++ b/com.oracle.truffle.r.native/library/stats/Makefile
@@ -21,4 +21,12 @@
# questions.
#
+LIB_PKG_PRE = lib/fft.o
+GNUR_FFT := $(GNUR_DIR)/src/library/stats/src/fft.c
+
include ../lib.mk
+
+C_SOURCES := $(C_SOURCES) src/fft.c
+
+src/fft.c: $(GNUR_FFT) src/ed_fft
+ ed $(GNUR_FFT) < src/ed_fft
diff --git a/com.oracle.truffle.r.native/library/stats/src/ed_fft b/com.oracle.truffle.r.native/library/stats/src/ed_fft
new file mode 100644
index 0000000000..6ac2d1f3d4
--- /dev/null
+++ b/com.oracle.truffle.r.native/library/stats/src/ed_fft
@@ -0,0 +1,44 @@
+/HAVE_CONFIG.H/
+s!^!//!
++1
+s!^!//!
++1
+s!^!//!
+/include <Rmath.h>/
+s!^!//!
++1
+s!^!//!
+a
+#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
+.
+/^void fft_factor(int n/
+s!int n!int *nptr!
+i
+// signature modification (pointers to scalars) required to enable JNR call
+.
+/int j/
+i
+ int n = *nptr;
+.
+/fft_work(double/
+d
+i
+// 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,
+.
+/int nf/
+i
+ int nseg = *nsegptr;
+ int n = *nptr;
+ int nspn = *nspnptr;
+ int isn = *isnptr;
+ double *b=&(a[1]);
+.
+w src/fft.c
diff --git a/com.oracle.truffle.r.runtime.ffi/src/com/oracle/truffle/r/runtime/ffi/jnr/JNR_RFFIFactory.java b/com.oracle.truffle.r.runtime.ffi/src/com/oracle/truffle/r/runtime/ffi/jnr/JNR_RFFIFactory.java
index d386cc3f33..795edd3d01 100644
--- a/com.oracle.truffle.r.runtime.ffi/src/com/oracle/truffle/r/runtime/ffi/jnr/JNR_RFFIFactory.java
+++ b/com.oracle.truffle.r.runtime.ffi/src/com/oracle/truffle/r/runtime/ffi/jnr/JNR_RFFIFactory.java
@@ -33,6 +33,7 @@ import jnr.posix.*;
import com.oracle.truffle.api.CompilerDirectives.TruffleBoundary;
import com.oracle.truffle.r.runtime.*;
import com.oracle.truffle.r.runtime.ffi.*;
+import com.oracle.truffle.r.runtime.ffi.DLL.DLLInfo;
/**
* JNR-based factory.
@@ -389,7 +390,9 @@ public class JNR_RFFIFactory extends RFFIFactory implements RFFI, BaseRFFI, Stat
@TruffleBoundary
private static Stats createAndLoadLib() {
- return LibraryLoader.create(Stats.class).load("appl");
+ // fft is in the stats package .so
+ DLLInfo dllInfo = DLL.findLibraryContainingSymbol("fft");
+ return LibraryLoader.create(Stats.class).load(dllInfo.path);
}
static Stats fft() {
diff --git a/com.oracle.truffle.r.test/rpackages/rffi/testrffi/src/Makefile b/com.oracle.truffle.r.test/rpackages/rffi/testrffi/src/Makefile
index 1686772388..73ab7c69aa 100644
--- a/com.oracle.truffle.r.test/rpackages/rffi/testrffi/src/Makefile
+++ b/com.oracle.truffle.r.test/rpackages/rffi/testrffi/src/Makefile
@@ -1,5 +1,5 @@
#
-# Copyright (c) 2014, Oracle and/or its affiliates. All rights reserved.
+# Copyright (c) 2014, 2015, 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
--
GitLab