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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (show annotations) (download) (as text)
Fri Oct 29 20:54:12 2004 UTC (17 years, 6 months ago) by aw0a
File MIME type: text/x-csrc
File size: 21615 byte(s)
Setting up web subdirectory in repository
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