Re: Proposal for 'math.e' (2007-08-23)
Juergen Luethje wrote:
>
> This is my current (2007-08-23) proposal for a "math.e" standard
> include file, according to the discussion here on EUforum.
>
> Changes compared to the previous version:
> - Changed function sum(), so that it can't only add numbers,
> but also sequences.
> - Fixed a glitch in function rect_to_polar().
>
> }}}
<eucode>
> global constant
> E = 2.7182818284590452,
> PI = 3.1415926535897932,
> EULER_GAMMA = 0.57721566490153286, -- Euler-Mascheroni constant
> LN2 = log(2),
> LN10 = log(10),
> SQRT2 = sqrt(2),
> HALF_SQRT2 = SQRT2/2.0,
> HALF_PI = PI/2.0,
> QUARTER_PI = PI/4.0,
> TWO_PI = PI*2.0,
> EULER_NORMAL = 1/sqrt(TWO_PI)
>
>
> global function log2 (object x)
> -- logarithm base 2
> -- in : (sequence of) real number(s) > 0
> -- out: (sequence of) real number(s)
> -- Note: This function returns _exact_ results for all integral
> -- powers of 2 in the half-closed interval ]0,#FFFFFFFF]
>
> if atom(x) then
> if x = #20000000 then
> return 29 -- log(x)/LN2 is imprecise in this case
> elsif x = #80000000 then
> return 31 -- log(x)/LN2 is imprecise in this case
> else
> return log(x)/LN2
> end if
> end if
>
> for i = 1 to length(x) do
> x[i] = log2(x[i])
> end for
> return x
> end function
>
> global function log10 (object x)
> -- logarithm base 10
> -- in : (sequence of) real number(s) > 0
> -- out: (sequence of) real number(s)
>
> return log(x)/LN10
> end function
>
> type positive_not_1 (object x)
> if atom(x) and x > 0 then
> return x != 1
> end if
> return 0
> end type
>
> global function logx (object x, positive_not_1 base)
> -- general logarithm function
> -- in : x : (sequence of) atom(s) > 0
> -- base: atom > 0 and != 1
> -- out : (sequence of) atom(s)
> -- Note: If x = 1 then the function returns 0 for any possible base.
>
> return log(x)/log(base)
> end function
>
> global function exp (object x)
> return power(E, x)
> end function
>
> global function sinh (object x)
> return (exp(x) - exp(-x)) / 2
> end function
>
> global function cosh (object x)
> return (exp(x) + exp(-x)) / 2
> end function
>
> global function tanh (object x)
> return sinh(x) / cosh(x)
> end function
>
> global function arcsinh (object x)
> return log(x + sqrt(x*x+1))
> end function
>
> type not_below_1 (object x)
> if atom(x) then
> return x >= 1.0
> end if
>
> for i = 1 to length(x) do
> if not not_below_1(x[i]) then
> return 0
> end if
> end for
> return 1
> end type
>
> global function arccosh (not_below_1 x)
> return log(x + sqrt(x*x-1))
> end function
>
> type abs_below_1 (object x)
> if atom(x) then
> return x > -1.0 and x < 1.0
> end if
>
> for i = 1 to length(x) do
> if not abs_below_1(x[i]) then
> return 0
> end if
> end for
> return 1
> end type
>
> global function arctanh (abs_below_1 x)
> return log((1+x)/(1-x)) / 2
> end function
>
> global function abs (object x)
> -- return the absolute value of (all elements of) x
>
> if atom(x) then
> if x < 0 then
> return -x
> else
> return x
> end if
> end if
>
> for i = 1 to length(x) do
> x[i] = abs(x[i])
> end for
> return x
> end function
>
> global function sign (object x)
> -- x < 0 ==> -1
> -- x = 0 ==> 0
> -- x > 0 ==> +1
>
> if atom(x) then
> return compare(x, 0)
> end if
>
> for i = 1 to length(x) do
> x[i] = sign(x[i])
> end for
> return x
> end function
>
> global function ceil (object x)
> -- the opposite of floor()
> -- Examples: ? ceil(3.2) --> 4
> -- ? ceil({-3.2,7,1.6}) --> {-3,7,2}
>
> return -floor(-x)
> end function
>
> type sequence_of_a_xor_s (object x)
> -- A sequence whose top-level elements are either only atoms or only
> -- sequences (or which is empty).
> integer object_type
>
> if atom(x) then
> return 0
> end if
>
> if length(x) = 0 then
> return 1
> end if
>
> object_type = atom(x[1])
> for i = 2 to length(x) do
> if object_type != atom(x[i]) then
> return 0
> end if
> end for
>
> return 1
> end type
>
As I had mentioned earlier, this type doesn't check for equal lengths of added
items. As a result, sum() may crash in the middle of its operation, and the user
has to figure out what happened. See my replacement proposal (search for the
"sequence_of_addables" keyword). The name of the type is mysterious too.
> global function sum (sequence_of_a_xor_s list)
> -- Return the sum of all elements in 'list'.
> -- Note: This does not do a recursive sum of sub-sequences.
> object ret
>
> if length(list) = 0 then
> return 0
> end if
>
> ret = list[1]
> for i = 2 to length(list) do
> ret += list[i]
> end for
> return ret
> end function
>
> constant RADIANS_TO_DEGREES = 180.0/PI
>
> global function radians_to_degrees (object x)
> -- in : (sequence of) angle(s) in radians
> -- out: (sequence of) angle(s) in degrees
>
> return x * RADIANS_TO_DEGREES
> end function
>
> constant DEGREES_TO_RADIANS = PI/180.0
>
> global function degrees_to_radians (object x)
> -- in : (sequence of) angle(s) in degrees
> -- out: (sequence of) angle(s) in radians
>
> return x * DEGREES_TO_RADIANS
> end function
>
> type trig_range (object x)
> -- values passed to arccos and arcsin must be [-1,+1]
> if atom(x) then
> return x >= -1 and x <= 1
> end if
>
> for i = 1 to length(x) do
> if not trig_range(x[i]) then
> return 0
> end if
> end for
> return 1
> end type
>
> global function arcsin (trig_range x)
> -- returns angle in radians
> return 2 * arctan(x / (1.0 + sqrt(1.0 - x * x)))
> end function
>
> global function arccos (trig_range x)
> -- returns angle in radians
> return HALF_PI - 2 * arctan(x / (1.0 + sqrt(1.0 - x * x)))
> end function
>
> type point_pol (object x)
> if sequence(x) and (length(x) = 2)
> and atom(x[1]) and (x[1] >= 0)
> and atom(x[2]) then
> return 1
> end if
> return 0
> end type
>
> global function polar_to_rect (point_pol p)
> -- convert polar coordinates to rectangular coordinates
> -- in : sequence of two atoms: {distance, angle};
> -- 'distance' must be >= 0, 'angle' is in radians
> -- out: sequence of two atoms: {x, y}
> atom distance, angle, x, y
>
> distance = p[1]
> angle = p[2]
> x = distance*cos(angle)
> y = distance*sin(angle)
> return {x, y}
> end function
>
> type point_xy (object x)
> if sequence(x) and (length(x) = 2)
> and atom(x[1])
> and atom(x[2]) then
> return 1
> end if
> return 0
> end type
>
> global function rect_to_polar (point_xy p)
> -- convert rectangular coordinates to polar coordinates
> -- in : sequence of two atoms: {x, y}
> -- out: sequence of two atoms: {distance, angle}
> -- - 'distance' is always >= 0
> -- - 'angle' is an atom that expresses radians,
> -- and is in the half-closed interval ]-PI,+PI].
> -- If 'distance' equals 0, then 'angle' is undefined.
> object angle
> atom distance, x, y
>
> x = p[1]
> y = p[2]
> distance = sqrt(x*x + y*y)
> if x > 0 then
> angle = arctan(y/x)
> elsif x < 0 then
> if y < 0 then
> angle = arctan(y/x) - PI
> else
> angle = arctan(y/x) + PI
> end if
> else
> if y < 0 then
> angle = -HALF_PI
> elsif y > 0 then
> angle = HALF_PI
> else
> angle = 0 -- The angle is undefined in this case.
> end if
> end if
> return {distance, angle}
> end function
>
> global function find_min (sequence list, integer start)
> -- Search for the minimum value in 'list', beginning at index 'start'.
> -- Return the index of that element.
> -- Notes: This does not do a recursive compare on sub-sequences.
> -- An empty sequence will cause a runtime error.
> object temp
> integer ret
>
> temp = list[start]
> ret = start
> for i = start+1 to length(list) do
> if compare(temp, list[i]) = 1 then
> temp = list[i]
> ret = i
> end if
> end for
> return ret
> end function
>
> global function find_max (sequence list, integer start)
> -- Search for the maximum value in 'list', beginning at index 'start'.
> -- Return the index of that element.
> -- Notes: This does not do a recursive compare on sub-sequences.
> -- An empty sequence will cause a runtime error.
> object temp
> integer ret
>
> temp = list[start]
> ret = start
> for i = start+1 to length(list) do
> if compare(temp, list[i]) = -1 then
> temp = list[i]
> ret = i
> end if
> end for
> return ret
> end function
>
> global function min (sequence list, integer start)
> -- Search for the minimum value in 'list', beginning at index 'start'.
> -- Return the value of that element.
> -- Notes: This does not do a recursive compare on sub-sequences.
> -- An empty sequence will cause a runtime error.
> object ret
>
> ret = list[start]
> for i = start+1 to length(list) do
> if compare(ret, list[i]) = 1 then
> ret = list[i]
> end if
> end for
> return ret
> end function
>
> global function max (sequence list, integer start)
> -- Search for the maximum value in 'list', beginning at index 'start'.
> -- Return the value of that element.
> -- Notes: This does not do a recursive compare on sub-sequences.
> -- An empty sequence will cause a runtime error.
> object ret
>
> ret = list[start]
> for i = start+1 to length(list) do
> if compare(ret, list[i]) = -1 then
> ret = list[i]
> end if
> end for
> return ret
> end function
>
> global function lesser (object x1, object x2)
> -- Return the argument with the lesser value.
> -- Note: This does not do a recursive compare on sub-sequences.
>
> if compare(x1, x2) <= 0 then
> return x1
<snip>
Also suggested, and not rejected AFAIK, were:
* a gamma(a) function - standard generalisation of the integer factorial;
* an erf(x) function (distribution function at x of the centered reduced Gauss
probability density);
* an incomplete gamma function, part_gamma(a,x) - integrate the integrand of
gamma(a) over [0,x] and normalise by gamma(a).
There is ample supply of free C and Fortran code around to compute them, and
porting is not difficult, just tedious when the format of constants is not a
little endian sequence of bytes.
By the way, if I port (a small part of) the C library I have in mind, I'll have
to port a polynom evaluator (uses Horner's method). Could be nice to make it
global, since it would be here anyway.
CChris
|
Not Categorized, Please Help
|
|
