library(Bessel)

### Replicate some testing "ideas" from  zqcai.f (TOMS 644 test program)

if(!require("sfsmisc")) # from >>> sfsmisc/R/relErr.R <<<
## Componentwise aka "Vectorized" relative error:
## Must not be NA/NaN unless one of the components is  ==> deal with {0, Inf, NA}
relErrV <- function(target, current, eps0 = .Machine$double.xmin) {
    n <- length(target <- as.vector(target))
    ## assert( <length current> is multiple of <length target>) :
    lc <- length(current)
    if(!n) {
	if(!lc) return(numeric()) # everything length 0
	else stop("length(target) == 0 differing from length(current)")
    } else if(!lc)
	stop("length(current) == 0 differing from length(target)")
    ## else n, lc  > 0
    if(lc %% n)
	stop("length(current) must be a multiple of length(target)")
    recycle <- (lc != n) # explicitly recycle
    R <- if(recycle)
	     target[rep(seq_len(n), length.out=lc)]
	 else
	     target # (possibly "mpfr")
    R[] <- 0
    ## use *absolute* error when target is zero {and deal with NAs}:
    t0 <- abs(target) < eps0 & !(na.t <- is.na(target))
    R[t0] <- current[t0]
    ## absolute error also when it is infinite, as (-Inf, Inf) would give NaN:
    dInf <- is.infinite(E <- current - target)
    R[dInf] <- E[dInf]
    useRE <- !dInf & !t0 & (na.t | is.na(current) | (current != target))
    R[useRE] <- (current/target)[useRE] - 1
    ## preserve {dim, dimnames, names}  from 'current' :
    if(!is.null(d <- dim(current)))
	array(R, dim=d, dimnames=dimnames(current))
    else if(!is.null(nm <- names(current)) && is.null(names(R))) # not needed for mpfr
	`names<-`(R, nm)
    else R
}
##  Generates airy functions and their derivatives from zairy
##  and zbiry and checks them against the wronskian evaluation in the
##  region -pi/3 <= arg(z) <= pi/3:
##
##              Ai(z)*Bi'(z)-Ai'(z)*Bi(z) = 1/pi.
##
##  in the remainder of the cut plane, the identities
##
##    Ai(z)  = sqrt(-z)*( J(-1/3,zr) + J(1/3,zr) )/3
##
##    Ai'(z) =        z*( J(-2/3,zr) - J(2/3,zr) )/3
##
##    Bi(z)  = i* sqrt(-z/3) *( c1* H(1/3,1,zr) - c2* H(1/3,2,zr) )/2
##
##    Bi'(z) = i*(-z)/sqrt(3)*( c2* H(2/3,1,zr) - c1* H(2/3,2,zr) )/2
##
##  are checked where  zr = (2/3)(-z)^(3/2)  with
##  c1 = exp(pi*i/6),
##  c2 = conjg(c1)    and i^2 = -1.


all.equal12  <- function(a,b) all.equal(a,b, tolerance = 1e-12)
all.equal9e9 <- function(a,b) all.equal(a,b, tolerance = 9e-9)
op <- options(digits = 4, nwarnings = 1e4)

### Wronskian  Ai(z)*Bi'(z) - Ai'(z)*Bi(z) == 1/pi --- *only* for pi/3 <= Arg(z) <= pi/3
N <- 100
I.pi <- rep.int(1/pi + 0*1i, N)
c1 <- exp(pi * 1i/6) ## = sqrt(3)/2 + i/2
c2 <- Conj(c1)
verbose <- TRUE # <-- activate by commenting the next line:
verbose <- FALSE

set.seed(101)
for(n in 1:250) { # use longtailed different random samples
    cat(".")         ## limit to not more extreme than Cauchy (df = 1)
    z <- complex(real      = rt(N, df = max(1, 1/rexp(1))),
                 imaginary = rt(N, df = max(1, 1/rexp(1))))
    ai  <- AiryA(z)
    dai <- AiryA(z, deriv=1)
    bi  <- AiryB(z)
    dbi <- AiryB(z, deriv=1)

    ## First identity: only for  |Arg(z)| <= pi/3   <==>  z[in1] :
    Lz <- abs(Arg(z)) > pi/3
    in1 <- !Lz
    print(table(in1))
    addb <- ai * dbi - dai * bi
    rE1 <- relErrV(I.pi[in1], addb[in1])
    cat("summary(|relE.id.1|):\n"); print(summary(Mod(rE1))) # sometimes shows NAs

    ## next  stopifnot() fails for large |Mod(z[in1])| where addb[] is NaN;;
    ## partly Ai, Bi is already NaN, but less extreme:  {Ai = 0,  Bi = Inf} ==> addb = NaN (0*Inf or Inf-Inf!)
    if(verbose)
      print(data.frame(z=z, Ai=ai, "Ai'"=dai, Bi=bi, "Bi'"=dbi, "|z|"= Mod(z), "|addb|"=Mod(addb))[in1,])
    if(anyNA(addb[in1])) {
        cat("==> NAs in addb[in1]  at [", deparse1(which(is.na(addb[in1]))), "]\n")
        in1 <- in1 & !is.na(addb)
    }
    stopifnot(all.equal(addb[in1], I.pi[in1], tolerance = 1e-13))

    ## The remaining checks are only valid in this z-plane "sector":
    z   <-   z[Lz]
    ai  <-  ai[Lz]
    dai <- dai[Lz]
    bi  <-  bi[Lz]
    dbi <- dbi[Lz]

    zr <- 2/3 * (-z)^(3/2)
    if(any(Lrg <- abs(Im(zr)) > 700.921)) {
        cat("n = ", n, "; Lrg zr[i=", deparse1(which(Lrg), control={}),"]\n", sep="")
        ok <- which(!Lrg)
        z  <- z [ok]
        zr <- zr[ok]
        ai  <-  ai[ok]
        dai <- dai[ok]
        bi  <-  bi[ok]
        dbi <- dbi[ok]
    }
    if(any(isInf <- !is.finite(ai) | !is.finite(bi) | !is.finite(dai) | !is.finite(dbi))) {
        io <- which(!isInf)
        ai <- ai [io]
        bi <- bi [io]
        dai <-dai[io]
        dbi <-dbi[io]
        z <-    z[io]
        zr <-  zr[io]
    }
    stopifnot(exprs = {
        ##    Ai(z)  = sqrt(-z)*( J(-1/3,zr) + J(1/3,zr) )/3
        ##    Ai'(z) =        z*( J(-2/3,zr) - J(2/3,zr) )/3
        ## needs BesselJ() for *NEGATIVE* nu :
        all.equal9e9(ai, sqrt(-z)*(BesselJ(zr, -1/3) + BesselJ(zr, 1/3))/3)
        all.equal9e9(dai,       z*(BesselJ(zr, -2/3) - BesselJ(zr, 2/3))/3)
        ##    Bi(z)  = i* sqrt(-z/3) *( c1* H(1/3,1,zr) - c2* H(1/3,2,zr) )/2
        ##    Bi'(z) = i*(-z)/sqrt(3)*( c2* H(2/3,1,zr) - c1* H(2/3,2,zr) )/2
        all.equal9e9(bi, 1i*sqrt(-z/3)*
                        (c1*BesselH(1, zr, 1/3) - c2*BesselH(2, zr, 1/3))/2)
        all.equal9e9(dbi, -1i*z/sqrt(3)*
                         (c2*BesselH(1, zr, 2/3) - c1*BesselH(2, zr, 2/3))/2)
    })
}; cat("\n")
warnings() # .... large arguments --> precision loss (of at least half machine accuracy)