/[ascend]/trunk/ascend4/solver/calc.c
ViewVC logotype

Annotation of /trunk/ascend4/solver/calc.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (hide annotations) (download) (as text)
Fri Oct 29 20:54:12 2004 UTC (17 years, 8 months ago) by aw0a
File MIME type: text/x-csrc
File size: 21615 byte(s)
Setting up web subdirectory in repository
1 aw0a 1 /*
2     * SLV: Ascend Numeric Solver
3     * by Karl Michael Westerberg
4     * Created: 2/6/90
5     * Version: $Revision: 1.8 $
6     * Version control file: $RCSfile: calc.c,v $
7     * Date last modified: $Date: 1998/02/05 15:59:18 $
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     *
16     * The SLV solver is free software; you can redistribute
17     * it and/or modify it under the terms of the GNU General Public License as
18     * published by the Free Software Foundation; either version 2 of the
19     * License, or (at your option) any later version.
20     *
21     * The SLV solver is distributed in hope that it will be
22     * useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
23     * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
24     * General Public License for more details.
25     *
26     * You should have received a copy of the GNU General Public License
27     * along with the program; if not, write to the Free Software Foundation,
28     * Inc., 675 Mass Ave, Cambridge, MA 02139 USA. Check the file named
29     * COPYING. COPYING is found in ../compiler.
30     *
31     */
32    
33     #include <math.h>
34     #include <stdarg.h>
35     #include "utilities/ascConfig.h"
36     #include "utilities/ascMalloc.h"
37     #include "utilities/mem.h"
38     #include "compiler/fractions.h"
39     #include "utilities/ascPanic.h"
40     #include "compiler/functype.h"
41     #include "compiler/dimen.h"
42     #include "compiler/func.h"
43     #include "solver/calc.h"
44    
45     #define BIGNUM 1.0e8
46    
47     int32 calc_ok = TRUE;
48     boolean calc_print_errors = TRUE;
49    
50     static real64 rec_1_p_sqr_Dn(x,n)
51     real64 x;
52     int n;
53     /**
54     *** Calculates n-th derivative of 1/(1+x^2).
55     **/
56     {
57     /* g[n](x) = -1/(1+x^2) * (n*(n-1)*g[n-2](x) + 2*n*x*g[n-1](x))
58     where g[n] := n-th derivative of 1/(1+x^2). */
59     int k;
60     real64 prev[3],val; /* prev[j] := g[k-j](x), val = g[0](x) */
61    
62     prev[2] = val = 1.0/(1.0 + calc_sqr_D0(x));
63     if( n==0 )
64     return(val);
65     prev[1] = -2.0*x*val*val; /* first derivative */
66     for( k=2 ; k<=n ; ++k ) {
67     prev[0] = -val*(calc_mul_D0((double)(k*(k-1)),prev[2]) +
68     calc_mul_D0(2*k*x,prev[1]));
69     prev[2] = prev[1];
70     prev[1] = prev[0];
71     }
72    
73     return(prev[1]);
74     }
75    
76     static real64 *alloc_poly(order)
77     int order;
78     /**
79     *** Allocates a polynominal of given order and returns it. The
80     *** polynominal need not be freed, but this function should not be
81     *** called again until the old polynominal is not needed anymore.
82     **/
83     {
84     static real64 *poly = NULL;
85     static int poly_cap = 0;
86    
87     if( order + 1 > poly_cap ) {
88     poly_cap = order+1;
89     if( poly != NULL )
90     ascfree( (POINTER)poly );
91     poly = (real64 *)ascmalloc( poly_cap * sizeof(real64) );
92     }
93     return(poly);
94     }
95    
96     static real64 calc_poly(poly,order,x,x2)
97     real64 *poly;
98     int order;
99     real64 x,x2; /* x2 = x^2 */
100     /**
101     *** Calculates the value of the given polynomial, where only the
102     *** (order%1 ? odd : even) degree terms are used.
103     **/
104     {
105     real64 val;
106    
107     for( val=poly[order] ; (order -= 2) >= 0 ; )
108     val = calc_mul_D0(val,x2) + poly[order];
109     return( order==-1 ? calc_mul_D0(val,x) : val );
110     }
111    
112     static real64 exp_msqr_Dn(x,n)
113     real64 x;
114     int n;
115     /**
116     *** Computes n-th derivative of exp(-x^2).
117     **/
118     {
119     /**
120     *** n-th derivative of exp(-x^2) = f[n](x)*exp(-x^2), where f[n] is an
121     *** n-th degree polynominal of definite parity satisfying f[0](x) = 1
122     *** & f[n+1](x) = -2x*f[n](x) + f[n]'(x).
123     **/
124    
125     real64 x2 = calc_sqr_D0(x);
126     int k,r;
127     real64 *poly;
128     poly = alloc_poly(n);
129    
130     poly[0] = exp(-x2);
131     for( k=0 ; k<n ; ++k ) { /* Calculate f[k+1] from f[k] */
132     poly[k+1] = 0.0;
133     for( r=k ; r >= 1 ; r -= 2 )
134     poly[r-1] = (double)r * poly[r];
135     for( r=k ; r >= 0 ; r -= 2 )
136     poly[r+1] += -2.0*poly[r];
137     }
138     return( calc_poly(poly,n,x,x2) );
139     }
140    
141     static real64 sqrt_rec_1_m_sqr_Dn(x,n)
142     real64 x;
143     int n;
144     /**
145     *** Calculates n-th derivative of (1-x^2)^-.5
146     **/
147     {
148     /**
149     *** n-th derivative of (1-x^2)^-.5 = f[n](x) * (1-x^2)^-(n+.5), where
150     *** f[n] is an n-degree polynominal of definite parity satisfying
151     *** f[0] = 1 and f[n+1](x) = f[n]'(x)*(1+x^2) + (2n+1)*x*f[n](x).
152     **/
153    
154     int k,r;
155     real64 x2;
156     real64 *poly;
157     x2 = calc_sqr_D0(x);
158     poly = alloc_poly(n);
159    
160     poly[0] = calc_rec(calc_upower(calc_sqrt_D0(1.0-x2),2*n+1));
161     for( k=0 ; k<n ; ++k ) { /* Calculate f[k+1] from f[k] */
162     poly[k+1] = 0.0;
163     for( r=k ; r >= 1 ; r -= 2 )
164     poly[r-1] = (double)r * poly[r];
165     for( r=k ; r >= 0 ; r -= 2 )
166     poly[r+1] += (double)(r+2*k+1) * poly[r];
167     }
168     return( calc_poly(poly,n,x,x2) );
169     }
170    
171    
172     real64 calc_upower(real64 x, unsigned n)
173     {
174     double y = 1.0;
175     for( ; n-- > 0 ; y = calc_mul_D0(y,x) )
176     ;
177     return(y);
178     }
179    
180     real64 calc_factorial(unsigned n)
181     {
182     double x,y;
183     for( x = (double)n , y = 1.0 ; x >= 0.5 ; y = calc_mul_D0(y,x--) )
184     ;
185     return(y);
186     }
187    
188     real64 calc_rec(real64 x)
189     {
190     if( x == 0.0 ) {
191     real64 bogus = BIGNUM;
192     if( calc_print_errors ) {
193     FPRINTF(stderr,"ERROR: (calc) calc_rec\n");
194     FPRINTF(stderr," Divide by zero undefined.\n");
195     FPRINTF(stderr," Returning %g.\n",bogus);
196     }
197     calc_ok = FALSE;
198     return( bogus );
199     }
200     return( 1.0/x );
201     }
202    
203     #ifndef INV_CUBRTHUGE
204     #define INV_CUBRTHUGE 1.7718548704178434e-103
205     /* smallest cubeable number in an 8 byte ieee double */
206     #endif
207     #ifndef CUBRTHUGE
208     #define CUBRTHUGE 5.6438030941223618e102
209     /* largest cubeable number in an 8 byte ieee double */
210     #endif
211     real64 calc_cube(x)
212     register real64 x;
213     {
214     if( fabs(x) > CUBRTHUGE || fabs(x) < INV_CUBRTHUGE ) {
215     real64 bogus;
216     if ( fabs(x) > CUBRTHUGE) {
217     bogus=BIGNUM;
218     if( calc_print_errors ) {
219     FPRINTF(stderr,"ERROR: (calc) calc_cube\n");
220     FPRINTF(stderr," Overflow calculation requested \n");
221     FPRINTF(stderr," Returning %g.\n",bogus);
222     }
223     } else {
224     bogus=0.0;
225     if( calc_print_errors ) {
226     FPRINTF(stderr,"ERROR: (calc) calc_cube\n");
227     FPRINTF(stderr," Underflow calculation requested \n");
228     FPRINTF(stderr," Returning %g.\n",bogus);
229     }
230     }
231     calc_ok = FALSE;
232     return( bogus );
233     }
234     return ( x*x*x );
235     }
236    
237     real64 calc_Dn_rec(x,n)
238     real64 x;
239     int n;
240     {
241     return( -calc_upower(-calc_rec(x),n+1) * calc_factorial(n) );
242     }
243    
244     #ifdef HAVE_ERF
245     real64 calc_erf_inv(x)
246     real64 x;
247     {
248     #define CONV (1e-7)
249     real64 y,dy,sign;
250    
251     sign = (x<0.0) ? -1.0 : 1.0;
252     x *= sign;
253     if( x >= 1.0 ) {
254     real64 bogus = sign*BIGNUM;
255     if( calc_print_errors ) {
256     FPRINTF(stderr,"ERROR: (calc) calc_erf_inv\n");
257     FPRINTF(stderr," Inverse erf is undefined at %g.\n",x);
258     FPRINTF(stderr," Returning %g.\n",bogus);
259     }
260     calc_ok = FALSE;
261     return(bogus);
262     }
263    
264     y = 0.0;
265     do {
266     dy = (x-calc_erf_D0(y))*calc_exp_D0(calc_mul_D0(y,y))/calc_ERF_COEF;
267     y += dy;
268     /**
269     *** Since erf is concave and x >= erf(y0), dy should
270     *** always be positive (?).
271     **/
272     if( dy < 0.0 && calc_print_errors ) {
273     FPRINTF(stderr,"ERROR: (calc) calc_erf_inv\n");
274     FPRINTF(stderr," Found negative slope of %g on erf.\n",dy);
275     }
276     } while( calc_ok && dy > CONV );
277     return( y * sign );
278     #undef CONV
279     }
280     #endif /* HAVE_ERF */
281    
282     real64 calc_lnm_inv(x)
283     real64 x;
284     {
285     register double eps=FuncGetLnmEpsilon();
286     if( x > (MAXDOUBLE > 1.0e308 ? 709.196209 : log(MAXDOUBLE)) ) {
287     real64 bogus = BIGNUM;
288     if( calc_print_errors ) {
289     FPRINTF(stderr,"ERROR: (calc) calc_lnm_inv\n");
290     FPRINTF(stderr," Argument %g too large.\n",x);
291     FPRINTF(stderr," Returning %g.\n",bogus);
292     }
293     calc_ok = FALSE;
294     return(bogus);
295     }
296     /* since lnm(eps)= log(eps), use cheaper one. eps known >0 */
297     return( (x >= log(eps)) ? calc_exp_D0(x) : eps*(x+1.0-log(eps)) );
298     }
299    
300     int calc_is_int(x)
301     real64 x;
302     {
303     #define TOLERANCE (1e-4)
304     unsigned n;
305    
306     if( x<0.0 )
307     x = -x;
308    
309     /**
310     *** Theorem: Define P(n) :<==> x-TOLERANCE <= n <= x+TOLERANCE. Then
311     *** P(n) for some integer n <==> P(floor(x+TOLERANCE))
312     *** <==> x-TOLERANCE <= floor(x+TOLERANCE).
313     ***
314     *** Proof: Since floor(x+TOLERANCE) <= x+TOLERANCE, the second
315     *** equivalence holds. The backward direction of the first
316     *** equivalence is a tautology. If P(n) holds for some n, then
317     *** n <= floor(x+TOLERANCE) <= x+TOLERANCE, the second inequality
318     *** follows since floor(z) <= z, and the first inequality follows
319     *** since floor(z) is the largest integer with that property.
320     **/
321     if( x-TOLERANCE <= (double)(n = (unsigned)(x+TOLERANCE)) )
322     return( n&01 );
323     else
324     return( -1 );
325     #undef TOLERANCE
326     }
327    
328    
329     real64 calc_tan_D0(x)
330     real64 x;
331     {
332     return( calc_div_D0( sin(x) , cos(x) ) );
333     }
334    
335     real64 calc_arcsin_D0(x)
336     real64 x;
337     {
338     if( x < -1.0 || 1.0 < x ) {
339     real64 bogus = 0.0;
340     if( calc_print_errors ) {
341     FPRINTF(stderr,"ERROR: (calc) calc_arcsin_D0\n");
342     FPRINTF(stderr," Function arcsin is undefined at %g.\n",x);
343     FPRINTF(stderr," Returning %g.\n",bogus);
344     }
345     calc_ok = FALSE;
346     return(bogus);
347     }
348     return( asin(x) );
349     }
350    
351     real64 calc_arccos_D0(x)
352     real64 x;
353     {
354     if( x < -1.0 || 1.0 < x ) {
355     real64 bogus = 0.0;
356     if( calc_print_errors ) {
357     FPRINTF(stderr,"ERROR: (calc) calc_arccos_D0\n");
358     FPRINTF(stderr," Function arccos is undefined at %g.\n",x);
359     FPRINTF(stderr," Returning %g.\n",bogus);
360     }
361     calc_ok = FALSE;
362     return(bogus);
363     }
364     return( acos(x) );
365     }
366    
367     real64 calc_arccosh_D0(x)
368     real64 x;
369     {
370     if( x < 1.0 ) {
371     real64 bogus = 0.0;
372     if( calc_print_errors ) {
373     FPRINTF(stderr,"ERROR: (calc) calc_arccosh_D0\n");
374     FPRINTF(stderr," Function arccosh is undefined at %g.\n",x);
375     FPRINTF(stderr," Returning %g.\n",bogus);
376     }
377     calc_ok = FALSE;
378     return(bogus);
379     }
380    
381     return( arccosh(x) );
382     }
383    
384     real64 calc_arctanh_D0(x)
385     real64 x;
386     {
387     if( x < -1.0 || 1.0 < x ) {
388     real64 bogus = 0.0;
389     if( calc_print_errors ) {
390     FPRINTF(stderr,"ERROR: (calc) calc_arctanh_D0\n");
391     FPRINTF(stderr," Function arctanh is undefined at %g.\n",x);
392     FPRINTF(stderr," Returning %g.\n",bogus);
393     }
394     calc_ok = FALSE;
395     return(bogus);
396     }
397     return( arctanh(x) );
398     }
399    
400     real64 calc_exp_D0(x)
401     real64 x;
402     {
403     if( x > (MAXDOUBLE > 1.0e308 ? 709.196209 : log(MAXDOUBLE)) ) {
404     real64 bogus = BIGNUM;
405     if( calc_print_errors ) {
406     FPRINTF(stderr,"ERROR: (calc) calc_exp_D0\n");
407     FPRINTF(stderr," Argument %g too large.\n",x);
408     FPRINTF(stderr," Returning %g.\n",bogus);
409     }
410     calc_ok = FALSE;
411     return(bogus);
412     }
413     return( exp(x) );
414     }
415    
416     real64 calc_ln_D0(x)
417     real64 x;
418     {
419     if( x <= 0.0 ) {
420     real64 bogus = -BIGNUM;
421     if( calc_print_errors ) {
422     FPRINTF(stderr,"ERROR: (calc) calc_ln_D0\n");
423     FPRINTF(stderr," Natural log undefined at %g.\n",x);
424     FPRINTF(stderr," Returning %g.\n",bogus);
425     }
426     calc_ok = FALSE;
427     return(bogus);
428     }
429     return( log(x) );
430     }
431    
432     real64 calc_log_D0(x)
433     real64 x;
434     {
435     return( calc_LOG10_COEF * calc_ln_D0(x) );
436     }
437    
438     real64 calc_sqrt_D0(x)
439     real64 x;
440     {
441     if( x < 0.0 ) {
442     real64 bogus;
443     bogus = -sqrt(-x);
444     if( calc_print_errors ) {
445     FPRINTF(stderr,"ERROR: (calc) calc_sqrt_D0\n");
446     FPRINTF(stderr," Square root undefined at %g.\n",x);
447     FPRINTF(stderr," Returning %g.\n",bogus);
448     }
449     calc_ok = FALSE;
450     return(bogus);
451     }
452     return(sqrt(x));
453     }
454    
455     real64 calc_tan_D1(x)
456     register real64 x;
457     {
458     x=cos(x);
459     return( calc_rec(calc_sqr_D0(x)) );
460     }
461    
462     real64 calc_arcsin_D1(x)
463     real64 x;
464     {
465     return( calc_rec(calc_sqrt_D0(1.0 - calc_sqr_D0(x))) );
466     }
467    
468     real64 calc_arccos_D1(x)
469     real64 x;
470     {
471     return( -calc_arcsin_D1(x) );
472     }
473    
474     real64 calc_arccosh_D1(x)
475     real64 x;
476     {
477     return( calc_rec(calc_sqrt_D0(calc_sqr_D0(x) - 1.0)) );
478     }
479    
480     real64 calc_arctanh_D1(x)
481     real64 x;
482     {
483     return( calc_rec(1.0 - calc_sqr_D0(x)) );
484     }
485    
486     #ifdef HAVE_ERF
487     real64 calc_erf_D1(x)
488     real64 x;
489     {
490     return( calc_ERF_COEF * calc_exp_D0(-calc_sqr_D0(x)) );
491     }
492     #endif /* HAVE_ERF */
493    
494     real64 calc_log_D1(x)
495     real64 x;
496     {
497     return( calc_LOG10_COEF * calc_ln_D1(x) );
498     }
499    
500     real64 calc_sqrt_D1(x)
501     real64 x;
502     {
503     return( 0.5 * calc_rec(calc_sqrt_D0(x)) );
504     }
505    
506     real64 calc_cbrt_D1(x)
507     real64 x;
508     {
509     return( 0.3333333333333333*calc_rec(calc_sqr_D0(calc_cbrt_D0(x))) );
510     }
511     /* KHACK */
512     real64 calc_fabs_D1(x)
513     real64 x;
514     {
515     (void)x;
516     return( 1.0 );
517     }
518    
519     real64 calc_tan_D2(x)
520     real64 x;
521     {
522     if( fabs(cos(x)) == 0.0 ) {
523     real64 bogus = BIGNUM;
524     if( calc_print_errors ) {
525     FPRINTF(stderr,"ERROR: (calc) calc_tan_D2\n");
526     FPRINTF(stderr," Tan derivative infinite at %g.\n",x);
527     FPRINTF(stderr," Returning %g.\n",bogus);
528     }
529     calc_ok = FALSE;
530     return(bogus);
531     }
532     return dtan2(x);
533     }
534    
535     real64 calc_arcsin_D2(x)
536     real64 x;
537     {
538     register double t;
539     t = calc_rec(1.0 - calc_sqr_D0(x));
540     return( calc_mul_D0( calc_mul_D0(x,t) , calc_mul_D0(t,t) ) );
541     }
542    
543     real64 calc_arccos_D2(x)
544     real64 x;
545     {
546     return( -calc_arcsin_D2(x) );
547     }
548    
549     real64 calc_arccosh_D2(x)
550     real64 x;
551     {
552     if( fabs(x) <= 1.0 ) {
553     real64 bogus = -BIGNUM;
554     if( calc_print_errors ) {
555     FPRINTF(stderr,"ERROR: (calc) calc_arccosh_D2\n");
556     FPRINTF(stderr," Arccosh undefined at %g.\n",x);
557     FPRINTF(stderr," Returning %g.\n",bogus);
558     }
559     calc_ok = FALSE;
560     return(bogus);
561     }
562     return darccosh2(x);
563     }
564    
565     real64 calc_arctanh_D2(x)
566     real64 x;
567     {
568     if( fabs(x) == 1.0 ) {
569     real64 bogus = BIGNUM;
570     if( calc_print_errors ) {
571     FPRINTF(stderr,"ERROR: (calc) calc_arctanh_D2\n");
572     FPRINTF(stderr," Arctanh undefined at %g.\n",x);
573     FPRINTF(stderr," Returning %g.\n",bogus);
574     }
575     calc_ok = FALSE;
576     return(bogus);
577     }
578     return( darctanh2(x) );
579     }
580    
581     #ifdef HAVE_ERF
582     real64 calc_erf_D2(x)
583     real64 x;
584     {
585     return( -ldexp(calc_ERF_COEF * calc_mul_D0(x,calc_exp_D0(-calc_sqr_D0(x))),
586     1) );
587     }
588     #endif /* HAVE_ERF */
589    
590     real64 calc_ln_D2(x)
591     real64 x;
592     {
593     return( -calc_rec(calc_sqr_D0(x)) );
594     }
595    
596     real64 calc_log_D2(x)
597     real64 x;
598     {
599     return( calc_LOG10_COEF * calc_ln_D2(x) );
600     }
601    
602     real64 calc_sqrt_D2(x)
603     real64 x;
604     {
605     return( -ldexp( calc_rec( calc_mul_D0(x,calc_sqrt_D0(x))), -2) );
606     }
607    
608     real64 calc_cbrt_D2(x)
609     real64 x;
610     {
611     if( fabs(x) == 0.0 ) {
612     real64 bogus = BIGNUM;
613     if( calc_print_errors ) {
614     FPRINTF(stderr,"ERROR: (calc) calc_cbrt_D2\n");
615     FPRINTF(stderr," Cbrt undefined at %g.\n",x);
616     FPRINTF(stderr," Returning %g.\n",bogus);
617     }
618     calc_ok = FALSE;
619     return(bogus);
620     }
621     return dcbrt2(x);
622     }
623     /* KHACK */
624     real64 calc_fabs_D2(x)
625     real64 x;
626     {
627     (void)x;
628     return( 0.0 );
629     }
630    
631     real64 calc_arcsin_Dn(x,n)
632     real64 x;
633     int n;
634     {
635     return( n==0 ? calc_arcsin_D0(x) : sqrt_rec_1_m_sqr_Dn(x,n-1) );
636     }
637    
638     real64 calc_arccos_Dn(x,n)
639     real64 x;
640     int n;
641     {
642     return( n==0 ? calc_arccos_D0(x) : -sqrt_rec_1_m_sqr_Dn(x,n-1) );
643     }
644    
645     real64 calc_arctan_Dn(x,n)
646     real64 x;
647     int n;
648     {
649     return( n==0 ? calc_arctan_D0(x) : rec_1_p_sqr_Dn(x,n-1) );
650     }
651    
652     real64 calc_cos_Dn(x,n)
653     real64 x;
654     int n;
655     {
656     switch( n%4 ) {
657     case 0:
658     return( cos(x) );
659     case 1:
660     return( -sin(x) );
661     case 2:
662     return( -cos(x) );
663     case 3:
664     return( sin(x) );
665     }
666     if( calc_print_errors ) {
667     FPRINTF(stderr,"ERROR: (calc) calc_cos_Dn\n");
668     FPRINTF(stderr," Unreachable point reached.\n");
669     }
670     calc_ok = FALSE;
671     return 0.0;
672     }
673    
674     real64 calc_sin_Dn(x,n)
675     real64 x;
676     int n;
677     {
678     return( calc_cos_Dn(x,n+3) );
679     }
680    
681     #ifdef HAVE_ERF
682     real64 calc_erf_Dn(x,n)
683     real64 x;
684     int n;
685     {
686     return( n==0 ? calc_erf_D0(x) : calc_ERF_COEF*exp_msqr_Dn(x,n-1) );
687     }
688     #endif /* HAVE_ERF */
689    
690     real64 calc_exp_Dn(x,n)
691     real64 x;
692     int n;
693     {
694     (void)n;
695     return( calc_exp_D0(x) );
696     }
697    
698     real64 calc_ln_Dn(x,n)
699     real64 x;
700     int n;
701     {
702     return( n==0 ? calc_ln_D0(x) : calc_Dn_rec(x,n-1) );
703     }
704    
705     real64 calc_log_Dn(x,n)
706     real64 x;
707     int n;
708     {
709     return( calc_LOG10_COEF * calc_ln_Dn(x,n) );
710     }
711    
712     real64 calc_sqr_Dn(x,n)
713     real64 x;
714     int n;
715     {
716     switch(n) {
717     case 0:
718     return( calc_sqr_D0(x) );
719     case 1:
720     return( 2.0*x );
721     case 2:
722     return(2.0);
723     default:
724     return(0.0);
725     }
726     }
727    
728     real64 calc_sqrt_Dn(x,n)
729     real64 x;
730     int n;
731     {
732     double a,b;
733    
734     for( b = 1.0 , a = 1.5-n ; a < 1.0 ; ++a )
735     b *= a;
736     return( b * calc_sqrt_D0(x) * calc_rec(calc_upower(x,n)) );
737     }
738    
739     /*
740     * nth derivative of tan^k evaluated at x, where u=tan(x)
741     * k >= 0
742     */
743     static double tan_k_Dn(double u, int k, int n)
744     {
745     /* (tan^k)' = k * (tan^(k+1) + tan^(k-1)) */
746     if( n == 0 )
747     return( calc_upower(u,k) );
748     else if( k == 0 )
749     return( 0.0 );
750     else
751     return( k * (tan_k_Dn(u,k+1,n-1) + tan_k_Dn(u,k-1,n-1)) );
752     }
753    
754     extern real64 calc_tan_Dn(real64 x, int n)
755     {
756     return( tan_k_Dn(calc_tan_D0(x),1,n) );
757     }
758    
759     /*
760     * Computes x^y, x>0
761     */
762     static double cheap_pow(double x, double y)
763     {
764     return( calc_exp_D0(calc_mul_D0(y,calc_ln_D0(x))) );
765     }
766    
767     extern real64 calc_pow_D0(real64 x, real64 y)
768     {
769     if( x > 0.0 ) {
770     return( cheap_pow(x,y) );
771     }
772     else if( x == 0.0 ) {
773     if( y <= 0.0 ) {
774     real64 bogus;
775     bogus = y < 0.0 ? BIGNUM : 0.0;
776     if( calc_print_errors ) {
777     FPRINTF(stderr,"ERROR: (calc) calc_pow_D0\n");
778     FPRINTF(stderr," %g raised to %g power undefined.\n",x,y);
779     FPRINTF(stderr," Returning %g.\n",bogus);
780     }
781     calc_ok = FALSE;
782     return(bogus);
783     }
784     return( 0.0 );
785     }
786    
787     /* x < 0.0 */
788     switch( calc_is_int(y) ) {
789     case 0:
790     return( cheap_pow(-x,y) );
791     case 1:
792     return( -cheap_pow(-x,y) );
793     default: {
794     real64 bogus = cheap_pow(-x,y);
795     if( calc_print_errors ) {
796     FPRINTF(stderr,"ERROR: (calc) calc_pow_D0\n");
797     FPRINTF(stderr," %g raised to %g power undefined.\n",x,y);
798     FPRINTF(stderr," Returning %g.\n",bogus);
799     }
800     calc_ok = FALSE;
801     return(bogus);
802     }
803     }
804     }
805    
806     extern real64 calc_div_D1(real64 x, real64 y, int wrt)
807     {
808     y = calc_rec(y);
809     switch( wrt ) {
810     case 0:
811     return(y);
812     case 1:
813     return( -calc_mul_D0(calc_sqr_D0(y),x) );
814     default:
815     Asc_Panic(2, "calc_div_D1",
816     "calc_div_D1: Passed %d as value for \'wrt\'", wrt);
817     exit(2);/* Needed to keep gcc from whining */
818     }
819     }
820    
821     extern real64 calc_pow_D1(real64 x, real64 y, int wrt)
822     {
823     switch( wrt ) {
824     case 0:
825     return( calc_mul_D0(y,calc_pow_D0(x,y-1.0)) );
826     case 1:
827     return( calc_mul_D0( calc_ln_D0(x) , calc_pow_D0(x,y) ) );
828     default:
829     Asc_Panic(2, "calc_pow_D1",
830     "calc_pow_D1: Passed %d as value for \'wrt\'", wrt);
831     exit(2);/* Needed to keep gcc from whining */
832     }
833     }
834    
835     extern real64 calc_div_D2(real64 x, real64 y, int wrt1, int wrt2)
836     {
837     (void)x;
838     switch( wrt1+wrt2 ) {
839     case 0:
840     return(0.0);
841     case 1:
842     return( -calc_rec(calc_sqr_D0(y)) );
843     case 2:
844     return( calc_rec(0.5*calc_mul_D0(y,calc_sqr_D0(y))) );
845     default:
846     Asc_Panic(2, "calc_div_D2",
847     "calc_div_D2: Passed bad values for \'wrt1\' and \'wrt2\'");
848     exit(2);/* Needed to keep gcc from whining */
849     }
850     }
851    
852     real64 calc_pow_D2(real64 x, real64 y, int wrt1, int wrt2)
853     {
854     double lnx;
855     switch( wrt1+wrt2 ) {
856     case 0:
857     return( calc_mul_D0( calc_mul_D0(y,y-1.0),calc_pow_D0(x,y-2.0) ) );
858     case 1:
859     lnx = calc_ln_D0(x);
860     return( calc_mul_D0(calc_rec(x)+calc_sqr_D0(lnx), calc_pow_D0(x,y)) );
861     case 2:
862     lnx = calc_ln_D0(x);
863     return( calc_mul_D0(calc_sqr_D0(lnx) , calc_pow_D0(x,y)) );
864     default:
865     Asc_Panic(2, "calc_pow_D2",
866     "calc_pow_D2: Passed bad values for \'wrt1\' and \'wrt2\'");
867     exit(2);/* Needed to keep gcc from whining */
868     }
869     }
870    
871     real64 calc_add_Dn(real64 x, real64 y, int nwrt0, int nwrt1)
872     {
873     switch( nwrt0+nwrt1 ) {
874     case 0:
875     return( x+y );
876     case 1:
877     return( 1.0 );
878     default:
879     return( 0.0 );
880     }
881     }
882    
883     real64 calc_sub_Dn(real64 x, real64 y, int nwrt0, int nwrt1)
884     {
885     switch( nwrt0+nwrt1 ) {
886     case 0:
887     return( x-y );
888     case 1:
889     return( nwrt1==0 ? 1.0 : -1.0 );
890     default:
891     return( 0.0 );
892     }
893     }
894    
895     real64 calc_mul_Dn(real64 x, real64 y, int nwrt0, int nwrt1)
896     {
897     switch( nwrt0 ) {
898     case 0:
899     break;
900     case 1:
901     x = 1.0;
902     break;
903     default:
904     return(0.0);
905     }
906    
907     switch( nwrt1 ) {
908     case 0:
909     return( calc_mul_D0(x,y) );
910     case 1:
911     return(x);
912     default:
913     return(0.0);
914     }
915     }
916    
917     real64 calc_div_Dn(real64 x, real64 y, int nwrt0, int nwrt1)
918     {
919     switch( nwrt0 ) {
920     case 1:
921     x = 1.0;
922     /*** FALL THROUGH ***/
923     case 0:
924     return( calc_mul_D0(x,calc_Dn_rec(y,nwrt1)) );
925     default:
926     return(0.0);
927     }
928     }
929    
930    
931     /*
932     * The n-th derivative wrt x of (lnx)^m * x^y (the m-th derivative wrt y
933     * of x^y) = x^(y-n) * P[n](lnx), where P[0](z) = z^m and
934     * P[n+1](z) = (y-n)*P[n](z) + P[n]'(z). (n==nwrt0,m==nwrt1)
935     */
936     real64 calc_pow_Dn(real64 x, real64 y, int nwrt0, int nwrt1)
937     {
938     real64 *poly;
939     real64 lnx,lnx2; /* To be calculated later if necessary */
940     int k,r;
941     poly = alloc_poly(nwrt1);
942    
943     mem_zero_byte_cast(poly,0,nwrt1*sizeof(real64));
944     poly[nwrt1] = calc_pow_D0(x,y-(double)nwrt0);
945     for( k=0 ; k<nwrt0 ; ++k ) {
946     /* calculate P[k+1] from P[k] */
947     for( r=0 ; r<nwrt1 ; ++r ) {
948     poly[r] = (y-k)*poly[r] + (r+1)*poly[r+1];
949     }
950     poly[nwrt1] *= y-k;
951     }
952    
953     if( nwrt1 == 0 ) {
954     /* polynominal is just a constant */
955     return( poly[0] );
956     }
957    
958     lnx2 = calc_sqr_D0( lnx = calc_ln_D0(x) );
959     return( calc_poly(poly,nwrt1,lnx,lnx2)+calc_poly(poly,nwrt1-1,lnx,lnx2) );
960     }

Properties

Name Value
svn:executable *

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