
/* Some little inline code examples (Pure 0.45). */

/* The default language for inline code is C. This requires that you have an
   LLVM C compiler installed. Either clang, llvm-gcc or gcc with the dragonegg
   plugin (gcc -fplugin=dragonegg) will do. By default, the interpreter uses
   llvm-gcc; if you want clang or dragonegg instead, just set the PURE_CC
   environment variable accordingly, see the Pure manual for details. */

%<
int mygcd(int x, int y)
{
  if (y == 0)
    return x;
  else
    return mygcd(y, x%y);
}
%>

map (mygcd 25) (30..35);

/****************************************************************************/

/* C++ is supported as well, but this requires that you call into the C++ code
   using functions with C binding (extern "C"). The following example shows
   how to work with STL maps from Pure. (The same compilers are supported as
   with C, but use the PURE_CXX environment variable instead.) */

%< -*- C++ -*-

#include <pure/runtime.h>
#include <string>
#include <map>

// An STL map mapping strings to Pure expressions.

using namespace std;
typedef map<string,pure_expr*> exprmap;

// Since we can't directly deal with C++ classes in Pure, provide some C
// functions to create, destroy and manipulate these objects.

extern "C" exprmap *map_create()
{
  return new exprmap;
}

extern "C" void map_add(exprmap *m, const char *key, pure_expr *x)
{
  exprmap::iterator it = m->find(string(key));
  if (it != m->end()) pure_free(it->second);
  (*m)[key] = pure_new(x);
}

extern "C" void map_del(exprmap *m, const char *key)
{
  exprmap::iterator it = m->find(key);
  if (it != m->end()) {
    pure_free(it->second);
    m->erase(it);
  }
}

extern "C" pure_expr *map_get(exprmap *m, const char *key)
{
  exprmap::iterator it = m->find(key);
  return (it != m->end())?it->second:0;
}

extern "C" pure_expr *map_keys(exprmap *m)
{
  size_t i = 0, n = m->size();
  pure_expr **xs = new pure_expr*[n];
  for (exprmap::iterator it = m->begin(); it != m->end(); ++it)
    xs[i++] = pure_string_dup(it->first.c_str());
  pure_expr *x = pure_listv(n, xs);
  delete[] xs;
  return x;
}

extern "C" void map_destroy(exprmap *m)
{
  for (exprmap::iterator it = m->begin(); it != m->end(); ++it)
    pure_free(it->second);
  delete m;
}

%>

public bar; // constructor symbol used below

// Create the STL map and add a sentry so that it garbage-collects itself.
let m = sentry map_destroy map_create;

// Populate the map with some arbitrary Pure data.
do (\(x=>y) -> map_add m x y) ["foo"=>99, "bar"=>bar 4711L, "baz"=>1..5];

// Query the map.
map_keys m; // => ["bar","baz","foo"]
map (map_get m) (map_keys m); // => [bar 4711L,[1,2,3,4,5],99]

// Delete an element.
map_del m "foo";
map_keys m; // => ["bar","baz"]
map (map_get m) (map_keys m); // => [bar 4711L,[1,2,3,4,5]]

/****************************************************************************/

/* Fortran. This requires either llvm-gfortran (default) or gfortran
   -fplugin=dragonegg (set with the PURE_FC environment variable). */

%< -*- Fortran90 -*-
function fact(n) result(p)
  integer n, p
  p = 1
  do i = 1, n
     p = p*i
  end do
end function fact
%>

/* Note that Fortran function names have a trailing underscore, and parameters
   are always passed by reference. (That's true for llvm-gfortran; gfortran
   supports the 'value' attribute on arguments, though.) Thus the function
   above is actually to be called as 'fact_(int*)', you can verify this with
   'show fact_'. The following little wrapper function takes care of this. */

fact n::int = fact_ {n};
map fact (1..10);

/* Another Fortran example, using matrices (Gaussian elimination). This is
   basically the same as in examples/fortran.f90. */

%< -*- Fortran90 -*-
! Matrix example (Gaussian elimination). This brings the given mxn matrix A
! into (column) echelon form. The algorithm uses partial pivoting. The
! permutation of the columns is returned in the index array.

! Note that this algorithm is prepared to work on transposed matrices, as Pure
! matrices are stored in row-major order. Hence it computes a column echelon
! form. In Pure land this becomes a row echelon form which is what we want.

subroutine gauss(m, n, A, index)
  implicit none
  integer i, j, k, p, q, n, m
  double precision A(m,n), pivot, x, y
  integer index(n)
  do i = 1, n
     index(i) = i
  end do
  do i = 1, n
     ! partial pivoting
     k = 0; pivot = 0.0
     do j = i, n
        x = A(i, index(j))
        if (abs(x) > abs(pivot)) then
           k = j; pivot = x
        end if
     end do
     x = pivot
     if (abs(x) == 0.0) exit ! zero pivot, bail out
     ! the pivot column
     p = index(k)
     if (i /= k) then
        index(k) = index(i); index(i) = p
     end if
     ! normalize the pivot column
     A(:, p) = A(:, p) / x
     ! subtract multiples of the pivot column from the remaining columns
     do k = i+1, n
        q = index(k); y = A(i, q)
        A(:, q) = A(:, q) - y*A(:, p)
     end do
  end do
end subroutine gauss
%>

/* Pure matrices are passed as void* parameters here (meaning that the raw data
   pointer will be passed), enabling the Fortran routine to modify the matrix
   in-place. The following wrapper function takes care of this. */

gauss A::matrix = index, A when
  // Get the dimensions and pack the matrix into contiguous storage. This also
  // allocates a new matrix which can be modified in-place.
  n,m = dim A; A = pack A;
  // Create the index vector to be passed to the Fortran routine.
  index = imatrix n;
  // Call the Fortran routine to do all the hard work. Note that in this
  // example the Fortran routine is prepared to work on the columns of the
  // input matrix, so that we don't have to transpose inputs and outputs here.
  gauss_ {m} {n} A index;
  // Translate to zero-based indices.
  index = map pred index;
  // Apply the row permutation to the result matrix, in order to bring the
  // matrix into triangle form. (You may want to comment this out to just
  // return the result as is.)
  A = rowcat [row A (index!i) | i = 0..n-1];
end if dmatrixp A;

/* For convenience, print a double matrix in "short" format a la Octave. */

using system;
__show__ x::matrix
= strcat [printd j (x!(i,j))|i=0..n-1; j=0..m-1] + "\n"
with printd 0 = sprintf "\n%10.5f"; printd _ = sprintf "%10.5f" end
when n,m = dim x end if dmatrixp x;

/* Example: */

let x = dmatrix {2,1,-1,8; -3,-1,2,-11; -2,1,2,-3};
x; gauss x;
