*! version 1.0.0  07jul2005 

mata

real colvector rank1_prod(real matrix x)
{
	real scalar i, n
	real colvector p

	p = x[.,1]
	n = cols(x)
	for (i=2; i<=n; i++) p = p :* x[.,i]

	return(p)
}

real colvector rank1_Korobov(real matrix x)
{
	return(rank1_prod((1:-2:*x):^2))
}

real colvector rank1_F2(real matrix x)
{
	external f1, f2 

	return(rank1_prod(1:+f1:*(x:*(x:-1.0):+f2)))
}

real colvector rank1_F4(real matrix x)
{
	external f1, f2

	return(rank1_prod(1:+f1:*(1.0:-f2:*x:^2:*(1.0:-x):^2)))
}

real matrix rank1_mod1(real matrix a) 
{
	return(a:-floor(a))
}

real scalar rank1_Qf(real scalar z, real scalar s, real scalar p,
		pointer(real scalar function) Fa ) 
{
  	real scalar i, k1, k2, result, p2, nb
	real matrix zk
	real scalar RANK1_BLOCKSZ

	RANK1_BLOCKSZ = 50000

	zk = J(1,s,1)
	for (i=2; i<=s; i++) {
		zk[i] = mod(zk[i-1]*z,p)
	}
	p2 = floor(p/2)
	nb = floor(p2/RANK1_BLOCKSZ)+1
	result = 0
	k2 = 0
	for (i=1; i<=nb; i++) {
		k1 = k2 + 1
		k2 = min((k2+RANK1_BLOCKSZ,p2))
		result = result + sum((*Fa)(rank1_mod1(((k1::k2):/p)*zk)))
	}
	return(result)
}

real rowvector rank1coef(real scalar smax, real scalar p, 
		| real scalar alpha) 
{
	real scalar z, a, min, h, p2, result
	pointer(real scalar function) fa
	external real scalar f1, f2

	p2 = floor(p/2)
	if (missing(alpha)) {
		/* Fa using Korobov worst function */
		fa = &rank1_Korobov()
	}
	else if (alpha==2) {
		/* Fa using order 2 Bernoulli polynomial */
		f1 = 2*st_numscalar("c(pi)")^2
		f2 = 1/6
		fa = &rank1_F2()
	}
	else {
		if (alpha != 4) {
			errprintf("alpha must be 2 or 4")
			exit(198)
		}
		/* Fa using order 4 Bernoulli polynomial */
		f1 = st_numscalar("c(pi)")^4/45.0
		f2 = 30.0
		fa = &rank1_F4()
	}
	result = J(1,smax,.)
	result[1] = p
	printf("\ng_p[++i] = %g\ng_eta[i] = &(", p)
	for (s=2; s<=smax; s++) { 
		min = st_numscalar("c(maxdouble)")
  		/* Find minimum of Qf */
 		for(z=2; z<=p2; z++) {
    			h = rank1_Qf(z,s,p,fa)
    			if(h < min) {
     	 			min = h
      				a = z
    			}
  		}

  		/* Write coefficients to stdout */
		if (s<smax) printf("%g,",a)
		else printf("%g)\n",a)
		displayflush()
		result[s] = a
  	}

  	return(result)
}
 
end