1 /****************************************************************************
2 * chi2.cpp
3 *
4 * This module contains the function for the chi square distribution.
5 *
6 * All functions in this module are taken from the Cephes math library.
7 *
8 * Cephes Math Library Release 2.0: April, 1987
9 * Copyright 1984, 1987 by Stephen L. Moshier
10 * Direct inquiries to 30 Frost Street, Cambridge, MA 02140
11 *
12 * from Persistence of Vision(tm) Ray Tracer version 3.6.
13 * Copyright 1991-2003 Persistence of Vision Team
14 * Copyright 2003-2004 Persistence of Vision Raytracer Pty. Ltd.
15 *---------------------------------------------------------------------------
16 * NOTICE: This source code file is provided so that users may experiment
17 * with enhancements to POV-Ray and to port the software to platforms other
18 * than those supported by the POV-Ray developers. There are strict rules
19 * regarding how you are permitted to use this file. These rules are contained
20 * in the distribution and derivative versions licenses which should have been
21 * provided with this file.
22 *
23 * These licences may be found online, linked from the end-user license
24 * agreement that is located at http://www.povray.org/povlegal.html
25 *---------------------------------------------------------------------------
26 * This program is based on the popular DKB raytracer version 2.12.
27 * DKBTrace was originally written by David K. Buck.
28 * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
29 *---------------------------------------------------------------------------
30 *
31 *===========================================================================
32 * This file is part of MegaPOV, a modified and unofficial version of POV-Ray
33 * For more information on MegaPOV visit our website:
34 * http://megapov.inetart.net/
35 *===========================================================================
36 *
37 * $RCSfile: chi2.cpp,v $
38 * $Revision: 1.11 $
39 * $Author: chris $
40 *
41 *****************************************************************************/
42
43 #include "frame.h"
44 #include "userio.h"
45 #include "chi2.h"
46
47 BEGIN_POV_NAMESPACE
48
49 /*
50 Cephes Math Library Release 2.0: April, 1987
51 Copyright 1984, 1987 by Stephen L. Moshier
52 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
53 */
54
55 /*****************************************************************************
56 * Local preprocessor defines
57 ******************************************************************************/
58
59 const DBL MAXLGM = 2.556348e305;
60 const DBL BIG = 1.44115188075855872E+17;
61 const DBL MACHEP = 1.38777878078144567553E-17; /* 2**-56 */
62 const DBL MAXLOG = 8.8029691931113054295988E1; /* log(2**127) */
63 const DBL MAXNUM = 1.701411834604692317316873e38; /* 2**127 */
64 const DBL LOGPI = 1.14472988584940017414;
65
66
67
68 /*****************************************************************************
69 * Local typedefs
70 ******************************************************************************/
71
72
73
74 /*****************************************************************************
75 * Local variables
76 ******************************************************************************/
77
78 static int sgngam = 0; // GLOBAL VARIABLE
79
80 /*
81 * A[]: Stirling's formula expansion of log gamma
82 * B[], C[]: log gamma function between 2 and 3
83 */
84
85 const DBL A[] =
86 {
87 8.11614167470508450300E-4,
88 -5.95061904284301438324E-4,
89 7.93650340457716943945E-4,
90 -2.77777777730099687205E-3,
91 8.33333333333331927722E-2
92 };
93
94 const DBL B[] =
95 {
96 -1.37825152569120859100E3,
97 -3.88016315134637840924E4,
98 -3.31612992738871184744E5,
99 -1.16237097492762307383E6,
100 -1.72173700820839662146E6,
101 -8.53555664245765465627E5
102 };
103
104 const DBL C[] =
105 {
106 1.00000000000000000000E0,
107 -3.51815701436523470549E2,
108 -1.70642106651881159223E4,
109 -2.20528590553854454839E5,
110 -1.13933444367982507207E6,
111 -2.53252307177582951285E6,
112 -2.01889141433532773231E6
113 };
114
115 /* log(sqrt(2pi)) */
116
117 const DBL LS2PI = 0.91893853320467274178;
118
119 /* sqrt(2pi) */
120
121 const DBL s2pi = 2.50662827463100050242E0;
122
123 /* approximation for 0 <= |y - 0.5| <= 3/8 */
124
125 const DBL P0[5] =
126 {
127 -5.99633501014107895267E1,
128 9.80010754185999661536E1,
129 -5.66762857469070293439E1,
130 1.39312609387279679503E1,
131 -1.23916583867381258016E0,
132 };
133
134 const DBL Q0[8] =
135 {
136 /* 1.00000000000000000000E0,*/
137 1.95448858338141759834E0,
138 4.67627912898881538453E0,
139 8.63602421390890590575E1,
140 -2.25462687854119370527E2,
141 2.00260212380060660359E2,
142 -8.20372256168333339912E1,
143 1.59056225126211695515E1,
144 -1.18331621121330003142E0,
145 };
146
147 /*
148 * Approximation for interval z = sqrt(-2 log y ) between 2 and 8
149 * i.e., y between exp(-2) = .135 and exp(-32) = 1.27e-14.
150 */
151
152 const DBL P1[9] =
153 {
154 4.05544892305962419923E0,
155 3.15251094599893866154E1,
156 5.71628192246421288162E1,
157 4.40805073893200834700E1,
158 1.46849561928858024014E1,
159 2.18663306850790267539E0,
160 -1.40256079171354495875E-1,
161 -3.50424626827848203418E-2,
162 -8.57456785154685413611E-4,
163 };
164
165 const DBL Q1[8] =
166 {
167 /* 1.00000000000000000000E0,*/
168 1.57799883256466749731E1,
169 4.53907635128879210584E1,
170 4.13172038254672030440E1,
171 1.50425385692907503408E1,
172 2.50464946208309415979E0,
173 -1.42182922854787788574E-1,
174 -3.80806407691578277194E-2,
175 -9.33259480895457427372E-4,
176 };
177
178 /*
179 * Approximation for interval z = sqrt(-2 log y ) between 8 and 64
180 * i.e., y between exp(-32) = 1.27e-14 and exp(-2048) = 3.67e-890.
181 */
182
183 const DBL P2[9] =
184 {
185 3.23774891776946035970E0,
186 6.91522889068984211695E0,
187 3.93881025292474443415E0,
188 1.33303460815807542389E0,
189 2.01485389549179081538E-1,
190 1.23716634817820021358E-2,
191 3.01581553508235416007E-4,
192 2.65806974686737550832E-6,
193 6.23974539184983293730E-9,
194 };
195
196 const DBL Q2[8] =
197 {
198 /* 1.00000000000000000000E0,*/
199 6.02427039364742014255E0,
200 3.67983563856160859403E0,
201 1.37702099489081330271E0,
202 2.16236993594496635890E-1,
203 1.34204006088543189037E-2,
204 3.28014464682127739104E-4,
205 2.89247864745380683936E-6,
206 6.79019408009981274425E-9,
207 };
208
209
210 /*****************************************************************************
211 * Static functions
212 ******************************************************************************/
213
214 static DBL igami (DBL a, DBL y0);
215 static DBL lgam (DBL x);
216 static DBL polevl (DBL x, const DBL * coef, int N);
217 static DBL p1evl (DBL x, const DBL * coef, int N);
218 static DBL igamc (DBL a, DBL x);
219 static DBL igam (DBL a, DBL x);
220 static DBL ndtri (DBL y0);
221
222
223
224 /* chdtri()
225 *
226 * Inverse of complemented Chi-square distribution
227 *
228 *
229 *
230 * SYNOPSIS:
231 *
232 * DBL df, x, y, chdtri();
233 *
234 * x = chdtri( df, y );
235 *
236 *
237 *
238 *
239 * DESCRIPTION:
240 *
241 * Finds the Chi-square argument x such that the integral
242 * from x to infinity of the Chi-square density is equal
243 * to the given cumulative probability y.
244 *
245 * This is accomplished using the inverse gamma integral
246 * function and the relation
247 *
248 * x/2 = igami( df/2, y );
249 *
250 *
251 *
252 *
253 * ACCURACY:
254 *
255 * See igami.c.
256 *
257 * ERROR MESSAGES:
258 *
259 * message condition value returned
260 * chdtri domain y < 0 or y > 1 0.0
261 * v < 1
262 *
263 */
264
chdtri(DBL df,DBL y)265 DBL chdtri(DBL df, DBL y)
266 {
267 DBL x;
268
269 if ((y < 0.0) || (y > 1.0) || (df < 1.0))
270 {
271 Error("Illegal values fd=%f and y=%f in chdtri().", df, y);
272
273 return (0.0);
274 }
275
276 x = igami(0.5 * df, y);
277
278 return (2.0 * x);
279 }
280
281
282
283 /* lgam()
284 *
285 * Natural logarithm of gamma function
286 *
287 *
288 *
289 * SYNOPSIS:
290 *
291 * DBL x, y, lgam();
292 * extern int sgngam;
293 *
294 * y = lgam( x );
295 *
296 *
297 *
298 * DESCRIPTION:
299 *
300 * Returns the base e (2.718...) logarithm of the absolute
301 * value of the gamma function of the argument.
302 * The sign (+1 or -1) of the gamma function is returned in a
303 * global (extern) variable named sgngam.
304 *
305 * For arguments greater than 13, the logarithm of the gamma
306 * function is approximated by the logarithmic version of
307 * Stirling's formula using a polynomial approximation of
308 * degree 4. Arguments between -33 and +33 are reduced by
309 * recurrence to the interval [2,3] of a rational approximation.
310 * The cosecant reflection formula is employed for arguments
311 * less than -33.
312 *
313 * Arguments greater than MAXLGM return MAXNUM and an error
314 * message. MAXLGM = 2.035093e36 for DEC
315 * arithmetic or 2.556348e305 for IEEE arithmetic.
316 *
317 *
318 *
319 * ACCURACY:
320 *
321 *
322 * arithmetic domain # trials peak rms
323 * DEC 0, 3 7000 5.2e-17 1.3e-17
324 * DEC 2.718, 2.035e36 5000 3.9e-17 9.9e-18
325 * IEEE 0, 3 28000 5.4e-16 1.1e-16
326 * IEEE 2.718, 2.556e305 40000 3.5e-16 8.3e-17
327 * The error criterion was relative when the function magnitude
328 * was greater than one but absolute when it was less than one.
329 *
330 * The following test used the relative error criterion, though
331 * at certain points the relative error could be much higher than
332 * indicated.
333 * IEEE -200, -4 10000 4.8e-16 1.3e-16
334 *
335 */
336
lgam(DBL x)337 static DBL lgam(DBL x)
338 {
339 DBL p, q, w, z;
340 int i;
341
342 sgngam = 1;
343
344 if (x < -34.0)
345 {
346 q = -x;
347 w = lgam(q); /* note this modifies sgngam! */
348 p = floor(q);
349
350 if (p == q)
351 {
352 goto loverf;
353 }
354
355 #ifdef AVOID_TYPE_CONVERSION_WARNINGS_PATCH
356 i = (int)p;
357 #else
358 i = p;
359 #endif
360
361 if ((i & 1) == 0)
362 {
363 sgngam = -1;
364 }
365 else
366 {
367 sgngam = 1;
368 }
369
370 z = q - p;
371
372 if (z > 0.5)
373 {
374 p += 1.0;
375
376 z = p - q;
377 }
378
379 z = q * sin(M_PI * z);
380
381 if (z == 0.0)
382 {
383 goto loverf;
384 }
385
386 /* z = log(M_PI) - log( z ) - w;*/
387
388 z = LOGPI - log(z) - w;
389
390 return (z);
391 }
392
393 if (x < 13.0)
394 {
395 z = 1.0;
396
397 while (x >= 3.0)
398 {
399 x -= 1.0;
400
401 z *= x;
402 }
403
404 while (x < 2.0)
405 {
406 if (x == 0.0)
407 {
408 goto loverf;
409 }
410
411 z /= x;
412
413 x += 1.0;
414 }
415
416 if (z < 0.0)
417 {
418 sgngam = -1;
419
420 z = -z;
421 }
422 else
423 {
424 sgngam = 1;
425 }
426
427 if (x == 2.0)
428 {
429 return (log(z));
430 }
431
432 x -= 2.0;
433
434 p = x * polevl(x, B, 5) / p1evl(x, C, 6);
435
436 return (log(z) + p);
437 }
438
439 if (x > MAXLGM)
440 {
441 loverf:
442 /*
443 mtherr("lgam", OVERFLOW);
444 */
445 return (sgngam * MAXNUM);
446 }
447
448 q = (x - 0.5) * log(x) - x + LS2PI;
449
450 if (x > 1.0e8)
451 {
452 return (q);
453 }
454
455 p = 1.0 / (x * x);
456
457 if (x >= 1000.0)
458 {
459 q += ((7.9365079365079365079365e-4 * p -
460 2.7777777777777777777778e-3) * p +
461 0.0833333333333333333333) / x;
462 }
463 else
464 {
465 q += polevl(p, A, 4) / x;
466 }
467
468 return (q);
469 }
470
471
472
473 /* igamc()
474 *
475 * Complemented incomplete gamma integral
476 *
477 *
478 *
479 * SYNOPSIS:
480 *
481 * DBL a, x, y, igamc();
482 *
483 * y = igamc( a, x );
484 *
485 *
486 *
487 * DESCRIPTION:
488 *
489 * The function is defined by
490 *
491 *
492 * igamc(a,x) = 1 - igam(a,x)
493 *
494 * inf.
495 * -
496 * 1 | | -t a-1
497 * = ----- | e t dt.
498 * - | |
499 * | (a) -
500 * x
501 *
502 *
503 * In this implementation both arguments must be positive.
504 * The integral is evaluated by either a power series or
505 * continued fraction expansion, depending on the relative
506 * values of a and x.
507 *
508 *
509 *
510 * ACCURACY:
511 *
512 * Relative error:
513 * arithmetic domain # trials peak rms
514 * DEC 0,30 2000 2.7e-15 4.0e-16
515 * IEEE 0,30 60000 1.4e-12 6.3e-15
516 *
517 */
518
igamc(DBL a,DBL x)519 static DBL igamc(DBL a, DBL x)
520 {
521 DBL ans, c, yc, ax, y, z;
522 DBL pk, pkm1, pkm2, qk, qkm1, qkm2;
523 DBL r, t;
524
525 if ((x <= 0) || (a <= 0))
526 {
527 return (1.0);
528 }
529
530 if ((x < 1.0) || (x < a))
531 {
532 return (1.0 - igam(a, x));
533 }
534
535 ax = a * log(x) - x - lgam(a);
536
537 if (ax < -MAXLOG)
538 {
539 /*
540 mtherr("igamc", UNDERFLOW);
541 */
542 return (0.0);
543 }
544
545 ax = exp(ax);
546
547 /* continued fraction */
548
549 y = 1.0 - a;
550 z = x + y + 1.0;
551 c = 0.0;
552
553 pkm2 = 1.0;
554 qkm2 = x;
555 pkm1 = x + 1.0;
556 qkm1 = z * x;
557
558 ans = pkm1 / qkm1;
559
560 do
561 {
562 c += 1.0;
563 y += 1.0;
564 z += 2.0;
565
566 yc = y * c;
567
568 pk = pkm1 * z - pkm2 * yc;
569 qk = qkm1 * z - qkm2 * yc;
570
571 if (qk != 0)
572 {
573 r = pk / qk;
574 t = fabs((ans - r) / r);
575 ans = r;
576 }
577 else
578 {
579 t = 1.0;
580 }
581
582 pkm2 = pkm1;
583 pkm1 = pk;
584 qkm2 = qkm1;
585 qkm1 = qk;
586
587 if (fabs(pk) > BIG)
588 {
589 pkm2 /= BIG;
590 pkm1 /= BIG;
591 qkm2 /= BIG;
592 qkm1 /= BIG;
593 }
594 }
595 while (t > MACHEP);
596
597 return (ans * ax);
598 }
599
600
601
602 /* igam.c
603 *
604 * Incomplete gamma integral
605 *
606 *
607 *
608 * SYNOPSIS:
609 *
610 * DBL a, x, y, igam();
611 *
612 * y = igam( a, x );
613 *
614 *
615 *
616 * DESCRIPTION:
617 *
618 * The function is defined by
619 *
620 * x
621 * -
622 * 1 | | -t a-1
623 * igam(a,x) = ----- | e t dt.
624 * - | |
625 * | (a) -
626 * 0
627 *
628 *
629 * In this implementation both arguments must be positive.
630 * The integral is evaluated by either a power series or
631 * continued fraction expansion, depending on the relative
632 * values of a and x.
633 *
634 *
635 *
636 * ACCURACY:
637 *
638 * Relative error:
639 * arithmetic domain # trials peak rms
640 * DEC 0,30 4000 4.4e-15 6.3e-16
641 * IEEE 0,30 10000 3.6e-14 5.1e-15
642 *
643 */
644
645 /* left tail of incomplete gamma function:
646 *
647 * inf. k
648 * a -x - x
649 * x e > ----------
650 * - -
651 * k=0 | (a+k+1)
652 *
653 */
654
igam(DBL a,DBL x)655 static DBL igam(DBL a, DBL x)
656 {
657 DBL ans, ax, c, r;
658
659 if ((x <= 0) || (a <= 0))
660 {
661 return (0.0);
662 }
663
664 if ((x > 1.0) && (x > a))
665 {
666 return (1.0 - igamc(a, x));
667 }
668
669 /* Compute x**a * exp(-x) / gamma(a) */
670 ax = a * log(x) - x - lgam(a);
671
672 if (ax < -MAXLOG)
673 {
674 /*
675 mtherr("igam", UNDERFLOW);
676 */
677 return (0.0);
678 }
679
680 ax = exp(ax);
681
682 /* power series */
683 r = a;
684 c = 1.0;
685 ans = 1.0;
686
687 do
688 {
689 r += 1.0;
690 c *= x / r;
691 ans += c;
692 }
693 while (c / ans > MACHEP);
694
695 return (ans * ax / a);
696 }
697
698
699
700 /* igami()
701 *
702 * Inverse of complemented imcomplete gamma integral
703 *
704 *
705 *
706 * SYNOPSIS:
707 *
708 * DBL a, x, y, igami();
709 *
710 * x = igami( a, y );
711 *
712 *
713 *
714 * DESCRIPTION:
715 *
716 * Given y, the function finds x such that
717 *
718 * igamc( a, x ) = y.
719 *
720 * Starting with the approximate value
721 *
722 * 3
723 * x = a t
724 *
725 * where
726 *
727 * t = 1 - d - ndtri(y) sqrt(d)
728 *
729 * and
730 *
731 * d = 1/9a,
732 *
733 * the routine performs up to 10 Newton iterations to find the
734 * root of igamc(a,x) - y = 0.
735 *
736 *
737 * ACCURACY:
738 *
739 * Tested for a ranging from 0.5 to 30 and x from 0 to 0.5.
740 *
741 * Relative error:
742 * arithmetic domain # trials peak rms
743 * DEC 0,0.5 3400 8.8e-16 1.3e-16
744 * IEEE 0,0.5 10000 1.1e-14 1.0e-15
745 *
746 */
747
igami(DBL a,DBL y0)748 static DBL igami(DBL a, DBL y0)
749 {
750 DBL d, y, x0, lgm;
751 int i;
752
753 /* approximation to inverse function */
754 d = 1.0 / (9.0 * a);
755 y = (1.0 - d - ndtri(y0) * sqrt(d));
756
757 x0 = a * y * y * y;
758
759 lgm = lgam(a);
760
761 for (i = 0; i < 10; i++)
762 {
763 if (x0 <= 0.0)
764 {
765 /*
766 mtherr("igami", UNDERFLOW);
767 */
768 return (0.0);
769 }
770
771 y = igamc(a, x0);
772
773 /* compute the derivative of the function at this point */
774 d = (a - 1.0) * log(x0) - x0 - lgm;
775
776 if (d < -MAXLOG)
777 {
778 /*
779 mtherr("igami", UNDERFLOW);
780 */
781 goto done;
782 }
783
784 d = -exp(d);
785
786 /* compute the step to the next approximation of x */
787 if (d == 0.0)
788 {
789 goto done;
790 }
791
792 d = (y - y0) / d;
793
794 x0 = x0 - d;
795
796 if (i < 3)
797 {
798 continue;
799 }
800
801 if (fabs(d / x0) < 2.0 * MACHEP)
802 {
803 goto done;
804 }
805 }
806
807 done:
808
809 return (x0);
810 }
811
812
813
814 /* ndtri.c
815 *
816 * Inverse of Normal distribution function
817 *
818 *
819 *
820 * SYNOPSIS:
821 *
822 * DBL x, y, ndtri();
823 *
824 * x = ndtri( y );
825 *
826 *
827 *
828 * DESCRIPTION:
829 *
830 * Returns the argument, x, for which the area under the
831 * Gaussian probability density function (integrated from
832 * minus infinity to x) is equal to y.
833 *
834 *
835 * For small arguments 0 < y < exp(-2), the program computes
836 * z = sqrt( -2.0 * log(y) ); then the approximation is
837 * x = z - log(z)/z - (1/z) P(1/z) / Q(1/z).
838 * There are two rational functions P/Q, one for 0 < y < exp(-32)
839 * and the other for y up to exp(-2). For larger arguments,
840 * w = y - 0.5, and x/sqrt(2pi) = w + w**3 R(w**2)/S(w**2)).
841 *
842 *
843 * ACCURACY:
844 *
845 * Relative error:
846 * arithmetic domain # trials peak rms
847 * DEC 0.125, 1 5500 9.5e-17 2.1e-17
848 * DEC 6e-39, 0.135 3500 5.7e-17 1.3e-17
849 * IEEE 0.125, 1 20000 7.2e-16 1.3e-16
850 * IEEE 3e-308, 0.135 50000 4.6e-16 9.8e-17
851 *
852 *
853 * ERROR MESSAGES:
854 *
855 * message condition value returned
856 * ndtri domain x <= 0 -MAXNUM
857 * ndtri domain x >= 1 MAXNUM
858 *
859 */
860
ndtri(DBL y0)861 static DBL ndtri(DBL y0)
862 {
863 DBL x, y, z, y2, x0, x1;
864 int code;
865
866 if (y0 <= 0.0)
867 {
868 /*
869 mtherr("ndtri", DOMAIN);
870 */
871 return (-MAXNUM);
872 }
873
874 if (y0 >= 1.0)
875 {
876 /*
877 mtherr("ndtri", DOMAIN);
878 */
879 return (MAXNUM);
880 }
881
882 code = 1;
883
884 y = y0;
885
886 if (y > (1.0 - 0.13533528323661269189)) /* 0.135... = exp(-2) */
887 {
888 y = 1.0 - y;
889 code = 0;
890 }
891
892 if (y > 0.13533528323661269189)
893 {
894 y = y - 0.5;
895 y2 = y * y;
896 x = y + y * (y2 * polevl(y2, P0, 4) / p1evl(y2, Q0, 8));
897 x = x * s2pi;
898
899 return (x);
900 }
901
902 x = sqrt(-2.0 * log(y));
903 x0 = x - log(x) / x;
904
905 z = 1.0 / x;
906
907 if (x < 8.0) /* y > exp(-32) = 1.2664165549e-14 */
908 {
909 x1 = z * polevl(z, P1, 8) / p1evl(z, Q1, 8);
910 }
911 else
912 {
913 x1 = z * polevl(z, P2, 8) / p1evl(z, Q2, 8);
914 }
915
916 x = x0 - x1;
917
918 if (code != 0)
919 {
920 x = -x;
921 }
922
923 return (x);
924 }
925
926
927
928 /* polevl.c
929 * p1evl.c
930 *
931 * Evaluate polynomial
932 *
933 *
934 *
935 * SYNOPSIS:
936 *
937 * int N;
938 * DBL x, y, coef[N+1], polevl[];
939 *
940 * y = polevl( x, coef, N );
941 *
942 *
943 *
944 * DESCRIPTION:
945 *
946 * Evaluates polynomial of degree N:
947 *
948 * 2 N
949 * y = C + C x + C x +...+ C x
950 * 0 1 2 N
951 *
952 * Coefficients are stored in reverse order:
953 *
954 * coef[0] = C , ..., coef[N] = C .
955 * N 0
956 *
957 * The function p1evl() assumes that coef[N] = 1.0 and is
958 * omitted from the array. Its calling arguments are
959 * otherwise the same as polevl().
960 *
961 *
962 * SPEED:
963 *
964 * In the interest of speed, there are no checks for out
965 * of bounds arithmetic. This routine is used by most of
966 * the functions in the library. Depending on available
967 * equipment features, the user may wish to rewrite the
968 * program in microcode or assembly language.
969 *
970 */
971
polevl(DBL x,const DBL coef[],int N)972 static DBL polevl(DBL x, const DBL coef[], int N)
973 {
974 DBL ans;
975 int i;
976 DBL const *p;
977
978 p = coef;
979 ans = *p++;
980 i = N;
981
982 do
983 {
984 ans = ans * x + *p++;
985 }
986 while (--i);
987
988 return (ans);
989 }
990
991 /* p1evl() */
992 /* N
993 * Evaluate polynomial when coefficient of x is 1.0.
994 * Otherwise same as polevl.
995 */
996
p1evl(DBL x,const DBL coef[],int N)997 static DBL p1evl(DBL x, const DBL coef[], int N)
998 {
999 DBL ans;
1000 DBL const *p;
1001 int i;
1002
1003 p = coef;
1004 ans = x + *p++;
1005 i = N - 1;
1006
1007 do
1008 {
1009 ans = ans * x + *p++;
1010 }
1011 while (--i);
1012
1013 return (ans);
1014 }
1015
1016 END_POV_NAMESPACE
1017