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

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