/[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 123 - (hide annotations) (download) (as text)
Mon Dec 19 06:59:40 2005 UTC (14 years, 7 months ago) by johnpye
File MIME type: text/x-csrc
File size: 27844 byte(s)
Changing 'log' function to 'log10'
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 jds 114 static double *alloc_poly(int order)
104 aw0a 1 /**
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 johnpye 89 #ifdef HAVE_ERF
141 aw0a 1 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 johnpye 89 #endif
168 aw0a 1
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 johnpye 89 (void)safe;
187    
188 aw0a 1 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 johnpye 89 double y = 1.0;
203     (void)safe;
204    
205     for( ; n-- > 0 ; y = safe_mul_D0(y,x,safe) );
206     return(y);
207 aw0a 1 }
208    
209     double safe_factorial(unsigned n, enum safe_err *safe)
210     {
211     double x,y;
212 johnpye 89 (void)safe;
213    
214 aw0a 1 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 johnpye 91 error_reporter(ASC_USER_ERROR,NULL,0,"Divide by zero in 1/x expr: returning %g",bogus);
225 aw0a 1 }
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 johnpye 91 error_reporter(ASC_USER_ERROR,NULL,0,"Overflow in x^3 expr: returning %g.",bogus);
248 aw0a 1 }
249     } else {
250     bogus=0.0;
251     if( safe_print_errors ) {
252 johnpye 91 error_reporter(ASC_USER_ERROR,NULL,0,"Underflow in x^3 expr: returning %g.",bogus);
253 aw0a 1 }
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 johnpye 91 error_reporter(ASC_USER_ERROR,NULL,0,"inv_erf undefined at %g: returning %g",x,bogus);
278 aw0a 1 }
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 johnpye 91 error_reporter(ASC_USER_ERROR,NULL,0,"Found negative slope of %g on erf_inv.",dy);
294 aw0a 1 }
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 johnpye 89
305 aw0a 1 if( x > (DBL_MAX > 1.0e308 ? 709.196209 : log(DBL_MAX)) ) {
306     double bogus = BIGNUM;
307     if( safe_print_errors ) {
308 johnpye 91 error_reporter(ASC_USER_ERROR,NULL,0,"Argument %g too large in lnm_inv: returning %g.",x,bogus);
309 aw0a 1 }
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 johnpye 89 (void)safe;
322 aw0a 1
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 johnpye 91 error_reporter(ASC_USER_ERROR,NULL,0,"arcsin_D0: arcsin undefined at %g: Returning %g.",x,bogus);
388 aw0a 1 }
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 johnpye 91 error_reporter(ASC_USER_ERROR,NULL,0,"arccos_D0: arccos undefined at %g: returning %g.",x,bogus);
407 aw0a 1 }
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 johnpye 91 error_reporter(ASC_USER_ERROR,NULL,0,"arccosh_D0: undefined at %g: returning %g.",x,bogus);
420 aw0a 1 }
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 johnpye 91 error_reporter(ASC_USER_ERROR,NULL,0,"arctanh_D0: undefined at %g: returning %g.",x,bogus);
440 aw0a 1 }
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 johnpye 91 error_reporter(ASC_USER_ERROR,NULL,0,"exp_D0: Argument %g too large: returning %g.",x,bogus);
461 aw0a 1 }
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 johnpye 91 error_reporter(ASC_USER_ERROR,NULL,0,"ln_D0: undefined at %g: Returning %g.",x,bogus);
474 aw0a 1 }
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 johnpye 123 double safe_log10_D0(double x,enum safe_err *safe)
488 aw0a 1 {
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 johnpye 91 error_reporter(ASC_USER_ERROR,NULL,0,"sqrt_D0: undefined at %g: returning %g.",x,bogus);
505 aw0a 1 }
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 johnpye 123 double safe_log10_D1(double x,enum safe_err *safe)
631 aw0a 1 {
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 johnpye 91 error_reporter(ASC_USER_ERROR,NULL,0,"tan_D2: Infinite at %g: returning %g.",x,bogus);
701 aw0a 1 }
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 johnpye 91 error_reporter(ASC_USER_ERROR,NULL,0,"arccosh_D2: Undefined at %g: returning %g",x,bogus);
738 aw0a 1 }
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 johnpye 91 error_reporter(ASC_USER_ERROR,NULL,0,"arctanh_D2: undefined at %g: returning %g.",x,bogus);
757 aw0a 1 }
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 johnpye 123 double safe_log10_D2(double x,enum safe_err *safe)
789 aw0a 1 {
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 johnpye 89 (void)safe;
807 aw0a 1 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 johnpye 91 error_reporter(ASC_USER_ERROR,NULL,0,"cbrt_D2: undefined at %g: returning %g.",x,bogus);
816 aw0a 1 }
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 johnpye 91 error_reporter(ASC_USER_ERROR,NULL,0,"cos_D: Unreachable point reached?!?");
859 aw0a 1 }
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 johnpye 123 double safe_log10_Dn(double x,int n,enum safe_err *safe)
890 aw0a 1 {
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 johnpye 91 error_reporter(ASC_USER_ERROR,NULL,0,"pow_D0: %g raised to power %g undefined: returning %g",x,y,bogus);
955 aw0a 1 }
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 johnpye 91 error_reporter(ASC_USER_ERROR,NULL,0,"pow_D0: %g raised to power %g undefined: returning %g",x,y,bogus);
972 johnpye 62
973 aw0a 1 }
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 johnpye 89 (void)safe;
997    
998 aw0a 1 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 johnpye 89 (void)safe;
1112    
1113 aw0a 1 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