## -*-SM-*- Set SM-mode in emacs
#
# Written by Daniel Pfenniger (PFENNIGER@obs.unige.ch)
# `A&S' is Abramowitz & Stegun
#
J0 1 # Bessel func. J0(|x|), A&S, 9.4.1,.3
  SET _x=ABS($1) SET _t=_x/3 SET _u=(_t<1) ? _t*_t : 1/(_t+1e-9)
  SET _v=_u/3
  SET $0=(_t<1) ? 1 + \
    _u*(-2.2499997+_u*( 1.2656208+_u*(-0.3163866+ \
    _u*( 0.0444479+_u*(-0.0039444+_u*0.00021))))) : \
    (0.79788456+_u*(-0.00000077+_u*(-0.0055274 +_u*(-0.00009512+ \
                _u*( 0.00137237+_u*(-0.00072805+_u*0.00014476))))))* \
    cos(_x-0.78539816+_u*(-0.04166397+_u*(-0.00003954+_u*(0.00262573+ \
      _u*(-0.00054125+_u*(-0.00029333+_u*0.00013558))))))*SQRT(_v)
  DELETE _t DELETE _u DELETE _v DELETE _x
Y0 1 # Bessel func. Y0(|x|), |x|>0, A&S, 9.4.2-3
  SET _x=ABS($1) SET _t=_x/3 SET _u=(_t<1) ? _t*_t : 1/_t
  SET $0=(_t<1) ? \
    (1+_u*(-2.2499997+_u*( 1.2656208+_u*(-0.3163866+ \
       _u*( 0.0444479 +_u*(-0.0039444 +_u*0.00021))))))*2/PI*LN(.5*_x)+ \
            0.36746691+_u*( 0.60559366+_u*(-0.74350384+_u*(0.25300117+ \
       _u*(-0.04261214+_u*( 0.00427916-_u*0.00024846))))) : \
    (0.79788456+_u*(-0.00000077+_u*(-0.0055274 +_u*(-0.00009512+ \
                _u*( 0.00137237+_u*(-0.00072805+_u*0.00014476))))))* \
    sin(_x-0.78539816+_u*(-0.04166397+_u*(-0.00003954+_u*(0.00262573+ \
      _u*(-0.00054125+_u*(-0.00029333+_u*( 0.00013558)))))))/SQRT(_x)
  DELETE _t DELETE _u DELETE _x
J1 1 # Bessel func. J1(|x|), A&S, 9.4.4,.6
  SET _x=ABS($1) SET _t=_x/3 SET _u=(_t<1) ? _t*_t : 1/(_t+1e-9)
  SET _v=_u/3
  SET $0=(_t<1) ? \
    _x*( 0.5       +_u*(-0.56249985+_u*( 0.21093573+_u*(-0.03954289+ \
    _u*( 0.00443319+_u*(-0.00031761+_u*0.00001109)))))) : \
       ( 0.79788456+_u*( 0.00000156+_u*( 0.01659667+_u*(0.00017105+ \
    _u*(-0.00249511+_u*( 0.00113653-_u*0.00020033))))))* \
   cos(_x-2.35619449+_u*(0.12499612+_u*(0.0000565+_u*(-0.00637879+ \
      _u*(0.00074348+_u*(0.00079824-_u*0.00029166))))))*SQRT(_v)
  DELETE _t DELETE _u DELETE _v DELETE _x
Y1 1 # Bessel func. Y1(|x|), |x|>0, A&S, 9.4.5-6
  SET _x=ABS($1) SET _t=_x/3 SET _u=(_t<1) ? _t*_t : 1/_t
  SET $0=(_t<1) ? \
    _x*( 0.5       +_u*(-0.56249985+_u*(0.21093573+_u*(-0.03954289+ \
    _u*( 0.00443319+_u*(-0.00031761+_u*0.00001109))))))*2/PI*LN(.5*_x)+ \
       (-0.6366198 +_u*( 0.2212091 +_u*(2.1682709 +_u*(-1.3164827 + \
    _u*( 0.3123951 +_u*(-0.0400976 +_u*0.0027873)))))) / _x : \
       ( 0.79788456+_u*( 0.00000156+_u*(0.01659667+_u*( 0.00017105+ \
    _u*(-0.00249511+_u*( 0.00113653-_u*0.00020033))))))* \
  sin(_x-2.35619449+_u*( 0.12499612+_u*(0.0000565 +_u*(-0.00637879+ \
    _u*( 0.00074348+_u*( 0.00079824-_u*0.00029166))))))/SQRT(_x)
  DELETE _t DELETE _u DELETE _x
