raushon
需要在 VC++ 中调用R的函数。 按照R Manual --- 《Writing R Extensions》的第6部分 ：
The R API:entry points for Ccode 的说明， 编译总是通过不了。VC++代码如下：
<br />
#include <R.h><br />
<br />
void main()<br />
{<br />
double *X;<br />
X = (double *) R_alloc(10, sizeof(double));<br />
}<br />
编译出错：
--------------------Configuration: callR - Win32 Debug--------------------
Compiling...
callR.c
Linking...
callR.obj : error LNK2001: unresolved external symbol _R_alloc
Debug/callR.exe : fatal error LNK1120: 1 unresolved externals
Error executing link.exe.
Creating browse info file...
callR.exe - 2 error(s), 0 warning(s)
VC++6.0的编译环境 和 链接库 (Tools | Options | Directories | include...)的设置似乎是没有问题的。
这应该是动态链接库的问题，但是我实在是 不知道该怎么解决，把R.dll拷贝到Debug目录下也是不行。
我的R是2.50RC版本的，各位大虾帮帮忙看如何解决C++调用R的问题，谢谢了！
karlqi
这个我已经实现了，部分代码如下，明眼人一看就清楚了。当然了，提取时不要忘记提取底层函数和一些个全局变量定义。
p.s.:怎么找相关函数？在R下用UltraEdit搜一下，轻松搞定！代码比较长，就隐藏了，想看的回复下就好了，不然这么多行的代码很晕的呢，隐了吧。
//////////////////////////RNG.h/////////////////////////////////////
#defined(RNG_EXPORTS)
#define RNG_API __declspec(dllexport)
#else
#define RNG_API __declspec(dllimport)
#endif
#ifdef __cplusplus
extern "C" {
#endif
RNG_API void __stdcall rmultinom(long* n, double* prob, long* size, long* K, long* randomNum);
#ifdef __cplusplus
}
#endif
//////////////////////////RNG.cpp/////////////////////////////////////
#include "RNG.h"
...
void __stdcall rmultinom(long* n, double* prob, long* size, long* K, long* randomNum)
{
long* rN = randomNum;
for(long i = 0; i < n[0]; i++)
{
__rmultinom(size[0], prob, K[0], rN);
for(long j = 0; j < K[0]; j++){
randomNum[i*K[0]+j] = rN[j];
}
}
return;
}
///以下来自R package：
void __rmultinom(long n, double* prob, long K, long* rN)
{
long k;
double pp, p_tot = 0.;
for(k = 0; k < K; k++) {
pp = prob[k];
p_tot += pp;
rN[k] = 0;
}
if (n == 0) return;
if (K == 1 && p_tot == 0.) return;
for(k = 0; k < K-1; k++) {
if(prob[k]) {
pp = prob[k] / p_tot;
rN[k] = ((pp < 1.) ? (long) __rbinom((double) n, pp) : n);
n -= rN[k];
}
else rN[k] = 0;
if(n <= 0) return;
p_tot -= prob[k];
}
rN[K-1] = n;
return;
}
///以下代码来自R
long __rbinom(double nin, double pp)
{
static double c, fm, npq, p1, p2, p3, p4, qn;
static double xl, xll, xlr, xm, xr;
static double psave = -1.0;
static long nsave = -1;
static long m;
double f, f1, f2, u, v, w, w2, x, x1, x2, z, z2;
double p, q, np, g, r, al, alv, amaxp, ffm, ynorm;
long i,ix,k,n;
n = (long)floor(nin + 0.5);
if (n == 0 || pp == 0.) return 0;
if (pp == 1.) return n;
p = fmin2(pp, 1. - pp);
q = 1. - p;
np = n * p;
r = p / q;
g = r * (n + 1);
/* Setup, perform only when parameters change [using static (globals): */
/* FIXING: Want this thread safe
-- use as little (thread globals) as possible
*/
if (pp != psave || n != nsave) {
psave = pp;
nsave = n;
if (np < 30.0) {
/* inverse cdf logic for mean less than 30 */
qn = pow(q, (double) n);
goto L_np_small;
} else {
ffm = np + p;
m = (long)ffm;
fm = m;
npq = np * q;
p1 = (long)(2.195 * sqrt(npq) - 4.6 * q) + 0.5;
xm = fm + 0.5;
xl = xm - p1;
xr = xm + p1;
c = 0.134 + 20.5 / (15.3 + fm);
al = (ffm - xl) / (ffm - xl * p);
xll = al * (1.0 + 0.5 * al);
al = (xr - ffm) / (xr * q);
xlr = al * (1.0 + 0.5 * al);
p2 = p1 * (1.0 + c + c);
p3 = p2 + c / xll;
p4 = p3 + c / xlr;
}
} else if (n == nsave) {
if (np < 30.0)
goto L_np_small;
}
/*-------------------------- np = n*p >= 30 : ------------------- */
repeat {
u = unif_rand() * p4;
v = unif_rand();
/* triangular region */
if (u <= p1) {
ix = (long)(xm - p1 * v + u);
goto finis;
}
/* parallelogram region */
if (u <= p2) {
x = xl + (u - p1) / c;
v = v * c + 1.0 - fabs(xm - x) / p1;
if (v > 1.0 || v <= 0.)
continue;
ix = (long)x;
} else {
if (u > p3) { /* right tail */
ix = (long)(xr - log(v) / xlr);
if (ix > n)
continue;
v = v * (u - p3) * xlr;
} else {/* left tail */
ix = (long)(xl + log(v) / xll);
if (ix < 0)
continue;
v = v * (u - p2) * xll;
}
}
/* determine appropriate way to perform accept/reject test */
k = abs(ix - m);
if (k <= 20 || k >= npq / 2 - 1) {
/* explicit evaluation */
f = 1.0;
if (m < ix) {
for (i = m + 1; i <= ix; i++)
f *= (g / i - r);
} else if (m != ix) {
for (i = ix + 1; i <= m; i++)
f /= (g / i - r);
}
if (v <= f)
goto finis;
} else {
/* squeezing using upper and lower bounds on log(f(x)) */
amaxp = (k / npq) * ((k * (k / 3. + 0.625) + 0.1666666666666) / npq + 0.5);
ynorm = -k * k / (2.0 * npq);
alv = log(v);
if (alv < ynorm - amaxp)
goto finis;
if (alv <= ynorm + amaxp) {
/* stirling's formula to machine accuracy */
/* for the final acceptance/rejection test */
x1 = ix + 1;
f1 = fm + 1.0;
z = n + 1 - fm;
w = n - ix + 1.0;
z2 = z * z;
x2 = x1 * x1;
f2 = f1 * f1;
w2 = w * w;
if (alv <= xm * log(f1 / x1) + (n - m + 0.5) * log(z / w) + (ix - m) * log(w * p / (x1 * q)) + (13860.0 - (462.0 - (132.0 - (99.0 - 140.0 / f2) / f2) / f2) / f2) / f1 / 166320.0 + (13860.0 - (462.0 - (132.0 - (99.0 - 140.0 / z2) / z2) / z2) / z2) / z / 166320.0 + (13860.0 - (462.0 - (132.0 - (99.0 - 140.0 / x2) / x2) / x2) / x2) / x1 / 166320.0 + (13860.0 - (462.0 - (132.0 - (99.0 - 140.0 / w2) / w2) / w2) / w2) / w / 166320.)
goto finis;
}
}
}
L_np_small:
/*---------------------- np = n*p < 30 : ------------------------- */
repeat {
ix = 0;
f = qn;
u = unif_rand();
repeat {
if (u < f)
goto finis;
if (ix > 110)
break;
u -= f;
ix++;
f *= (g / ix - r);
}
}
finis:
if (psave > 0.5)
ix = n - ix;
return ix;
}
加油！加油！
