 |
 |
 |
 |
 |
 |
 |
 |
When you use the PROTOTYPE comment to generate C code wrappers, you cannot take advantage of many of the features of the yorick language. In this example, there is no easy way to pass an indeterminate number of arguments to your function. When you write your own C code wrapper, all such limitations disappear. The price is that you need to use the yorick-specific C functions declared in ‘yapi.h’ to interact with the intepreter.
Since you no longer need an interpreted wrapper (although you might still choose to have one – some tasks are easier, occasionally dramatically easier, in interpreted code than in compiled code), your inteface definition file ‘rtleg.i’ becomes trivial:
plug_in, "rtleg"; extern rtleg; |
But your C source becomes much more complicated:
#include "yapi.h"
#include <math.h>
extern void Y_rtleg(int argc);
void
Y_rtleg(int argc)
{
int iarg;
long i, n, dimsr[Y_DIMSIZE], dims[Y_DIMSIZE];
double *r, *rr, **px, *x, *a, xval;
int *flag, flagr;
if (argc < 2) y_error("rtleg needs at least two arguments");
/* get some scratch space to deal with the indeterminate
* argument count, this increases iarg for arguments by 1
*/
px = ypush_scratch((argc-1)*(sizeof(double)+sizeof(int)), 0);
flag = (int *)&px[argc];
/* get r argument as double */
r = ygeta_d(argc, 0, dimsr);
/* first pass through remaining arguments checks they are
* conformable, and builds the result dimension list
*/
for (flagr=0, iarg=argc-1 ; iarg>0 ; iarg--) {
flagr |= flag[iarg-1] = yarg_dims(iarg, dims, dimsr);
if (flag[iarg-1] < 0)
y_errorn("rtleg argument %ld not conformable", argc-iarg);
}
/* broadcast r argument to result shape if necessary */
if (flagr & 1) yarg_bcast(argc, dimsr);
/* second pass broadcasts all arguments to result shape
* and converts to type double
*/
for (iarg=argc-1 ; iarg>0 ; iarg--) {
if (flag[iarg-1] & 2) yarg_bcast(iarg, dimsr);
px[iarg-1] = ygeta_d(iarg, &n, 0);
flag[iarg-1] = yarg_scratch(iarg);
if (flag[iarg-1]) flagr = iarg;
}
/* push new array and fill with reciprocal of r */
rr = ypush_d(dimsr);
for (i=0 ; i<n ; i++)
if (r[i] > 0.) rr[i] = 1./r[i];
else if (r[i] < 0.) rr[i] = -1./r[i];
else rr[i] = -1.;
if (flagr) { /* use scratch x as accumulator */
a = px[flagr-1];
for (i=0 ; i<n ; i++) {
xval = a[i]*rr[i];
a[i] = -xval*xval;
}
} else { /* need to create accumulator on stack */
a = ypush_d(dimsr);
}
/* compute 1 - sum((x/r)^2) */
for (i=0 ; i<n ; i++) a[i] += (rr[i]<0.)? 0. : 1.;
for (iarg=0 ; iarg<argc-1 ; iarg++) {
x = px[iarg];
if (x == a) continue;
for (i=0 ; i<n ; i++) {
xval = x[i]*rr[i];
a[i] -= xval*xval;
}
}
/* compute final result */
for (i=0 ; i<n ; i++) a[i] = r[i]*sqrt(a[i]);
/* discard overlying stack elements to leave result on top */
if (flagr) yarg_drop(flagr+1);
}
|
The ypush_scratch function enables you to create scratch space that is safe against asynchronous interrupt. You have to write every compiled function assuming that the user will type Control-C and interrupt your function. If you used malloc to obtain scratch space, and your function were asynchronously interrupted before the matching free, your routine would leak memory. In this case, the scratch space is only required to cope with the fact that the function accepts an indeterminate number of arguments.
The implementation unavoidably requires a temporary array to hold the reciprocal of the first argument, which is used to normalize the remaining arguments to prevent overflows or erronious underflows when the other arguments are squared. However, there is an opportunity to reuse an argument array which happens to be a temporary (either because it is the result of an expression, or because it had to be type converted or broadcast) as the result. Taking advantage of that optimization also costs a few additional lines of source code towards the end of the function.
What you get for this modest complexity is an rtleg function that is equal in performance and features to the built in abs function.
|
|
|
|
|
|
|
|
