1 /****************************************************************
2 Copyright (C) 1997-1998 Lucent Technologies
3 All Rights Reserved
4
5 Permission to use, copy, modify, and distribute this software and
6 its documentation for any purpose and without fee is hereby
7 granted, provided that the above copyright notice appear in all
8 copies and that both that the copyright notice and this
9 permission notice and warranty disclaimer appear in supporting
10 documentation, and that the name of Lucent or any of its entities
11 not be used in advertising or publicity pertaining to
12 distribution of the software without specific, written prior
13 permission.
14
15 LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
16 INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
17 IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
18 SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
19 WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
20 IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
21 ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
22 THIS SOFTWARE.
23 ****************************************************************/
24
25 /* sample funcadd */
26
27 #include "math.h" /* for sqrt */
28 #include "funcadd.h" /* includes "stdio1.h" */
29
30 char *ix_details_ASL[] = {0}; /* no -i command-line option */
31
32 static real
ginv(arglist * al)33 ginv(arglist *al) /* generalized inverse of a single argument */
34 {
35 real x = al->ra[0];
36 x = x ? 1./x : 0.;
37 if (al->derivs) {
38 *al->derivs = -x*x;
39 if (al->hes)
40 *al->hes = -2.*x * *al->derivs;
41 }
42 return x;
43 }
44
45 static real
myhypot(arglist * al)46 myhypot(arglist *al) /* sqrt(x*x + y*y) */
47 {
48 real *d, *h, rv, x, x0, y, y0;
49
50 x = x0 = al->ra[0];
51 y = y0 = al->ra[1];
52
53 if (x < 0.)
54 x = -x;
55 if (y < 0.)
56 y = -y;
57 rv = x;
58 if (x < y) {
59 rv = y;
60 y = x;
61 x = rv;
62 }
63 if (rv) {
64 y /= x;
65 rv = x * sqrt(1. + y*y);
66 if (d = al->derivs) {
67 d[0] = x0 / rv;
68 d[1] = y0 / rv;
69 if (h = al->hes) {
70 h[0] = d[1]*d[1] / rv;
71 h[1] = -d[0]*d[1] / rv;
72 h[2] = d[0]*d[0] / rv;
73 }
74 }
75 }
76 else if (d = al->derivs) {
77 d[0] = d[1] = 0;
78 if (h = al->hes)
79 h[0] = h[1] = h[2] = 0;
80 }
81 return rv;
82 }
83
84 static real
mean(register arglist * al)85 mean(register arglist *al) /* mean of arbitrarily many arguments */
86 {
87 real x, z;
88 real *d, *de, *ra;
89 int *at, i, j, n;
90 char *se;
91 const char *sym;
92 AmplExports *ae = al->AE; /* for fprintf and strtod */
93
94 if ((n = al->n) <= 0)
95 return 0;
96 at = al->at;
97 ra = al->ra;
98 d = de = al->derivs;
99 x = 0.;
100 for(i = 0; i < n;)
101 if ((j = at[i++]) >= 0) {
102 x += ra[j];
103 ++de;
104 }
105 else {
106 x += z = strtod(sym = al->sa[-(j+1)], &se);
107 if (*se) {
108 fprintf(Stderr,
109 "mean treating arg %d = \"%s\" as %.g\n",
110 i, sym, z);
111 }
112 }
113 if (d) {
114 z = 1. / n;
115 while(d < de)
116 *d++ = z;
117 /* The Hessian is == 0 and, if needed, has been */
118 /* initialized to 0. */
119 }
120 return x / n;
121 }
122
123 void
funcadd(AmplExports * ae)124 funcadd(AmplExports *ae){
125 /* Insert calls on addfunc here... */
126 /* Arg 3, called argtype, can be 0 or 1:
127 * 0 ==> force all arguments to be numeric
128 * 1 ==> pass both symbolic and numeric arguments.
129 *
130 * Arg 4, called nargs, is interpretted as follows:
131 * >= 0 ==> the function has exactly nargs arguments
132 * <= -1 ==> the function has >= -(nargs+1) arguments.
133 *
134 * Arg 5, funcinfo, is passed to the functions in struct arglist;
135 * it is not used in these examples, so we just pass 0.
136 */
137
138 addfunc("ginv", (ufunc*)ginv, 0, 1, 0);
139 addfunc("hypot", (ufunc*)myhypot, 0, 2, 0);
140 addfunc("mean", (ufunc*)mean, 1, -1, 0);
141 }
142