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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 100 by aw0a, Sat Nov 13 16:45:56 2004 UTC revision 101 by jds, Sat Dec 10 04:22:07 2005 UTC
# Line 47  Line 47 
47  int32 calc_ok = TRUE;  int32 calc_ok = TRUE;
48  boolean calc_print_errors = TRUE;  boolean calc_print_errors = TRUE;
49    
50  static real64 rec_1_p_sqr_Dn(x,n)  static real64 rec_1_p_sqr_Dn(real64 x, int n)
 real64 x;  
 int n;  
51  /**  /**
52   ***  Calculates n-th derivative of 1/(1+x^2).   ***  Calculates n-th derivative of 1/(1+x^2).
53   **/   **/
# Line 73  int n; Line 71  int n;
71     return(prev[1]);     return(prev[1]);
72  }  }
73    
74  static real64 *alloc_poly(order)  static real64 *alloc_poly(int order)
 int order;  
75  /**  /**
76   ***  Allocates a polynominal of given order and returns it.  The   ***  Allocates a polynominal of given order and returns it.  The
77   ***  polynominal need not be freed, but this function should not be   ***  polynominal need not be freed, but this function should not be
# Line 93  int order; Line 90  int order;
90     return(poly);     return(poly);
91  }  }
92    
93  static real64 calc_poly(poly,order,x,x2)  static real64 calc_poly(real64 *poly, int order, real64 x, real64 x2)
 real64 *poly;  
 int order;  
 real64 x,x2;   /* x2 = x^2 */  
94  /**  /**
95   ***  Calculates the value of the given polynomial, where only the   ***  Calculates the value of the given polynomial, where only the
96   ***  (order%1 ? odd : even) degree terms are used.   ***  (order%1 ? odd : even) degree terms are used.  x2 = x^2
97   **/   **/
98  {  {
99     real64 val;     real64 val;
# Line 109  real64 x,x2;   /* x2 = x^2 */ Line 103  real64 x,x2;   /* x2 = x^2 */
103     return( order==-1 ? calc_mul_D0(val,x) : val );     return( order==-1 ? calc_mul_D0(val,x) : val );
104  }  }
105    
106  static real64 exp_msqr_Dn(x,n)  static real64 exp_msqr_Dn(real64 x, int n)
 real64 x;  
 int n;  
107  /**  /**
108   ***  Computes n-th derivative of exp(-x^2).   ***  Computes n-th derivative of exp(-x^2).
109   **/   **/
# Line 138  int n; Line 130  int n;
130     return( calc_poly(poly,n,x,x2) );     return( calc_poly(poly,n,x,x2) );
131  }  }
132    
133  static real64 sqrt_rec_1_m_sqr_Dn(x,n)  static real64 sqrt_rec_1_m_sqr_Dn(real64 x, int n)
 real64 x;  
 int n;  
134  /**  /**
135   ***  Calculates n-th derivative of (1-x^2)^-.5   ***  Calculates n-th derivative of (1-x^2)^-.5
136   **/   **/
# Line 208  real64 calc_rec(real64 x) Line 198  real64 calc_rec(real64 x)
198  #define CUBRTHUGE      5.6438030941223618e102  #define CUBRTHUGE      5.6438030941223618e102
199  /* largest cubeable number in an 8 byte ieee double */  /* largest cubeable number in an 8 byte ieee double */
200  #endif  #endif
201  real64 calc_cube(x)  real64 calc_cube(register real64 x)
 register real64 x;  
