# math functions
{
	eps_1.e-7
	x_c(1,2,3,4,-1,-2)
	y_as.numeric(c(-4.5,6.03,-7.7,"2.00",3))
	T
}
is.numeric(x)
mode(y)=="numeric"
abs(4.7)==abs(4.7)
abs(-6.3)==abs(6.3)
all(approx(c(1,2,3,4),c(1,2,3.001,4),xout=c(1.5,2.5,3.5),method="linear")$x-
c(1.5,2.5,3.5))<eps
all(floor(x)==ceiling(x))
all(ceiling(x+eps)-floor(x+eps)==1)
all(trunc(x*(1+eps))==x)
all(cumsum(c(x[1],diff(x))==x))
all(range(y)==c(min(y),max(y)))
min(x,y,runif(77))==min(min(x),min(y))
max(x,y,runif(77))==max(max(x),max(y))
all(pmax(1:6,0:5,-3:2)==1:6)
all(pmin(1:6+eps,0:5+eps,-3:2+eps)==-3:2+eps)
{foo_pmin(c(1,NA,3),2); is.na(foo[2]) && all(foo[-2]==1:2)}
median(y)==quantile(y,c(.5))
mean(x,trim=.5)==median(x)
rank(x[6])==1
all(round(x*100,digits=-2)==x*100)
all(round(y)-as.integer(y)<eps)
all(signif(y,2)-y<eps)
sum(y)-max(cumsum(y))<eps
prod(x)==x[1]*x[2]*x[3]*x[4]*x[5]*x[6]
is.numeric(x_abs(x))
mode(as.numeric(y_abs(y)))=="numeric"
all(3^log(x,base=3)-x<eps)
all(10**log10(y)-y<eps)
all(sqrt(y)-(10^(1/2*log10(y)))<eps)
exp(1)-(1/exp(-1))<eps
all(exp(cumsum(log(1:10)))-gamma(2:11)<eps)
{x_runif(500);T}
all(lgamma(x)-log(gamma(x))<10*eps)
{
	cis <<- function(theta) exp(1i*theta)
	norm <<- function(z) sqrt(Re(sum(z*Conj(z)))/length(z))
	tsign <- function(z) Mod(sign(z)) - 1
	texp <- function(z) c(log(Mod(exp(z)))-Re(z),Arg(exp(z))-Arg(cis(Im(z))))
	tlog <- function(z) c(log(Mod(z))-Re(log(z)),Arg(z)-Im(log(z)))
	tpow <- function(z,w) z^w-exp(w*log(z))
	tsin <- function(z) sin(z)-(exp(1i*z)-exp(-1i*z))/2i
	tcos <- function(z) cos(z)-(exp(1i*z)+exp(-1i*z))/2
	ttan <- function(z) tan(z)-sin(z)/cos(z)
	tasin <- function(z) asin(z)+1i*log(1i*z+sqrt(1-z*z))
	tacos <- function(z) acos(z)+1i*log(z+1i*sqrt(1-z*z))
	tatan <- function(z) atan(z)+1i*log((1+(1i*z))*sqrt(1/(1+z*z)))
	tsinh <- function(z) sinh(z)-(exp(z)-exp(-z))/2
	tcosh <- function(z) cosh(z)-(exp(z)+exp(-z))/2
	ttanh <- function(z) tanh(z)-sinh(z)/cosh(z)
	tasinh <- function(z) asinh(z)-log(z+sqrt(1+z*z))
	tacosh <- function(z) acosh(z)-log(z+(z+1)*sqrt((z-1)/(z+1)))
	tatanh <- function(z) atanh(z)-log((1+z)*sqrt(1/(1-z*z)))
	n <- 500
	z1 <- rnorm(n)/runif(n) + 1i*rnorm(n)/runif(n)
	z2 <- runif(n,-20,20) + 1i*runif(n,-20,20)
	z3 <- rnorm(n) + 1i*rnorm(n)
	w3 <- rnorm(n) + 1i*rnorm(n)
	r <- tan(runif(n,0,pi/2)) + 0i
	e2 <- exp(z2)
	l1 <- log(z1)
	b1 <<- c(1+r, -1-r)
	b2 <- c(-1+0i, 1+0i)
	eps <- 5e-8
	T
}
{current.warn <- options("warn")[[1]];options(warn=-1);T}
is.finite(log(0+0i))==F
all(is.finite(atan(1i*b2))==F)
all(is.na(atanh(b2)))
{options(warn=current.warn);T}
# This is split up into separate expressions below:
# all(c(Arg(0)==0,norm(Arg(-r)-pi)<1e-12,Arg(z1)>=-pi,Arg(z1)<=pi,
# 	norm(Mod(z1)*cis(Arg(z1))-z1)<eps))
vectors.equal(Arg(0), 0)	# Arg(0)==0
vectors.equal(norm(Arg(-r)-pi), 0, tol=1e-12)	# norm(Arg(-r)-pi)<1e-12
{
	# elementwise min/max
	"Min"<- function(x, y) ifelse(x < y, x, y)
	"Max"<- function(x, y) ifelse(x > y, x, y)
	T
}
vectors.equal(Min(Arg(z1), -pi), -pi, tol=1e-15)	# all(Arg(z1)>=-pi)
vectors.equal(Max(Arg(z1),  pi),  pi, tol=1e-15)	# all(Arg(z1)<=pi)
vectors.equal(norm(Mod(z1)*cis(Arg(z1))-z1), 0, tol=eps) 
	# norm(Mod(z1)*cis(Arg(z1))-z1)<eps
all(norm(tsign(z1))<eps)
all(norm(texp(z2)<eps))
all(c(norm(tlog(z1)),norm(tlog(-r)))<eps)
all(c(norm(tpow(z3,w3))<eps,norm(0i^r)==0,0i^0i==1))
all(c(norm(tsin(z2)),norm(tcos(z2)),norm(ttan(z2)))<eps)
all(c(norm(tsin(z3)),norm(tcos(z3)),norm(ttan(z3)))<eps)
all(c(norm(tasin(z1)),norm(tacos(z1)),norm(tatan(z1)))<eps)
all(c(norm(xx<<-tasin(b1)),norm(yy<<-tacos(b1)),norm(zz<<-tatan(b1)))<eps)
all(c(norm(tasin(b2)),norm(tacos(b2)))<eps)
all(c(norm(tsinh(z2)),norm(tcosh(z2)),norm(ttanh(z2)))<eps)
all(c(norm(tsinh(z3)),norm(tcosh(z3)),norm(ttanh(z3)))<eps)
all(c(norm(tasinh(z1)),norm(tacosh(z1)),norm(tatanh(z1)))<eps)
all(c(norm(tasinh(1i*b1)),norm(tacosh(1-r)),norm(tatanh(b1)))<eps)
all(c(norm(tasinh(1i*b2)),norm(tacosh(1)))<eps)
{
	tfft <- function(z)fft(fft(z),inverse=T)/length(z)-z
	x <- runif(1000)
	z <- runif(990) + 1i*runif(990)
	a <- array(z,9:11)
	T
}
all(norm(tfft(x))<eps)
all(norm(tfft(z))<eps)
	{ta_tfft(a);T}
all(dim(a)==dim(ta))
all(norm(ta)<eps)