NUMERIC DIGITS 10

n = 1
nums. = 0
sum = 0
DO WHILE LINES() <> 0
    PARSE LINEIN L
    IF L <> "" THEN DO
        INTERPRET "sum = sum + " L
        INTERPRET "nums." || n || "=" || L
        n = n + 1
    END
END

n = n - 1
nums.0 = n
mean = sum/n
average_deviation = 0
standard_deviation = 0
variance = 0
skew = 0
kurtosis = 0
DO I = 1 TO nums.0
    deviation = nums.I - mean
    average_deviation = average_deviation + ABS(deviation)
    variance = variance + deviation ** 2
    skew = skew + deviation ** 3
    kurtosis = kurtosis + deviation ** 4
END
average_deviation = average_deviation / n
variance = variance / (n-1)
standard_deviation = SQRT(variance)

IF variance <> 0 THEN DO
    skew = skew / (n * variance * standard_deviation)
    kurtosis = kurtosis / (n * variance * variance) - 3
END

call qqsort 1, n

mid = n%2+1
if n//2 <> 0 then do
    median = nums.mid
end 
else do
    pmid = mid - 1
    median = (nums.mid + nums.pmid) / 2
end

SAY sprintf("n:                  %d", n)
SAY sprintf("median:             %-.6f", median)
SAY sprintf("mean:               %-.6f", mean)
SAY sprintf("average_deviation:  %-.6f", average_deviation)
SAY sprintf("standard_deviation: %-.6f", standard_deviation)
SAY sprintf("variance:           %-.6f", variance)
SAY sprintf("skew:               %-.6f", skew)
SAY sprintf("kurtosis:           %-.6f", kurtosis)

EXIT

SQRT: PROCEDURE
  ARG n
        
        
  ans = n / 2                   
  prevans = n                   
  do until prevans = ans        
     prevans = ans              
     ans = ( prevans + ( n / prevans ) ) / 2
  end
return ans




SQRTOLD:
    PROCEDURE
    PARSE ARG n
    parse value format(n,,,,0) with mant "E" exp
    if exp = "" then exp = 0
    if exp//2 < > 0 then do
       mant = mant * 10
       exp = exp -1
    end
    root = 0
    do 10
       do digit = 9 by -1 to 0,
          while,
          (root + digit) ** 2 > mant
       end
       root = root + digit
       if root**2 = mant then leave
       root = root * 10
       mant = mant * 100
       exp = exp -2
    end
return root * 10**(exp/2)
    
    









































































sprintf: procedure
  argno = 1                     
  string = ""
  start = 1                     
  len = length(arg(1))

  do until(p >= len)
    s = ""
    argno = argno + 1
    p = pos("%", arg(1), start)
    if p = 0 then
    do
      p = len + 1
    end
    if substr(arg(1), p, 1) == "%" then
    do
      s = "%"
    end
    string = string || substr(arg(1), start, p - start)
    start = p + 1
    p = verify(arg(1), "%cdfsx", "M", start)
    if p = 0 then
      leave
    spec = substr(arg(1), start, p - start + 1)
    start = p + 1
    r = right(spec, 1)
    spec = delstr(spec, length(spec), 1)
    if left(spec,1) == "-" then
    do                          
      left = 1
      spec = substr(spec, 2)
    end
    else
    do
      left = 0
      spec = substr(spec, 1)
    end
    if spec \== "" then                 
      parse var spec width "." prec
    else
    do
      width = 0
      prec = 0
    end
    if \datatype(width, "W") then
      width = 0
    if \datatype(prec, "W") then
      prec = 0
    pad = " "

    select

      when r == "s" then
      do
        if width = 0 then
          width = length(arg(argno))
        if prec \= 0 then
          s = left(arg(argno), prec)     
        else
          s = arg(argno)
      end

      when r == "d" then
      do
        if width = 0 then
          width = length(arg(argno))
        s = format(arg(argno), length(arg(argno)), 0)
      end

      when r == "f" then
      do
        if arg(argno) > -1 & arg(argno) < 1 then
          pad = "0"
        parse value arg(argno) with int "." frac
        if width = 0 & prec = 0 then
        do
          d = 1
          if arg(argno) < 0 then d = 2
          width = digits() + d
          prec = digits() - (length(int)) + d - 1
        end
        if width = 0 then
          width = len - prec
        s = format(arg(argno), width, prec, 0)
      end

      when r == "x" then
      do
        if width = 0 then
          width = length(arg(argno))
        s = d2x(arg(argno))
        if prec \= 0 then
          s = left(s, prec)     
      end

      when r == "%" then
      do
        argno = argno - 1
      end

      otherwise
        nop

    end 

    if r \== "%" then
    do
      if left then
        s = left(strip(s), width, pad)      
      else
        s = right(strip(s), width, pad)
    end
    string = string || s
  end 
return string




Fast Quick sort



/*
This is a fast quick sort routine. 

Author: Ruediger Wilke 
*/
 











qqsort: procedure expose nums.

  arg lf, re

  if re -lf < 9 then
    do lf = lf to re -1

      m = lf

      do j = lf +1 to re
        if nums.j < nums.m then
          m = j
      end 

      t = nums.m; nums.m = nums.lf; nums.lf = t

    end 
    else
    do
      i = lf
      j = re
      k = (lf + re)%2
      t = nums.k

      do until i > j

        do while nums.i < t
          i = i + 1
        end 

        do while nums.j > t
          j = j - 1
        end 

        if i <= j then
        do
          xchg = nums.i
          nums.i = nums.j
          nums.j = xchg
          i = i + 1
          j = j - 1
        end 

      end 

      call qqsort lf, j
      call qqsort i, re
    end 

return