-- FILE:     MATH.E
-- PURPOSE:  Math routines for Euphoria 3.1.1.
-- AUTHOR:   Shian Lee
-- VERSION:  1.02  September/21/2019
-- LANGUAGE: Euphoria version 3.1.1 (http://www.RapidEuphoria.com)
-- PLATFORM: Any
-- LICENCE:  Free. Use at your own risk.
-- NOTE:     * See the attached reference manuals MATH.DOC and MACHINE2.DOC,
--             and also Euphoria 3.1.1 LIBRARY.DOC, for details.
--           * Created using QE 2.3.9 & EDU 2.33 (RapidEuphoria311.com archive).
-- HISTORY:  v1.00 - initial version
--           v1.01 - new function shuffle().
--           v1.02 - new function percent().
-----------------------------------------------------------------------------



--bookmark: Math Constants;

-- CREDIT: Euphoria 4.1.0 STD, http://www.openeuphoria.org/.
-- NOTE:   PI is already defined as a global constant in misc.e.

constant
    PI = 3.14159265358979323846


global constant
    -- PINF - the positive infinity value inf (1E308 * 1000)
    PINF      = 1E308 * 1000,

    -- MINF - the negative infinity value -inf (-1E308 * 1000)
    MINF      = -PINF,

    -- MNAN - the negative not-a-number value -nan (PINF / MINF)
 -- MNAN      = PINF / MINF,                -- (tested only on DOS platform)

    -- PNAN - the positive not-a-number value nan (-(PINF / MINF))
 -- PNAN      = -MNAN,                      -- (tested only on DOS platform)

    -- E - the base of the natural logarithm Euler (e)
    E         = 2.71828182845904523536,

    -- PHI - the golden ratio phi (1+sqrt(5))/2
    PHI       = 1.61803398874989484820,

    -- EULER_GAMMA - the value Euler Gamma
    EULER_GAMMA = 0.57721566490153286061,

    -- HALFSQRT2 - the value sqrt(2)/2
    HALFSQRT2 = 0.70710678118654752440,

    -- SQRT2 - the value sqrt(2)
    SQRT2     = 1.41421356237309504880,

    -- SQRT3 - the value sqrt(3)
    SQRT3     = 1.73205080756887729353,

    -- SQRT5 - the value sqrt(5)
    SQRT5     = 2.23606797749978969641,

    -- SQRTE - the value sqrt(E)
    SQRTE     = 1.64872127070012814684,

    -- LN2 - the value ln(2) (log(2))
    LN2       = 0.69314718055994530941,

    -- LN10 - the value ln(10) (log(10))
    LN10      = 2.30258509299404568401,

    -- INVLN2 - the value 1/ln(2) (1/log(2))
    INVLN2    = 1.44269504088896340736,

    -- INVLN10 - the value 1/ln(10) (1/log(10))
    INVLN10   = 0.43429448190325182765,

    -- QUARTPI - the value PI/4
    QUARTPI   = 0.78539816339744830962,

    -- HALFPI - the value PI/2
    HALFPI    = 1.57079632679489661923,

    -- TWOPI - the value PI*2
    TWOPI     = 6.28318530717958647692,

    -- PISQR - the value power(PI,2)
    PISQR     = 9.86960440108935861883,

    -- INVSQ2PI - the value 1/sqrt(PI*2)
    INVSQ2PI  = 0.39894228040143367794,

    -- DEG_TO_RAD - conversion factor degrees to radians PI/180
    DEG_TO_RAD  = 0.01745329251994329576,

    -- RAD_TO_DEG - conversion factor radians to degrees 180/PI
    RAD_TO_DEG  = 57.29577951308232090712



--bookmark: Local Types;

-- PURPOSE: maximum long integer value with 15 decimal digits of accuracy.
--          used by prime_range() routine.
constant
    MAX_LONG_INTEGER = 999999999999999


-- PURPOSE: test if an object is <= MAX_LONG_INTEGER
type max_long_integer(object x)
    return equal(x <= MAX_LONG_INTEGER, x or 1)
end type


-- PURPOSE: test if an object does not contain a fractional part
type long_integers(object x)
    return equal(remainder(x, 1), x and 0)
end type


-- PURPOSE: test if an object is >= 0
type greater_or_equal_0(object x)
    return equal(x >= 0, x or 1)
end type


-- PURPOSE: test if an object is >= 1
type greater_or_equal_1(object x)
    return equal(x >= 1, x or 1)
end type


-- PURPOSE: test if an object != 0
type not_equal_zero(object x)
    return equal(x != 0, x or 1)
end type


-- PURPOSE: test if an object is <= -1 or >= 1
type acsc_asec_domain(object x)
    return equal(x <= -1 or x >= 1, x or 1)
end type


-- PURPOSE: test if an object is > -1 and < 1
type atanh_domain(object x)
    return equal(x > -1 and x < 1, x or 1)
end type


-- PURPOSE: test if an object is > 0 and <= 1
type asech_domain(object x)
    return equal(x > 0 and x <= 1, x or 1)
end type


-- PURPOSE: test if an object is < -1 or > 1
type acoth_domain(object x)
    return equal(x < -1 or x > 1, x or 1)
end type


-- PURPOSE: test if integer is within range 1 to 3
type choice123(integer i)
    return i >= 1 and i <= 3
end type



--bookmark: Comparisons;

-- NOTE: sign(), psign(), nsign(), zsign() are in machine2.e.


-- PURPOSE: return a value indicating equality of numbers
-- RETURN:  object; 0 if x1 = x2; 1 if x1 > x2; -1 if x1 < x2.
-- EXAMPLE: x = comp({4.2, -8.5, 3}, {4, -8.5, 11})  -- x is {1, 0, -1}
global function comp(object x1, object x2)
    return (x1 > x2) - (x1 < x2)
end function



-- PURPOSE: return a value indicating approximate equality of numbers
-- DETAILS: x1 and x2 are two objects to compare;
--          range defines amount of inequality allowed (>= 0)
-- RETURN:  object; 0 if x1 = x2 (approximately); 1 if x1 > x2; -1 if x1 < x2.
-- EXAMPLE: x = comp_approx({5.9, 8, 3.2}, {5, 7.6, 4}, 0.5)  -- x is {1,0,-1}
global function comp_approx(object x1, object x2, greater_or_equal_0 range)
    return (x1 > (x2 + range)) - (x1 < (x2 - range))
end function



-- PURPOSE: return true if a number is an integer
-- RETURN:  object;
-- EXAMPLE: x = is_int({{1,3.05},{{4,5},6.1},-9}) -- x is {{1,0},{{1,1},0},1}
global function is_int(object x)
    return not remainder(x, 1)
end function



-- PURPOSE: return true if a number is an even integer
-- RETURN:  object;
-- EXAMPLE: x = is_even({{1,-2},{{4,5},6.1},9}) -- x is {{0,1},{{1,0},0},0}
global function is_even(object x)
    return not remainder(x, 2)
end function



-- PURPOSE: return true if a number is an odd integer
-- RETURN:  object;
-- EXAMPLE: x = is_odd({{-1,2},{{3,5.7},6.1},9})  -- x is {{1,0},{{1,0},0},1}
global function is_odd(object x)
    return remainder(x, 2) and not remainder(x, 1)
end function



-- NOTE: is_prime() function is in bookmark: Misc.



-- PURPOSE: return the greatest common divisor of two numbers
-- RETURN:  object; positive integer larger then 0
-- EXAMPLE: x = gcd({-3, 15, -18.5, 0, 16}, 9) -- x is {3, 3, 9, 9, 1}
global function gcd(object x1, object x2)
    atom a
    sequence r

    if atom(x1) and atom(x2) then
	-- floor(abs(x1)); floor(abs(x2)) + avoid all conflicts with zero
	x1 = floor(x1 * ((x1 > 0) - (x1 < 0)))
	x2 = floor(x2 * ((x2 > 0) - (x2 < 0))) + ((x2 = 0) * (x1 + (x1 = 0)))
	a = 1
	while a do                          -- calculate gcd
	    a = remainder(x1, x2)
	    x1 = x2
	    x2 = a
	end while
	return x1
    else -- [nested] sequences
	-- make sure that all objects are sequences of same length
	r = x1 and x2 and 0
	x1 += r
	x2 += r
	for i = 1 to length(r) do
	    r[i] = gcd(x1[i], x2[i])
	end for
	return r
    end if
end function



-- PURPOSE: return the least common multiple of two numbers
-- RETURN:  object; positive integer, or 0
-- EXAMPLE: x = lcm({-3, 15, -18.3, 0, 16.7}, 9) -- x is {9, 45, 18, 0, 144}
global function lcm(object x1, object x2)
    return                                  -- lcm = x2 / gcd(x1, x2) * x1
    floor(x2 * ((x2 > 0) - (x2 < 0))) /     -- floor(abs(x2))
    gcd(x1, x2) *                           -- (gcd floor(abs(x1,x2)) as well)
    floor(x1 * ((x1 > 0) - (x1 < 0)))       -- floor(abs(x1))
end function



--bookmark: Accumulation;

-- NOTE: sum() is in machine2.e.


-- PURPOSE: return the minimum value in a sequence
-- RETURN:  atom;
-- EXAMPLE: a = min({{7},9,{{3,4,{5,6}}}})  --> a is 3
global function min(sequence s)
    atom a, a2

    a = PINF
    for i = 1 to length(s) do
	if atom(s[i]) then
	    if s[i] < a then a = s[i] end if
	else
	    a2 = min(s[i])
	    if a2 < a then a = a2 end if
	end if
    end for
    return a
end function



-- PURPOSE: return the maximum value in a sequence
-- RETURN:  atom;
-- EXAMPLE: a = max({{7},9,{{3,4,{5,6}}}})  --> a is 9
global function max(sequence s)
    atom a, a2

    a = MINF
    for i = 1 to length(s) do
	if atom(s[i]) then
	    if s[i] > a then a = s[i] end if
	else
	    a2 = max(s[i])
	    if a2 > a then a = a2 end if
	end if
    end for
    return a
end function



-- PURPOSE: calculate the product of a sequence
-- RETURN:  atom;
-- EXAMPLE: a = product({{1},2,{{3,4,{5,6}}}})   -- a is 720 (1*2*3*4*5*6)
global function product(sequence s)
    atom a

    a = (length(s) > 0)
    for i = 1 to length(s) do
	if atom(s[i]) then
	    a *= s[i]
	elsif length(s[i]) then
	    a *= product(s[i])
	end if
    end for
    return a
end function



-- PURPOSE: calculate the average of a sequence
-- RETURN:  sequence; {sum, count, average}
-- EXAMPLE: s = average({{1},2,{{3,4,{5,6}}}})
--          -- s is {21, 6, 3.5} -- i.e. {1+2+3+4+5+6, 6, 21 / 6}
global function average(sequence s)
    atom total, count
    sequence r

    total = 0
    count = 0
    for i = 1 to length(s) do
	if atom(s[i]) then
	    total += s[i]
	    count += 1
	else
	    r = average(s[i])
	    total += r[1]
	    count += r[2]
	end if
    end for
    return {total, count, total / (count + (count = 0))}
end function



--bookmark: Rounding and Remainders;

-- NOTE: abs() and fix() are in machine2.e.


-- PURPOSE: round up to the nearest integer
-- RETURN:  object;
-- EXAMPLE: x = ceil({-15.5, 15.5})  -- x is {-15, 16}
global function ceil(object x)
    return -floor(-x)
end function



-- PURPOSE: round down or up to the nearest integer
-- NOTE:    rounds half-way up (the common method)
-- RETURN:  object;
-- EXAMPLE: x = round({-1.5, 1.5, -1.1, 1.1, -1.9, 1.9})
--           -- x is {-1, 2, -1, 1, -2, 2}
global function round(object x)
    return floor(x + 0.5)
end function



-- PURPOSE: round down or up to some precision
-- NOTE:    rounds half-way up (the common method)
-- DETAILS: d is number of fractional part digits: 0, 1, 2, 3, ...
-- RETURN:  object;
-- EXAMPLE: x = round_to(9.12645, {0, 1, 2, 3, 4})
--           -- x is {9, 9.1, 9.13, 9.126, 9.1265}
--          x = round_to(-9.12645, 4)  -- x is -9.1264
global function round_to(object x, greater_or_equal_0 d)
    return floor(x * power(10, d) + 0.5) / power(10, d)
end function



-- PURPOSE: round towards infinity to the nearest integer
-- NOTE:    to round towards zero use fix(), in machine2.e
-- RETURN:  object;
-- EXAMPLE: x = fixup({-15.5, 15.5})  -- x is {-16, 16}
global function fixup(object x)
    return floor(x) + (remainder(x, 1) > 0)
end function



-- PURPOSE: round towards infinity to the next even integer
-- RETURN:  object;
-- EXAMPLE: x = even({2.3, 2, 0, -0.5})  -- x is {4, 2, 0, -2}
global function even(object x)
    x = floor(x) + (remainder(x, 1) > 0)    -- fixup(x)
    return x + remainder(x, 2)
end function



-- PURPOSE: round towards infinity to the next odd integer
-- RETURN:  object;
-- EXAMPLE: x = odd({2, 1.2, 1, 0, -3, -3.1, -4})  -- x is {3,3,1,1,-3,-5,-5}
global function odd(object x)
    x = floor(x) + (remainder(x, 1) > 0)    -- fixup(x)
    return x + remainder(x + 1, 2)
end function



-- PURPOSE: calculate the remainder when a number is divided by another,
--          using floored division (modulus division)
-- NOTE:    modulus division example:
--          floor( 9 /  4) =  2, remainder is  1; since  9 =  4 * 2 + 1
--          floor(-9 /  4) = -3, remainder is  3; since -9 =  4 * (-3) + 3
--          floor( 9 / -4) = -3, remainder is -3; since  9 = -4 * (-3) + (-3)
--          floor(-9 / -4) =  2, remainder is -1; since -9 = -4 * 2 + (-1)
-- RETURN:  object;
-- EXAMPLE: x = mod({81, -3.5, -9, 5.5}, {8, -1.7, 2, -4})
--           -- x is {1, -0.1, 1, -2.5}
global function mod(object x1, object x2)
    return x1 - floor(x1 / x2) * x2
end function



-- PURPOSE: return the fractional part of a number
-- RETURN:  object;
-- EXAMPLE: x = frac({-15.5, 15.5})  -- x is {-0.5, 0.5}
global function frac(object x)
    return remainder(x, 1)
end function



--bookmark: Logarithms and Powers;

-- PURPOSE: calculate the base 10 logarithm
-- NOTE:    x must be greater than zero
-- RETURN:  object;
-- EXAMPLE: x = log10(12)  -- x is 1.079181246
global function log10(object x)
    return log(x) * INVLN10
end function



-- PURPOSE: calculate E (Euler (e)) raised to a specified power
-- RETURN:  object;
-- EXAMPLE: x = exp(5.4)  -- x is 221.4064162 (the exponential of 5.4)
global function exp(atom x)
    return power(E, x)
end function



-- PURPOSE: calculate the nth Fibonacci Number
-- CREDIT:  Euphoria 4.1.0 STD, http://www.openeuphoria.org/.
-- RETURN:  object;
-- EXAMPLE: x = fib({6, 9}) -- x is {8, 34} (the Fibonacci Number of {6, 9})
global function fib(long_integers x)
    return floor((power(PHI, x) / SQRT5) + 0.5)
end function



--bookmark: Random Numbers;

-- PURPOSE: generate random integer from a specified inclusive range
-- RETURN:  object;
-- EXAMPLE: x = rand_range(18, 24) -- x may be 18, 19, 20, 21, 22, 23 or 24
global function rand_range(greater_or_equal_1 min, greater_or_equal_1 max)
    return rand(max - min + 1) + min - 1
end function



-- PURPOSE: generate random floating point number (> 0 and <= 1)
-- DETAILS: x is 1 to the largest positive integer (1073741823),
--          as x is larger the resolution is higher
-- RETURN:  object;
-- EXAMPLE: rnd({10, 20, 30}) -- x may be {0.2, 0.05555555556, 0.05263157895}
global function rnd(greater_or_equal_1 x)
    return 1 / rand(x)
end function



--bookmark: Trigonometry;

-- NOTE: arcsin() and arccos() are in misc.e.


-- PURPOSE: calculate the cosecant of an angle
-- FORMULA: y = csc(x) = 1 / sin(x)
-- DOMAIN:  x >= -HALFPI and x <= HALFPI and x != 0;  (y <= -1 or y >= 1).
-- ERROR:   if x is 0 or a multiple of PI then y is infinity
-- RETURN:  object;
global function csc(object x)
    return 1 / sin(x)
end function



-- PURPOSE: calculate the secant of an angle
-- FORMULA: y = sec(x) = 1 / cos(x)
-- DOMAIN:  x >= 0 and x <= PI and x != HALFPI;  (y <= -1 or y >= 1).
-- ERROR:   if x is an odd multiple of HALFPI then y is infinity
-- RETURN:  object;
global function sec(object x)
    return 1 / cos(x)
end function



-- PURPOSE: calculate the cotangent of an angle
-- FORMULA: y = cot(x) = 1 / tan(x)
-- DOMAIN:  x > 0 and x < PI;  (y > -inf and y < +inf).
-- ERROR:   if x is 0 or a multiple of PI then y is infinity
-- RETURN:  object;
global function cot(object x)
    return 1 / tan(x)
end function



-- PURPOSE: calculate the arctangent of a ratio
-- DETAILS: x1 is numerator of the ratio; x2 is denominator of the ratio
-- RETURN:  object; which is equal to arctan(x1 / x2), except that it can
--                  handle zero denominator and is more accurate
-- EXAMPLE: x = atan2({10.5, -10.5, 2}, {3.1, 3.1, 0})
--           -- x is {1.283713958, -1.283713958, 1.570796327}
global function atan2(object x1, object x2)
    sequence r

    if atom(x1) and atom(x2) then
	if x2 > 0 then
	    return arctan(x1 / x2)
	elsif x2 < 0 then
	    if x1 < 0 then
		return arctan(x1 / x2) - PI
	    else -- x1 >= 0
		return arctan(x1 / x2) + PI
	    end if
	elsif x1 > 0 then
	    return HALFPI
	elsif x1 < 0 then
	    return -(HALFPI)
	else -- x1 = 0
	    return 0
	end if
    else -- [nested] sequence(s)
	-- make sure that all objects are sequences of same length
	r = x1 and x2 and 0
	x1 += r
	x2 += r
	for i = 1 to length(r) do
	    r[i] = atan2(x1[i], x2[i])
	end for
	return r
    end if
end function



-- PURPOSE: calculate the angle with a given cosecant
--          (y = arccsc(x) <--> x = csc(y)) -- but x should be within domain.
-- FORMULA: y = arccsc(x) = arcsin(1 / x)
-- DOMAIN:  x <= -1 or x >= 1;  (y >= -HALFPI and y <= HALFPI and y != 0).
-- ERROR:   if x > -1 and x < 1 then y is undefined
-- RETURN:  object;
global function arccsc(acsc_asec_domain x)
    object z

    z = 1 / x                               -- return arcsin(1 / x), (misc.e)
    return 2 * arctan(z / (1.0 + sqrt(1.0 - z * z)))
end function



-- PURPOSE: calculate the angle with a given secant
--          (y = arcsec(x) <--> x = sec(y)) -- but x should be within domain.
-- FORMULA: y = arcsec(x) = arccos(1 / x)
-- DOMAIN:  x <= -1 or x >= 1;  (y >= 0 and y <= PI and y != HALFPI).
-- ERROR:   if x > -1 and x < 1 then y is undefined
-- RETURN:  object;
global function arcsec(acsc_asec_domain x)
    object z

    z = 1 / x                               -- return arccos(1 / x), (misc.e)
    return HALFPI - 2 * arctan(z / (1.0 + sqrt(1.0 - z * z)))
end function



-- PURPOSE: calculate the angle with a given cotangent
--          (y = arccot(x) <--> x = cot(y)) -- but x should be within domain.
-- FORMULA: y = arccot(x) = arctan(1 / x)
-- DOMAIN:  x > -inf and x < +inf and x != 0;  (y > 0 and y < PI).
-- ERROR:   if x is 0 then y is infinity
-- RETURN:  object;
global function arccot(not_equal_zero x)
    return arctan(1 / x)
end function



-- PURPOSE: convert an angle measured in degrees to an angle measured in radians
-- RETURN:  object;
-- EXAMPLE: x = deg_to_rad(194)  -- x is 3.385938749
global function deg_to_rad(object x)
   return x * DEG_TO_RAD                    -- x * (PI / 180)
end function



-- PURPOSE: convert an angle measured in radians to an angle measured in degrees
-- RETURN:  object;
-- EXAMPLE: x = rad_to_deg(3.385938749)  -- x is 194
global function rad_to_deg(object x)
   return x * RAD_TO_DEG                    -- x * (180 / PI)
end function



--bookmark: Hyperbolic Trigonometry;

-- PURPOSE: calculate the hyperbolic sine of a number
-- FORMULA: sinh(x) = (exp(x) - exp(-x)) / 2
-- DOMAIN:  (may overflow if x is too large)
-- RETURN:  object;
global function sinh(object x)
    return (power(E, x) - power(E, -x)) / 2
end function



-- PURPOSE: calculate the hyperbolic cosine of a number
-- FORMULA: cosh(x) = (exp(x) + exp(-x)) / 2
-- DOMAIN:  (may overflow if x is too large)
-- RETURN:  object;
global function cosh(object x)
    return (power(E, x) + power(E, -x)) / 2
end function



-- PURPOSE: calculate the hyperbolic tangent of a number
-- FORMULA: tanh(x) = (exp(x) - exp(-x)) / (exp(x) + exp(-x))
--          tanh(x) = (sinh(x) / cosh(x))
-- DOMAIN:  (may overflow if x is too large)
-- RETURN:  object;
global function tanh(object x)
    return (power(E, x) - power(E, -x)) / (power(E, x) + power(E, -x))
end function



-- PURPOSE: calculate the hyperbolic cosecant of a number
-- FORMULA: csch(x) = 2 / (exp(x) - exp(-x)) = (1 / sinh(x))
-- DOMAIN:  x != 0; (may overflow if x is too large)
-- RETURN:  object;
global function csch(not_equal_zero x)
    return 2 / (power(E, x) - power(E, -x))
end function



-- PURPOSE: calculate the hyperbolic secant of a number
-- FORMULA: sech(x) = 2 / (exp(x) + exp(-x)) = (1 / cosh(x))
-- DOMAIN:  (may overflow if x is too large)
-- RETURN:  object;
global function sech(object x)
    return 2 / (power(E, x) + power(E, -x))
end function



-- PURPOSE: calculate the hyperbolic cotangent of a number
-- FORMULA: coth(x) = (exp(x) + exp(-x)) / (exp(x) - exp(-x))
--          coth(x) = cosh(x) / sinh(x) = (1 / tanh(x))
-- DOMAIN:  x != 0; (may overflow if x is too large)
-- RETURN:  object;
global function coth(not_equal_zero x)
    return (power(E, x) + power(E, -x)) / (power(E, x) - power(E, -x))
end function



-- PURPOSE: calculate the inverse hyperbolic sine of a number
--          (y = arcsinh(x) <--> x = sinh(y))
-- FORMULA: arcsinh(x) = log(x + sqrt(x * x + 1))
-- DOMAIN:  (non)
-- RETURN:  object;
global function arcsinh(object x)
    return log(x + sqrt(x * x + 1))
end function



-- PURPOSE: calculate the inverse hyperbolic cosine of a number
--          (y = arccosh(x) <--> x = cosh(y)) -- but x should be within domain.
-- FORMULA: arccosh(x) = log(x + sqrt(x * x - 1))
-- DOMAIN:  x >= 1
-- RETURN:  object;
global function arccosh(greater_or_equal_1 x)
    return log(x + sqrt(x * x - 1))
end function



-- PURPOSE: calculate the inverse hyperbolic tangent of a number
--          (y = arctanh(x) <--> x = tanh(y)) -- but x should be within domain.
-- FORMULA: arctanh(x) = log((1 + x) / (1 - x)) / 2
-- DOMAIN:  x > -1 and x < 1
-- RETURN:  object;
global function arctanh(atanh_domain x)
    return log((1 + x) / (1 - x)) / 2
end function



-- PURPOSE: calculate the inverse hyperbolic cosecant of a number
--          (y = arccsch(x) <--> x = csch(y)) -- but x should be within domain.
-- FORMULA: arccsch(x) = arcsinh(1 / x)
-- DOMAIN:  x != 0
-- RETURN:  object;
global function arccsch(not_equal_zero x)
    object z

    z = 1 / x
    return log(z + sqrt(z * z + 1))         -- return arcsinh(1 / x)
end function



-- PURPOSE: calculate the inverse hyperbolic secant of a number
--          (y = arcsech(x) <--> x = sech(y)) -- but x should be within domain.
-- FORMULA: arcsech(x) = arccosh(1 / x)
-- DOMAIN:  x > 0 and x <= 1
-- RETURN:  object;
global function arcsech(asech_domain x)
    object z

    z = 1 / x
    return log(z + sqrt(z * z - 1))         -- return arccosh(1 / x)
end function



-- PURPOSE: calculate the inverse hyperbolic cotangent of a number
--          (y = arccoth(x) <--> x = coth(y)) -- but x should be within domain.
-- FORMULA: arccoth(x) = arctanh(1 / x)
-- DOMAIN:  x < -1 or x > 1
-- RETURN:  object;
global function arccoth(acoth_domain x)
    object z

    z = 1 / x
    return log((1 + z) / (1 - z)) / 2       -- return arctanh(1 / x)
end function



--bookmark: Misc;

-- PURPOSE: calculate the percentages of an object (in 3 modes)
-- DETAILS: if i is 1: x1 = %   , x2 = total, return part
--          if i is 2: x1 = part, x2 = total, return %
--          if i is 3: x1 = part, x2 = %    , return total
-- RETURN:  object;
-- EXAMPLE: x = percent(19, {80, 20, -200}, PCT_PART)  -- x is {15.2, 3.8, -38}
global constant
	PCT_PART    = 1,
	PCT_PERCENT = 2,
	PCT_TOTAL   = 3
global function percent(object x1, object x2, choice123 i)
    if i = PCT_PART then
	return x1 * (x2 / 100)  -- part
    elsif i = PCT_PERCENT then
	return (x1 / x2) * 100  -- percentages
    elsif i = PCT_TOTAL then
	return x1 * (100 / x2)  -- total
    end if
end function



-- PURPOSE: shuffle the elements of a sequence into random order
-- RETURN:  sequence;
-- EXAMPLE: s = shuffle({1, -5.9, {19, 4}}) -- s may be {-5.9, {19, 4}, 1}
global function shuffle(sequence s)
    integer r, len
    object swap

    len = length(s)

    for i = 1 to len do
	r = rand(len)
	swap = s[i]
	s[i] = s[r]
	s[r] = swap
    end for

    return s
end function



-- PURPOSE: return a sequence of serial numbers
-- RETURN:  sequence;
-- EXAMPLE: s = serial_numbers(1, 5, {1, -0.5})
--          -- s is {{1, 2, 3, 4, 5}, {1, 0.5, 0, -0.5, -1}}
global function serial_numbers(object num, greater_or_equal_1 len, object inc)
    sequence r

    if atom(num) and atom(len) and atom(inc) then
	r = repeat(num, len)
	for i = 2 to len do
	    r[i] = r[i - 1] + inc
	end for
	return r
    else -- [nested] sequences
	-- make sure that all objects are sequences of same length
	r = num and len and inc and 0
	num += r
	len += r
	inc += r
	for i = 1 to length(r) do
	    r[i] = serial_numbers(num[i], len[i], inc[i])
	end for
	return r
    end if
end function



-- PURPOSE: return a sequence of prime numbers by specified range
-- CREDIT:  http://www.theologyweb.com/campus/showthread.php?
--           4854-A-QBasic-Program-that-calculates-a-range-of-prime-numbers
-- DETAILS: num,max are lower,upper limits for search (positive integers)
-- RETURN:  sequence;
-- EXAMPLE: s = prime_range(740, 805)
--            -- s is {743, 751, 757, 761, 769, 773, 787, 797}
--          s = prime_range({0, 100, 1000}, {0, 100, 1000} + 10)
--            -- s is {{2,3,5,7}, {101,103,107,109}, {1009}}
global function prime_range(object num, max_long_integer max)
    sequence list, r
    integer by6, by2, by2_end
    atom c, a1, a2

    -- [nested] sequences?
    if not (atom(num) and atom(max)) then
	-- make sure that all objects are sequences of same length
	r = num and max and 0
	num += r
	max += r
	for i = 1 to length(r) do
	    r[i] = prime_range(num[i], max[i])
	end for
	return r
    end if

    -- verify prime numbers range, 2..MAX_LONG_INTEGER
    num = floor(num)
    max = floor(max)
    if num < 2 then
       num = 2
    end if
    if max > MAX_LONG_INTEGER then
       max = MAX_LONG_INTEGER
    end if

    -- initialize list with 2 and 3 prime numbers
       if max < num then           return {}
    elsif num = 2 and max = 2 then return {2}
    elsif num = 2 and max > 2 then list = {2, 3}
    elsif num = 3 then             list = {3}
    else                           list = {}
    end if

    -- search for prime numbers above 3
    by6 = 1
    by2 = 1
    by2_end = 0
    while num <= max do
	-- initially search one by one, then search faster by 6 numbers
	if by6 = 1 and not remainder((num + 1) / 6, 1) then
	    by6 = 6
	    by2 = 2
	    by2_end = 2
	end if

	-- initially search for one number, then search faster for 2 numbers
	for n = num to num + by2_end by by2 do
	    c = (n + 1) / 6
	    if not remainder(c, 1) then
		for b = 1 to c do
		    a1 = (c + b) / (6 * b + 1)
		    a2 = (c - b) / (6 * b - 1)
		    if b > a1 and b > a2 then
			list = append(list, n)
			exit
		    elsif not remainder(a1, 1) or not remainder(a2, 1) then
			exit                -- not prime
		    end if
		end for
	    else
		c = (n - 1) / 6
		if not remainder(c, 1) then
		    for b = 1 to c do
			a1 = (c - b) / (6 * b + 1)
			a2 = (c + b) / (6 * b - 1)
			if b > a1 and b > a2 then
			    list = append(list, n)
			    exit
			elsif not remainder(a1, 1) or not remainder(a2, 1) then
			    exit            -- not prime
			end if
		    end for
		end if
	    end if
	end for

	num += by6                          -- main loop increment
    end while

    if length(list) and list[$] > max then  -- return list <= max limit
	return list[1..$ - 1]
    else
	return list
    end if
end function



-- PURPOSE: return true if a number is a prime integer
-- RETURN:  object;
-- EXAMPLE: x = is_prime({1,7,11.2,{1627},18,19}) -- x is {0,1,0,{1},0,1}
global function is_prime(max_long_integer x)
    sequence r

    if atom(x) and x >= 2 and not remainder(x, 1) then
	return length(prime_range(x, x))    -- 1 = prime number, else 0
    elsif atom(x) then
	return 0                            -- negative or float number
    else -- [nested] sequence(x)
	r = repeat(0, length(x))
	for i = 1 to length(x) do
	    r[i] = is_prime(x[i])
	end for
	return r
    end if
end function



-- PURPOSE: return a list of prime numbers
-- NOTE:    * fast search algorithm is variation of the Sieve of Eratosthenes
--          * prime number is a positive integer >= 2
--          * high upper limit may run your program out of memory
-- DETAILS: max is the upper limit for search (positive integer)
-- RETURN:  string;
-- EXAMPLE: s = prime_list(40)  -- s is {2,3,5,7,11,13,17,19,23,29,31,37}
global function prime_list(integer max)
    sequence mark, primes

    mark = repeat(0, max)                   -- allocate memory for search
    primes = repeat(2, max >= 2)

    for n = 3 to max by 2 do
	if not mark[n] then
	    primes = append(primes, n)
	    for i = 3 * n to max by 2 * n do
		mark[i] = 1
	    end for
	end if
    end for
    return primes
end function



-- End of file.
