static char sccsid[]="%Z% %M% %I% %E% %U%";
/************************************************/
/*												*/
/*	akxgmath.c									*/
/*												*/
/*		  coded by A.Kobayashi 2010.05.27		*/
/*												*/
/************************************************/
#include "akxcommon.h"

#define XDEBUG(x)

int akxgisqrt(b)
int b;
{
	int x0,d,y;

	x0=b;
	if (x0 > 0) {
		x0=x0/2;
		if (x0 > 1) {
			while (d=((y=x0*x0)-b)/(2*x0)) {
				x0-=d;
			}
			if (y>b) x0--;
		}
		else x0=1;
	}
	return x0;
}

#define _DIFF(a,b) (a>b?a-b:b-a)
#define _ABS(a) (a>=0?a:-a)

static int _b_search(sosu,max_size,l)
int  *sosu,max_size;
int l;
{
	int max=max_size;
	int ll=l;
	int i,m,i1,i2,m1;

	if (!sosu) return -1;
	else if (ll==(m1=sosu[0]) || ll==(m=sosu[max-1])) return 1;
	else if (ll>m1 && ll<m) {
		i1=0;
		i2=max-1;
		while (i2-i1>1) {
			i = (i1+i2)/2;
			m = sosu[i];

XDEBUG(XDEBUGOUTL5(0,"_b_search: ll=%d i1=%d i2=%d i=%d m=%d",ll,i1,i2,i,m);)

			if (ll == m) return 1;
			else if (ll < m) i2 = i;
			else i1 = i;
		}
		return 0;
	}
	return -2;
}

static int _h_search(sosu,max_size,l)
int  *sosu,max_size;
int l;
{
	int max=max_size;
	int ll=l;
	int i,m,i1,i2,m1,m2,d,dy,dy1,a,dx1,dx1dy,dy2;

	if (!sosu) return -1;
	else if (ll==(m1=sosu[0]) || ll==(m2=sosu[max-1])) return 1;
	else if (ll>m1 && ll<m2) {
		i1=0;
		i2=max-1;

XDEBUG(XDEBUGOUTL5(0,"_h_search: ll=%d",ll,0,0,0,0);)

		for (;;) {
			dx1 = i2 - i1;
			dy  = ll - m2;

XDEBUG(XDEBUGOUTL5(0,"    dx1=%d dy=%d",dx1,dy,0,0,0);)

			if ((_ABS(dx1)<=1) && (dy*(ll-m1))<=0) break;
			dy1 = m2 - m1;
			dx1dy = dx1 * dy;
			if (!(d=dx1dy/dy1)) {
				if (dx1dy*dy1>=0) d = 1;
				else d = - 1;
			}
			i = d + i2;
			m = sosu[i];

XDEBUG(XDEBUGOUTL5(0,"    dx1dy=%d dy1=%d d=%d P1(%d,%d)",dx1dy,dy1,d,i1,m1);)
XDEBUG(XDEBUGOUTL5(0,"            P2(%d,%d) P(%d,%d)",i2,m2,i,m,0);)

			if (ll == m) return 1;
			dy2 = m2 - m;
			if (((dy1>0) && (dy2>(-dy))) ||
			    ((dy1<0) && ((-dy2)>dy))) {
				m1 = m2;
				i1 = i2;
			}
			m2 = m;
			i2 = i;
		}
		return 0;
	}
	return -2;
}

static int  gSosuTblSize=0;
static int *gSosuTbl=NULL;

int akxg_sosu_chk_by_tbl_opt(sosu,max_size,l,opt)
int  *sosu,max_size,opt;
int l;
{
	int i,max=max_size,*ip,m;
	int ll,l1,l2;

	if (!(ip=sosu)) return -1;
	ll = l;

	if (ll < 2) return 0;
	else if (ll == 2) return 1;
	else if (ll <= 7) {
		if (ll&0x01) return 1;
		return 0;
	}

	/* search sosu table */
	if (ll==sosu[0] || ll==(m=sosu[max-1])) return 1;
	else if (ll < m) {
		if (opt) {
			if ((i=_h_search(sosu,max,ll)) >= 0) return i;
		}
		else {
			if ((i=_b_search(sosu,max,ll)) >= 0) return i;
		}
	}

	l2 = akxgisqrt(ll);
	for (i=0;i<max;i++,ip++) {
		if (!(ll % (m=*ip))) break;
		if (m > l2) break;
	}
	if (i>=max || m>l2) return 1;
	return 0;
}

int akxg_sosu_chk_by_tbl(sosu,max_size,l)
int  *sosu,max_size;
int l;
{
	return akxg_sosu_chk_by_tbl_opt(sosu,max_size,l,0);
}

#define MAX_SOSU	32
int akxg_sosu_tbl(pi)
int **pi;
{
	static int sosu[MAX_SOSU]={
	  2,  3,  5,  7, 11, 13, 17, 19, 23, 29,
	 31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
	 73, 79, 83, 89, 97,101,103,107,109,113,
	127,131};
	int *piw,size;

	if (!pi) return -1;
	if (gSosuTbl) {
		size = gSosuTblSize;
		piw = gSosuTbl;
	}
	else {
		size = MAX_SOSU;
		piw = sosu;
	}
	*pi = piw;
	return size;
}

int akxg_sosu_chk(l)
int l;
{
	int  max_size,*sosu,i;

	max_size = akxg_sosu_tbl(&sosu);
	if ((i=akxg_sosu_chk_by_tbl(sosu,max_size,l)) < 0) i = 0;
	return i;
}