I0 1 # Modified Bessel func. I0(|x|), A&S, 9.8.1-2
  SET _x=ABS($1) SET _t=_x/3.75 SET _u=(_t<1) ? _t*_t : 1/(_t+1e-9)
  SET _v=_u/3.75
  SET $0=(_t<1) ? 1 + \
    _u*(3.5156229 +_u*( 3.0899424 +_u*( 1.2067492+ \
    _u*(0.2659732 +_u*( 0.0360768 +_u*0.0045813))))) : \
       (0.39894228+_u*( 0.01328592+_u*( 0.00225319+_u*(-0.00157565+ \
    _u*(0.00916281+_u*(-0.02057706+_u*( 0.02635537+_u*(-0.01647633+ \
    _u*0.00392377))))))))*EXP(_x)*SQRT(_v)
  DELETE _t DELETE _u DELETE _v DELETE _x
I1 1 # Modified Bessel func. I1(|x|), A&S, 9.8.3-4
  SET _x=ABS($1) SET _t=_x/3.75 SET _u=(_t<1) ? _t*_t : 1/(_t+1e-9)
  SET _v=_u/3.75
  SET $0=(_t<1) ? \
    _x*( 0.5       +_u*( 0.87890594+_u*( 0.51498869+_u*(0.15084934+ \
    _u*( 0.02658733+_u*( 0.00301532+_u*0.00032411)))))) : \
       ( 0.39894228+_u*(-0.03988024+_u*(-0.00362018+_u*(0.00163801+ \
    _u*(-0.01031555+_u*( 0.02282967+_u*(-0.02895312+_u*(0.01787654- \
    _u*0.00420059))))))))*EXP(_x)*SQRT(_v)
  DELETE _t DELETE _u DELETE _v DELETE _x
K0 1 # Modified Bessel func. K0(|x|), |x|>0, A&S, 9.8.5-6
  SET _x=ABS($1) SET _t=0.5*_x SET _u=(_t<1) ? _t*_t : 1/_t
  SET _v=_u*0.28444444
  SET $0=(_t<1) ? -LN(_t)*( 1 + \
    _v*(3.5156229 +_v*( 3.0899424 +_v*(1.2067492+ \
    _v*(0.2659732 +_v*( 0.0360768 +_v*0.0045813))))))-0.57721566 + \
    _u*(0.4227842 +_u*( 0.23069756+_u*(0.0348859+ \
    _u*(0.00262698+_u*( 0.0001075 +_u*0.0000074))))) : \
       (1.25331414+_u*(-0.07832358+_u*(0.02189568+_u*(-0.01062446+ \
    _u*(0.00587872+_u*(-0.0025154 +_u*0.00053208))))))/(EXP(_x)*SQRT(_x))
  DELETE _t DELETE _u DELETE _v DELETE _x
K1 1 # Modified Bessel func. K1(|x|), |x|>0, A&S, 9.8.7-8
  SET _x=ABS($1) SET _t=0.5*_x SET _u=(_t<1) ? _t*_t : 1/_t
  SET _v=_u*0.28444444
  SET $0=(_t<1) ? LN(_t)*( \
    _x*(0.5        +_v*( 0.87890594+_v*( 0.51498869+_v*(0.15084934+ \
    _v*( 0.02658733+_v*( 0.00301532+_v*0.00032411))))))) + \
    (1+_u*( 0.15443144+_u*(-0.67278579+_u*(-0.18156897+ \
       _u*(-0.01919402+_u*(-0.00110404-_u*0.00004686))))))/_x : \
    (1.25331414+_u*(0.23498619+_u*(-0.0365562+_u*(0.01504268+_u*(-0.00780353+ \
    _u*(0.00325614-_u*0.00068245))))))/(EXP(_x)*SQRT(_x))
  DELETE _t DELETE _u DELETE _v DELETE _x
K 1 # Complete elliptic integral K(m), 0<=m<1, A&S 17.3.34
  SET _n=1-$1
  SET $0=1.38629436+ \
    _n*(0.09666344+_n*(0.03590092+_n*(0.03742564+_n*0.01451196)))- \
    LN(_n)*(0.5+_n*(0.12498594+_n*(0.06880249+_n*(0.03328355+_n*0.00441787))))
  DELETE _n
E 1 # Complete elliptic integral E(m), 0<=m<=1, A&S 17.3.36
  SET _n=1-$1
  SET $0=(_n<1e-9) ? 1+0*_n : \
    1+_n*(0.44325141+_n*(0.06260601+_n*(0.04757384+_n*0.01736506)))- \
    LN(_n+1e-9)*(_n*(0.24998368+_n*(0.0920018+_n*(0.04069698+_n*0.0052645))))
  DELETE _n
