-- include graphics.e

-- This include file extends the math functions available with the following :

-- PI        This is the value 3.1415926
-- e         This is the constaint 'e'
-- deg
-- rad_deg
-- abs
-- ceiling
-- min
-- max
-- fact
-- frac 
-- dec_frac
-- fib
-- perm_n_by_r
-- mod
-- log_10
-- exp
-- sinh
-- cosh
-- tanh
-- matrix_add  This function will add or subtract to n x n matrixes
-- matrix_mult Thsi function will multiply to n x n matrixes
-----------------------------------------------------------------------------

-- and then the following (DAM)
-- sum
-- mean
-- variance
-- std

---------------------

global function sum (sequence a)
  -- returns the sum of the elements in a
  atom s
  s = 0
  for i = 1 to length(a) do
    s = s + a[i]
  end for
  return s
end function


global function mean (sequence a)
  -- returns the average of the elements in a
  return sum(a)/length(a)
end function


global function variance (sequence a)
  -- returns the sample variance of a
  atom m, s

  if length(a) = 1 then 
    return 0 -- avoid divide by zero
  end if

  s = 0
  m = mean(a)
  for i=1 to length(a) do
    s = s + power(a[i] - m, 2)
  end for
  return s / (length(a)-1)
end function


global function std (sequence a)
  -- returns the standard deviation of a 
  return sqrt(variance(a))
end function

---

global constant pi = 3.141592654
global constant e = 2.718281828


global function deg(atom x)
--syntax sin(deg(x))           
	   x = x*.0174533
	   return x
       end function
       
       
global function rad_deg(atom x)
	   x = x * 57.29578
	   return x
       end function


global function log_10(atom x)
	   x = log(x) * .4342945
	   return x
       end function


global function abs(atom x)
	   if x<0 then
	     x = (-x)
	   end if
	   return x
	end function


global function frac(atom x)
	  if abs((floor(x))) = x then 
	    return 0
	  else               
	    return x - floor(x)
	  end if
       end function


global function fact(atom x)
	if x = 1 then
	  return 1
	else
	  return ( x * fact(x -1))  
	end if
       end function


global function fib (atom n)
atom previous,current,new

    if n = 0 then
	return 0
    elsif n = 1 then
	return 1
    else 
	previous = 0
	current =1
	for i = 2 to n do
	    new = previous + current
	    previous = current
	    current = new
	end for
	return new
    end if
end function


global function min (sequence s)
    object current
    current = s[1]
    for i = 2 to length(s) do
	if current > s[i] then
	    current =s[i]
	end if
    end for
    return current
end function


global function max (sequence s)
    object current
    current = s[1]
    for i = 2 to length(s) do
	if current < s[i] then
	    current =s[i]
	end if
    end for
    return current
end function


global function dec_frac(atom x1)
atom numerator, denominator
numerator = 1
denominator = 2
    for q = denominator to 1000 do
	for r = numerator to q - 1 do
	  if abs(x1 - (r/q)) < .00003  then
	    return {r,q}
	  end if
	end for 
    end for    
end function


global function mod(atom x,atom y)
	 x = abs(x -y * abs(x/y))
	 return x
       end function
       
       
global function ceiling(atom x1,)
    if (floor(x1)< x1) then
	return (floor(x1)+1)
	else
	return x1   
    end if
end function


global function exp(atom x)
	 x = power (e,x)
	 return x
       end function
       

global function perm_n_by_r(atom n, atom r)
	  n = (fact(n))/(fact(n-r))
	  return n
       end function
	  

global procedure pascal_triangle(atom rows)
       atom c
	 for n = 0  to rows-1  do
	   for r = 0  to n   do
	     c = 1
	     for x = n  to n-r+1  by -1  do
	       c=c*x/(n-x+1)
	     end for
	     printf(1,"%d ",c)
	    end for
	    puts(1,"\n")
	  end for
	end procedure
	  

global function sinh(atom x)
	 x = .5 * (exp(x)-exp(-x))
	 return x
       end function
       
       
global function cosh(atom x)
	 x = .5 * (exp(x)+exp(-x))
	 return x
       end function
       
       
global function tanh(atom x)
	 x = 1-2*exp(-x)/(exp(x)+exp(-x))
	 return x
       end function
       

sequence c

global function matrix_add(sequence a,sequence b,atom row,atom column,atom s)
     c= repeat (repeat({0},row),column)
     for j = 1  to row   do
	for i = 1 to column do
	    if s = 2 then
		b[j][i] = -b[j][i] 
	    end if
	    c[j][i] = { a[j][i] + b[j][i]}
	end for
    end for
    return c
end function           
	   
global function matrix_mult(sequence a,sequence b,atom row,atom column)
     c= repeat (repeat({0},row),column)
     for j = 1  to row   do
	for i = 1 to column do
	    --if s = 2 then
		--b[j][i] = -b[j][i] 
	    --end if
	    c[j][i] = { a[j][i] * b[j][i]}
	end for
    end for
    return c
end function