/[ascend]/trunk/ascend/compiler/func.c
ViewVC logotype

Contents of /trunk/ascend/compiler/func.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2497 - (show annotations) (download) (as text)
Wed Aug 17 09:34:25 2011 UTC (10 years, 10 months ago) by jpye
File MIME type: text/x-csrc
File size: 22742 byte(s)
Fixes bug 509 (thanks again, Ujjaval!)
1 /*
2 * Function Module
3 * by Tom Epperly
4 * Created: 8/11/1990
5 * Version: $Revision: 1.18 $
6 * Version control file: $RCSfile: func.c,v $
7 * Date last modified: $Date: 2001/01/31 22:23:53 $
8 * Last modified by: $Author: ballan $
9 *
10 * This file is part of the Ascend Language Interpreter.
11 *
12 * Copyright (C) 1990, 1993, 1994 Thomas Guthrie Epperly
13 *
14 * The Ascend Language Interpreter is free software; you can redistribute
15 * it and/or modify it under the terms of the GNU General Public License as
16 * published by the Free Software Foundation; either version 2 of the
17 * License, or (at your option) any later version.
18 *
19 * The Ascend Language Interpreter is distributed in hope that it will be
20 * useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
21 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
22 * General Public License for more details.
23 *
24 * You should have received a copy of the GNU General Public License
25 * along with the program; if not, write to the Free Software Foundation,
26 * Inc., 675 Mass Ave, Cambridge, MA 02139 USA. Check the file named
27 * COPYING.
28 *
29 */
30
31 #include "func.h"
32
33 #include<math.h>
34 #include "safe.h"
35
36 #ifndef M_PI
37 #define M_PI F_PI
38 #endif
39 #ifndef M_LOG10E
40 #define M_LOG10E F_LOG10_COEF
41 #endif
42
43 #ifndef NULL
44 #define NULL 0L
45 #endif
46
47 double g_lnm_epsilon = 1.0e-8;
48
49
50 #if defined(__WIN32__) || defined(_WIN32) || defined(WIN32)
51 double cbrt(register double d)
52 {
53 return pow(d,(double)0.3333333333333333333333);
54 }
55 #endif
56
57 int ascnintF(register double d)
58 {
59 return ((d)>=0.0 ? (int)floor((d) + 0.5) : -(int)floor(0.5 - (d)));
60 }
61
62
63 double dln(register double d)
64 {
65 return 1.0/d;
66 }
67
68 double dln2(register double d)
69 {
70 return -1.0/(d*d);
71 }
72
73
74 double lnm(register double d)
75 {
76 return (d>g_lnm_epsilon?log(d):d/g_lnm_epsilon + (log(g_lnm_epsilon) -1));
77
78 }
79
80 double dlnm(register double d)
81 {
82 return ( d>g_lnm_epsilon ? (double)1.0/d : 1/g_lnm_epsilon);
83 }
84
85 double dlnm2(register double d)
86 {
87 return (d>g_lnm_epsilon ? (double)-1.0/(d*d) : (double)0.0);
88 }
89
90 double dlog10(register double d)
91 {
92 return M_LOG10E/d;
93 }
94
95 double dlog102(register double d)
96 {
97 return -M_LOG10E/(d*d);
98 }
99
100 double dcos(register double d)
101 {
102 return -sin(d);
103 }
104
105 double dcos2(register double d)
106 {
107 return -cos(d);
108 }
109
110 double dtan(register double d)
111 {
112 register double t;
113 t=cos(d);
114 return 1.0/(t*t);
115 }
116
117 double dtan2(register double d)
118 {
119 register double t;
120 t=cos(d);
121 return ldexp(tan(d)/(t*t),1);
122 }
123
124 double sqr(register double d)
125 {
126 return d*d;
127 }
128
129 double dsqr(register double d)
130 {
131 return ldexp(d,1);
132 }
133
134 double dsqr2(register double d)
135 {
136 (void)d; /* stop gcc whine about unused parameter */
137 return 2.0;
138 }
139
140 double hold(double d)
141 {
142 return d;
143 }
144
145 double dsqrt(register double d)
146 {
147 return 1.0/(ldexp(sqrt(d),1));
148 }
149
150 double dsqrt2(register double d)
151 {
152 return -1.0/ldexp(sqrt(d)*d,2);
153 }
154
155 double dfabs(register double d)
156 {
157 return ((d > 0.0) ? 1.0 : ((d<0.0 ) ? -1 : 0));
158 }
159
160 double dfabs2(register double d)
161 {
162 (void)d; /* stop gcc whine about unused parameter */
163 return 0.0;
164 }
165
166 double dhold(double d)
167 {
168 (void)d; /* stop gcc whine about unused parameter */
169 return 0;
170 }
171
172 /* The next 4 are new */
173 double asc_ipow(register double d, int i) {
174 unsigned negative = 0;
175 negative = (i<0);
176 if (negative) i = (-i);
177 if (d==0 && i!=0) return 0.0;
178 switch (i) {
179 case 0: return 1.0; /* a^0 = 1, for a==0 pow is undefined. */
180 case 1: break;
181 case 2: d *= d; break;
182 case 3: d = d*d*d; break;
183 case 4: d = d*d*d*d;break;
184 case 5: d = d*d*d*d*d; break;
185 case 6: d = d*d*d*d*d*d; break;
186 case 7: d = d*d*d*d*d*d*d; break;
187 case 8: d = d*d*d*d*d*d*d*d; break;
188 case 9: d = d*d*d*d*d*d*d*d*d; break;
189 default:
190 {
191 register double res;
192 res = d;
193 for (--i; i > 0; i--) res *= d;
194 d = res;
195 }
196 break;
197 }
198 return (!negative ? d : 1.0/d);
199 }
200
201 /*
202 * Note that the following derivative functions do not
203 * set calc_ok to FALSE in the event of errors. This
204 * checking is done in the solver so we are baisicaly
205 * double checking now -> this should be fixed
206 */
207
208 double asc_d1ipow(double d, int i) {
209 if (d == 0 && i <= 1) {
210 FPRINTF(stderr,"ERROR:\t(calc) calc_ipow_D1\n");
211 FPRINTF(stderr,
212 "\t1st derivative, %g raised to %d <= 1 power undefined.\n",
213 d,i);
214 FPRINTF(stderr,"\tReturning %g.\n",0.0);
215 return(0.0);
216 }
217 return( i * asc_ipow(d,i-1));
218 }
219
220 double asc_d2ipow(double d, int i) {
221 if (d == 0 && i <= 2) {
222 FPRINTF(stderr,"ERROR:\t(calc) calc_ipow_D2\n");
223 FPRINTF(stderr,
224 "\t2nd derivative, %g raised to %d <= 2 power undefined.\n",
225 d,i);
226 FPRINTF(stderr,"\tReturning %g.\n",0.0);
227 return(0.0);
228 }
229 return( i * (i - 1) * asc_ipow(d,i-2));
230 }
231
232
233 double cube(register double d)
234 {
235 return d*d*d;
236 }
237 double dcube(register double d)
238 {
239 return 3.0*d*d;
240 }
241 double dcube2(register double d)
242 {
243 return 6.0*d;
244 }
245
246 double dcbrt(register double d)
247 {
248 register double c;
249 c=cbrt(d);
250 return (double)0.3333333333333333/(c*c);
251 }
252
253 double dcbrt2(register double d)
254 {
255 register double c;
256 c=cbrt(d);
257 return (double)-0.2222222222222222/pow(c,5.0);
258 }
259
260 double dasin(register double d)
261 {
262 return 1.0/sqrt(1.0-d*d);
263 }
264
265 double dasin2(register double d)
266 {
267 register double c;
268 c=1.0-d*d;
269 return d/(c*sqrt(c));
270 }
271
272 double dacos(register double d)
273 {
274 return -1.0/sqrt(1-d*d);
275 }
276
277 double dacos2(register double d)
278 {
279 register double c;
280 c=1.0-d*d;
281 return -d/(c*sqrt(c));
282 }
283
284 double datan(register double d)
285 {
286 return 1.0/(1.0+d*d);
287 }
288
289 double datan2(register double d)
290 {
291 return -2*d/sqr(1+d*d);
292 }
293
294 #ifdef HAVE_ERF
295 double derf(register double d)
296 {
297 return ldexp(exp(-(d*d))/sqrt(M_PI),1);
298 }
299
300 double derf2(register double d)
301 {
302 return -ldexp(d*exp(-(d*d))/sqrt(M_PI),2);
303 }
304 #endif /* HAVE_ERF */
305
306 double dtanh(register double d)
307 {
308 register double c;
309 c = cosh(d);
310 c = 1/(c*c);
311 return c;
312 }
313
314 double dtanh2(register double d)
315 {
316 register double c;
317 c = cosh(d);
318 return -ldexp(tanh(d),1)/(c*c);
319 }
320
321 double arcsinh(register double d)
322 {
323 return log(d+sqrt(d*d+1.0));
324 }
325
326 double darcsinh(register double d)
327 {
328 return 1.0/sqrt(d*d+1.0);
329 }
330
331 double darcsinh2(register double d)
332 {
333 register double c;
334 c=d*d+1.0;
335 return -d/sqrt(c*c*c);
336 }
337
338 double arccosh(register double d)
339 {
340 return log(d+sqrt(d*d-1.0));
341 }
342
343 double darccosh(register double d)
344 {
345 return 1.0/sqrt(d*d-1.0);
346 }
347
348 double darccosh2(register double d)
349 {
350 register double c;
351 c=d*d-1.0;
352 return -d/sqrt(c*c*c);
353 }
354
355 double arctanh(register double d)
356 {
357 return ldexp( log((d+1.0)/(1.0-d)) ,-1);
358 /* an alternative, more expensive but perhaps less exception prone
359 * coding of arctanh is:
360 * return log(sqrt((d+1.0)/(1.0-d)));
361 * which for d near -1 will be less likely to underflow and send 0
362 * to the log function. Until otherwise noted we are running the
363 * cheap version.
364 */
365 }
366
367 double darctanh(register double d)
368 {
369 return 1.0/(1-d*d);
370 }
371
372 double darctanh2(register double d)
373 {
374 register double c;
375 c=1.0-d*d;
376 return ldexp( d/(c*c) ,1);
377 }
378
379 #ifdef CHRIS_FUNC
380 void ExpSlope(unsigned long int nvar,
381 struct Interval *center, struct Interval *range,
382 struct Interval *slope)
383 {
384 *center = ExpInterval(*center);
385 *range = ExpInterval(*range);
386 while (nvar--){
387 *slope = MulIntervals(*range,*slope);
388 slope++;
389 }
390 }
391
392 void LnSlope(unsigned long int nvar,
393 struct Interval *center, struct Interval *range,
394 struct Interval *slope)
395 {
396 while(nvar--){
397 *slope = DivIntervals(*slope,*range);
398 slope++;
399 }
400 *center = LnInterval(*center);
401 *range = LnInterval(*range);
402 }
403
404 void LogSlope(unsigned long int nvar,
405 struct Interval *center, struct Interval *range,
406 struct Interval *slope)
407 {
408 struct Interval temp;
409 temp = MulIntervals(CreateThin(M_LN10),*range);
410 while(nvar--){
411 *slope = DivIntervals(*slope,temp);
412 slope++;
413 }
414 *center = LogInterval(*center);
415 *range = LogInterval(*range);
416 }
417
418 void SqrSlope(unsigned long int nvar,
419 struct Interval *center, struct Interval *range,
420 struct Interval *slope)
421 {
422 struct Interval temp;
423 temp = AddIntervals(*center,*range);
424 while(nvar--){
425 *slope = MulIntervals(temp,*slope);
426 slope++;
427 }
428 *center = SqrInterval(*center);
429 *range = SqrInterval(*range);
430 }
431
432 void SqrtSlope(unsigned long int nvar,
433 struct Interval *center, struct Interval *range,
434 struct Interval *slope)
435 {
436 struct Interval temp;
437 *center = SqrtInterval(*center);
438 *range = SqrtInterval(*range);
439 temp = AddIntervals(*center,*range);
440 while(nvar--){
441 *slope = DivIntervals(*slope,temp);
442 slope++;
443 }
444 }
445
446 #ifdef HAVE_ERF
447 void ErfSlope(unsigned long int nvar,
448 struct Interval *center, struct Interval *range,
449 struct Interval *slope)
450 {
451 struct Interval temp;
452 temp =
453 DivIntervals(MulIntervals(CreateThinInteger(2l),
454 ExpInterval(NegInterval(SqrInterval(*range)))),
455 SqrtInterval(CreateThin(M_PI)));
456 while(nvar--){
457 *slope =
458 MulIntervals(temp,*slope);
459 slope++;
460 }
461 *center = ErfInterval(*center);
462 *range = ErfInterval(*range);
463 }
464
465 struct Interval ErfDeriv(struct Interval i)
466 {
467 return
468 DivIntervals(MulIntervals(CreateThinInteger(2l),
469 ExpInterval(NegInterval(SqrInterval(i)))),
470 SqrtInterval(CreateThin(M_PI)));
471 }
472 #endif /* HAVE_ERF */
473
474 struct Interval LnDeriv(struct Interval i)
475 {
476 return DivIntervals(CreateThinInteger(1L),i);
477 }
478
479 struct Interval LogDeriv(struct Interval i)
480 {
481 return DivIntervals(CreateThin(M_LOG10E),i);
482 }
483
484 struct Interval SqrDeriv(struct Interval i)
485 {
486 return MulIntervals(CreateThinInteger(2L),i);
487 }
488
489 struct Interval SqrtDeriv(struct Interval i)
490 {
491 return DivIntervals(CreateThinInteger(1L),
492 MulIntervals(CreateThinInteger(2L),
493 SqrtInterval(i)));
494 }
495
496 double MinOfRange(double lower, double upper)
497 {
498 return lower;
499 }
500
501 double MaxOfRange(double lower, double upper)
502 {
503 return upper;
504 }
505
506 double ArgMinSqr(double lower, double upper)
507 {
508 if (upper < 0.0) return upper;
509 if (lower > 0.0) return lower;
510 return 0.0;
511 }
512
513 double ArgMaxSqr(double lower, double upper)
514 {
515 return (ABS(lower)>ABS(upper))?lower:upper;
516 }
517
518 double ConvexOfSqr(double x, double lower, double upper,
519 double (*value) (/* ??? */))
520 {
521 return x*x;
522 }
523
524 double ConvexDOfSqr(double x, double lower, double upper,
525 double (*value) (/* ??? */))
526 {
527 return 2*x;
528 }
529
530 double ConcaveOfSqr(double x, double lower, double upper,
531 double (*value) (/* ??? */))
532 {
533 return (lower+upper)*x-lower*upper;
534 }
535
536 double ConcaveDOfSqr(double x, double lower, double upper,
537 double (*value) (/* ??? */))
538 {
539 return lower+upper;
540 }
541
542 double ConvexOfExp(double x, double lower, double upper,
543 double (*value) (/* ??? */))
544 {
545 return exp(x);
546 }
547
548 double ConcaveOfLn(double x, double lower, double upper,
549 double (*value) (/* ??? */))
550 {
551 return log(x);
552 }
553
554 double ConcaveDOfLn(double x, double lower, double upper,
555 double (*value) (/* ??? */))
556 {
557 return 1.0/x;
558 }
559
560 double ConcaveOfLog(double x, double lower, double upper,
561 double (*value) (/* ??? */))
562 {
563 return log10(x);
564 }
565
566 double ConcaveDOfLog(double x, double lower, double upper,
567 double (*value) (/* ??? */))
568 {
569 return M_LOG10E/x;
570 }
571
572 double ConcaveOfSqrt(double x, double lower, double upper,
573 double (*value) (/* ??? */))
574 {
575 return sqrt(x);
576 }
577
578 double ConcaveDOfSqrt(double x, double lower, double upper,
579 double (*value) (/* ??? */))
580 {
581 return 0.5/sqrt(x);
582 }
583
584 double Interpolate(double x, double lower, double upper,
585 double (*value) (/* ??? */))
586 {
587 register double vl,vu;
588 vl = (*value)(lower);
589 vu = (*value)(upper);
590 return ((vu-vl)*x+upper*vl-lower*vu)/(upper-lower);
591 }
592
593 double InterpolateD(double x, double lower, double upper,
594 double (*value) (/* ??? */))
595 {
596 return ((*value)(upper)-(*value)(lower))/(upper-lower);
597 }
598
599 #endif
600
601 struct Func g_exp_f = {
602 "exp",
603 "exp",
604 "Exp",
605 "exp",
606 "exp",
607 F_EXP,
608 exp,
609 exp,
610 exp,
611 safe_exp_D0,
612 safe_exp_D1,
613 safe_exp_D2,
614 #ifdef CHRIS_FUNC
615 ExpInterval,
616 ExpSlope,
617 ExpInterval,
618 MinOfRange,
619 MaxOfRange,
620 ConvexOfExp,
621 ConvexOfExp,
622 Interpolate,
623 InterpolateD
624 #endif
625 };
626
627 struct Func g_ln_f = {
628 "ln",
629 "log",
630 "Ln",
631 "dln",
632 "dln2",
633 F_LN,
634 log,
635 dln,
636 dln2,
637 safe_ln_D0,
638 safe_ln_D1,
639 safe_ln_D2,
640 #ifdef CHRIS_FUNC
641 LnInterval,
642 LnSlope,
643 LnDeriv,
644 MinOfRange,
645 MaxOfRange,
646 Interpolate,
647 InterpolateD,
648 ConcaveOfLn,
649 ConcaveDOfLn
650 #endif
651 };
652
653 struct Func g_lnm_f = {
654 "lnm",
655 "lnm",
656 "Lnm", /** FIXME Yacas does not have an equivalent for Lnm??? */
657 "dlnm",
658 "dlnm2",
659 F_LNM,
660 lnm,
661 dlnm,
662 dlnm2,
663 safe_lnm_D0,
664 safe_lnm_D1,
665 safe_lnm_D2,
666 #ifdef CHRIS_FUNC
667 NULL,
668 NULL,
669 NULL,
670 NULL,
671 NULL,
672 NULL,
673 NULL,
674 NULL,
675 NULL
676 #endif
677 };
678
679 struct Func g_log10_f = {
680 "log10",
681 "log10",
682 "Log10", /** FIXME Yacas does not have an equivalent for Log10??? */
683 "dlog10",
684 "dlog102",
685 F_LOG10,
686 log10,
687 dlog10,
688 dlog102,
689 safe_log10_D0,
690 safe_log10_D1,
691 safe_log10_D2,
692 #ifdef CHRIS_FUNC
693 Log10Interval,
694 Log10Slope,
695 Log10Deriv,
696 MinOfRange,
697 MaxOfRange,
698 Interpolate,
699 InterpolateD,
700 ConcaveOfLog10,
701 ConcaveDOfLog10
702 #endif
703 };
704
705 struct Func g_sin_f = {
706 "sin",
707 "sin",
708 "Sin",
709 "cos",
710 "dcos",
711 F_SIN,
712 sin,
713 cos,
714 dcos,
715 safe_sin_D0,
716 safe_sin_D1,
717 safe_sin_D2,
718 #ifdef CHRIS_FUNC
719 NULL,
720 NULL,
721 NULL,
722 NULL,
723 NULL,
724 NULL,
725 NULL,
726 NULL,
727 NULL
728 #endif
729 };
730
731 struct Func g_cos_f = {
732 "cos",
733 "cos",
734 "Cos",
735 "dcos",
736 "dcos2",
737 F_COS,
738 cos,
739 dcos,
740 dcos2,
741 safe_cos_D0,
742 safe_cos_D1,
743 safe_cos_D2,
744 #ifdef CHRIS_FUNC
745 NULL,
746 NULL,
747 NULL,
748 NULL,
749 NULL,
750 NULL,
751 NULL,
752 NULL,
753 NULL
754 #endif
755 };
756
757 struct Func g_tan_f = {
758 "tan",
759 "tan",
760 "Tan",
761 "dtan",
762 "dtan2",
763 F_TAN,
764 tan,
765 dtan,
766 dtan2,
767 safe_tan_D0,
768 safe_tan_D1,
769 safe_tan_D2,
770 #ifdef CHRIS_FUNC
771 NULL,
772 NULL,
773 NULL,
774 NULL,
775 NULL,
776 NULL,
777 NULL,
778 NULL,
779 NULL
780 #endif
781 };
782
783 struct Func g_sqr_f = {
784 "sqr",
785 "sqr",
786 "Sqr", /** FIXME Yacas does not have an equivalent for Sqr??? */
787 "dsqr",
788 "dsqr2",
789 F_SQR,
790 sqr,
791 dsqr,
792 dsqr2,
793 safe_sqr_D0,
794 safe_sqr_D1,
795 safe_sqr_D2,
796 #ifdef CHRIS_FUNC
797 SqrInterval,
798 SqrSlope,
799 SqrDeriv,
800 ArgMinSqr,
801 ArgMaxSqr,
802 ConvexOfSqr,
803 ConvexDOfSqr,
804 ConcaveOfSqr,
805 ConcaveDOfSqr
806 #endif
807 };
808
809 struct Func g_sqrt_f = {
810 "sqrt",
811 "sqrt",
812 "Sqrt",
813 "dsqrt",
814 "dsqrt2",
815 F_SQRT,
816 sqrt,
817 dsqrt,
818 dsqrt2,
819 safe_sqrt_D0,
820 safe_sqrt_D1,
821 safe_sqrt_D2,
822 #ifdef CHRIS_FUNC
823 SqrtInterval,
824 SqrtSlope,
825 SqrtDeriv,
826 MinOfRange,
827 MaxOfRange,
828 Interpolate,
829 InterpolateD,
830 ConcaveOfSqrt,
831 ConcaveDOfSqrt
832 #endif
833 };
834
835 struct Func g_abs_f = {
836 "abs",
837 "fabs",
838 "Abs",
839 "dfabs",
840 "dfabs2",
841 F_ABS,
842 fabs,
843 dfabs,
844 dfabs2,
845 safe_fabs_D0,
846 safe_fabs_D1,
847 safe_fabs_D2,
848 #ifdef CHRIS_FUNC
849 NULL,
850 NULL,
851 NULL,
852 NULL,
853 NULL,
854 NULL,
855 NULL,
856 NULL,
857 NULL
858 #endif
859 };
860
861 struct Func g_hold_f = {
862 "hold",
863 "hold",
864 "Hold", /** FIXME Yacas does not have an equivalent for Hold??? */
865 "dhold",
866 "dhold2",
867 F_HOLD,
868 hold,
869 dhold,
870 dhold2,
871 safe_hold_D0,
872 safe_hold_D1,
873 safe_hold_D2,
874 #ifdef CHRIS_FUNC
875 NULL,
876 NULL,
877 NULL,
878 NULL,
879 NULL,
880 NULL,
881 NULL,
882 NULL,
883 NULL
884 #endif
885 };
886
887 struct Func g_arcsin_f = {
888 "arcsin",
889 "asin",
890 "ArcSin",
891 "dasin",
892 "dasin2",
893 F_ARCSIN,
894 asin,
895 dasin,
896 dasin2,
897 safe_arcsin_D0,
898 safe_arcsin_D1,
899 safe_arcsin_D2,
900 #ifdef CHRIS_FUNC
901 NULL,
902 NULL,
903 NULL,
904 NULL,
905 NULL,
906 NULL,
907 NULL,
908 NULL,
909 NULL
910 #endif
911 };
912
913 struct Func g_arccos_f = {
914 "arccos",
915 "acos",
916 "ArcCos",
917 "dacos",
918 "dacos2",
919 F_ARCCOS,
920 acos,
921 dacos,
922 dacos2,
923 safe_arccos_D0,
924 safe_arccos_D1,
925 safe_arccos_D2,
926 #ifdef CHRIS_FUNC
927 NULL,
928 NULL,
929 NULL,
930 NULL,
931 NULL,
932 NULL,
933 NULL,
934 NULL,
935 NULL
936 #endif
937 };
938
939 struct Func g_arctan_f = {
940 "arctan",
941 "atan",
942 "ArcTan",
943 "datan",
944 "datan2",
945 F_ARCTAN,
946 atan,
947 datan,
948 datan2,
949 safe_arctan_D0,
950 safe_arctan_D1,
951 safe_arctan_D2,
952 #ifdef CHRIS_FUNC
953 NULL,
954 NULL,
955 NULL,
956 NULL,
957 NULL,
958 NULL,
959 NULL,
960 NULL,
961 NULL
962 #endif
963 };
964
965 #ifdef HAVE_ERF
966 struct Func g_erf_f = {
967 "erf",
968 "erf",
969 "Erf", /** FIXME Yacas does not have an equivalent for Erf??? */
970 "derf",
971 "derf2",
972 F_ERF,
973 erf,
974 derf,
975 derf2,
976 safe_erf_D0,
977 safe_erf_D1,
978 safe_erf_D2,
979 #ifdef CHRIS_FUNC
980 ErfInterval,
981 ErfSlope,
982 ErfDeriv
983 #endif
984 };
985 #endif /* HAVE_ERF */
986
987 struct Func g_sinh_f = {
988 "sinh",
989 "sinh",
990 "Sinh", /** FIXME Yacas does not have an equivalent for Sinh??? */
991 "cosh",
992 "sinh",
993 F_SINH,
994 sinh,
995 cosh,
996 sinh,
997 safe_sinh_D0,
998 safe_sinh_D1,
999 safe_sinh_D2,
1000 #ifdef CHRIS_FUNC
1001 NULL,
1002 NULL,
1003 NULL,
1004 NULL,
1005 NULL,
1006 NULL,
1007 NULL,
1008 NULL,
1009 NULL
1010 #endif
1011 };
1012
1013 struct Func g_cosh_f = {
1014 "cosh",
1015 "cosh",
1016 "Cosh", /** FIXME Yacas does not have an equivalent for Cosh??? */
1017 "sinh",
1018 "cosh",
1019 F_COSH,
1020 cosh,
1021 sinh,
1022 cosh,
1023 safe_cosh_D0,
1024 safe_cosh_D1,
1025 safe_cosh_D2,
1026 #ifdef CHRIS_FUNC
1027 NULL,
1028 NULL,
1029 NULL,
1030 NULL,
1031 NULL,
1032 NULL,
1033 NULL,
1034 NULL,
1035 NULL
1036 #endif
1037 };
1038
1039 struct Func g_tanh_f = {
1040 "tanh",
1041 "tanh",
1042 "Tanh", /** FIXME Yacas does not have an equivalent for Tanh??? */
1043 "dtanh",
1044 "dtanh2",
1045 F_TANH,
1046 tanh,
1047 dtanh,
1048 dtanh2,
1049 safe_tanh_D0,
1050 safe_tanh_D1,
1051 safe_tanh_D2,
1052 #ifdef CHRIS_FUNC
1053 NULL,
1054 NULL,
1055 NULL,
1056 NULL,
1057 NULL,
1058 NULL,
1059 NULL,
1060 NULL,
1061 NULL
1062 #endif
1063 };
1064
1065 struct Func g_arcsinh_f = {
1066 "arcsinh",
1067 "arcsinh",
1068 "ArcSinh", /** FIXME Yacas does not have an equivalent for ArcSinh??? */
1069 "darcsinh",
1070 "darcsinh2",
1071 F_ARCSINH,
1072 arcsinh,
1073 darcsinh,
1074 darcsinh2,
1075 safe_arcsinh_D0,
1076 safe_arcsinh_D1,
1077 safe_arcsinh_D2,
1078 #ifdef CHRIS_FUNC
1079 NULL,
1080 NULL,
1081 NULL,
1082 NULL,
1083 NULL,
1084 NULL,
1085 NULL,
1086 NULL,
1087 NULL
1088 #endif
1089 };
1090
1091 struct Func g_arccosh_f = {
1092 "arccosh",
1093 "arccosh",
1094 "ArcCosh", /** FIXME Yacas does not have an equivalent for ArcCosh??? */
1095 "darccosh",
1096 "darccosh2",
1097 F_ARCCOSH,
1098 arccosh,
1099 darccosh,
1100 darccosh2,
1101 safe_arccosh_D0,
1102 safe_arccosh_D1,
1103 safe_arccosh_D2,
1104 #ifdef CHRIS_FUNC
1105 NULL,
1106 NULL,
1107 NULL,
1108 NULL,
1109 NULL,
1110 NULL,
1111 NULL,
1112 NULL,
1113 NULL
1114 #endif
1115 };
1116
1117 struct Func g_arctanh_f = {
1118 "arctanh",
1119 "arctanh",
1120 "ArcTanh", /** FIXME Yacas does not have an equivalent for ArcTanh??? */
1121 "darctanh",
1122 "darctanh2",
1123 F_ARCTANH,
1124 arctanh,
1125 darctanh,
1126 darctanh2,
1127 safe_arctanh_D0,
1128 safe_arctanh_D1,
1129 safe_arctanh_D2,
1130 #ifdef CHRIS_FUNC
1131 NULL,
1132 NULL,
1133 NULL,
1134 NULL,
1135 NULL,
1136 NULL,
1137 NULL,
1138 NULL,
1139 NULL
1140 #endif
1141 };
1142
1143 struct Func g_cube_f = {
1144 "cube",
1145 "cube",
1146 "Cube", /** FIXME Yacas does not have an equivalent for Cube??? */
1147 "dcube",
1148 "dcube2",
1149 F_CUBE,
1150 cube,
1151 dcube,
1152 dcube2,
1153 safe_cube,
1154 safe_cube_D1,
1155 safe_cube_D2,
1156 #ifdef CHRIS_FUNC
1157 NULL,
1158 NULL,
1159 NULL,
1160 NULL,
1161 NULL,
1162 NULL,
1163 NULL,
1164 NULL,
1165 NULL
1166 #endif
1167 };
1168
1169 struct Func g_cbrt_f = {
1170 "cbrt",
1171 "cbrt",
1172 "Cbrt", /** FIXME Yacas does not have an equivalent for Cbrt??? */
1173 "dcbrt",
1174 "dcbrt2",
1175 F_CBRT,
1176 cbrt,
1177 dcbrt,
1178 dcbrt2,
1179 safe_cbrt_D0,
1180 safe_cbrt_D1,
1181 safe_cbrt_D2,
1182 #ifdef CHRIS_FUNC
1183 NULL,
1184 NULL,
1185 NULL,
1186 NULL,
1187 NULL,
1188 NULL,
1189 NULL,
1190 NULL,
1191 NULL
1192 #endif
1193 };
1194
1195
1196 struct Func *g_func_list[]={
1197 &g_log10_f,
1198 &g_ln_f,
1199 &g_exp_f,
1200 &g_sin_f,
1201 &g_cos_f,
1202 &g_tan_f,
1203 &g_sqr_f,
1204 &g_sqrt_f,
1205 &g_arcsin_f,
1206 &g_arccos_f,
1207 &g_arctan_f,
1208 #ifdef HAVE_ERF
1209 &g_erf_f,
1210 #endif /* HAVE_ERF */
1211 &g_lnm_f,
1212 &g_sinh_f,
1213 &g_cosh_f,
1214 &g_tanh_f,
1215 &g_arcsinh_f,
1216 &g_arccosh_f,
1217 &g_arctanh_f,
1218 &g_cube_f,
1219 &g_cbrt_f,
1220 &g_abs_f,
1221 &g_hold_f,
1222 NULL /* must be last */
1223 };
1224
1225 CONST struct Func *LookupFunc(CONST char *name)
1226 {
1227 unsigned f=0;
1228 while(g_func_list[f]!=NULL){
1229 if(strcmp(g_func_list[f]->name,name)==0)
1230 return g_func_list[f];
1231 f++;
1232 }
1233 return NULL;
1234 }
1235
1236 CONST struct Func *LookupFuncById(enum Func_enum id)
1237 {
1238 unsigned f=0;
1239 while(g_func_list[f]!=NULL){
1240 if (g_func_list[f]->id==id)
1241 return g_func_list[f];
1242 f++;
1243 }
1244 return NULL;
1245 }
1246
1247 CONST char *FuncName(CONST struct Func *f)
1248 {
1249 return f->name;
1250 }
1251
1252 CONST char *FuncCName(CONST struct Func *f)
1253 {
1254 return f->cname;
1255 }
1256
1257 CONST char *FuncYName(CONST struct Func *f)
1258 {
1259 return f->yname;
1260 }
1261
1262 CONST char *FuncDeriv1CName(CONST struct Func *f)
1263 {
1264 return f->deriv1cname;
1265 }
1266
1267 CONST char *FuncDeriv2CName(CONST struct Func *f)
1268 {
1269 return f->deriv2cname;
1270 }
1271
1272 enum Func_enum FuncId(CONST struct Func *f)
1273 {
1274 return f->id;
1275 }
1276
1277 CONST dim_type *FuncDimens(CONST struct Func *f)
1278 {
1279 if (!f) return Dimensionless();
1280 switch (FuncId(f)) {
1281 case F_LOG10:
1282 case F_LN:
1283 case F_EXP:
1284 #ifdef HAVE_ERF
1285 case F_ERF:
1286 #endif /* HAVE_ERF */
1287 case F_LNM:
1288 case F_ARCSIN:
1289 case F_ARCCOS:
1290 case F_ARCTAN:
1291 case F_SINH:
1292 case F_COSH:
1293 case F_TANH:
1294 case F_ARCSINH:
1295 case F_ARCCOSH:
1296 case F_ARCTANH:
1297 return Dimensionless();
1298 case F_SQR:
1299 case F_SQRT:
1300 case F_CUBE:
1301 case F_CBRT:
1302 case F_ABS:
1303 case F_HOLD:
1304 return WildDimension();
1305 case F_SIN:
1306 case F_COS:
1307 case F_TAN:
1308 return TrigDimension();
1309 default: return Dimensionless();
1310 }
1311 }
1312
1313 double FuncEval(CONST struct Func *f, double d)
1314 {
1315 return (*(f->value))(d);
1316 }
1317
1318 double FuncEvalSafe(CONST struct Func *f, double d,enum safe_err *not_safe)
1319 {
1320 return (*(f->safevalue))(d,not_safe);
1321 }
1322
1323 double FuncDeriv(CONST struct Func *f, double d)
1324 {
1325 return (*(f->deriv))(d);
1326 }
1327
1328 double FuncDerivSafe(CONST struct Func *f, double d,enum safe_err *not_safe)
1329 {
1330 return (*(f->safederiv))(d,not_safe);
1331 }
1332
1333 double FuncDeriv2(CONST struct Func *f, double d)
1334 {
1335 return (*(f->deriv2))(d);
1336 }
1337
1338 double FuncDeriv2Safe(CONST struct Func *f, double d,enum safe_err *not_safe)
1339 {
1340 return (*(f->safederiv2))(d,not_safe);
1341 }
1342
1343 #ifdef CHRIS_FUNC
1344 struct Interval FuncRange(CONST struct Func *f, struct Interval u)
1345 {
1346 return (*(f->ivalue))(u);
1347 }
1348
1349 void FuncSlope(CONST struct Func *f, unsigned long int nvar,
1350 struct Interval *center, struct Interval *range,
1351 struct Interval *slope)
1352 {
1353 (*(f->slope))(nvar,center,range,slope);
1354 }
1355
1356 struct Interval FuncIDeriv(CONST struct Func *f, struct Interval i)
1357 {
1358 return (*(f->ideriv))(i);
1359 }
1360
1361 double ArgMin(CONST struct Func *f, double lower, double upper)
1362 {
1363 return (*(f->tmin))(lower,upper);
1364 }
1365
1366 double ArgMax(CONST struct Func *f, double lower, double upper)
1367 {
1368 return (*(f->tmax))(lower,upper);
1369 }
1370
1371 double ConvexEnv(CONST struct Func *f, double x, double lower, double upper)
1372 {
1373 assert((x>=lower)&&(x<=upper));
1374 return (*(f->e))(x,lower,upper,f->value);
1375 }
1376
1377 double ConvexEnvDeriv(CONST struct Func *f, double x,
1378 double lower, double upper)
1379 {
1380 assert((x>=lower)&&(x<=upper));
1381 return (*(f->ed))(x,lower,upper,f->value);
1382 }
1383
1384 double ConcaveEnv(CONST struct Func *f, double x, double lower, double upper)
1385 {
1386 assert((x>=lower)&&(x<=upper));
1387 return (*(f->E))(x,lower,upper,f->value);
1388 }
1389
1390 double ConcaveEnvDeriv(CONST struct Func *f, double x,
1391 double lower, double upper)
1392 {
1393 assert((x>=lower)&&(x<=upper));
1394 return (*(f->Ed))(x,lower,upper,f->value);
1395 }
1396 #endif

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22