fresnel	1	# Calculate the Fresnel integrals C, S and F (fresnel)
  # Using the same approximations as in NAG libraries (Mark 14-16)
  # fresnel(x) returns the value of the full fresnel integral,
  # while C and S returns the sin and cos integrals respectively,
  # written by Paul Bartholdi
  set Abs_$1  = abs($1)
  set Xvsmall = 2.8e-13
  set _N_  = 0,dimen($1)-1
  # the next 3 lines select the for each element the approximation used
  # In that way, approximations are calculated only for the relevant 
  # values.
  # index 0 for very small values
  #       1     values < 3
  #       2     the larger ones
  set _N0_ = _N_ if( Abs_$1 < Xvsmall )
  set _N1_ = _N_ if( Abs_$1 >= Xvsmall && Abs_$1 < 3 )
  set _N2_ = _N_ if( Abs_$1 >= 3.0 )
  set S = 0 * $1
  set C = 0 * $1
  # for $1 very small :
  if( dimen(_N0_) > 0 ) {
     set $1_0 = $1[_N0_]
     set S[_N0_] = pi * $1_0 ** 3 / 6.0
     set C[_N0_] = $1_0 
  }
  # for small $1 < 3
  if( dimen(_N1_) > 0 ) {
     set $1_1 = $1[_N1_]
     set T = $1_1 / 3.0
     set T = T * T
     set T = 2.0 * T * T - 1.0
     set S[_N1_] = (((((((((((((((+6.39492767e-6)*T-5.03118359e-5 )     \
         *T+3.16611738e-4)*T-1.76549525e-3)*T+8.12986510e-3 )   \
         *T-2.99521169e-2)*T+8.52522455e-2)*T-1.76968509e-1 )   \
         *T+2.39974031e-1)*T-1.52262948e-1)*T-6.68325797e-2 )   \
         *T+1.58081431e-1)*T-3.37698367e-2)*T-4.96904567e-2 )*T \
         +3.79136372e-2) * $1_1 ** 3
     set C[_N1_] = ((((((((((((((((-1.30342972e-5)*T+1.09422578e-4)     \
         *T-7.41035352e-4)*T+4.45420608e-3)*T-2.22215620e-2)    \
         *T+8.92428847e-2)*T-2.78502318e-1)*T+6.36119616e-1)    \
         *T-9.40751488e-1)*T+5.74957819e-1)*T+6.50111442e-1)    \
         *T-1.37383189e+0)*T+4.46656174e-1)*T+4.96540249e-1)    \
         *T-2.53584714e-1)*T + 1.73361153e-1) * $1_1
  }
  # for larger $1 >= 3
  if( dimen(_N2_) > 0 ) {
     set $1_2 = $1[_N2_]
     set Pix2 = 0.5 * pi * $1_2 * $1_2
     set Cos = cos(Pix2)
     set Sin = sin(Pix2)
     set T = 3.0 / $1_2
     set T = T * T
     set T = 2.0 * T * T - 1.0
     set S05 = $1_2 > 0 ? 0.5 : -0.5
     set _F = (((((-5.93583900e-9)*T+4.72717282e-8)*T-5.42175825e-7)    \
         *T+1.10965120e-5)*T-5.73211031e-4)*T + 3.17724983e-1
     set _G = ((((((+6.24327451e-9)*T-3.29585381e-8)*T+2.12006168e-7)   \
         *T-2.00950571e-6)*T+2.97191981e-5)*T-8.84035394e-4)*T  \
         +1.00405168e-1
     set S[_N2_] = S05 - _F * Cos / $1_2 - Sin * _G / $1_2**3
     set C[_N2_] = S05 + _F * Sin / $1_2 - Cos * _G / $1_2**3
  }
  set $0 = 0.5 * ( ( 0.5 + S ) ** 2 + ( 0.5 + C ) ** 2 )
  foreach 2 ( T Abs_$1 _F _G Pix2 S05 Xvsmall Sin Cos ) { delete $2 }
acosh	1	# return the positive inverse cosh of $1 (>= 0)
		set $0 = ln($1 + sqrt($1**2 - 1))
asinh	1	# return the inverse sinh of $1
		set $0 = ($1 >= 0) ? \
		ln($1 + sqrt(1 + $1**2)) : -ln(-$1 + sqrt(1 + $1**2))
cosh	1	# return the cosh of $1
		set $0 = 0.5*(exp($1) + exp(-$1))
sinh	1	# return the sinh of $1
		set $0 = 0.5*(exp($1) - exp(-$1))

