home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / cmath / mthutl.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  14.2 KB  |  571 lines
  1. /*                            mconf.h
  2.  *
  3.  *    Common include file for math routines
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * #include "mconf.h"
  10.  *
  11.  *
  12.  *
  13.  * DESCRIPTION:
  14.  *
  15.  * This file contains definitions for error codes that are
  16.  * passed to the common error handling routine mtherr()
  17.  * (which see).
  18.  *
  19.  * The file also includes a conditional assembly definition
  20.  * for the type of computer arithmetic (IEEE, DEC, Motorola
  21.  * IEEE, or UNKnown).
  22.  *
  23.  * For Digital Equipment PDP-11 and VAX computers, certain
  24.  * IBM systems, and others that use numbers with a 56-bit
  25.  * significand, the symbol DEC should be defined.  In this
  26.  * mode, most floating point constants are given as arrays
  27.  * of octal integers to eliminate decimal to binary conversion
  28.  * errors that might be introduced by the compiler.
  29.  *
  30.  * For computers, such as IBM PC, that follow the IEEE 
  31.  * Standard for Binary Floating Point Arithmetic (ANSI/IEEE
  32.  * Std 754-1985), the symbol IBMPC should be defined.  These
  33.  * numbers have 53-bit significands.  In this mode, constants
  34.  * are provided as arrays of hexadecimal 16 bit integers.
  35.  *
  36.  * To accommodate other types of computer arithmetic, all
  37.  * constants are also provided in a normal decimal radix
  38.  * which one can hope are correctly converted to a suitable
  39.  * format by the available C language compiler.  To invoke
  40.  * this mode, the symbol UNK is defined.
  41.  *
  42.  * An important difference among these modes is a predefined
  43.  * set of machine arithmetic constants for each.  The numbers
  44.  * MACHEP (the machine roundoff error), MAXNUM (largest number
  45.  * represented), and several other parameters are preset by
  46.  * the configuration symbol.  Check the file const.c to
  47.  * ensure that these values are correct for your computer.
  48.  *
  49.  */
  50.  
  51. /*
  52. Cephes Math Library Release 2.0:  April, 1987
  53. Copyright 1984, 1987 by Stephen L. Moshier
  54. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  55. */
  56.  
  57.  
  58. /* Constant definitions for math error conditions
  59.  */
  60.  
  61. #define DOMAIN        1    /* argument domain error */
  62. #define SING        2    /* argument singularity */
  63. #define OVERFLOW    3    /* overflow range error */
  64. #define UNDERFLOW    4    /* underflow range error */
  65. #define TLOSS        5    /* total loss of precision */
  66. #define PLOSS        6    /* partial loss of precision */
  67.  
  68. #define EDOM        33
  69. #define ERANGE        34
  70.  
  71. typedef struct
  72.     {
  73.     double r;
  74.     double i;
  75.     }cmplx;
  76.  
  77. /* Type of computer arithmetic */
  78.  
  79. /* PDP-11, Pro350, VAX:
  80.  */
  81. /*define DEC 1*/
  82.  
  83. /* Intel IEEE, low order words come first:
  84.  */
  85. #define IBMPC 1
  86.  
  87. /* Motorola IEEE, high order words come first
  88.  * (Sun workstation):
  89.  */
  90. /*define MIEEE 1*/
  91.  
  92. /* UNKnown arithmetic, invokes coefficients given in
  93.  * normal decimal format.  Beware of range boundary
  94.  * problems (MACHEP, MAXLOG, etc. in const.c) and
  95.  * roundoff problems in pow.c:
  96.  /*define UNK 1*/
  97.  
  98. /*                            const.c
  99.  *
  100.  *    Globally declared constants
  101.  *
  102.  *
  103.  *
  104.  * SYNOPSIS:
  105.  *
  106.  * extern double nameofconstant;
  107.  *
  108.  *
  109.  *
  110.  *
  111.  * DESCRIPTION:
  112.  *
  113.  * This file contains a number of mathematical constants and
  114.  * also some needed size parameters of the computer arithmetic.
  115.  * The values are supplied as arrays of hexadecimal integers
  116.  * for IEEE arithmetic; arrays of octal constants for DEC
  117.  * arithmetic; and in a normal decimal scientific notation for
  118.  * other machines.  The particular notation used is determined
  119.  * by a symbol (DEC, IBMPC, or UNK) defined in the include file
  120.  * mconf.h.
  121.  *
  122.  * The default size parameters are as follows.
  123.  *
  124.  * For DEC and UNK modes:
  125.  * MACHEP =  1.38777878078144567553E-17       2**-56
  126.  * MAXLOG =  8.8029691931113054295988E1       log(2**127)
  127.  * MINLOG = -8.872283911167299960540E1        log(2**-128)
  128.  * MAXNUM =  1.701411834604692317316873e38    2**127
  129.  *
  130.  * For IEEE arithmetic (IBMPC):
  131.  * MACHEP =  1.11022302462515654042E-16       2**-53
  132.  * MAXLOG =  7.09782712893383996843E2         log(2**1024)
  133.  * MINLOG = -7.08396418532264106224E2         log(2**-1022)
  134.  * MAXNUM =  1.7976931348623158E308           2**1024
  135.  *
  136.  * The global symbols for mathematical constants are
  137.  * PI     =  3.14159265358979323846           pi
  138.  * PIO2   =  1.57079632679489661923           pi/2
  139.  * PIO4   =  7.85398163397448309616E-1        pi/4
  140.  * SQRT2  =  1.41421356237309504880           sqrt(2)
  141.  * SQRTH  =  7.07106781186547524401E-1        sqrt(2)/2
  142.  * LOG2E  =  1.4426950408889634073599         1/log(2)
  143.  * SQ2OPI =  7.9788456080286535587989E-1      sqrt( 2/pi )
  144.  * LOGE2  =  6.93147180559945309417E-1        log(2)
  145.  * LOGSQ2 =  3.46573590279972654709E-1        log(2)/2
  146.  * THPIO4 =  2.35619449019234492885           3*pi/4
  147.  * TWOOPI =  6.36619772367581343075535E-1     2/pi
  148.  *
  149.  * These lists are subject to change.
  150.  */
  151.  
  152. /*                            const.c */
  153.  
  154. /*
  155. Cephes Math Library Release 2.0:  April, 1987
  156. Copyright 1984, 1987 by Stephen L. Moshier
  157. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  158. */
  159.  
  160. #include "mconf.h"
  161.  
  162. #ifdef UNK
  163. double MACHEP =  1.38777878078144567553E-17;   /* 2**-56 */
  164. double MAXLOG =  8.8029691931113054295988E1;   /* log(2**127) */
  165. double MINLOG = -8.872283911167299960540E1;    /* log(2**-128) */
  166. double MAXNUM =  1.701411834604692317316873e38; /* 2**127 */
  167. double PI     =  3.14159265358979323846;       /* pi */
  168. double PIO2   =  1.57079632679489661923;       /* pi/2 */
  169. double PIO4   =  7.85398163397448309616E-1;    /* pi/4 */
  170. double SQRT2  =  1.41421356237309504880;       /* sqrt(2) */
  171. double SQRTH  =  7.07106781186547524401E-1;    /* sqrt(2)/2 */
  172. double LOG2E  =  1.4426950408889634073599;     /* 1/log(2) */
  173. double SQ2OPI =  7.9788456080286535587989E-1;  /* sqrt( 2/pi ) */
  174. double LOGE2  =  6.93147180559945309417E-1;    /* log(2) */
  175. double LOGSQ2 =  3.46573590279972654709E-1;    /* log(2)/2 */
  176. double THPIO4 =  2.35619449019234492885;       /* 3*pi/4 */
  177. double TWOOPI =  6.36619772367581343075535E-1; /* 2/pi */
  178. #endif
  179.  
  180. #ifdef IBMPC
  181.             /* 2**-53 =  1.11022302462515654042E-16 */
  182. short MACHEP[4] = {0x0000,0x0000,0x0000,0x3ca0};
  183.             /* log(2**1024) =   7.09782712893383996843E2 */
  184. short MAXLOG[4] = {0x39ef,0xfefa,0x2e42,0x4086};
  185.             /* log(2**-1022) = - 7.08396418532264106224E2 */
  186. short MINLOG[4] = {0xbcd2,0xdd7a,0x232b,0xc086};
  187.             /* 2**1024*(1-MACHEP) =  1.7976931348623158E308 */
  188. short MAXNUM[4] = {0xffff,0xffff,0xffff,0x7fef};
  189. short PI[4]     = {0x2d18,0x5444,0x21fb,0x4009};
  190. short PIO2[4]   = {0x2d18,0x5444,0x21fb,0x3ff9};
  191. short PIO4[4]   = {0x2d18,0x5444,0x21fb,0x3fe9};
  192. short SQRT2[4]  = {0x3bcd,0x667f,0xa09e,0x3ff6};
  193. short SQRTH[4]  = {0x3bcd,0x667f,0xa09e,0x3fe6};
  194. short LOG2E[4]  = {0x82fe,0x652b,0x1547,0x3ff7};
  195. short SQ2OPI[4] = {0x3651,0x33d4,0x8845,0x3fe9};
  196. short LOGE2[4]  = {0x39ef,0xfefa,0x2e42,0x3fe6};
  197. short LOGSQ2[4] = {0x39ef,0xfefa,0x2e42,0x3fd6};
  198. short THPIO4[4] = {0x21d2,0x7f33,0xd97c,0x4002};
  199. short TWOOPI[4] = {0xc883,0x6dc9,0x5f30,0x3fe4};
  200. #endif
  201.  
  202. #ifdef MIEEE
  203.             /* 2**-53 =  1.11022302462515654042E-16 */
  204. short MACHEP[4] = {
  205. 0x3ca0,0x0000,0x0000,0x0000
  206. };
  207.             /* log(2**1024) =   7.09782712893383996843E2 */
  208. short MAXLOG[4] = {
  209. 0x4086,0x2e42,0xfefa,0x39ef
  210. };
  211.             /* log(2**-1022) = - 7.08396418532264106224E2 */
  212. short MINLOG[4] = {
  213. 0xc086,0x232b,0xdd7a,0xbcd2
  214. };
  215.             /* 2**1024*(1-MACHEP) =  1.7976931348623158E308 */
  216. short MAXNUM[4] = {
  217. 0x7fef,0xffff,0xffff,0xffff
  218. };
  219. short PI[4]     = {
  220. 0x4009,0x21fb,0x5444,0x2d18
  221. };
  222. short PIO2[4]   = {
  223. 0x3ff9,0x21fb,0x5444,0x2d18
  224. };
  225. short PIO4[4]   = {
  226. 0x3fe9,0x21fb,0x5444,0x2d18
  227. };
  228. short SQRT2[4]  = {
  229. 0x3ff6,0xa09e,0x667f,0x3bcd
  230. };
  231. short SQRTH[4]  = {
  232. 0x3fe6,0xa09e,0x667f,0x3bcd
  233. };
  234. short LOG2E[4]  = {
  235. 0x3ff7,0x1547,0x652b,0x82fe
  236. };
  237. short SQ2OPI[4] = {
  238. 0x3fe9,0x8845,0x33d4,0x3651
  239. };
  240. short LOGE2[4]  = {
  241. 0x3fe6,0x2e42,0xfefa,0x39ef
  242. };
  243. short LOGSQ2[4] = {
  244. 0x3fd6,0x2e42,0xfefa,0x39ef
  245. };
  246. short THPIO4[4] = {
  247. 0x4002,0xd97c,0x7f33,0x21d2
  248. };
  249. short TWOOPI[4] = {
  250. 0x3fe4,0x5f30,0x6dc9,0xc883
  251. };
  252. #endif
  253.  
  254. #ifdef DEC
  255.             /* 2**-56 =  1.38777878078144567553E-17 */
  256. short MACHEP[4] = {0022200,0000000,0000000,0000000};
  257.             /* log 2**127 = 88.029691931113054295988 */
  258. short MAXLOG[4] = {041660,007463,0143742,025733,};
  259.             /* log 2**-128 = -88.72283911167299960540 */
  260. short MINLOG[4] = {0141661,071027,0173721,0147572,};
  261.             /* 2**127 = 1.701411834604692317316873e38 */
  262. short MAXNUM[4] = {077777,0177777,0177777,0177777,};
  263. short PI[4]     = {040511,007732,0121041,064302,};
  264. short PIO2[4]   = {040311,007732,0121041,064302,};
  265. short PIO4[4]   = {040111,007732,0121041,064302,};
  266. short SQRT2[4]  = {040265,002363,031771,0157145,};
  267. short SQRTH[4]  = {040065,002363,031771,0157144,};
  268. short LOG2E[4]  = {040270,0125073,024534,013761,};
  269. short SQ2OPI[4] = {040114,041051,0117241,0131204,};
  270. short LOGE2[4]  = {040061,071027,0173721,0147572,};
  271. short LOGSQ2[4] = {037661,071027,0173721,0147572,};
  272. short THPIO4[4] = {040426,0145743,0174631,007222,};
  273. short TWOOPI[4] = {040042,0174603,067116,042025,};
  274. #endif
  275.  
  276. extern short MACHEP[];
  277. extern short MAXLOG[];
  278. extern short MINLOG[];
  279. extern short MAXNUM[];
  280. extern short PI[];
  281. extern short PIO2[];
  282. extern short PIO4[];
  283. extern short SQRT2[];
  284. extern short SQRTH[];
  285. extern short LOG2E[];
  286. extern short SQ2OPI[];
  287. extern short LOGE2[];
  288. extern short LOGSQ2[];
  289. extern short THPIO4[];
  290. extern short TWOOPI[];
  291. static char x[] =
  292. "Cephes Math Library Copyright 1987 by Stephen L. Moshier";
  293.  
  294. /*                            polevl.c
  295.  *                            p1evl.c
  296.  *
  297.  *    Evaluate polynomial
  298.  *
  299.  *
  300.  *
  301.  * SYNOPSIS:
  302.  *
  303.  * int N;
  304.  * double x, y, coef[N+1], polevl[];
  305.  *
  306.  * y = polevl( C, x, N );
  307.  *
  308.  *
  309.  *
  310.  * DESCRIPTION:
  311.  *
  312.  * Evaluates polynomial of degree N:
  313.  *
  314.  *                     2          N
  315.  * y  =  C  + C x + C x  +...+ C x
  316.  *        0    1     2          N
  317.  *
  318.  * Coefficients are stored in reverse order:
  319.  *
  320.  * coef[0] = C  , ..., coef[N] = C  .
  321.  *            N                   0
  322.  *
  323.  *  The function p1evl() assumes that coef[N] = 1.0 and is
  324.  * omitted from the array.  Its calling arguments are
  325.  * otherwise the same as polevl().
  326.  *
  327.  *
  328.  * SPEED:
  329.  *
  330.  * In the interest of speed, there are no checks for out
  331.  * of bounds arithmetic.  This routine is used by most of
  332.  * the functions in the library.  Depending on available
  333.  * equipment features, the user may wish to rewrite the
  334.  * program in microcode or assembly language.
  335.  *
  336.  */
  337.  
  338. /*
  339. Cephes Math Library Release 2.0:  April, 1987
  340. Copyright 1984, 1987 by Stephen L. Moshier
  341. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  342. */
  343.  
  344.  
  345. double polevl( x, coef, N )
  346. double x;
  347. double coef[];
  348. int N;
  349. {
  350. double ans;
  351. short i;
  352. register double *p;
  353.  
  354. p = coef;
  355. ans = *p++;
  356. i = N;
  357.  
  358. do
  359.     ans = ans * x  +  *p++;
  360. while( --i );
  361.  
  362. return( ans );
  363. }
  364.  
  365. /*                            p1evl()    */
  366. /*                                          N
  367.  * Evaluate polynomial when coefficient of x  is 1.0.
  368.  * Otherwise same as polevl.
  369.  */
  370.  
  371. double p1evl( x, coef, N )
  372. double x;
  373. double coef[];
  374. int N;
  375. {
  376. double ans;
  377. register double *p;
  378. short i;
  379.  
  380. p = coef;
  381. ans = x + *p++;
  382. i = N-1;
  383.  
  384. do
  385.     ans = ans * x  + *p++;
  386. while( --i );
  387.  
  388. return( ans );
  389. }
  390.  
  391.  
  392. /*                            chbevl.c
  393.  *
  394.  *    Evaluate Chebyshev series
  395.  *
  396.  *
  397.  *
  398.  * SYNOPSIS:
  399.  *
  400.  * int N;
  401.  * double x, y, coef[N], chebevl();
  402.  *
  403.  * y = chbevl( x, coef, N );
  404.  *
  405.  *
  406.  *
  407.  * DESCRIPTION:
  408.  *
  409.  * Evaluates the series
  410.  *
  411.  *        N-1
  412.  *         - '
  413.  *  y  =   >   coef[i] T (x/2)
  414.  *         -            i
  415.  *        i=0
  416.  *
  417.  * of Chebyshev polynomials Ti at argument x/2.
  418.  *
  419.  * Coefficients are stored in reverse order, i.e. the zero
  420.  * order term is last in the array.  Note N is the number of
  421.  * coefficients, not the order.
  422.  *
  423.  * If coefficients are for the interval a to b, x must
  424.  * have been transformed to x -> 2(2x - b - a)/(b-a) before
  425.  * entering the routine.  This maps x from (a, b) to (-1, 1),
  426.  * over which the Chebyshev polynomials are defined.
  427.  *
  428.  * If the coefficients are for the inverted interval, in
  429.  * which (a, b) is mapped to (1/b, 1/a), the transformation
  430.  * required is x -> 2(2ab/x - b - a)/(b-a).  If b is infinity,
  431.  * this becomes x -> 4a/x - 1.
  432.  *
  433.  *
  434.  *
  435.  * SPEED:
  436.  *
  437.  * Taking advantage of the recurrence properties of the
  438.  * Chebyshev polynomials, the routine requires one more
  439.  * addition per loop than evaluating a nested polynomial of
  440.  * the same degree.
  441.  *
  442.  */
  443. /*                            chbevl.c    */
  444.  
  445. /*
  446. Cephes Math Library Release 2.0:  April, 1987
  447. Copyright 1985, 1987 by Stephen L. Moshier
  448. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  449. */
  450.  
  451. double chbevl( x, array, n )
  452. double x;
  453. double array[];
  454. int n;
  455. {
  456. double b0, b1, b2, *p;
  457. int i;
  458.  
  459. p = array;
  460. b0 = *p++;
  461. b1 = 0.0;
  462. i = n - 1;
  463.  
  464. do
  465.     {
  466.     b2 = b1;
  467.     b1 = b0;
  468.     b0 = x * b1  -  b2  + *p++;
  469.     }
  470. while( --i );
  471.  
  472. return( 0.5*(b0-b2) );
  473. }
  474.  
  475. /*                            mtherr.c
  476.  *
  477.  *    Library common error handling routine
  478.  *
  479.  *
  480.  *
  481.  * SYNOPSIS:
  482.  *
  483.  * char *fctnam;
  484.  * int code;
  485.  * void mtherr();
  486.  *
  487.  * mtherr( fctnam, code );
  488.  *
  489.  *
  490.  *
  491.  * DESCRIPTION:
  492.  *
  493.  * This routine may be called to report one of the following
  494.  * error conditions (in the include file mconf.h).
  495.  *  
  496.  *   Mnemonic        Value          Significance
  497.  *
  498.  *    DOMAIN            1       argument domain error
  499.  *    SING              2       function singularity
  500.  *    OVERFLOW          3       overflow range error
  501.  *    UNDERFLOW         4       underflow range error
  502.  *    TLOSS             5       total loss of precision
  503.  *    PLOSS             6       partial loss of precision
  504.  *    EDOM             33       Unix domain error code
  505.  *    ERANGE           34       Unix range error code
  506.  *
  507.  * The default version of the file prints the function name,
  508.  * passed to it by the pointer fctnam, followed by the
  509.  * error condition.  The display is directed to the standard
  510.  * output device.  The routine then returns to the calling
  511.  * program.  Users may wish to modify the program to abort by
  512.  * calling exit() under severe error conditions such as domain
  513.  * errors.
  514.  *
  515.  * Since all error conditions pass control to this function,
  516.  * the display may be easily changed, eliminated, or directed
  517.  * to an error logging device.
  518.  *
  519.  * SEE ALSO:
  520.  *
  521.  * mconf.h
  522.  *
  523.  */
  524.  
  525. /*
  526. Cephes Math Library Release 2.0:  April, 1987
  527. Copyright 1984, 1987 by Stephen L. Moshier
  528. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  529. */
  530.  
  531. #include "mconf.h"
  532.  
  533. /* Notice: the order of appearance of the following
  534.  * messages is bound to the error codes defined
  535.  * in mconf.h.
  536.  */
  537. static char *ermsg[7] = {
  538. "unknown",      /* error code 0 */
  539. "domain",       /* error code 1 */
  540. "singularity",  /* et seq.      */
  541. "overflow",
  542. "underflow",
  543. "total loss of precision",
  544. "partial loss of precision"
  545. };
  546.  
  547.  
  548.  
  549. mtherr( name, code )
  550. char *name;
  551. int code;
  552. {
  553.  
  554. /* Display string passed by calling program,
  555.  * which is supposed to be the name of the
  556.  * function in which the error occurred:
  557.  */
  558. printf( "\n%s ", name );
  559.  
  560. /* Display error message defined
  561.  * by the code argument.
  562.  */
  563. if( (code <= 0) || (code >= 6) )
  564.     code = 0;
  565. printf( "%s error\n", ermsg[code] );
  566.  
  567. /* Return to calling
  568.  * program
  569.  */
  570. }
  571.