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

Properties

Name Value
svn:executable *

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