brentMin 55                             
		# one-dimensional search for a function minimum over the given range
		#
		# The position of the minimum is determined with an accuracy of
		#	3*SQRT_EPSILON*abs(x) + tol.
		# The function always obtains a local minimum which coincides with
		# the global one only if a function under investigation being unimodular.
		#
		# Algorithm
		#	G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical
		#	computations. M., Mir, 1980, p.202 of the Russian edition
		#
		#	The function makes use of the "gold section" procedure combined with
		#	the parabolic interpolation.
		#	At every step program operates three abscissae - x,v, and w.
		#	x - the last and the best approximation to the minimum location,
		#	    i.e. f(x) <= f(a) or/and f(x) <= f(b)
		#	    (if the function f has a local minimum in (a,b), then the both
		#	    conditions are fulfiled after one or two steps).
		#	v,w are previous approximations to the minimum location. They may
		#	coincide with a, b, or x (although the algorithm tries to make all
		#	u, v, and w distinct). Points x, v, and w are used to construct
		#	interpolating parabola whose minimum will be treated as a new
		#	approximation to the minimum location if the former falls within
		#	[a,b] and reduces the range enveloping minimum more efficient than
		#	the gold section procedure. 
		#	When f(x) has a second derivative positive at the minimum location
		#	(not coinciding with a or b) the procedure converges superlinearly
		#	at a rate order about 1.324
		#
		local set a = $1        # Left border | of the range
		local set b = $2        # Right border| the min is sought
		local define func $3    # Function under investigation
		local set niter = $4    # max number of iterations
		local set tol = $5      # Acceptable tolerance	for minimum, must be > 0
		
		set x local  set v local  set w local # Abscissae, descr. see above
		set fx local            # func(x)
		set fv local            # func(v)
		set fw local            # func(w)
		local set r = (3 - sqrt(5))/2 # Gold section ratio
		
		if (tol <= 0) {
		   user abort "tol == $!(tol) is <= 0"
		}
		if (b < a) {
		   user abort "(b = $!(b)) < (a == $!(a))"
		   
		}
		
		set v = a + r*(b-a) set  fv = $func(v) # First step - always gold section
		set x = v  set w = v
		set fx = fv  set fw = fv
		
		local define SQRT_EPSILON 1.49012e-08
		
		set rng local   set middle_rng local
		set tol_act local  set new_step local
		set p local             # Interpolation step is calculated as p/q;
		set q local             # division is delayed until the last moment
		set t local  set ft local
		
		define i local
		do i=0, niter - 1 {     # Main iteration loop
		   set rng = b - a      # Range in which the minimum is sought
		   set middle_rng = (a + b)/2
		   set tol_act = $SQRT_EPSILON*abs(x) + tol/3 # Actual tolerance
		   
		   if(abs(x - middle_rng) + rng/2 <= 2*tol_act) {
		      set $0 = x        # Acceptable approx. is found
		      return
		   }
		   
		   # Obtain the gold section step
		   set new_step = r*(x < middle_rng ? b - x : a - x)
		   
		   # Decide if the interpolation can be tried
		   if(abs(x-w) >= tol_act) {	# If x and w are distinct, try interpolation
		      set t = (x - w)*(fx - fv)
		      set q = (x - v)*(fx - fw)
		      set p = (x - v)*q - (x - w)*t
		      set q = 2*(q - t)
		      
		      if(q > 0) {       # make q positive and assign
		         set p = -p     # possible minus to p
		      } else {				
		         set q = -q
		      }
		      #
		      # If x+p/q falls in [a,b] not too close to a and b, and isn't too large,
		      # it is accepted.
		      #
		      # If p/q is too large then the gold section procedure can reduce [a,b]
		      # range to more extent
		      if(abs(p) < abs(new_step*q) && p > q*(a - x + 2*tol_act) && p < q*(b - x - 2*tol_act)) {
		         set new_step = p/q
		      }
		   }
		   
		   if(abs(new_step) < tol_act) {	# Adjust the step to be not less than tolerance
		      set new_step = (new_step > 0) ? tol_act : -tol_act
		   }
		   #      
		   # Obtain the next approximation to min and reduce the enveloping range
		   #
		   set t = x + new_step # guess for the min
		   set ft = $func(t)
		   if(ft <= fx) {       # t is a better approximation
		      if(t < x) {       # Reduce the range so that
		         set b = x      # t falls within it
		      } else {
		         set a = x
		      }
		      
		      set v = w   set w = x    set x = t # Assign the best approx to x
		      set fv = fw set fw = fx  set fx = ft
		   } else {             # x remains the better approx
		      if(t < x) {       # Reduce the range enclosing x
		         set a = t
		      } else {
		         set b = t
		      }
		      
		      if(ft <= fw || w == x) {
		         set v = w    set w = t
		         set fv = fw  set fw = ft
		      } else { if(ft <= fv || v == x || v == w) {
		         set v = t
		         set fv = ft
		      }}
		   }
		}
		
		set $0 = x
		#
test_brent	
		macro test 11 {
		   local set x = $1
		   set $0 = -cos(x) + (x - 1)**2
		}
   
		set x=0,2,0.001
		set f = test(x)
		lim x f box con x f 
		ctype red
		rel $(brentMin(x[0], x[dimen(x)-1], test, 100, 1e-5)) $fy1 draw $uxp $fy2
		ctype 0

		macro test DELETE
