Dear All,
could someone please shed some light on the use of the .C or .Fortran function:
I am trying load and running on R the following function
// psi.cpp -- psi function for real arguments.
// Algorithms and coefficient values from "Computation of Special
// Functions", Zhang and Jin, John Wiley and Sons, 1996.
//
// (C) 2003, C. Bond. All rights reserved.
//
// Returns psi function for real argument 'x'.
// NOTE: Returns 1e308 for argument 0 or a negative integer.
//
#include <iostream.h> // or <stdio.h>?
#include <math.h>
#define el 0.5772156649015329
double psi(double x)
{
double s,ps,xa,x2;
int n,k;
static double a[] = {
-0.8333333333333e-01,
0.83333333333333333e-02,
-0.39682539682539683e-02,
0.41666666666666667e-02,
-0.75757575757575758e-02,
0.21092796092796093e-01,
-0.83333333333333333e-01,
0.4432598039215686};
xa = fabs(x);
s = 0.0;
if ((x == (int)x) && (x <= 0.0)) {
ps = 1e308;
return ps;
}
if (xa == (int)xa) {
n = xa;
for (k=1;k<n;k++) {
s += 1.0/k;
}
ps = s-el;
}
else if ((xa+0.5) == ((int)(xa+0.5))) {
n = xa-0.5;
for (k=1;k<=n;k++) {
s += 1.0/(2.0*k-1.0);
}
ps = 2.0*s-el-1.386294361119891;
}
else {
if (xa < 10.0) {
n = 10-(int)xa;
for (k=0;k<n;k++) {
s += 1.0/(xa+k);
}
xa += n;
}
x2 = 1.0/(xa*xa);
ps = log(xa)-0.5/xa+x2*(((((((a[7]*x2+a[6])*x2+a[5])*x2+
a[4])*x2+a[3])*x2+a[2])*x2+a[1])*x2+a[0]);
ps -= s;
}
if (x < 0.0)
ps = ps - M_PI*cos(M_PI*x)/sin(M_PI*x)-1.0/x;
return ps;
}
However, when applicated the codes
>digamma(-0.9)
[1] -9.312644 OK!
But when
> dyn.load("psi.so")
> out<-.C("psi", as.double(-0.9))
Erro en .C("psi", as.double(-0.9)) : C symbol name "psi" not in load table
> out
Error: objeto "out" no encontrado
More information on OS:
> version
_
platform i486-pc-linux-gnu
arch i486
os linux-gnu
system i486, linux-gnu
status
major 2
minor 4.1
year 2006
month 12
day 18
svn rev 40228
language R
version.string R version 2.4.1 (2006-12-18)
Many thanks for any insight or comments.
Bernardo Lagos A.
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
