/[ascend]/trunk/base/generic/solver/calc.c
ViewVC logotype

Contents of /trunk/base/generic/solver/calc.c

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

Name Value
svn:executable *

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