202  {  {
203     if( fabs(x) > CUBRTHUGE || fabs(x) < INV_CUBRTHUGE ) {     if( fabs(x) > CUBRTHUGE || fabs(x) < INV_CUBRTHUGE ) {
204        real64 bogus;        real64 bogus;
# Line 234  register real64 x; Line 223  register real64 x;
223     return ( x*x*x );     return ( x*x*x );
224  }  }
225    
226  real64 calc_Dn_rec(x,n)  real64 calc_Dn_rec(real64 x, int n)
 real64 x;  
 int n;  
227  {  {
228     return( -calc_upower(-calc_rec(x),n+1) * calc_factorial(n) );     return( -calc_upower(-calc_rec(x),n+1) * calc_factorial(n) );
229  }  }
230    
231  #ifdef HAVE_ERF  #ifdef HAVE_ERF
232  real64 calc_erf_inv(x)  real64 calc_erf_inv(real64 x)
 real64 x;  
233  {  {
234  #define CONV (1e-7)  #define CONV (1e-7)
235     real64 y,dy,sign;     real64 y,dy,sign;
# Line 279  real64 x; Line 265  real64 x;
265  }  }
266  #endif /* HAVE_ERF */  #endif /* HAVE_ERF */
267    
268  real64 calc_lnm_inv(x)  real64 calc_lnm_inv(real64 x)
 real64 x;  
269  {  {
270     register double eps=FuncGetLnmEpsilon();     register double eps=FuncGetLnmEpsilon();
271     if( x > (MAXDOUBLE > 1.0e308 ? 709.196209 : log(MAXDOUBLE)) ) {     if( x > (MAXDOUBLE > 1.0e308 ? 709.196209 : log(MAXDOUBLE)) ) {
# Line 297  real64 x; Line 282  real64 x;
282     return( (x >= log(eps)) ? calc_exp_D0(x) : eps*(x+1.0-log(eps)) );     return( (x >= log(eps)) ? calc_exp_D0(x) : eps*(x+1.0-log(eps)) );
283  }  }
284    
285  int calc_is_int(x)  int calc_is_int(real64 x)
 real64 x;  
286  {  {
287  #define TOLERANCE (1e-4)  #define TOLERANCE (1e-4)
288     unsigned n;     unsigned n;
# Line 326  real64 x; Line 310  real64 x;
310  }  }
311    
312    
313  real64 calc_tan_D0(x)  real64 calc_tan_D0(real64 x)
 real64 x;  
314  {  {
315     return( calc_div_D0( sin(x) , cos(x) ) );     return( calc_div_D0( sin(x) , cos(x) ) );
316  }  }
317    
318  real64 calc_arcsin_D0(x)  real64 calc_arcsin_D0(real64 x)
 real64 x;  
319  {  {
320     if( x < -1.0 || 1.0 < x ) {     if( x < -1.0 || 1.0 < x ) {
321        real64 bogus = 0.0;        real64 bogus = 0.0;
# Line 348  real64 x; Line 330  real64 x;
330     return( asin(x) );     return( asin(x) );
331  }  }
332    
333  real64 calc_arccos_D0(x)  real64 calc_arccos_D0(real64 x)
 real64 x;  
334  {  {
335     if( x < -1.0 || 1.0 < x ) {     if( x < -1.0 || 1.0 < x ) {
336        real64 bogus = 0.0;        real64 bogus = 0.0;
# Line 364  real64 x; Line 345  real64 x;
345     return( acos(x) );     return( acos(x) );
346  }  }
347    
348  real64 calc_arccosh_D0(x)  real64 calc_arccosh_D0(real64 x)
 real64 x;  
349  {  {
350     if( x < 1.0 ) {     if( x < 1.0 ) {
351        real64 bogus = 0.0;        real64 bogus = 0.0;
# Line 381  real64 x; Line 361  real64 x;
361     return( arccosh(x) );     return( arccosh(x) );
362  }  }
363    
364  real64 calc_arctanh_D0(x)  real64 calc_arctanh_D0(real64 x)
 real64 x;  
365  {  {
366     if( x < -1.0 || 1.0 < x ) {     if( x < -1.0 || 1.0 < x ) {
367        real64 bogus = 0.0;        real64 bogus = 0.0;
# Line 397  real64 x; Line 376  real64 x;
376     return( arctanh(x) );     return( arctanh(x) );
377  }  }
378    
379  real64 calc_exp_D0(x)  real64 calc_exp_D0(real64 x)
 real64 x;  
380  {  {
381     if( x > (MAXDOUBLE > 1.0e308 ? 709.196209 : log(MAXDOUBLE)) ) {     if( x > (MAXDOUBLE > 1.0e308 ? 709.196209 : log(MAXDOUBLE)) ) {
382        real64 bogus = BIGNUM;        real64 bogus = BIGNUM;
# Line 413  real64 x; Line 391  real64 x;
391     return( exp(x) );     return( exp(x) );
392  }  }
393    
394  real64 calc_ln_D0(x)  real64 calc_ln_D0(real64 x)
 real64 x;  
395  {  {
396     if( x <= 0.0 ) {     if( x <= 0.0 ) {
397        real64 bogus = -BIGNUM;        real64 bogus = -BIGNUM;
# Line 429  real64 x; Line 406  real64 x;
406     return( log(x) );     return( log(x) );
407  }  }
408    
409  real64 calc_log_D0(x)  real64 calc_log_D0(real64 x)
 real64 x;  
410  {  {
411     return( calc_LOG10_COEF * calc_ln_D0(x) );     return( calc_LOG10_COEF * calc_ln_D0(x) );
412  }  }
413    
414  real64 calc_sqrt_D0(x)  real64 calc_sqrt_D0(real64 x)
 real64 x;  
415  {  {
416     if( x < 0.0 ) {     if( x < 0.0 ) {
417        real64 bogus;        real64 bogus;
# Line 452  real64 x; Line 427  real64 x;
427     return(sqrt(x));     return(sqrt(x));
428  }  }
429    
430  real64 calc_tan_D1(x)  real64 calc_tan_D1(register real64 x)
 register real64 x;  
431  {  {
432     x=cos(x);     x=cos(x);
433     return( calc_rec(calc_sqr_D0(x)) );     return( calc_rec(calc_sqr_D0(x)) );
434  }  }
435    
436  real64 calc_arcsin_D1(x)  real64 calc_arcsin_D1(real64 x)
 real64 x;  
437  {  {
438     return( calc_rec(calc_sqrt_D0(1.0 - calc_sqr_D0(x))) );     return( calc_rec(calc_sqrt_D0(1.0 - calc_sqr_D0(x))) );
439  }  }
440    
441  real64 calc_arccos_D1(x)  real64 calc_arccos_D1(real64 x)
 real64 x;  
442  {  {
443     return( -calc_arcsin_D1(x) );     return( -calc_arcsin_D1(x) );
444  }  }
445    
446  real64 calc_arccosh_D1(x)  real64 calc_arccosh_D1(real64 x)
 real64 x;  
447  {  {
448     return( calc_rec(calc_sqrt_D0(calc_sqr_D0(x) - 1.0)) );     return( calc_rec(calc_sqrt_D0(calc_sqr_D0(x) - 1.0)) );
449  }  }
450    
451  real64 calc_arctanh_D1(x)  real64 calc_arctanh_D1(real64 x)
 real64 x;  
452  {  {
453     return( calc_rec(1.0 - calc_sqr_D0(x)) );     return( calc_rec(1.0 - calc_sqr_D0(x)) );
454  }  }
455    
456  #ifdef HAVE_ERF  #ifdef HAVE_ERF
457  real64 calc_erf_D1(x)  real64 calc_erf_D1(real64 x)
 real64 x;  
458  {  {
459     return( calc_ERF_COEF * calc_exp_D0(-calc_sqr_D0(x)) );     return( calc_ERF_COEF * calc_exp_D0(-calc_sqr_D0(x)) );
460  }  }
461  #endif /* HAVE_ERF */  #endif /* HAVE_ERF */
462    
463  real64 calc_log_D1(x)  real64 calc_log_D1(real64 x)
 real64 x;  
464  {  {
465     return( calc_LOG10_COEF * calc_ln_D1(x) );     return( calc_LOG10_COEF * calc_ln_D1(x) );
466  }  }
467    
468  real64 calc_sqrt_D1(x)  real64 calc_sqrt_D1(real64 x)
 real64 x;  
469  {  {
470     return( 0.5 * calc_rec(calc_sqrt_D0(x)) );     return( 0.5 * calc_rec(calc_sqrt_D0(x)) );
471  }  }
472    
473  real64 calc_cbrt_D1(x)  real64 calc_cbrt_D1(real64 x)
 real64 x;  
474  {  {
475     return( 0.3333333333333333*calc_rec(calc_sqr_D0(calc_cbrt_D0(x))) );     return( 0.3333333333333333*calc_rec(calc_sqr_D0(calc_cbrt_D0(x))) );
476  }  }
477  /* KHACK */  /* KHACK */
478  real64 calc_fabs_D1(x)  real64 calc_fabs_D1(real64 x)
 real64 x;  
479  {  {
480     (void)x;     (void)x;
481     return( 1.0 );     return( 1.0 );
482  }  }
483    
484  real64 calc_tan_D2(x)  real64 calc_tan_D2(real64 x)
 real64 x;  
485  {  {
486     if( fabs(cos(x)) == 0.0 ) {     if( fabs(cos(x)) == 0.0 ) {
487        real64 bogus = BIGNUM;        real64 bogus = BIGNUM;
# Line 532  real64 x; Line 496  real64 x;
496     return dtan2(x);     return dtan2(x);
497  }  }
498    
499  real64 calc_arcsin_D2(x)  real64 calc_arcsin_D2(real64 x)
 real64 x;  
500  {  {
501     register double t;     register double t;
502     t = calc_rec(1.0 - calc_sqr_D0(x));     t = calc_rec(1.0 - calc_sqr_D0(x));
503     return( calc_mul_D0( calc_mul_D0(x,t) , calc_mul_D0(t,t) ) );     return( calc_mul_D0( calc_mul_D0(x,t) , calc_mul_D0(t,t) ) );
504  }  }
505    
506  real64 calc_arccos_D2(x)  real64 calc_arccos_D2(real64 x)
 real64 x;  
507  {  {
508     return( -calc_arcsin_D2(x) );     return( -calc_arcsin_D2(x) );
509  }  }
510    
511  real64 calc_arccosh_D2(x)  real64 calc_arccosh_D2(real64 x)
 real64 x;  
512  {  {
513     if( fabs(x) <= 1.0 ) {     if( fabs(x) <= 1.0 ) {
514        real64 bogus = -BIGNUM;        real64 bogus = -BIGNUM;
# Line 562  real64 x; Line 523  real64 x;
523     return darccosh2(x);     return darccosh2(x);
524  }  }
525    
526  real64 calc_arctanh_D2(x)  real64 calc_arctanh_D2(real64 x)
 real64 x;  
527  {  {
528     if( fabs(x) == 1.0 ) {     if( fabs(x) == 1.0 ) {
529        real64 bogus = BIGNUM;        real64 bogus = BIGNUM;
# Line 579  real64 x; Line 539  real64 x;
539  }  }
540    
541  #ifdef HAVE_ERF  #ifdef HAVE_ERF
542  real64 calc_erf_D2(x)  real64 calc_erf_D2(real64 x)
 real64 x;  
543  {  {
544     return( -ldexp(calc_ERF_COEF * calc_mul_D0(x,calc_exp_D0(-calc_sqr_D0(x))),     return( -ldexp(calc_ERF_COEF * calc_mul_D0(x,calc_exp_D0(-calc_sqr_D0(x))),
545                    1) );                    1) );
546  }  }
547  #endif /* HAVE_ERF */  #endif /* HAVE_ERF */
548    
549  real64 calc_ln_D2(x)  real64 calc_ln_D2(real64 x)
 real64 x;  
550  {  {
551     return( -calc_rec(calc_sqr_D0(x)) );     return( -calc_rec(calc_sqr_D0(x)) );
552  }  }
553    
554  real64 calc_log_D2(x)  real64 calc_log_D2(real64 x)
 real64 x;  
555  {  {
556     return( calc_LOG10_COEF * calc_ln_D2(x) );     return( calc_LOG10_COEF * calc_ln_D2(x) );
557  }  }
558    
559  real64 calc_sqrt_D2(x)  real64 calc_sqrt_D2(real64 x)
 real64 x;  
560  {  {
561     return( -ldexp( calc_rec( calc_mul_D0(x,calc_sqrt_D0(x))), -2) );     return( -ldexp( calc_rec( calc_mul_D0(x,calc_sqrt_D0(x))), -2) );
562  }  }
563    
564  real64 calc_cbrt_D2(x)  real64 calc_cbrt_D2(real64 x)
 real64 x;  
565  {  {
566     if( fabs(x) == 0.0 ) {     if( fabs(x) == 0.0 ) {
567        real64 bogus = BIGNUM;        real64 bogus = BIGNUM;
# Line 621  real64 x; Line 576  real64 x;
576     return dcbrt2(x);     return dcbrt2(x);
577  }  }
578  /* KHACK */  /* KHACK */
579  real64 calc_fabs_D2(x)  real64 calc_fabs_D2(real64 x)
 real64 x;  
580  {  {
581     (void)x;     (void)x;
582     return( 0.0 );     return( 0.0 );
583  }  }
584    
585  real64 calc_arcsin_Dn(x,n)  real64 calc_arcsin_Dn(real64 x, int n)
 real64 x;  
 int n;  
586  {  {
587     return( n==0 ? calc_arcsin_D0(x) : sqrt_rec_1_m_sqr_Dn(x,n-1) );     return( n==0 ? calc_arcsin_D0(x) : sqrt_rec_1_m_sqr_Dn(x,n-1) );
588  }  }
589    
590  real64 calc_arccos_Dn(x,n)  real64 calc_arccos_Dn(real64 x, int n)
 real64 x;  
 int n;  
591  {  {
592     return( n==0 ? calc_arccos_D0(x) : -sqrt_rec_1_m_sqr_Dn(x,n-1) );     return( n==0 ? calc_arccos_D0(x) : -sqrt_rec_1_m_sqr_Dn(x,n-1) );
593  }  }
594    
595  real64 calc_arctan_Dn(x,n)  real64 calc_arctan_Dn(real64 x, int n)
 real64 x;  
 int n;  
596  {  {
597     return( n==0 ? calc_arctan_D0(x) : rec_1_p_sqr_Dn(x,n-1) );     return( n==0 ? calc_arctan_D0(x) : rec_1_p_sqr_Dn(x,n-1) );
598  }  }
599    
600  real64 calc_cos_Dn(x,n)  real64 calc_cos_Dn(real64 x, int n)
 real64 x;  
 int n;  
601  {  {
602     switch( n%4 ) {     switch( n%4 ) {
603        case 0:        case 0:
# Line 671  int n; Line 617  int n;
617     return 0.0;     return 0.0;
618  }  }
619    
620  real64 calc_sin_Dn(x,n)  real64 calc_sin_Dn(real64 x, int n)
 real64 x;  
 int n;  
621  {  {
622     return( calc_cos_Dn(x,n+3) );     return( calc_cos_Dn(x,n+3) );
623  }  }
624    
625  #ifdef HAVE_ERF  #ifdef HAVE_ERF
626  real64 calc_erf_Dn(x,n)  real64 calc_erf_Dn(real64 x, int n)
 real64 x;  
 int n;  
627  {  {
628     return( n==0 ? calc_erf_D0(x) : calc_ERF_COEF*exp_msqr_Dn(x,n-1) );     return( n==0 ? calc_erf_D0(x) : calc_ERF_COEF*exp_msqr_Dn(x,n-1) );
629  }  }
630  #endif /* HAVE_ERF */  #endif /* HAVE_ERF */
631    
632  real64 calc_exp_Dn(x,n)  real64 calc_exp_Dn(real64 x, int n)
 real64 x;  
 int n;  
633  {  {
634     (void)n;     (void)n;
635     return( calc_exp_D0(x) );     return( calc_exp_D0(x) );
636  }  }
637    
638  real64 calc_ln_Dn(x,n)  real64 calc_ln_Dn(real64 x, int n)
 real64 x;  
 int n;  
639  {  {
640     return( n==0 ? calc_ln_D0(x) : calc_Dn_rec(x,n-1) );     return( n==0 ? calc_ln_D0(x) : calc_Dn_rec(x,n-1) );
641  }  }
642    
643  real64 calc_log_Dn(x,n)  real64 calc_log_Dn(real64 x, int n)
 real64 x;  
 int n;  
644  {  {
645     return( calc_LOG10_COEF * calc_ln_Dn(x,n) );     return( calc_LOG10_COEF * calc_ln_Dn(x,n) );
646  }  }
647    
648  real64 calc_sqr_Dn(x,n)  real64 calc_sqr_Dn(real64 x, int n)
 real64 x;  
 int n;  
649  {  {
650     switch(n) {     switch(n) {
651        case 0:        case 0:
# Line 725  int n; Line 659  int n;
659     }     }
660  }  }
661    
662  real64 calc_sqrt_Dn(x,n)  real64 calc_sqrt_Dn(real64 x, int n)
 real64 x;  
 int n;  
663  {  {
664     double a,b;     double a,b;
665    

Legend:
Removed from v.100  
changed lines
  Added in v.101

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