Saturday, June 22, 2013

Maybe

There are different approaches to the issue of not having a value to return. One idiom to deal with this in C++ is the use of boost::optional<T> or std::pair<bool, T>.
class boost::optional<T> //Discriminated-union wrapper for values.

Maybe is a polymorphic sum type with two constructors : Nothing or Just a.
Here's how Maybe is defined in Haskell.

  {- The Maybe type encapsulates an optional value. A value of type
  Maybe a either contains a value of type a (represented as Just a), or
  it is empty (represented as Nothing). Using Maybe is a good way to
  deal with errors or exceptional cases without resorting to drastic
  measures such as error.

  The Maybe type is also a monad.
  It is a simple kind of error monad, where all errors are
  represented by Nothing. -}

  data Maybe a = Nothing | Just a

  {- The maybe function takes a default value, a function, and a Maybe
  value. If the Maybe value is Nothing, the function returns the default
  value. Otherwise, it applies the function to the value inside the Just
  and returns the result. -}

  maybe :: b -> (a -> b) -> Maybe a -> b
  maybe n _ Nothing  = n
  maybe _ f (Just x) = f x

I haven't tried to compile the following OCaml yet but I think it should be roughly OK.
 type 'a option = None | Some of 'a  ;;

  let maybe n f a =
    match a with
      | None -> n
      | Some x -> f x
      ;;

Here's another variant on the Maybe monad this time in Felix. It is applied to the problem of "safe arithmetic" i.e. the usual integer arithmetic but with guards against under/overflow and division by zero.

  union success[T] =
    | Success of T
    | Failure of string
    ;

  fun str[T] (x:success[T]) =>
    match x with
      | Success ?t => "Success " + str(t)
      | Failure ?s => "Failure " + s
    endmatch
    ;

  typedef fun Fallible (t:TYPE) : TYPE => success[t] ;

  instance Monad[Fallible]
  {
    fun bind[a, b] (x:Fallible a, f: a -> Fallible b) =>
      match x with
        | Success ?a => f a
        | Failure[a] ?s => Failure[b] s
      endmatch
      ;

    fun ret[a](x:a):Fallible a => Success x ;
  }

  //Safe arithmetic.

  const INT_MAX:int requires Cxx_headers::cstdlib ;
  const INT_MIN:int requires Cxx_headers::cstdlib ;

  fun madd (x:int) (y:int) : success[int] =>
    if x > 0 and y > (INT_MAX - x) then
        Failure[int] "overflow"
    else
      Success (y + x)
    endif
    ;

  fun msub (x:int) (y:int) : success[int] =>
    if x > 0 and y < (INT_MIN + x) then
      Failure[int] "underflow"
    else
      Success (y - x)
    endif
    ;

  fun mmul (x:int) (y:int) : success[int] =>
    if x != 0 and y > (INT_MAX / x) then
      Failure[int] "overflow"
    else
      Success (y * x)
    endif
    ;

  fun mdiv (x:int) (y:int) : success[int] =>
      if (x == 0) then
          Failure[int] "attempted division by zero"
      else
        Success (y / x)
      endif
      ;

  //--
  //
  //Test.

  open Monad[Fallible] ;

  //Evalue some simple expressions.

  val zero = ret 0 ;
  val zero_over_one = bind ((Success 0), (mdiv 1)) ;
  val undefined = bind ((Success 1),(mdiv 0)) ;
  val two = bind((ret 1), (madd 1)) ;
  val two_by_one_plus_one = bind (two , (mmul 2)) ;

  println$ "zero = " + str zero ;
  println$ "1 / 0 = " + str undefined ;
  println$ "0 / 1 = " + str zero_over_one ;
  println$ "1 + 1 = " + str two ;
  println$ "2 * (1 + 1) = " + str (bind (bind((ret 1), (madd 1)) , (mmul 2))) ;
  println$ "INT_MAX - 1 = " + str (bind ((ret INT_MAX), (msub 1))) ;
  println$ "INT_MAX + 1 = " + str (bind ((ret INT_MAX), (madd 1))) ;
  println$ "INT_MIN - 1 = " + str (bind ((ret INT_MIN), (msub 1))) ;
  println$ "INT_MIN + 1 = " + str (bind ((ret INT_MIN), (madd 1))) ;

  println$ "--" ;

  //We do it again, this time using the "traditional" rshift-assign
  //syntax.

  syntax monad //Override the right shift assignment operator.
  {
    x[ssetunion_pri] := x[ssetunion_pri] ">>=" x[>ssetunion_pri] =># "`(ast_apply ,_sr (bind (,_1 ,_3)))";
  }
  open syntax monad;

  println$ "zero = " + str (ret 0) ;
  println$ "1 / 0 = " + str (ret 1 >>= mdiv 0) ;
  println$ "0 / 1 = " + str (ret 0 >>= mdiv 1) ;
  println$ "1 + 1 = " + str (ret 1 >>= madd 1) ;
  println$ "2 * (1 + 1) = " + str (ret 1 >>= madd 1 >>= mmul 2) ;
  println$ "INT_MAX = " + str (INT_MAX) ;
  println$ "INT_MAX - 1 = " + str (ret INT_MAX >>= msub 1) ;
  println$ "INT_MAX + 1 = " + str (ret INT_MAX >>= madd 1) ;
  println$ "INT_MIN = " + str (INT_MIN) ;
  println$ "INT_MIN - 1 = " + str (ret INT_MIN >>= msub 1) ;
  println$ "INT_MIN + 1 = " + str (ret INT_MIN >>= madd 1) ;
  println$ "2 * (INT_MAX/2) = " + str (ret INT_MAX >>= mdiv 2 >>= mmul 2 >>= madd 1) ; //The last one since we know INT_MAX is odd and that division will truncate.
  println$ "2 * (INT_MAX/2 + 1) = " + str (ret INT_MAX >>= mdiv 2 >>= madd 1 >>= mmul 2) ;

  //--
