1 /****************************************************************
2 Copyright (C) 1997-1999, 2000-2001 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 #include "nlp2.h"
26 #include "errchk.h"
27
28 #ifdef __cplusplus
29 extern "C" {
30 #endif
31
32 #ifdef MULTIPLE_THREADS
33 #define MTNot_Used(x) Not_Used(x)
34 #else
35 #define MTNot_Used(x) /*nothing*/
36 #define asl cur_ASL
37 #endif
38
39 static void
40 #ifdef KR_headers
jmp_check(J,jv)41 jmp_check(J, jv) Jmp_buf *J; int jv;
42 #else
43 jmp_check(Jmp_buf *J, int jv)
44 #endif
45 {
46 if (J)
47 longjmp(J->jb, jv);
48 }
49
50 static void
51 #ifdef KR_headers
Errprint(msg)52 Errprint(msg) char *msg;
53 #else
54 Errprint(char *msg)
55 #endif
56 {
57 #ifndef NO_PERROR
58 if (errno)
59 fprintf(Stderr, "\n%s: %s.\n", msg, strerror(errno));
60 else
61 #endif
62 fprintf(Stderr, "%s.\n", msg);
63 fflush(Stderr);
64 }
65
66 static void
67 #ifdef KR_headers
68 introuble(who, a, jv K_ASL) char *who; real a; int jv; D_ASL
69 #else
70 introuble(char *who, real a, int jv A_ASL)
71 #endif
72 {
73 char buf[64];
74 jmp_check(err_jmp, jv);
75 report_where(asl);
76 snprintf(buf, sizeof(buf), "can't evaluate %s(%g)", who, a);
77 Errprint(buf);
78 jmp_check(err_jmp1, jv);
79 exit(1);
80 }
81
82 static void
83 #ifdef KR_headers
84 introuble2(who, a, b, jv K_ASL) char *who; real a; real b; int jv; D_ASL
85 #else
86 introuble2(char *who, real a, real b, int jv A_ASL)
87 #endif
88 {
89 char buf[96];
90 jmp_check(err_jmp, jv);
91 report_where(asl);
92 snprintf(buf, sizeof(buf), "can't evaluate %s(%g,%g)", who, a, b);
93 Errprint(buf);
94 jmp_check(err_jmp1, jv);
95 exit(1);
96 }
97
98 real
99 #ifdef KR_headers
100 f_OPPLUS(e K_ASL) expr *e; D_ASL
101 #else
102 f_OPPLUS(expr *e A_ASL)
103 #endif
104 {
105 expr *L;
106 /* e->dL = e->dR = 1.; */
107 L = e->L.e;
108 e = e->R.e;
109 return (*L->op)(L K_ASL) + (*e->op)(e K_ASL);
110 }
111
112 static real
113 #ifdef KR_headers
114 f_OPMINUS(e K_ASL) expr *e; D_ASL
115 #else
116 f_OPMINUS(expr *e A_ASL)
117 #endif
118 {
119 expr *L;
120 /* e->dL = 1.; */
121 /* e->dR = -1.; */
122 L = e->L.e;
123 e = e->R.e;
124 return (*L->op)(L K_ASL) - (*e->op)(e K_ASL);
125 }
126
127 static real
128 #ifdef KR_headers
129 f_OPMULT(e K_ASL) expr *e; D_ASL
130 #else
131 f_OPMULT(expr *e A_ASL)
132 #endif
133 {
134 expr *e1, *e2;
135 e1 = e->L.e;
136 e2 = e->R.e;
137 return (e->dR = (*e1->op)(e1 K_ASL)) * (e->dL = (*e2->op)(e2 K_ASL));
138 }
139
140 static void
141 #ifdef KR_headers
142 zero_div(L, op K_ASL) real L; char *op; D_ASL
143 #else
144 zero_div(real L, char *op A_ASL)
145 #endif
146 {
147 errno_set(EDOM);
148 jmp_check(err_jmp, 1);
149 report_where(asl);
150 fprintf(Stderr, "can't compute %g%s0.\n", L, op);
151 fflush(Stderr);
152 jmp_check(err_jmp1, 1);
153 exit(1);
154 }
155
156 static real
157 #ifdef KR_headers
158 f_OPDIV(e K_ASL) expr *e; D_ASL
159 #else
160 f_OPDIV(expr *e A_ASL)
161 #endif
162 {
163 expr *e1;
164 real L, R, dL, rv;
165 e1 = e->L.e;
166 L = (*e1->op)(e1 K_ASL);
167 e1 = e->R.e;
168 if (!(R = (*e1->op)(e1 K_ASL)))
169 zero_div(L, "/" K_ASL);
170 rv = L / R;
171 if (want_deriv) {
172 e->dR = -rv * (dL = e->dL = 1. / R);
173 /*e->dL2 = 0;*/ /* !! could special-case division to avoid */
174 /* !! this and two multiplies in hv_back */
175 e->dLR = -dL*dL;
176 dL *= e->dR;
177 e->dR2 = -(dL + dL);
178 }
179 return rv;
180 }
181
182 static real
183 #ifdef KR_headers
184 f_OPREM(e K_ASL) expr *e; D_ASL
185 #else
186 f_OPREM(expr *e A_ASL)
187 #endif
188 {
189 expr *e1;
190 real L, R, rv;
191 /* e->dL = 1.; */
192 e1 = e->L.e;
193 L = (*e1->op)(e1 K_ASL);
194 e1 = e->R.e;
195 R = (*e1->op)(e1 K_ASL);
196 rv = fmod(L,R);
197 if (errchk(rv))
198 introuble2("fmod",L,R,1 K_ASL);
199 else if (want_deriv) {
200 e->dR = (rv - L) / R;
201 e->dR2 = 0;
202 }
203 return rv;
204 }
205
206 real
207 #ifdef KR_headers
208 f_OPPOW(e K_ASL) expr *e; D_ASL
209 #else
210 f_OPPOW(expr *e A_ASL)
211 #endif
212 {
213 expr *e1;
214 real L, R, rv, xlog, xym1;
215 e1 = e->L.e;
216 L = (*e1->op)(e1 K_ASL);
217 e1 = e->R.e;
218 R = (*e1->op)(e1 K_ASL);
219 rv = mypow(L,R);
220 if (errchk(rv))
221 introuble2("pow",L,R,1 K_ASL);
222 if (want_deriv) {
223 if (L < 0.)
224 goto badpow;
225 if (L > 0.) {
226 xym1 = rv / L;
227 xlog = log(L);
228 e->dL = R*xym1;
229 e->dR = xlog * rv;
230 e->dL2 = (R-1.)*(e->dL/L);
231 e->dLR = xym1 * (1. + R*xlog);
232 e->dR2 = xlog * e->dR;
233 }
234 else if (R > 1.) {
235 e->dL = 0.;
236 goto dRzero;
237 }
238 else if (R == 1.) {
239 e->dL = 1.;
240 dRzero:
241 e->dR = e->dL2 = e->dLR = e->dR2 = 0.;
242 }
243 else
244 badpow:
245 introuble2("pow'",L,R,2 K_ASL);
246 }
247 return rv;
248 }
249 real
250 #ifdef KR_headers
251 f_OP1POW(e K_ASL) expr *e; D_ASL
252 #else
253 f_OP1POW(expr *e A_ASL)
254 #endif
255 /* f_OPPOW for R = numeric constant */
256 {
257 expr *e1;
258 real L, R, rv;
259 e1 = e->L.e;
260 L = (*e1->op)(e1 K_ASL);
261 R = ((expr_n *)e->R.e)->v;
262 rv = mypow(L,R);
263 if (errchk(rv))
264 introuble2("pow",L,R,1 K_ASL);
265 if (want_deriv) {
266 if (L) {
267 e->dL = R * (rv/L);
268 e->dL2 = (R-1.)*(e->dL/L);
269 }
270 else if (R > 1.)
271 e->dL = e->dL2 = 0.;
272 else
273 introuble2("pow'",L,R,2 K_ASL);
274 }
275 return rv;
276 }
277 real
278 #ifdef KR_headers
279 f_OP2POW(e K_ASL) expr *e; D_ASL
280 #else
281 f_OP2POW(expr *e A_ASL)
282 #endif
283 /* f_OPPOW for R = 2 */
284 {
285 expr *e1;
286 real L;
287 e1 = e->L.e;
288 L = (*e1->op)(e1 K_ASL);
289 e->dL = L + L;
290 return L*L;
291 }
292
293 real
294 #ifdef KR_headers
295 f_OPCPOW(e K_ASL) expr *e; D_ASL
296 #else
297 f_OPCPOW(expr *e A_ASL)
298 #endif
299 /* f_OPPOW for L = numeric constant */
300 {
301 expr *e1;
302 real L, R, rv;
303 e1 = e->R.e;
304 rv = mypow(L = e->L.en->v, R = (*e1->op)(e1 K_ASL));
305 if (errchk(rv))
306 introuble2("pow",L,R,1 K_ASL);
307 if (want_deriv) {
308 if (L > 0.) {
309 if (e->dL == 1)
310 e->dL = log(L); /* cache value */
311 e->dR = e->dL * rv;
312 e->dR2 = e->dR * e->dL;
313 }
314 else if (L == 0 && R >= 1.)
315 e->dR = e->dR2 = 0.;
316 else
317 introuble2("pow'",L,R,2 K_ASL);
318 }
319 return rv;
320 }
321
322 static real
323 #ifdef KR_headers
324 f_OPLESS(e K_ASL) expr *e; D_ASL
325 #else
326 f_OPLESS(expr *e A_ASL)
327 #endif
328 {
329 expr *e1;
330 real L;
331 e1 = e->L.e;
332 L = (*e1->op)(e1 K_ASL);
333 e1 = e->R.e;
334 L -= (*e1->op)(e1 K_ASL);
335 if (L < 0.)
336 return e->dL = e->dR = 0;
337 e->dL = 1.;
338 e->dR = -1.;
339 return L;
340 }
341
342 static real
343 #ifdef KR_headers
344 f_MINLIST(e0 K_ASL) expr *e0; D_ASL
345 #else
346 f_MINLIST(expr *e0 A_ASL)
347 #endif
348 {
349 de *d, *d1;
350 real t, rv;
351 expr *e1, *e2;
352 derp *D;
353 expr_va *e = (expr_va *)e0;
354
355 d = e->L.d;
356 e1 = e2 = d->e;
357 rv = (*e1->op)(e1 K_ASL);
358 for(d1 = d++; e1 = d->e; d++) {
359 t = (*e1->op)(e1 K_ASL);
360 if (rv > t) {
361 rv = t;
362 d1 = d;
363 e2 = e1;
364 }
365 }
366 if (D = e->R.D) {
367 D->a.rp = d1->dv.rp;
368 D->next = d1->d;
369 }
370 e->val = e2;
371 e->vale = d1->ee;
372 e->valf = d1->ef;
373 return rv;
374 }
375
376 static real
377 #ifdef KR_headers
378 f_MAXLIST(e0 K_ASL) expr *e0; D_ASL
379 #else
380 f_MAXLIST(expr *e0 A_ASL)
381 #endif
382 {
383 de *d, *d1;
384 real t, rv;
385 expr *e1, *e2;
386 derp *D;
387 expr_va *e = (expr_va *)e0;
388
389 d = e->L.d;
390 e1 = e2 = d->e;
391 rv = (*e1->op)(e1 K_ASL);
392 for(d1 = d++; e1 = d->e; d++) {
393 t = (*e1->op)(e1 K_ASL);
394 if (rv < t) {
395 rv = t;
396 d1 = d;
397 }
398 }
399 if (D = e->R.D) {
400 D->a.rp = d1->dv.rp;
401 D->next = d1->d;
402 }
403 e->val = e2;
404 e->vale = d1->ee;
405 e->valf = d1->ef;
406 return rv;
407 }
408
409 static real
410 #ifdef KR_headers
411 f_FLOOR(e K_ASL) expr *e; D_ASL
412 #else
413 f_FLOOR(expr *e A_ASL)
414 #endif
415 {
416 /* e->dL = 0.; */
417 e = e->L.e;
418 return floor((*e->op)(e K_ASL));
419 }
420
421 static real
422 #ifdef KR_headers
423 f_CEIL(e K_ASL) expr *e; D_ASL
424 #else
425 f_CEIL(expr *e A_ASL)
426 #endif
427 {
428 /* e->dL = 0.; */
429 e = e->L.e;
430 return ceil((*e->op)(e K_ASL));
431 }
432
433 static real
434 #ifdef KR_headers
435 f_ABS(e K_ASL) expr *e; D_ASL
436 #else
437 f_ABS(expr *e A_ASL)
438 #endif
439 {
440 real rv;
441 expr *e1;
442 e1 = e->L.e;
443 rv = (*e1->op)(e1 K_ASL);
444 if (rv < 0.) {
445 e->dL = -1.;
446 return -rv;
447 }
448 e->dL = 1.;
449 return rv;
450 }
451
452 static real
453 #ifdef KR_headers
454 f_OPUMINUS(e K_ASL) expr *e; D_ASL
455 #else
456 f_OPUMINUS(expr *e A_ASL)
457 #endif
458 {
459 /* e->dL = -1.; */
460 e = e->L.e;
461 return -(*e->op)(e K_ASL);
462 }
463
464 static real
465 #ifdef KR_headers
466 f_OP_tanh(e K_ASL) expr *e; D_ASL
467 #else
468 f_OP_tanh(expr *e A_ASL)
469 #endif
470 {
471 real rv, t, t1;
472 expr *e1;
473 e1 = e->L.e;
474 rv = tanh(t = (*e1->op)(e1 K_ASL));
475 if (errchk(rv))
476 introuble("tanh",t,1 K_ASL);
477 if (want_deriv) {
478 t1 = cosh(t);
479 if (errchk(t1))
480 introuble("tanh'", t, 2 K_ASL);
481 t1 = 1. / t1;
482 e->dL2 = -(rv+rv) * (e->dL = t1*t1);
483 }
484 return rv;
485 }
486
487 static real
488 #ifdef KR_headers
489 f_OP_tan(e K_ASL) expr *e; D_ASL
490 #else
491 f_OP_tan(expr *e A_ASL)
492 #endif
493 {
494 real rv, t, t1;
495 expr *e1;
496 e1 = e->L.e;
497 rv = tan(t = (*e1->op)(e1 K_ASL));
498 if (errchk(rv))
499 introuble("tan",t,1 K_ASL);
500 if (want_deriv) {
501 t1 = cos(t);
502 if (errchk(t1) || !t1)
503 introuble("tan'",t,2 K_ASL);
504 t1 = 1. / t1;
505 e->dL2 = (rv+rv) * (e->dL = t1*t1);
506 }
507 return rv;
508 }
509
510 static real
511 #ifdef KR_headers
512 f_OP_sqrt(e K_ASL) expr *e; D_ASL
513 #else
514 f_OP_sqrt(expr *e A_ASL)
515 #endif
516 {
517 real t, rv;
518 expr *e1;
519 e1 = e->L.e;
520 t = (*e1->op)(e1 K_ASL);
521 if (t < 0. || (rv = sqrt(t), errchk(rv)))
522 introuble("sqrt",t,1 K_ASL);
523 if (want_deriv) {
524 if (rv <= 0.)
525 introuble("sqrt'",t,2 K_ASL);
526 e->dL = 0.5 / rv;
527 e->dL2 = -0.5 * e->dL / t;
528 }
529 return rv;
530 }
531
532 static real
533 #ifdef KR_headers
534 f_OP_sinh(e K_ASL) expr *e; D_ASL
535 #else
536 f_OP_sinh(expr *e A_ASL)
537 #endif
538 {
539 real t, rv;
540 expr *e1;
541 e1 = e->L.e;
542 rv = sinh(t = (*e1->op)(e1 K_ASL));
543 if (errchk(rv))
544 introuble("sinh",t,1 K_ASL);
545 if (want_deriv) {
546 e->dL = cosh(t);
547 if (errchk(e->dL))
548 introuble("sinh'",t,2 K_ASL);
549 e->dL2 = rv;
550 }
551 return rv;
552 }
553
554 static real
555 #ifdef KR_headers
556 f_OP_sin(e K_ASL) expr *e; D_ASL
557 #else
558 f_OP_sin(expr *e A_ASL)
559 #endif
560 {
561 real t, rv;
562 expr *e1;
563 e1 = e->L.e;
564 rv = sin(t = (*e1->op)(e1 K_ASL));
565 if (errchk(rv))
566 introuble("sin",t,1 K_ASL);
567 if (want_deriv) {
568 e->dL = cos(t);
569 if (errchk(e->dL))
570 introuble("sin'",t,2 K_ASL);
571 e->dL2 = -rv;
572 }
573 return rv;
574 }
575
576 static real
577 #ifdef KR_headers
578 f_OP_log10(e K_ASL) expr *e; D_ASL
579 #else
580 f_OP_log10(expr *e A_ASL)
581 #endif
582 {
583 real t, rv;
584 expr *e1;
585 static real Le10;
586
587 e1 = e->L.e;
588 rv = log10(t = (*e1->op)(e1 K_ASL));
589 if (errchk(rv))
590 introuble("log10",t,1 K_ASL);
591 if (want_deriv) {
592 if (!Le10)
593 Le10 = 1. / log(10.);
594 e->dL = Le10 / t;
595 e->dL2 = -e->dL / t;
596 }
597 return rv;
598 }
599
600 static real
601 #ifdef KR_headers
602 f_OP_log(e K_ASL) expr *e; D_ASL
603 #else
604 f_OP_log(expr *e A_ASL)
605 #endif
606 {
607 real t, rv;
608 expr *e1;
609 e1 = e->L.e;
610 rv = log(t = (*e1->op)(e1 K_ASL));
611 if (errchk(rv))
612 introuble("log",t,1 K_ASL);
613 if (want_deriv) {
614 t = e->dL = 1. / t;
615 e->dL2 = -t * t;
616 }
617 return rv;
618 }
619
620 static real
621 #ifdef KR_headers
622 f_OP_exp(e K_ASL) expr *e; D_ASL
623 #else
624 f_OP_exp(expr *e A_ASL)
625 #endif
626 {
627 real t, rv;
628 expr *e1;
629 e1 = e->L.e;
630 rv = e->dL = e->dL2 = exp(t = (*e1->op)(e1 K_ASL));
631 if (errchk(rv))
632 if (t >= 0.)
633 introuble("exp",t,1 K_ASL);
634 else {
635 errno_set(0);
636 rv = 0.;
637 }
638 return rv;
639 }
640
641 static real
642 #ifdef KR_headers
643 f_OP_cosh(e K_ASL) expr *e; D_ASL
644 #else
645 f_OP_cosh(expr *e A_ASL)
646 #endif
647 {
648 real t, rv;
649 expr *e1;
650 e1 = e->L.e;
651 rv = cosh(t = (*e1->op)(e1 K_ASL));
652 if (errchk(rv))
653 introuble("cosh",t,1 K_ASL);
654 if (want_deriv) {
655 e->dL = sinh(t);
656 if (errchk(e->dL))
657 introuble("cosh'",t,2 K_ASL);
658 e->dL2 = rv;
659 }
660 return rv;
661 }
662
663 static real
664 #ifdef KR_headers
665 f_OP_cos(e K_ASL) expr *e; D_ASL
666 #else
667 f_OP_cos(expr *e A_ASL)
668 #endif
669 {
670 real t, rv;
671 expr *e1;
672 e1 = e->L.e;
673 rv = cos(t = (*e1->op)(e1 K_ASL));
674 if (errchk(rv))
675 introuble("cos",t,1 K_ASL);
676 if (want_deriv) {
677 e->dL = -sin(t);
678 if (errchk(e->dL))
679 introuble("cos'",t,2 K_ASL);
680 e->dL2 = -rv;
681 }
682 return rv;
683 }
684
685 static real
686 #ifdef KR_headers
687 f_OP_atanh(e K_ASL) expr *e; D_ASL
688 #else
689 f_OP_atanh(expr *e A_ASL)
690 #endif
691 {
692 real rv, t, t1;
693 expr *e1;
694
695 e1 = e->L.e;
696 t = (*e1->op)(e1 K_ASL);
697 if (t <= -1. || t >= 1.) {
698 errno_set(EDOM);
699 bad_atanh:
700 rv = 0.;
701 introuble("atanh",t,1 K_ASL);
702 }
703 else {
704 rv = 0.5*log((1. + t) / (1. - t));
705 if (errchk(rv))
706 goto bad_atanh;
707 }
708 if (want_deriv) {
709 e->dL = t1 = 1. / (1. - t*t);
710 e->dL2 = (t+t)*t1*t1;
711 }
712 return rv;
713 }
714
715 static real
716 #ifdef KR_headers
717 f_OP_atan2(e K_ASL) expr *e; D_ASL
718 #else
719 f_OP_atan2(expr *e A_ASL)
720 #endif
721 {
722 real L, R, rv, t, t1;
723 expr *e1;
724
725 e1 = e->L.e;
726 L = (*e1->op)(e1 K_ASL);
727 e1 = e->R.e;
728 R = (*e1->op)(e1 K_ASL);
729 rv = atan2(L,R);
730 if (errchk(rv))
731 introuble2("atan2",L,R,1 K_ASL);
732 if (want_deriv) {
733 /*
734 t = L / R;
735 t1 = 1. / (1. + t*t);
736 e->dL = t1 /= R;
737 e->dR = -t*t1;
738 */
739
740 t = 1. / (L*L + R*R);
741 e->dL = t * R;
742 e->dR = -t * L;
743 t *= t;
744 t1 = L*R;
745 e->dL2 = -(e->dR2 = t * (t1+t1));
746 e->dLR = t * (L*L - R*R);
747 }
748 return rv;
749 }
750
751 static real
752 #ifdef KR_headers
753 f_OP_atan(e K_ASL) expr *e; D_ASL
754 #else
755 f_OP_atan(expr *e A_ASL)
756 #endif
757 {
758 real rv, t, t1;
759 expr *e1;
760
761 e1 = e->L.e;
762 rv = atan(t = (*e1->op)(e1 K_ASL));
763 if (errchk(rv))
764 introuble("atan",t,1 K_ASL);
765 if (want_deriv) {
766 e->dL = t1 = 1. / (1. + t*t);
767 e->dL2 = -(t+t)*t1*t1;
768 }
769 return rv;
770 }
771
772 static real
773 #ifdef KR_headers
774 f_OP_asinh(e K_ASL) expr *e; D_ASL
775 #else
776 f_OP_asinh(expr *e A_ASL)
777 #endif
778 {
779 expr *e1;
780 real rv, s, t, t1, t2;
781
782 e1 = e->L.e;
783 t = (*e1->op)(e1 K_ASL);
784 s = 1.;
785 if (t < 0.)
786 s = -1.;
787 rv = log(s*t + (t1 = sqrt(t2 = t*t + 1.)));
788 if (errchk(rv))
789 introuble("asinh",t,1 K_ASL);
790 if (want_deriv)
791 e->dL2 = -(t/t2) * (e->dL = 1. / t1);
792 return s*rv;
793 }
794
795 static real
796 #ifdef KR_headers
797 f_OP_asin(e K_ASL) expr *e; D_ASL
798 #else
799 f_OP_asin(expr *e A_ASL)
800 #endif
801 {
802 expr *e1;
803 real rv, t, t1;
804
805 e1 = e->L.e;
806 rv = asin(t = (*e1->op)(e1 K_ASL));
807 if (errchk(rv))
808 introuble("asin",t,1 K_ASL);
809 if (want_deriv) {
810 if ((t1 = 1. - t*t) <= 0.)
811 introuble("asin'",t,2 K_ASL);
812 e->dL2 = t * (e->dL = 1. / sqrt(t1)) / t1;
813 }
814 return rv;
815 }
816
817 static real
818 #ifdef KR_headers
819 f_OP_acosh(e K_ASL) expr *e; D_ASL
820 #else
821 f_OP_acosh(expr *e A_ASL)
822 #endif
823 {
824 expr *e1;
825 real rv, t, t1, t2;
826
827 e1 = e->L.e;
828 t = (*e1->op)(e1 K_ASL);
829 if (t < 1.) {
830 errno_set(EDOM);
831 bad_acosh:
832 rv = 0.;
833 introuble("acosh",t,1 K_ASL);
834 }
835 else {
836 rv = log(t + (t1 = sqrt(t2 = t*t - 1.)));
837 if (errchk(rv))
838 goto bad_acosh;
839 }
840 if (want_deriv) {
841 if (t2 <= 0.)
842 introuble("acosh'",t,1 K_ASL);
843 e->dL2 = -t * (e->dL = 1. / t1) / t2;
844 }
845 return rv;
846 }
847
848 static real
849 #ifdef KR_headers
850 f_OP_acos(e K_ASL) expr *e; D_ASL
851 #else
852 f_OP_acos(expr *e A_ASL)
853 #endif
854 {
855 expr *e1;
856 real rv, t, t1;
857
858 e1 = e->L.e;
859 rv = acos(t = (*e1->op)(e1 K_ASL));
860 if (errchk(rv))
861 introuble("acos",t,1 K_ASL);
862 if (want_deriv) {
863 if ((t1 = 1. - t*t) <= 0.)
864 introuble("acos'",t,2 K_ASL);
865 e->dL2 = t * (e->dL = -1. / sqrt(1. - t*t)) / t1;
866 }
867 return rv;
868 }
869
870 static real
871 #ifdef KR_headers
872 f_OPIFnl(e K_ASL) expr *e; D_ASL
873 #else
874 f_OPIFnl(expr *e A_ASL)
875 #endif
876 {
877 expr_if *eif = (expr_if *)e;
878 derp *D;
879
880 e = eif->e;
881 if ((*e->op)(e K_ASL)) {
882 eif->val = e = eif->T;
883 eif->vale = eif->Te;
884 eif->valf = eif->Tf;
885 if (D = eif->D) {
886 D->a.rp = eif->Tv.rp;
887 D->next = eif->dT;
888 }
889 }
890 else {
891 eif->val = e = eif->F;
892 eif->vale = eif->Fe;
893 eif->valf = eif->Ff;
894 if (D = eif->D) {
895 D->a.rp = eif->Fv.rp;
896 D->next = eif->dF;
897 }
898 }
899 return (*e->op)(e K_ASL);
900 }
901
902 static real
903 #ifdef KR_headers
904 f_OPOR(e K_ASL) expr *e; D_ASL
905 #else
906 f_OPOR(expr *e A_ASL)
907 #endif
908 {
909 expr *e2;
910 e2 = e->R.e;
911 e = e->L.e;
912 return (*e->op)(e K_ASL) || (*e2->op)(e2 K_ASL) ? 1. : 0.;
913 }
914
915 static real
916 #ifdef KR_headers
917 f_OPAND(e K_ASL) expr *e; D_ASL
918 #else
919 f_OPAND(expr *e A_ASL)
920 #endif
921 {
922 expr *e2;
923 e2 = e->R.e;
924 e = e->L.e;
925 return (*e->op)(e K_ASL) && (*e2->op)(e2 K_ASL) ? 1. : 0.;
926 }
927
928 static real
929 #ifdef KR_headers
930 f_LT(e K_ASL) expr *e; D_ASL
931 #else
932 f_LT(expr *e A_ASL)
933 #endif
934 {
935 expr *e2;
936 e2 = e->R.e;
937 e = e->L.e;
938 return (*e->op)(e K_ASL) < (*e2->op)(e2 K_ASL) ? 1. : 0;
939 }
940
941 static real
942 #ifdef KR_headers
943 f_LE(e K_ASL) expr *e; D_ASL
944 #else
945 f_LE(expr *e A_ASL)
946 #endif
947 {
948 expr *e2;
949 e2 = e->R.e;
950 e = e->L.e;
951 return (*e->op)(e K_ASL) <= (*e2->op)(e2 K_ASL) ? 1. : 0;
952 }
953
954 static real
955 #ifdef KR_headers
956 f_EQ(e K_ASL) expr *e; D_ASL
957 #else
958 f_EQ(expr *e A_ASL)
959 #endif
960 {
961 expr *e2;
962 e2 = e->R.e;
963 e = e->L.e;
964 return (*e->op)(e K_ASL) == (*e2->op)(e2 K_ASL) ? 1. : 0;
965 }
966
967 static real
968 #ifdef KR_headers
969 f_GE(e K_ASL) expr *e; D_ASL
970 #else
971 f_GE(expr *e A_ASL)
972 #endif
973 {
974 expr *e2;
975 e2 = e->R.e;
976 e = e->L.e;
977 return (*e->op)(e K_ASL) >= (*e2->op)(e2 K_ASL) ? 1. : 0;
978 }
979
980 static real
981 #ifdef KR_headers
982 f_GT(e K_ASL) expr *e; D_ASL
983 #else
984 f_GT(expr *e A_ASL)
985 #endif
986 {
987 expr *e2;
988 e2 = e->R.e;
989 e = e->L.e;
990 return (*e->op)(e K_ASL) > (*e2->op)(e2 K_ASL) ? 1. : 0;
991 }
992
993 static real
994 #ifdef KR_headers
995 f_NE(e K_ASL) expr *e; D_ASL
996 #else
997 f_NE(expr *e A_ASL)
998 #endif
999 {
1000 expr *e2;
1001 e2 = e->R.e;
1002 e = e->L.e;
1003 return (*e->op)(e K_ASL) != (*e2->op)(e2 K_ASL) ? 1. : 0;
1004 }
1005
1006 static real
1007 #ifdef KR_headers
1008 f_OPNOT(e K_ASL) expr *e; D_ASL
1009 #else
1010 f_OPNOT(expr *e A_ASL)
1011 #endif
1012 {
1013 e = e->L.e;
1014 return (*e->op)(e K_ASL) ? 0. : 1.;
1015 }
1016
1017 static real
1018 #ifdef KR_headers
1019 f_ANDLIST(e K_ASL) expr *e; D_ASL
1020 #else
1021 f_ANDLIST(expr *e A_ASL)
1022 #endif
1023 {
1024 expr **ep, **epe;
1025
1026 ep = e->L.ep;
1027 epe = e->R.ep;
1028 do {
1029 e = *ep++;
1030 if ((*e->op)(e K_ASL) == 0.)
1031 return 0.;
1032 }
1033 while(ep < epe);
1034 return 1.;
1035 }
1036
1037 static real
1038 #ifdef KR_headers
1039 f_ORLIST(e K_ASL) expr *e; D_ASL
1040 #else
1041 f_ORLIST(expr *e A_ASL)
1042 #endif
1043 {
1044 expr **ep, **epe;
1045
1046 ep = e->L.ep;
1047 epe = e->R.ep;
1048 do {
1049 e = *ep++;
1050 if ((*e->op)(e K_ASL) != 0.)
1051 return 1.;
1052 }
1053 while(ep < epe);
1054 return 0.;
1055 }
1056
1057 #define f_OPIMPELSE f_OPIMPELSE_rops2_ASL
1058 real
1059 #ifdef KR_headers
1060 f_OPIMPELSE(e K_ASL) expr *e; D_ASL
1061 #else
1062 f_OPIMPELSE(expr *e A_ASL)
1063 #endif
1064 {
1065 expr_if *eif = (expr_if *)e;
1066
1067 e = eif->e;
1068 e = ((*e->op)(e K_ASL) != 0.) ? eif->T : eif->F;
1069 return (*e->op)(e K_ASL);
1070 }
1071
1072 static real
1073 #ifdef KR_headers
1074 f_OP_IFF(e K_ASL) expr *e; D_ASL
1075 #else
1076 f_OP_IFF(expr *e A_ASL)
1077 #endif
1078 {
1079 expr *e1;
1080 int a, b;
1081
1082 e1 = e->L.e;
1083 a = (*e1->op)(e1 K_ASL) != 0.;
1084 e1 = e->R.e;
1085 b = (*e1->op)(e1 K_ASL) != 0.;
1086 return a == b ? 1. : 0.;
1087 }
1088
1089 real
1090 #ifdef KR_headers
1091 f_OPSUMLIST(e K_ASL) expr *e; D_ASL
1092 #else
1093 f_OPSUMLIST(expr *e A_ASL)
1094 #endif
1095 {
1096 expr **ep, **epe;
1097 real x;
1098 ep = e->L.ep;
1099 epe = e->R.ep;
1100 e = *ep;
1101 x = (*e->op)(e K_ASL);
1102 while(++ep < epe) {
1103 e = *ep;
1104 x += (*e->op)(e K_ASL);
1105 }
1106 return x;
1107 }
1108
1109 real
1110 #ifdef KR_headers
1111 f_OPPLTERM(e K_ASL) expr *e; D_ASL
1112 #else
1113 f_OPPLTERM(expr *e A_ASL)
1114 #endif
1115 {
1116 plterm *p = e->L.p;
1117 real r, t;
1118 int n = p->n;
1119 real *bs = p->bs;
1120 MTNot_Used(asl);
1121
1122 r = ((expr_v *)e->R.e)->v;
1123 if (r >= 0) {
1124 /* should use binary search for large n */
1125 while(bs[1] <= 0.) {
1126 bs += 2;
1127 if (--n <= 1)
1128 return r*(e->dL = *bs);
1129 }
1130 if (r <= bs[1])
1131 return r*(e->dL = *bs);
1132 for(t = bs[0]*bs[1]; --n > 1 && r > bs[3]; bs += 2)
1133 t += (bs[3]-bs[1])*bs[2];
1134 return t + (r-bs[1])*(e->dL = bs[2]);
1135 }
1136 bs += 2*(n-1);
1137 while(bs[-1] >= 0.) {
1138 bs -= 2;
1139 if (--n <= 1)
1140 return r*(e->dL = *bs);
1141 }
1142 if (r >= bs[-1])
1143 return r*(e->dL = *bs);
1144 for(t = bs[0]*bs[-1]; --n > 1 && r < bs[-3]; bs -= 2)
1145 t += (bs[-3]-bs[-1])*bs[-2];
1146 return t + (r-bs[-1])*(e->dL = bs[-2]);
1147 }
1148
1149 char *
1150 #ifdef KR_headers
1151 f_OPHOL(e K_ASL) expr *e; D_ASL
1152 #else
1153 f_OPHOL(expr *e A_ASL)
1154 #endif
1155 {
1156 MTNot_Used(asl);
1157 return ((expr_h *)e)->sym;
1158 }
1159
1160 real
1161 #ifdef KR_headers
1162 f_OPVARVAL(e K_ASL) expr *e; D_ASL
1163 #else
1164 f_OPVARVAL(expr *e A_ASL)
1165 #endif
1166 {
1167 MTNot_Used(asl);
1168 return ((expr_v *)e)->v;
1169 }
1170
1171 char *
1172 #ifdef KR_headers
1173 f_OPIFSYM(e K_ASL) expr *e; D_ASL
1174 #else
1175 f_OPIFSYM(expr *e A_ASL)
1176 #endif
1177 {
1178 expr_if *eif = (expr_if *)e;
1179 e = eif->e;
1180 e = (*e->op)(e K_ASL) ? eif->T : eif->F;
1181 return (*(sfunc*)e->op)(e K_ASL);
1182 }
1183
1184 struct
1185 TMInfo {
1186 union {
1187 TMInfo *prev;
1188 double align;
1189 } u;
1190 };
1191
1192 real
1193 #ifdef KR_headers
1194 f_OPFUNCALL(e K_ASL) expr *e; D_ASL
1195 #else
1196 f_OPFUNCALL(expr *e A_ASL)
1197 #endif
1198 {
1199 TMInfo T, *T1, *T1prev;
1200 arglist *al;
1201 argpair *ap, *ape;
1202 char *s;
1203 expr_f *f = (expr_f *)e;
1204 func_info *fi = f->fi;
1205 int jv;
1206 real rv;
1207
1208 for(ap = f->ap, ape = f->ape; ap < ape; ap++) {
1209 e = ap->e;
1210 *ap->u.v = (*e->op)(e K_ASL);
1211 }
1212 for(ap = f->sap, ape = f->sape; ap < ape; ap++) {
1213 e = ap->e;
1214 *ap->u.s = (*(sfunc*)e->op)(e K_ASL);
1215 }
1216 T.u.prev = 0;
1217 al = f->al;
1218 al->TMI = &T;
1219 al->Errmsg = 0;
1220 rv = (*fi->funcp)(al);
1221 errno_set(0);
1222 if ((s = al->Errmsg) && !err_jmp) {
1223 report_where(asl);
1224 jv = 1;
1225 if (*s == '\'') {
1226 ++jv;
1227 ++s;
1228 }
1229 fprintf(Stderr, "Error in function %s:\n\t%s\n", fi->name, s);
1230 fflush(Stderr);
1231 }
1232 for(T1 = T.u.prev; T1; T1 = T1prev) {
1233 T1prev = T1->u.prev;
1234 free(T1);
1235 }
1236 if (s) {
1237 jmp_check(err_jmp, jv);
1238 jmp_check(err_jmp1,jv);
1239 exit(1);
1240 }
1241 return rv;
1242 }
1243
1244 static real
1245 #ifdef KR_headers
1246 f_OPintDIV(e K_ASL) expr *e; D_ASL
1247 #else
1248 f_OPintDIV(expr *e A_ASL)
1249 #endif
1250 {
1251 expr *e1;
1252 real L, R;
1253
1254 e1 = e->L.e;
1255 L = (*e1->op)(e1 K_ASL);
1256 e1 = e->R.e;
1257 if (!(R = (*e1->op)(e1 K_ASL)))
1258 zero_div(L, " div " K_ASL);
1259 return (L /= R) >= 0 ? floor(L) : ceil(L);
1260 }
1261
1262 static real
1263 #ifdef KR_headers
1264 f_OPprecision(e K_ASL) expr *e; D_ASL
1265 #else
1266 f_OPprecision(expr *e A_ASL)
1267 #endif
1268 {
1269 expr *e1;
1270 real L, R;
1271 char buf[32];
1272
1273 e1 = e->L.e;
1274 L = (*e1->op)(e1 K_ASL);
1275 e1 = e->R.e;
1276 R = (*e1->op)(e1 K_ASL);
1277 g_fmtp(buf, L, (int)R);
1278 return strtod(buf, (char **)0);
1279 }
1280
1281 static real
1282 #ifdef KR_headers
Round(x,prec)1283 Round(x, prec) real x; int prec;
1284 #else
1285 Round(real x, int prec)
1286 #endif
1287 {
1288 char *b, *s, *s0, *se;
1289 int decpt, L, sign;
1290 char buf[96];
1291
1292 if (!x)
1293 return x;
1294 s = dtoa(x, 3, prec, &decpt, &sign, &se);
1295 if (decpt == 9999) {
1296 zreturn:
1297 freedtoa(s);
1298 return x;
1299 }
1300 L = se - s;
1301 if (L <= 0) {
1302 x = 0.;
1303 goto zreturn;
1304 }
1305 if (L > 80)
1306 se = s + 80;
1307 s0 = s;
1308 b = buf;
1309 if (sign)
1310 *b++ = '-';
1311 *b++ = '.';
1312 while(s < se)
1313 *b++ = *s++;
1314 *b = 0;
1315 freedtoa(s0);
1316 if (decpt)
1317 snprintf(b, buf + sizeof(buf) - b, "e%d", decpt);
1318 return strtod(buf, (char **)0);
1319 }
1320
1321 static real
1322 #ifdef KR_headers
1323 f_OPround(e K_ASL) expr *e; D_ASL
1324 #else
1325 f_OPround(expr *e A_ASL)
1326 #endif
1327 {
1328 expr *e1;
1329 real L, R;
1330
1331 e1 = e->L.e;
1332 L = (*e1->op)(e1 K_ASL);
1333 e1 = e->R.e;
1334 R = (*e1->op)(e1 K_ASL);
1335 return Round(L, (int)R);
1336 }
1337
1338 static real
1339 #ifdef KR_headers
1340 f_OPtrunc(e K_ASL) expr *e; D_ASL
1341 #else
1342 f_OPtrunc(expr *e A_ASL)
1343 #endif
1344 {
1345 expr *e1;
1346 real L, R;
1347 int prec;
1348
1349 e1 = e->L.e;
1350 L = (*e1->op)(e1 K_ASL);
1351 e1 = e->R.e;
1352 if (!(R = (*e1->op)(e1 K_ASL)))
1353 return L >= 0. ? floor(L) : ceil(L);
1354 R = Round(L, prec = (int)R);
1355 if (R != L) {
1356 R = 0.5*mypow(10., (real)-prec);
1357 R = Round(L > 0 ? L - R : L + R, prec);
1358 }
1359 return R;
1360 }
1361
1362
1363 static real
1364 #ifdef KR_headers
1365 f_OPCOUNT(e K_ASL) expr *e; D_ASL
1366 #else
1367 f_OPCOUNT(expr *e A_ASL)
1368 #endif
1369 {
1370 expr **ep, **epe;
1371 real x;
1372 ep = e->L.ep;
1373 epe = e->R.ep;
1374 e = *ep++;
1375 if (x = (*e->op)(e K_ASL))
1376 x = 1.;
1377 do {
1378 e = *ep++;
1379 if ((*e->op)(e K_ASL))
1380 x++;
1381 }
1382 while(ep < epe);
1383 return x;
1384 }
1385
1386 typedef struct
1387 jb_st {
1388 jmp_buf jb;
1389 } jb_st;
1390
1391 static int
1392 #ifdef KR_headers
rcompj(a,b,v)1393 rcompj(a, b, v) char *a, *b, *v;
1394 #else
1395 rcompj(const void *a, const void *b, void *v)
1396 #endif
1397 {
1398 jb_st *J;
1399 real t = *(real *)a - *(real *)b;
1400
1401 if (!t) {
1402 J = (jb_st*)v;
1403 longjmp(J->jb, 1);
1404 }
1405 return t < 0 ? -1 : 1;
1406 }
1407
1408 static real
1409 #ifdef KR_headers
1410 f_OPALLDIFF(e K_ASL) expr *e; D_ASL
1411 #else
1412 f_OPALLDIFF(expr *e A_ASL)
1413 #endif
1414 {
1415 int n;
1416 expr **ep, **epe;
1417 jb_st J;
1418 real *r, *r1, rbuf[128], t;
1419
1420 ep = e->L.ep;
1421 epe = e->R.ep;
1422 n = epe - ep;
1423 r = rbuf;
1424 if (n > sizeof(rbuf)/sizeof(real))
1425 r = (real*)Malloc(n*sizeof(real));
1426 r1 = r;
1427 while(ep < epe) {
1428 e = *ep++;
1429 *r1++ = (*e->op)(e K_ASL);
1430 }
1431 t = 1.;
1432 if (setjmp(J.jb)) {
1433 t = 0.;
1434 goto done;
1435 }
1436 qsortv(r, n, sizeof(real), rcompj, &J);
1437 done:
1438 if (r != rbuf)
1439 free(r);
1440 return t;
1441 }
1442
1443 static real
1444 #ifdef KR_headers
1445 f_OPNUMBEROF(e K_ASL) expr *e; D_ASL
1446 #else
1447 f_OPNUMBEROF(expr *e A_ASL)
1448 #endif
1449 {
1450 expr **ep, **epe;
1451 real n, t;
1452
1453 ep = e->L.ep;
1454 epe = e->R.ep;
1455 n = 0.;
1456 e = *ep++;
1457 t = (*e->op)(e K_ASL);
1458 while(ep < epe) {
1459 e = *ep++;
1460 if (t == (*e->op)(e K_ASL))
1461 n++;
1462 }
1463 return n;
1464 }
1465
1466 #define f_OPATLEAST f_LE
1467 #define f_OPATMOST f_GE
1468 #define f_OPEXACTLY f_EQ
1469 #define f_OPNOTATLEAST f_GT
1470 #define f_OPNOTATMOST f_LT
1471 #define f_OPNOTEXACTLY f_NE
1472 #define efunc efunc2
1473 #define f_OPNUMBEROFs f_OPNUMBEROF /*TEMPORARY*/
1474
1475 efunc *
1476 r_ops[] = {
1477 #include "r_op.hd"
1478 };
1479
1480 #ifdef __cplusplus
1481 }
1482 #endif
1483