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

Properties

Name Value
svn:executable *

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