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

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