That last block using the <<= syntax produces (in part) the following output (the last two print statments have been truncated away -- the very last one produces an expected overflow).

Tuesday, June 4, 2013

Powerset

This is an algorithm to compute a powerset. For example, given the set [1, 2, 3] the program should compute [[1, 2, 3], [1, 2], [1, 3], [2,3], [1], [2], [3], []]. The program is written in the Felix programming language.
fun find_seq[T] (k:size) (l:list[T]) : list[list[T]] =
  {
    return
      if k > len l
      then
        list[list[T]] () //There are no subsets of length k.
      elif k == 0uz
      then
        list[list[T]] (list[T] ()) //The empty set (the one subset of length zero).
      else
        match l with
          | Cons(?x, ?xs) => 
             join 
              (map (fun (l:list[T]):list[T]=>
                           (join (list[T] x) l)) (find_seq (k - 1) xs))
              (find_seq k xs)
        endmatch
      endif
    ;
  }
  
  fun power_set[T] (lst:list[T]):list[list[T]] =
  {
    fun loop[T] (k:size) (lst:list[T]) (acc:list[list[T]]):list[list[T]] =
    {
      return
        if k == size(-1) then acc
        else loop (k - 1) lst (join acc (find_seq k lst))
        endif
      ;  
    }
  
   return loop (len lst) lst (list[list[T]] ());
  }
  
  println$ str (power_set (list (1, 2, 3)));

Sunday, June 2, 2013

Taylor Polynomials

I find this little calculation a pleasing, succinct reminder of some the (many) things I like about functional programming. It concerns a formulation of the trigonometric sin function as the first few terms of its Taylor polynomial about x = 0:
The "definition" above can be minimally and naturally "coded up" from "first principles" (I mean, using nothing more than elementary arithmetic operations) as simply as this!:

 let rec power  = 
  fun u -> fun i ->
    if i < 0 then 1.0 / (power u (-i))
    else if i = 0 then 1.0 else u * (power u (i-1))

let rec fac = fun n -> if n = 0 then 1 else n * fac (n-1)
let pow = fun i -> fun u -> power u (2 * i + 1)
let coef = fun i -> (power (-1) i)/(fac (2 * i + 1))

let rec taylor = 
  fun i -> fun n -> fun u ->
    if i = n then (coef i) * (pow i u)
    else ((coef i) * (pow i u)) + (taylor (i + 1) n u)
  
let x = 0.78539816339744830961566084581988 (*pi/4*) in 

taylor 0 3 x 
We expect 0.70710678118654752440084436210485. I can get 0.70710646957517809, with the above which is correct to 6 decimal places. I came across this little example in Kluge's "Abstract Computing Machines" by the way, a title I'm happy to recommend to those who share my interest in symbolic computation!

Saturday, June 1, 2013

Building an OCaml package

Building an OCaml package

The details relating to OCaml package construction are provided in the the OCaml System Release 4.00 documentation and users manual. Here's how I go about it. The key to building an OCaml package is to first compile the package modules with the -for-pack flag and then combine the compiled modules with the -pack flag. The procedure can be nicely wrapped up in an OMake macro, presented in this OMakefile schematic. 

make_package(P, libname, S) =
  for_pack (P, f) = 
    OCAMLFLAGS = -for-pack $(P)
    $(f).cmx $(f)$(EXT_OBJ):
    $(f).cmo:
    $(f).cmi:
  foreach (f, $(S))
    for_pack ($(P), $(f))
  section
    .DEFAULT: $(OCamlPackage $(P), $(S))
  .DEFAULT: $(OCamlLibraryInstall _,$(BIN_DIR), $(libname), $(P))

S=$(rootname $(file $(basename $(glob *.ml)) ))

make_package(Foo, foo, $(S))