/[ascend]/trunk/base/generic/compiler/safe.c
ViewVC logotype

Contents of /trunk/base/generic/compiler/safe.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 909 - (show annotations) (download) (as text)
Thu Oct 26 12:44:41 2006 UTC (17 years, 9 months ago) by johnpye
File MIME type: text/x-csrc
File size: 27496 byte(s)
Added finite-difference evaluation of gradients in blackboxes.
Some work on J*v evaluation with IDA.
1 /*
2 * SLV: Ascend Numeric Solver
3 * by Karl Michael Westerberg
4 * Created: 2/6/90
5 * Version: $Revision: 1.16 $
6 * Version control file: $RCSfile: safe.c,v $
7 * Date last modified: $Date: 1998/06/11 15:28:35 $
8 * Last modified by: $Author: ballan $
9 *
10 * This file is part of the SLV solver.
11 *
12 * Copyright (C) 1990 Karl Michael Westerberg
13 * Copyright (C) 1993 Joseph Zaher
14 * Copyright (C) 1994 Joseph Zaher, Benjamin Andrew Allan
15 * Copyright (C) 1996 Benjamin Andrew Allan, Kenneth Tyner
16 *
17 * The SLV solver is free software; you can redistribute
18 * it and/or modify it under the terms of the GNU General Public License as
19 * published by the Free Software Foundation; either version 2 of the
20 * License, or (at your option) any later version.
21 *
22 * The SLV solver is distributed in hope that it will be
23 * useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
24 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
25 * General Public License for more details.
26 *
27 * You should have received a copy of the GNU General Public License
28 * along with the program; if not, write to the Free Software Foundation,
29 * Inc., 675 Mass Ave, Cambridge, MA 02139 USA. Check the file named
30 * COPYING. COPYING is found in ../compiler.
31 *
32 */
33
34 #include <math.h>
35 #include <stdarg.h>
36 #include <utilities/ascConfig.h>
37 #include <utilities/ascMalloc.h>
38 #include "compiler.h"
39 #include <utilities/ascPanic.h>
40 #include "fractions.h"
41 #include "dimen.h"
42 #include "functype.h"
43 #include "func.h"
44 #include "safe.h"
45
46 #define BIGNUM 1.0e8
47
48 /*boolean safe_ok = TRUE;*/
49 static int safe_print_errors = TRUE;
50
51 void safe_error_to_stderr(enum safe_err *not_safe){
52 if (!safe_print_errors) {
53 return;
54 }
55 switch(*not_safe)
56 {
57 case safe_problem: /* error msg should have already been printed */
58 case safe_ok:
59 return;
60 case safe_div_by_zero:
61 ERROR_REPORTER_NOLINE(ASC_PROG_WARNING,"Division by zero");
62 case safe_complex_result:
63 ERROR_REPORTER_NOLINE(ASC_PROG_WARNING,"Complex result");
64 break;
65 case safe_overflow:
66 ERROR_REPORTER_NOLINE(ASC_PROG_WARNING,"Overflow");
67 break;
68 case safe_underflow:
69 ERROR_REPORTER_NOLINE(ASC_PROG_WARNING,"Underflow");
70 break;
71 case safe_range_error:
72 ERROR_REPORTER_NOLINE(ASC_PROG_WARNING,"Function range error");
73 break;
74 }
75 }
76
77
78 static double rec_1_p_sqr_Dn(double x,int n, enum safe_err *safe)
79 /**
80 *** Calculates n-th derivative of 1/(1+x^2).
81 **/
82 {
83 /* g[n](x) = -1/(1+x^2) * (n*(n-1)*g[n-2](x) + 2*n*x*g[n-1](x))
84 where g[n] := n-th derivative of 1/(1+x^2). */
85 int k;
86 double prev[3],val; /* prev[j] := g[k-j](x), val = g[0](x) */
87
88 prev[2] = val = 1.0/(1.0 + safe_sqr_D0(x,safe));
89 if( n==0 )
90 return(val);
91 prev[1] = -2.0*x*val*val; /* first derivative */
92 for( k=2 ; k<=n ; ++k ) {
93 prev[0] = -val*(safe_mul_D0((double)(k*(k-1)),prev[2],safe) +
94 safe_mul_D0(2*k*x,prev[1],safe));
95 prev[2] = prev[1];
96 prev[1] = prev[0];
97 }
98
99 return(prev[1]);
100 }
101
102 static double *alloc_poly(int order)
103 /**
104 *** Allocates a polynominal of given order and returns it. The
105 *** polynominal need not be freed, but this function should not be
106 *** called again until the old polynominal is not needed anymore.
107 **/
108 {
109 static double *poly = NULL;
110 static int poly_cap = 0;
111
112 if( order + 1 > poly_cap ) {
113 poly_cap = order+1;
114 if( poly != NULL ) {
115 ascfree( poly );
116 }
117 poly = ASC_NEW_ARRAY(double,poly_cap );
118 }
119 return(poly);
120 }
121
122 static double safe_poly(double *poly, int order,
123 double x, double x2, enum safe_err *safe)
124 /* x2 = x^2 */
125 /**
126 *** Calculates the value of the given polynomial, where only the
127 *** (order%1 ? odd : even) degree terms are used.
128 **/
129 {
130 double val;
131 (void)safe;
132
133 for( val=poly[order] ; (order -= 2) >= 0 ; ) {
134 val = safe_mul_D0(val,x2,safe) + poly[order];
135 }
136 return( order==-1 ? safe_mul_D0(val,x,safe) : val );
137 }
138
139 #ifdef HAVE_ERF
140 static double exp_msqr_Dn(double x,int n,enum safe_err *safe)
141 /**
142 *** Computes n-th derivative of exp(-x^2).
143 **/
144 {
145 /**
146 *** n-th derivative of exp(-x^2) = f[n](x)*exp(-x^2), where f[n] is an
147 *** n-th degree polynominal of definite parity satisfying f[0](x) = 1
148 *** & f[n+1](x) = -2x*f[n](x) + f[n]'(x).
149 **/
150
151 double x2 = safe_sqr_D0(x,safe);
152 int k,r;
153 double *poly;
154 poly = alloc_poly(n);
155
156 poly[0] = exp(-x2);
157 for( k=0 ; k<n ; ++k ) { /* Calculate f[k+1] from f[k] */
158 poly[k+1] = 0.0;
159 for( r=k ; r >= 1 ; r -= 2 )
160 poly[r-1] = (double)r * poly[r];
161 for( r=k ; r >= 0 ; r -= 2 )
162 poly[r+1] += -2.0*poly[r];
163 }
164 return( safe_poly(poly,n,x,x2,safe) );
165 }
166 #endif
167
168 static double sqrt_rec_1_m_sqr_Dn(double x,int n,enum safe_err *safe)
169 /**
170 *** Calculates n-th derivative of (1-x^2)^-.5
171 **/
172 {
173 /**
174 *** n-th derivative of (1-x^2)^-.5 = f[n](x) * (1-x^2)^-(n+.5), where
175 *** f[n] is an n-degree polynominal of definite parity satisfying
176 *** f[0] = 1 and f[n+1](x) = f[n]'(x)*(1+x^2) + (2n+1)*x*f[n](x).
177 **/
178
179 int k,r;
180 double x2;
181 double *poly;
182 x2 = safe_sqr_D0(x,safe);
183 poly = alloc_poly(n);
184
185 (void)safe;
186
187 poly[0] = safe_rec(safe_upower(safe_sqrt_D0(1.0-x2,safe),2*n+1,safe),safe);
188 for( k=0 ; k<n ; ++k ) { /* Calculate f[k+1] from f[k] */
189 poly[k+1] = 0.0;
190 for( r=k ; r >= 1 ; r -= 2 )
191 poly[r-1] = (double)r * poly[r];
192 for( r=k ; r >= 0 ; r -= 2 )
193 poly[r+1] += (double)(r+2*k+1) * poly[r];
194 }
195 return( safe_poly(poly,n,x,x2,safe) );
196 }
197
198
199 double safe_upower(double x, unsigned n, enum safe_err *safe)
200 {
201 double y = 1.0;
202 (void)safe;
203
204 for( ; n-- > 0 ; y = safe_mul_D0(y,x,safe) );
205 return(y);
206 }
207
208 double safe_factorial(unsigned n, enum safe_err *safe)
209 {
210 double x,y;
211 (void)safe;
212
213 for( x = (double)n , y = 1.0 ; x >= 0.5 ; y = safe_mul_D0(y,x--,safe) )
214 ;
215 return(y);
216 }
217
218 double safe_rec(double x,enum safe_err *safe)
219 {
220 if( x == 0.0 ) {
221 double bogus = BIGNUM;
222 if( safe_print_errors ) {
223 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Divide by zero in 1/x expr: returning %g",bogus);
224 }
225 *safe = safe_div_by_zero;
226 return( bogus );
227 }
228 return( 1.0/x );
229 }
230
231 #ifndef INV_CUBRTHUGE
232 #define INV_CUBRTHUGE 1.7718548704178434e-103
233 /* smallest cubeable number in an 8 byte ieee double */
234 #endif
235 #ifndef CUBRTHUGE
236 #define CUBRTHUGE 5.6438030941223618e102
237 /* largest cubeable number in an 8 byte ieee double */
238 #endif
239 double safe_cube(register double x,enum safe_err *safe)
240 {
241 if( fabs(x) > CUBRTHUGE || fabs(x) < INV_CUBRTHUGE ) {
242 double bogus;
243 if ( fabs(x) > CUBRTHUGE) {
244 bogus=BIGNUM;
245 if( safe_print_errors ) {
246 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Overflow in x^3 expr: returning %g.",bogus);
247 }
248 } else {
249 bogus=0.0;
250 if( safe_print_errors ) {
251 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Underflow in x^3 expr: returning %g.",bogus);
252 }
253 }
254 *safe = safe_underflow;
255 return( bogus );
256 }
257 return ( x*x*x );
258 }
259
260 double safe_Dn_rec(double x,int n,enum safe_err *safe)
261 {
262 return( -safe_upower(-safe_rec(x,safe),n+1,safe) * safe_factorial(n,safe) );
263 }
264
265 #ifdef HAVE_ERF
266 double safe_erf_inv(double x,enum safe_err *safe)
267 {
268 #define CONV (1e-7)
269 double y,dy,sign;
270
271 sign = (x<0.0) ? -1.0 : 1.0;
272 x *= sign;
273 if( x >= 1.0 ) {
274 double bogus = sign*BIGNUM;
275 if( safe_print_errors ) {
276 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"inv_erf undefined at %g: returning %g",x,bogus);
277 }
278 *safe = safe_range_error;
279 return(bogus);
280 }
281
282 y = 0.0;
283 do {
284 dy = (x-safe_erf_D0(y,safe))*
285 safe_exp_D0(safe_mul_D0(y,y,safe),safe)/safe_ERF_COEF;
286 y += dy;
287 /**
288 *** Since erf is concave and x >= erf(y0), dy should
289 *** always be positive (?).
290 **/
291 if( dy < 0.0 && safe_print_errors ) {
292 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Found negative slope of %g on erf_inv.",dy);
293 }
294 } while( safe_ok && dy > CONV );
295 return( y * sign );
296 #undef CONV
297 }
298 #endif /* HAVE_ERF */
299
300 double safe_lnm_inv(double x,enum safe_err *safe)
301 {
302 register double eps=FuncGetLnmEpsilon();
303
304 if( x > (DBL_MAX > 1.0e308 ? 709.196209 : log(DBL_MAX)) ) {
305 double bogus = BIGNUM;
306 if( safe_print_errors ) {
307 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Argument %g too large in lnm_inv: returning %g.",x,bogus);
308 }
309 *safe = safe_range_error;
310 return(bogus);
311 }
312 /* since lnm(eps)= log(eps), use cheaper one. eps known >0 */
313 return( (x >= log(eps)) ? safe_exp_D0(x,safe) : eps*(x+1.0-log(eps)) );
314 }
315
316 int safe_is_int(double x,enum safe_err *safe)
317 {
318 #define TOLERANCE (1e-4)
319 unsigned n;
320 (void)safe;
321
322 if( x<0.0 ) {
323 x = -x;
324 }
325
326 /**
327 *** Theorem: Define P(n) :<==> x-TOLERANCE <= n <= x+TOLERANCE. Then
328 *** P(n) for some integer n <==> P(floor(x+TOLERANCE))
329 *** <==> x-TOLERANCE <= floor(x+TOLERANCE).
330 ***
331 *** Proof: Since floor(x+TOLERANCE) <= x+TOLERANCE, the second
332 *** equivalence holds. The backward direction of the first
333 *** equivalence is a tautology. If P(n) holds for some n, then
334 *** n <= floor(x+TOLERANCE) <= x+TOLERANCE, the second inequality
335 *** follows since floor(z) <= z, and the first inequality follows
336 *** since floor(z) is the largest integer with that property.
337 **/
338 if( x-TOLERANCE <= (double)(n = (unsigned)(x+TOLERANCE)) ) {
339 return( n&01 );
340 } else {
341 return( -1 );
342 }
343 #undef TOLERANCE
344 }
345
346 double safe_sin_D0(double x,enum safe_err *safe)
347 {
348 (void)safe;
349 return(sin(x));
350 }
351
352 double safe_sinh_D0(double x,enum safe_err *safe)
353 {
354 (void)safe;
355 return(sinh(x));
356 }
357
358 double safe_cos_D0(double x,enum safe_err *safe)
359 {
360 (void)safe;
361 return(cos(x));
362 }
363
364 double safe_cosh_D0(double x,enum safe_err *safe)
365 {
366 (void)safe;
367 return(cosh(x));
368 }
369
370 double safe_tan_D0(double x,enum safe_err *safe)
371 {
372 return( safe_div_D0( sin(x) , cos(x) , safe) );
373 }
374
375 double safe_tanh_D0(double x,enum safe_err *safe)
376 {
377 (void)safe;
378 return( tanh(x) );
379 }
380
381 double safe_arcsin_D0(double x,enum safe_err *safe)
382 {
383 if( x < -1.0 || 1.0 < x ) {
384 double bogus = 0.0;
385 if( safe_print_errors ) {
386 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"arcsin_D0: arcsin undefined at %g: Returning %g.",x,bogus);
387 }
388 *safe = safe_range_error;
389 return(bogus);
390 }
391 return( asin(x) );
392 }
393
394 double safe_arcsinh_D0(double x,enum safe_err *safe)
395 {
396 (void)safe;
397 return( arcsinh(x) );
398 }
399
400 double safe_arccos_D0(double x,enum safe_err *safe)
401 {
402 if( x < -1.0 || 1.0 < x ) {
403 double bogus = 0.0;
404 if( safe_print_errors ) {
405 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"arccos_D0: arccos undefined at %g: returning %g.",x,bogus);
406 }
407 *safe = safe_range_error;
408 return(bogus);
409 }
410 return( acos(x) );
411 }
412
413 double safe_arccosh_D0(double x,enum safe_err *safe)
414 {
415 if( x < 1.0 ) {
416 double bogus = 0.0;
417 if( safe_print_errors ) {
418 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"arccosh_D0: undefined at %g: returning %g.",x,bogus);
419 }
420 *safe = safe_range_error;
421 return(bogus);
422 }
423
424 return( arccosh(x) );
425 }
426
427 double safe_arctan_D0(double x,enum safe_err *safe)
428 {
429 (void)safe;
430 return( atan(x) );
431 }
432
433 double safe_arctanh_D0(double x,enum safe_err *safe)
434 {
435 if( x < -1.0 || 1.0 < x ) {
436 double bogus = 0.0;
437 if( safe_print_errors ) {
438 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"arctanh_D0: undefined at %g: returning %g.",x,bogus);
439 }
440 *safe = safe_range_error;
441 return(bogus);
442 }
443 return( arctanh(x) );
444 }
445
446 #ifdef HAVE_ERF
447 double safe_erf_D0(double x,enum safe_err *safe)
448 {
449 (void)safe;
450 return( erf(x) );
451 }
452 #endif /* HAVE_ERF */
453
454 double safe_exp_D0(double x,enum safe_err *safe)
455 {
456 if( x > (DBL_MAX > 1.0e308 ? 709.196209 : log(DBL_MAX)) ) {
457 double bogus = BIGNUM;
458 if( safe_print_errors ) {
459 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"exp_D0: Argument %g too large: returning %g.",x,bogus);
460 }
461 *safe = safe_range_error;
462 return(bogus);
463 }
464 return( exp(x) );
465 }
466
467 double safe_ln_D0(double x,enum safe_err *safe)
468 {
469 if( x <= 0.0 ) {
470 double bogus = -BIGNUM;
471 if( safe_print_errors ) {
472 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"ln_D0: undefined at %g: Returning %g.",x,bogus);
473 }
474 *safe = safe_range_error;
475 return(bogus);
476 }
477 return( log(x) );
478 }
479
480 double safe_lnm_D0(double x,enum safe_err *safe)
481 {
482 (void)safe;
483 return( lnm(x) );
484 }
485
486 double safe_log10_D0(double x,enum safe_err *safe)
487 {
488 return( safe_LOG10_COEF * safe_ln_D0(x,safe) );
489 }
490
491 double safe_sqr_D0(double x,enum safe_err *safe)
492 {
493 (void)safe;
494 return( sqr(x) );
495 }
496
497 double safe_sqrt_D0(double x,enum safe_err *safe)
498 {
499 if( x < 0.0 ) {
500 double bogus;
501 bogus = -sqrt(-x);
502 if( safe_print_errors ) {
503 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"sqrt_D0: undefined at %g: returning %g.",x,bogus);
504 }
505 *safe = safe_complex_result;
506 return(bogus);
507 }
508 return(sqrt(x));
509 }
510
511 double safe_cube_D0(double x,enum safe_err *safe)
512 {
513 (void)safe;
514 return( cbrt(x) );
515 }
516
517 double safe_cbrt_D0(double x,enum safe_err *safe)
518 {
519 (void)safe;
520 return( cbrt(x) );
521 }
522
523 double safe_fabs_D0(double x,enum safe_err *safe)
524 {
525 (void)safe;
526 return( fabs(x) );
527 }
528
529 double safe_hold_D0(double x,enum safe_err *safe)
530 {
531 (void)safe;
532 return x;
533 }
534
535 double safe_sin_D1(double x,enum safe_err *safe)
536 {
537 (void)safe;
538 return( cos(x) );
539 }
540
541 double safe_sinh_D1(double x,enum safe_err *safe)
542 {
543 (void)safe;
544 return( cosh(x) );
545 }
546
547 double safe_cos_D1(double x,enum safe_err *safe)
548 {
549 (void)safe;
550 return( sin(x) );
551 }
552
553 double safe_cosh_D1(double x,enum safe_err *safe)
554 {
555 (void)safe;
556 return( sinh(x) );
557 }
558
559 double safe_tan_D1(double x,enum safe_err *safe)
560 {
561 x = cos(x);
562 return( safe_rec(safe_sqr_D0(x,safe),safe) );
563 }
564
565 double safe_tanh_D1(double x,enum safe_err *safe)
566 {
567 (void)safe;
568 if (x < -200 || x > 200) {
569 return 0.0;
570 }
571 return( dtanh(x) );
572 }
573
574 double safe_arcsin_D1(double x,enum safe_err *safe)
575 {
576 return( safe_rec(safe_sqrt_D0(1.0 - safe_sqr_D0(x,safe),safe),safe) );
577 }
578
579 double safe_arcsinh_D1(double x,enum safe_err *safe)
580 {
581 (void)safe;
582 return( darcsinh(x) );
583 }
584
585 double safe_arccos_D1(double x,enum safe_err *safe)
586 {
587 return( -safe_arcsin_D1(x,safe) );
588 }
589
590 double safe_arccosh_D1(double x,enum safe_err *safe)
591 {
592 return( safe_rec(safe_sqrt_D0(safe_sqr_D0(x,safe) - 1.0,safe),safe) );
593 }
594
595 double safe_arctan_D1(double x,enum safe_err *safe)
596 {
597 (void)safe;
598 return( datan(x) );
599 }
600
601 double safe_arctanh_D1(double x,enum safe_err *safe)
602 {
603 return( safe_rec(1.0 - safe_sqr_D0(x,safe),safe) );
604 }
605
606 #ifdef HAVE_ERF
607 double safe_erf_D1(double x,enum safe_err *safe)
608 {
609 return( safe_ERF_COEF * safe_exp_D0(-safe_sqr_D0(x,safe),safe) );
610 }
611 #endif /* HAVE_ERF */
612
613 double safe_exp_D1(double x,enum safe_err *safe)
614 {
615 return( safe_exp_D0(x,safe) );
616 }
617
618 double safe_ln_D1(double x,enum safe_err *safe)
619 {
620 return( safe_rec(x,safe) );
621 }
622
623 double safe_lnm_D1(double x,enum safe_err *safe)
624 {
625 (void)safe;
626 return( dlnm(x) );
627 }
628
629 double safe_log10_D1(double x,enum safe_err *safe)
630 {
631 return( safe_LOG10_COEF * safe_ln_D1(x,safe) );
632 }
633
634 double safe_sqr_D1(double x,enum safe_err *safe)
635 {
636 (void)safe;
637 return( dsqr(x) );
638 }
639
640 double safe_sqrt_D1(double x,enum safe_err *safe)
641 {
642 return( 0.5 * safe_rec(safe_sqrt_D0(x,safe),safe) );
643 }
644
645 double safe_cube_D1(double x,enum safe_err *safe)
646 {
647 return( 3 * safe_sqr_D0(x,safe) );
648 }
649
650 double safe_cbrt_D1(double x,enum safe_err *safe)
651 {
652 return( 0.3333333333333333*
653 safe_rec(safe_sqr_D0(safe_cbrt_D0(x,safe),safe),safe) );
654 }
655
656 double safe_fabs_D1(double x,enum safe_err *safe)
657 {
658 (void)safe;
659 (void)x; /* stop gcc whine about unused parameter */
660 return( 1.0 );
661 }
662
663 double safe_hold_D1(double x,enum safe_err *safe)
664 {
665 (void)x; /* stop gcc whine about unused parameter */
666 (void)safe;
667 return( 0.0 );
668 }
669
670 double safe_sin_D2(double x,enum safe_err *safe)
671 {
672 (void)safe;
673 return( dcos(x) );
674 }
675
676 double safe_sinh_D2(double x,enum safe_err *safe)
677 {
678 (void)safe;
679 return( sinh(x) );
680 }
681
682 double safe_cos_D2(double x,enum safe_err *safe)
683 {
684 (void)safe;
685 return( dcos2(x) );
686 }
687
688 double safe_cosh_D2(double x,enum safe_err *safe)
689 {
690 (void)safe;
691 return( cosh(x) );
692 }
693
694 double safe_tan_D2(double x,enum safe_err *safe)
695 {
696 if( fabs(cos(x)) == 0.0 ) {
697 double bogus = BIGNUM;
698 if( safe_print_errors ) {
699 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"tan_D2: Infinite at %g: returning %g.",x,bogus);
700 }
701 *safe = safe_range_error;
702 return(bogus);
703 }
704 return dtan2(x);
705 }
706
707 double safe_tanh_D2(double x,enum safe_err *safe)
708 {
709 (void)safe;
710 return( dtanh2(x) );
711 }
712
713 double safe_arcsin_D2(double x,enum safe_err *safe)
714 {
715 register double t;
716 t = safe_rec(1.0 - safe_sqr_D0(x,safe),safe);
717 return( safe_mul_D0( safe_mul_D0(x,t,safe) , safe_mul_D0(t,t,safe) ,safe) );
718 }
719
720 double safe_arcsinh_D2(double x,enum safe_err *safe)
721 {
722 (void)safe;
723 return( darcsinh2(x) );
724 }
725
726 double safe_arccos_D2(double x,enum safe_err *safe)
727 {
728 return( -safe_arcsin_D2(x,safe) );
729 }
730
731 double safe_arccosh_D2(double x,enum safe_err *safe)
732 {
733 if( fabs(x) <= 1.0 ) {
734 double bogus = -BIGNUM;
735 if( safe_print_errors ) {
736 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"arccosh_D2: Undefined at %g: returning %g",x,bogus);
737 }
738 *safe = safe_range_error;
739 return(bogus);
740 }
741 return darccosh2(x);
742 }
743
744 double safe_arctan_D2(double x,enum safe_err *safe)
745 {
746 (void)safe;
747 return( datan2(x) );
748 }
749
750 double safe_arctanh_D2(double x,enum safe_err *safe)
751 {
752 if( fabs(x) == 1.0 ) {
753 double bogus = BIGNUM;
754 if( safe_print_errors ) {
755 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"arctanh_D2: undefined at %g: returning %g.",x,bogus);
756 }
757 *safe = safe_range_error;
758 return(bogus);
759 }
760 return( darctanh2(x) );
761 }
762
763 #ifdef HAVE_ERF
764 double safe_erf_D2(double x,enum safe_err *safe)
765 {
766 return( -ldexp(safe_ERF_COEF *
767 safe_mul_D0(x,safe_exp_D0(-safe_sqr_D0(x,safe),safe),safe),1) );
768 }
769 #endif /* HAVE_ERF */
770
771 double safe_exp_D2(double x,enum safe_err *safe)
772 {
773 return( safe_exp_D0(x,safe) );
774 }
775
776 double safe_ln_D2(double x,enum safe_err *safe)
777 {
778 return( -safe_rec(safe_sqr_D0(x,safe),safe) );
779 }
780
781 double safe_lnm_D2(double x,enum safe_err *safe)
782 {
783 (void)safe;
784 return( dlnm2(x) );
785 }
786
787 double safe_log10_D2(double x,enum safe_err *safe)
788 {
789 return( safe_LOG10_COEF * safe_ln_D2(x,safe) );
790 }
791
792 double safe_sqr_D2(double x,enum safe_err *safe)
793 {
794 (void)safe;
795 return( dsqr2(x) );
796 }
797
798 double safe_sqrt_D2(double x,enum safe_err *safe)
799 {
800 return(-ldexp(safe_rec(safe_mul_D0(x,safe_sqrt_D0( x,safe),safe),safe),-2));
801 }
802
803 double safe_cube_D2(double x,enum safe_err *safe)
804 {
805 (void)safe;
806 return( safe_mul_D0(6,x,safe) );
807 }
808
809 double safe_cbrt_D2(double x,enum safe_err *safe)
810 {
811 if( fabs(x) == 0.0 ) {
812 double bogus = BIGNUM;
813 if( safe_print_errors ) {
814 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"cbrt_D2: undefined at %g: returning %g.",x,bogus);
815 }
816 *safe = safe_range_error;
817 return(bogus);
818 }
819 return dcbrt2(x);
820 }
821
822 double safe_fabs_D2(double x,enum safe_err *safe)
823 {
824 (void)x; /* stop gcc whine about unused parameter */
825 (void)safe;
826 return( 0.0 );
827 }
828
829 double safe_arcsin_Dn(double x,int n,enum safe_err *safe)
830 {
831 return( n==0 ? safe_arcsin_D0(x,safe) : sqrt_rec_1_m_sqr_Dn(x,n-1,safe) );
832 }
833
834 double safe_arccos_Dn(double x,int n,enum safe_err *safe)
835 {
836 return( n==0 ? safe_arccos_D0(x,safe) : -sqrt_rec_1_m_sqr_Dn(x,n-1,safe) );
837 }
838
839 double safe_arctan_Dn(double x,int n,enum safe_err *safe)
840 {
841 return( n==0 ? safe_arctan_D0(x,safe) : rec_1_p_sqr_Dn(x,n-1,safe) );
842 }
843
844 double safe_cos_Dn(double x,int n,enum safe_err *safe)
845 {
846 switch( n%4 ) {
847 case 0:
848 return( cos(x) );
849 case 1:
850 return( -sin(x) );
851 case 2:
852 return( -cos(x) );
853 case 3:
854 return( sin(x) );
855 }
856 if( safe_print_errors ) {
857 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"cos_D: Unreachable point reached?!?");
858 }
859 *safe = safe_range_error;
860 return 0.0;
861 }
862
863 double safe_sin_Dn(double x,int n,enum safe_err *safe)
864 {
865 return( safe_cos_Dn(x,n+3,safe) );
866 }
867
868 #ifdef HAVE_ERF
869 double safe_erf_Dn(double x,int n,enum safe_err *safe)
870 {
871 return( (n==0)
872 ? (safe_erf_D0(x,safe))
873 : (safe_ERF_COEF*exp_msqr_Dn(x,n-1,safe)) );
874 }
875 #endif /* HAVE_ERF */
876
877 double safe_exp_Dn(double x,int n,enum safe_err *safe)
878 {
879 (void)n; /* stop gcc whine about unused parameter */
880 return( safe_exp_D0(x,safe) );
881 }
882
883 double safe_ln_Dn(double x,int n,enum safe_err *safe)
884 {
885 return( n==0 ? safe_ln_D0(x,safe) : safe_Dn_rec(x,n-1,safe) );
886 }
887
888 double safe_log10_Dn(double x,int n,enum safe_err *safe)
889 {
890 return( safe_LOG10_COEF * safe_ln_Dn(x,n,safe) );
891 }
892
893 double safe_sqr_Dn(double x,int n,enum safe_err *safe)
894 {
895 switch(n) {
896 case 0:
897 return( safe_sqr_D0(x,safe) );
898 case 1:
899 return( 2.0*x );
900 case 2:
901 return(2.0);
902 default:
903 return(0.0);
904 }
905 }
906
907 double safe_sqrt_Dn(double x,int n,enum safe_err *safe)
908 {
909 double a,b;
910
911 for( b = 1.0 , a = 1.5-n ; a < 1.0 ; ++a )
912 b *= a;
913 return( b * safe_sqrt_D0(x,safe) * safe_rec(safe_upower(x,n,safe),safe) );
914 }
915
916 static double tan_k_Dn(double u,int k,int n,enum safe_err *safe)
917 /* k >= 0 */
918 /**
919 *** nth derivative of tan^k evaluated at x, where u=tan(x)
920 **/
921 {
922 /* (tan^k)' = k * (tan^(k+1) + tan^(k-1)) */
923 if( n == 0 )
924 return( safe_upower(u,k,safe) );
925 else if( k == 0 )
926 return( 0.0 );
927 else
928 return( k * (tan_k_Dn(u,k+1,n-1,safe) + tan_k_Dn(u,k-1,n-1,safe)) );
929 }
930
931 double safe_tan_Dn(double x,int n,enum safe_err *safe)
932 {
933 return( tan_k_Dn(safe_tan_D0(x,safe),1,n,safe) );
934 }
935
936 static double cheap_pow(double x,double y,enum safe_err *safe)
937 /**
938 *** Computes x^y, x>0
939 **/
940 {
941 return( safe_exp_D0(safe_mul_D0(y,safe_ln_D0(x,safe),safe),safe) );
942 }
943
944 double safe_pow_D0(double x,double y,enum safe_err *safe)
945 {
946 if( x > 0.0 )
947 return( cheap_pow(x,y,safe) );
948 else if( x == 0.0 ) {
949 if( y <= 0.0 ) {
950 double bogus;
951 bogus = y < 0.0 ? BIGNUM : 0.0;
952 if( safe_print_errors ) {
953 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"pow_D0: %g raised to power %g undefined: returning %g",x,y,bogus);
954 }
955 *safe = safe_range_error;
956 return(bogus);
957 }
958 return( 0.0 );
959 }
960
961 /* x < 0.0 */
962 switch( safe_is_int(y,safe) ) {
963 case 0:
964 return( cheap_pow(-x,y,safe) );
965 case 1:
966 return( -cheap_pow(-x,y,safe) );
967 default: {
968 double bogus = cheap_pow(-x,y,safe);
969 if( safe_print_errors ) {
970 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"pow_D0: %g raised to power %g undefined: returning %g",x,y,bogus);
971
972 }
973 *safe = safe_range_error;
974 return(bogus);
975 }
976 }
977 }
978
979 double safe_div_D1(double x,double y,int wrt,enum safe_err *safe)
980 {
981 y = safe_rec(y,safe);
982 switch( wrt ) {
983 case 0:
984 return(y);
985 case 1:
986 return( -safe_mul_D0(safe_sqr_D0(y,safe),x,safe) );
987 default:
988 Asc_Panic(2, NULL, "'wrt' out of range!");
989
990 }
991 }
992
993 double safe_ipow_D1( double x,double y, int wrt,enum safe_err *safe)
994 {
995 (void)safe;
996
997 switch( wrt ) {
998 case 0:
999 return( safe_mul_D0(y,safe_ipow_D0(x,y-1,safe),safe) );
1000 case 1:
1001 return(0.0); /* d(x^i)/di where i is integer is 0.0 since we
1002 * assume integers are constant */
1003 default:
1004 Asc_Panic(2, NULL, "'wrt' out of range!");
1005
1006 }
1007 }
1008
1009 double safe_pow_D1(double x,double y,int wrt,enum safe_err *safe)
1010 {
1011 switch( wrt ) {
1012 case 0:
1013 return( safe_mul_D0(y,safe_pow_D0(x,y-1.0,safe),safe) );
1014 case 1:
1015 return( safe_mul_D0( safe_ln_D0(x,safe) , safe_pow_D0(x,y,safe) ,safe) );
1016 default:
1017 Asc_Panic(2, NULL, "'wrt' out of range!");
1018
1019 }
1020 }
1021
1022 double safe_div_D2(double x,double y,int wrt1,int wrt2,enum safe_err *safe)
1023 {
1024 (void)x; /* stop gcc whine about unused parameter */
1025 switch( wrt1+wrt2 ) {
1026 case 0:
1027 return(0.0);
1028 case 1:
1029 return( -safe_rec(safe_sqr_D0(y,safe),safe) );
1030 case 2:
1031 return( safe_rec(0.5*safe_mul_D0(y,safe_sqr_D0(y,safe),safe),safe) );
1032 default:
1033 Asc_Panic(2, NULL, "'wrt' out of range!");
1034
1035 }
1036 }
1037
1038 double safe_ipow_D2( double x,double y, int wrt1,int wrt2,enum safe_err *safe)
1039 {
1040 double lnx;
1041 int yint;
1042 yint = ascnint(y);
1043 switch( wrt1+wrt2 ) {
1044 case 0:
1045 return yint*(yint-1)*safe_ipow_D0(x,y-2,safe);
1046 case 1:
1047 lnx = safe_ln_D0(x,safe);
1048 return(
1049 safe_mul_D0(safe_rec(x,safe)+
1050 safe_sqr_D0(lnx,safe),safe_ipow_D0(x,y,safe),safe));
1051 case 2:
1052 lnx = safe_ln_D0(x,safe);
1053 return( safe_mul_D0(safe_sqr_D0(lnx,safe) ,
1054 safe_ipow_D0(x,y,safe),safe) );
1055 default:
1056 Asc_Panic(2, NULL, "'wrt' out of range!");
1057
1058 }
1059 }
1060
1061 double safe_pow_D2(double x,double y,int wrt1,int wrt2,enum safe_err *safe)
1062 {
1063 double lnx;
1064 switch( wrt1+wrt2 ) {
1065 case 0:
1066 return( safe_mul_D0( safe_mul_D0(y,y-1.0,safe),
1067 safe_pow_D0(x,y-2.0,safe),safe ) );
1068 case 1:
1069 lnx = safe_ln_D0(x,safe);
1070 return( safe_mul_D0(safe_rec(x,safe)+safe_sqr_D0(lnx,safe),
1071 safe_pow_D0(x,y,safe),safe) );
1072 case 2:
1073 lnx = safe_ln_D0(x,safe);
1074 return( safe_mul_D0(safe_sqr_D0(lnx,safe) ,
1075 safe_pow_D0(x,y,safe),safe) );
1076 default:
1077 Asc_Panic(2, NULL, "'wrt' out of range!");
1078
1079 }
1080 }
1081
1082 double safe_add_Dn(double x,double y,int nwrt0,int nwrt1,enum safe_err *safe)
1083 {
1084 (void)safe;
1085 switch( nwrt0+nwrt1 ) {
1086 case 0:
1087 return( x+y );
1088 case 1:
1089 return( 1.0 );
1090 default:
1091 return( 0.0 );
1092 }
1093 }
1094
1095 double safe_sub_Dn(double x,double y,int nwrt0,int nwrt1,enum safe_err *safe)
1096 {
1097 (void)safe;
1098 switch( nwrt0+nwrt1 ) {
1099 case 0:
1100 return( x-y );
1101 case 1:
1102 return( nwrt1==0 ? 1.0 : -1.0 );
1103 default:
1104 return( 0.0 );
1105 }
1106 }
1107
1108 double safe_mul_Dn(double x,double y,int nwrt0,int nwrt1,enum safe_err *safe)
1109 {
1110 (void)safe;
1111
1112 switch( nwrt0 ) {
1113 case 0:
1114 break;
1115 case 1:
1116 x = 1.0;
1117 break;
1118 default:
1119 return(0.0);
1120 }
1121
1122 switch( nwrt1 ) {
1123 case 0:
1124 return( safe_mul_D0(x,y,safe) );
1125 case 1:
1126 return(x);
1127 default:
1128 return(0.0);
1129 }
1130 }
1131
1132 double safe_div_Dn(double x,double y,int nwrt0,int nwrt1,enum safe_err *safe)
1133 {
1134 switch( nwrt0 ) {
1135 case 1:
1136 x = 1.0;
1137 case 0: /* FALL THROUGH */
1138 return( safe_mul_D0(x,safe_Dn_rec(y,nwrt1,safe),safe) );
1139 default:
1140 return(0.0);
1141 }
1142 }
1143
1144 double safe_pow_Dn(double x,double y,int nwrt0,int nwrt1,enum safe_err *safe)
1145 /**
1146 *** The n-th derivative wrt x of (lnx)^m * x^y (the m-th derivative wrt y
1147 *** of x^y) = x^(y-n) * P[n](lnx), where P[0](z) = z^m and
1148 *** P[n+1](z) = (y-n)*P[n](z) + P[n]'(z). (n==nwrt0,m==nwrt1)
1149 **/
1150 {
1151 double *poly;
1152 double lnx,lnx2; /* To be calculated later if necessary */
1153 int k,r;
1154 poly = alloc_poly(nwrt1);
1155
1156 /* mem_zero_byte_cast(poly,0,nwrt1*sizeof(double));*/
1157 ascbzero((void *)poly,nwrt1*sizeof(double));
1158 poly[nwrt1] = safe_pow_D0(x,y-(double)nwrt0,safe);
1159 for( k=0 ; k<nwrt0 ; ++k ) { /* calculate P[k+1] from P[k] */
1160 for( r=0 ; r<nwrt1 ; ++r )
1161 poly[r] = (y-k)*poly[r] + (r+1)*poly[r+1];
1162 poly[nwrt1] *= y-k;
1163 }
1164
1165 if( nwrt1 == 0 )
1166 return( poly[0] ); /* polynominal is just a constant */
1167
1168 lnx2 = safe_sqr_D0( lnx = safe_ln_D0(x,safe),safe );
1169 return( safe_poly(poly,nwrt1,lnx,lnx2,safe)+
1170 safe_poly(poly,nwrt1-1,lnx,lnx2,safe) );
1171 }

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