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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 91 - (hide annotations) (download) (as text)
Wed Dec 7 16:17:03 2005 UTC (18 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 aw0a 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 johnpye 91 error_reporter(ASC_PROG_WARNING,NULL,0,"Division by zero");
63 aw0a 1 case safe_complex_result:
64 johnpye 91 error_reporter(ASC_PROG_WARNING,NULL,0,"Complex result");
65 aw0a 1 break;
66     case safe_overflow:
67 johnpye 91 error_reporter(ASC_PROG_WARNING,NULL,0,"Overflow");
68 aw0a 1 break;
69     case safe_underflow:
70 johnpye 91 error_reporter(ASC_PROG_WARNING,NULL,0,"Underflow");
71 aw0a 1 break;
72     case safe_range_error:
73 johnpye 91 error_reporter(ASC_PROG_WARNING,NULL,0,"Function range error");
74 aw0a 1 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 johnpye 89 #ifdef HAVE_ERF
142 aw0a 1 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 johnpye 89 #endif
169 aw0a 1
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 johnpye 89 (void)safe;
188    
189 aw0a 1 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 johnpye 89 double y = 1.0;
204     (void)safe;
205    
206     for( ; n-- > 0 ; y = safe_mul_D0(y,x,safe) );
207     return(y);
208 aw0a 1 }
209    
210     double safe_factorial(unsigned n, enum safe_err *safe)
211     {
212     double x,y;
213 johnpye 89 (void)safe;
214    
215 aw0a 1 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 johnpye 91 error_reporter(ASC_USER_ERROR,NULL,0,"Divide by zero in 1/x expr: returning %g",bogus);
226 aw0a 1 }
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 johnpye 91 error_reporter(ASC_USER_ERROR,NULL,0,"Overflow in x^3 expr: returning %g.",bogus);
249 aw0a 1 }
250     } else {
251     bogus=0.0;
252     if( safe_print_errors ) {
253 johnpye 91 error_reporter(ASC_USER_ERROR,NULL,0,"Underflow in x^3 expr: returning %g.",bogus);
254 aw0a 1 }
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 johnpye 91 error_reporter(ASC_USER_ERROR,NULL,0,"inv_erf undefined at %g: returning %g",x,bogus);
279 aw0a 1 }
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 johnpye 91 error_reporter(ASC_USER_ERROR,NULL,0,"Found negative slope of %g on erf_inv.",dy);
295 aw0a 1 }
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 johnpye 89
306 aw0a 1 if( x > (DBL_MAX > 1.0e308 ? 709.196209 : log(DBL_MAX)) ) {
307     double bogus = BIGNUM;
308     if( safe_print_errors ) {
309 johnpye 91 error_reporter(ASC_USER_ERROR,NULL,0,"Argument %g too large in lnm_inv: returning %g.",x,bogus);
310 aw0a 1 }
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 johnpye 89 (void)safe;
323 aw0a 1
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 johnpye 91 error_reporter(ASC_USER_ERROR,NULL,0,"arcsin_D0: arcsin undefined at %g: Returning %g.",x,bogus);
389 aw0a 1 }
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 johnpye 91 error_reporter(ASC_USER_ERROR,NULL,0,"arccos_D0: arccos undefined at %g: returning %g.",x,bogus);
408 aw0a 1 }
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 johnpye 91 error_reporter(ASC_USER_ERROR,NULL,0,"arccosh_D0: undefined at %g: returning %g.",x,bogus);
421 aw0a 1 }
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 johnpye 91 error_reporter(ASC_USER_ERROR,NULL,0,"arctanh_D0: undefined at %g: returning %g.",x,bogus);
441 aw0a 1 }
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 johnpye 91 error_reporter(ASC_USER_ERROR,NULL,0,"exp_D0: Argument %g too large: returning %g.",x,bogus);
462 aw0a 1 }
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 johnpye 91 error_reporter(ASC_USER_ERROR,NULL,0,"ln_D0: undefined at %g: Returning %g.",x,bogus);
475 aw0a 1 }
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 johnpye 91 error_reporter(ASC_USER_ERROR,NULL,0,"sqrt_D0: undefined at %g: returning %g.",x,bogus);
506 aw0a 1 }
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 johnpye 91 error_reporter(ASC_USER_ERROR,NULL,0,"tan_D2: Infinite at %g: returning %g.",x,bogus);
702 aw0a 1 }
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 johnpye 91 error_reporter(ASC_USER_ERROR,NULL,0,"arccosh_D2: Undefined at %g: returning %g",x,bogus);
739 aw0a 1 }
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 johnpye 91 error_reporter(ASC_USER_ERROR,NULL,0,"arctanh_D2: undefined at %g: returning %g.",x,bogus);
758 aw0a 1 }
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 johnpye 89 (void)safe;
808 aw0a 1 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 johnpye 91 error_reporter(ASC_USER_ERROR,NULL,0,"cbrt_D2: undefined at %g: returning %g.",x,bogus);
817 aw0a 1 }
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 johnpye 91 error_reporter(ASC_USER_ERROR,NULL,0,"cos_D: Unreachable point reached?!?");
860 aw0a 1 }
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 johnpye 91 error_reporter(ASC_USER_ERROR,NULL,0,"pow_D0: %g raised to power %g undefined: returning %g",x,y,bogus);
956 aw0a 1 }
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 johnpye 91 error_reporter(ASC_USER_ERROR,NULL,0,"pow_D0: %g raised to power %g undefined: returning %g",x,y,bogus);
973 johnpye 62
974 aw0a 1 }
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 johnpye 89 (void)safe;
998    
999 aw0a 1 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 johnpye 89 (void)safe;
1113    
1114 aw0a 1 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