/***********************************************************************
[11-Jul-2026]
***********************************************************************/

#include <errno.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

typedef double fp_t;

#define iseven(n) ((n & 1) == 0)

fp_t  /* fp_t = float, double, long double, ... */
poch(fp_t x, int n)
{
    fp_t result;

    if (isnan(x))	/* generate run-time NaN, preserve payload */
    {
        errno = ERANGE;
	(void)nan("");
	result = x;
    }
    else if (n < 0)
    {
	errno = ERANGE;
	result = nan("0xfac50000");
    }
    else if (x == (fp_t)0.0)
	result = x;        /* NB: preserve sign of zero! */
    else if (n == 0)
	result = (fp_t)1.0;
    else if (isinf(x))
	result = (x > (fp_t)0.0) ? x : (iseven(n) ? -x : x);
    else
    {
	int k;

	result = x;

	for (k = 1; k < n; ++k)
	    result *= x + (fp_t)k;
    }

    return (result);
}

int
main(int argc, char* argv[])
{
    fp_t x;
    int n;
    
    x = (fp_t)1.0625;
    for (n = -2; n < 10; ++n)
	(void)printf("%7d\t% .17g\n", n, (fp_t)poch(x, n));

    return (EXIT_SUCCESS);
}
