#include <stdio.h>
#include <math.h>
#include "graphic.h"

#define	BIT15(x)	(int)((x)*(1<<15))/(float)(1<<15)

float pow15(float x, int n)
{
	int	i;

	BIT15(x);
	for (i=0; i<n; i++)
		x = BIT15(x*x);

	return (x);
}

/*
 * pow01()
 *
 * This function approximates the pow() function for the
 * range:   0 <= x <= 1  and 0 <= n <= PHONG_POWER_MAX
 * which is just fine for shading computations. The method
 * used is a Bernstein polynomial (think Bezier curves).
 *
 * Steve, Mon Jan 24 15:50:33 PST 1994
 */
#define PHONG_POWER_MAX		256
#define PHONG_POWER_MAX_LOG 	8
float
pow01(float x, int n)
{
    float t, t_2, t_3, t1, p3y;
    float s, y;
    int i;

    /* clamp inputs; be anal-retentive: */
    n = (n > PHONG_POWER_MAX) ? PHONG_POWER_MAX : n;
    n = (n < 0) ? 0 : n;

    /* simulate non-linear input domain. This is necessary
     * because the curve we are using to approximate the function
     * is parametized in "t" space, not in x space.
     * Need to find a way around this step...
     */
    s = 1.0;
    i = n;
    while ( i > 0) {
	s *= x;
	i -= 32;	/* no more than about 8 multiplies... */
    }

    /* compute Bernstein polynomial (partial): */
    p3y = (float) ((PHONG_POWER_MAX  - n) >> PHONG_POWER_MAX_LOG);
    t1 = 1.0 - s;
    t_2 = s * s;
    t_3 = t_2 * s;
    y = 3*t_2*t1*p3y + t_3;

    /* 5 multiplies, 3 adds, and a shift */

    return(y);
}

extern int en_powmul;
extern float pow_val;
extern float fog_val;

do_lighting(v, l)
struct Vertex	*v;
struct Light	*l;
{
	float	dc;
	float	sc;
	float	r[3];
	float	k;
	float	fog;

	/* diffuse cosine */
	dc = v->nx * l->nx + v->ny * l->ny + v->nz * l->nz;
	dc = clamp(0.0, dc, 1.0);

	/* specular cosine */
	sc = v->nx * l->bx + v->ny * l->by + v->nz * l->bz;
	sc = clamp(0.0, sc, 1.0);
	if (en_powmul)
#ifdef 0
		sc = pow15(sc, 6);	/* hardwire x**64 */
#endif
		sc = pow01(sc, (int) pow_val);
	else
		sc = powf(sc, (int) pow_val);

	v->r = (l->ra + dc * l->rd + sc * l->rs) * 255;
	v->g = (l->ga + dc * l->gd + sc * l->gs) * 255;
	v->b = (l->ba + dc * l->bd + sc * l->bs) * 255;
/*
	v->a = (l->aa + dc * l->ad + sc * l->as) * 255;
*/
	/* fog */
	if (fog_val > 0.0) {
	    fog = ((float)(v->z)/32767.0)*fog_val;
	    fog = 1.0/(float)exp(fog*fog);
	    v->r = fog*v->r + (1.0 - fog)*255;
	    v->g = fog*v->g + (1.0 - fog)*255;
	    v->b = fog*v->b + (1.0 - fog)*255;
	    v->a = fog*v->a + (1.0 - fog)*255;
	}

	v->r = clamp(0, v->r, 255);
	v->g = clamp(0, v->g, 255);
	v->b = clamp(0, v->b, 255);
	v->a = clamp(0, v->a, 255);

/*
printf("lx %4.2f ly %4.2f lz %4.2f vx %4.2f vy %4.2f vz %4.2f\n",
	l->nx, l->ny, l->nz, v->nx, v->ny, v->nz);
*/
}

#ifdef PWL	/* old piece wise linear code */
#define	MAX_PWL	17

int	pow_n;
float	pow_x[MAX_PWL];
float	pow_b[MAX_PWL];
float	pow_m[MAX_PWL];

float pwl_pow(float x, float n)
{
	int	i;
	float	y;

	for (i=0; i<pow_n; i++) {
		if (x >= pow_x[i] && x <= pow_x[i+1]) {
			y = pow_m[i] * (x - pow_x[i]) + pow_b[i];
		}
	}

	return (y);
}

pow_init(int p, int n, float x[])
{
	int	i;

	if (p > MAX_PWL)
		exit(EXIT_FAILURE);

	pow_n = p-1;

	for (i=0; i<p; i++) {
		pow_x[i] = x[i];
		pow_b[i] = powf(x[i], n);
	}

	for (i=0; i<pow_n; i++) {
		pow_m[i] = (pow_b[i+1] - pow_b[i]) / (pow_x[i+1] - pow_x[i]);
/*
printf("%6.4f %6.4f %6.4f\n", pow_x[i], pow_b[i], pow_m[i]);
*/
	}

}
#endif