determinant.c

// Generate fast C code for calculating the determinant of a matrix using the
// preprocessor and the Order library, get the dependecies from
// <http://sourceforge.net/projects/chaos-pp/>
// Note the code is BROKEN.
// Damiano Stoffie 2013, released into the public domain

#include <stdio.h>
#include "order/interpreter.h"

#define ORDER_PP_DEF_8lsquare ORDER_PP_CONST([)
#define ORDER_PP_DEF_8rsquare ORDER_PP_CONST(])

#define ORDER_PP_DEF_8seq_pair_with_index \
  ORDER_PP_FN( \
    8fn(8S, \
        8seq_pair_with(8fn(8L, 8R, 8seq(8L, 8R)), \
                       8S,\
                       8seq_iota(1, 8to_lit(8inc(8seq_size(8S)))))))

#define ORDER_PP_DEF_8seq_unpair \
  ORDER_PP_FN(8fn(8S, \
                  8seq_map(8fn(8E, 8seq_head(8E)), \
                           8S)))

// X: index,
// S: size of the matrix
// I: i-row
// J: j-column
// This is really badly written ^_^
#define ORDER_PP_DEF_8index_is_row_or_column \
  ORDER_PP_FN( \
    8fn(8X, 8S, 8I, 8J, \
        8or(8equal(8if(8is_0(8remainder(8X, 8S)), \
                       8div(8X, 8S), \
                       8inc(8div(8X, 8S))), \
                   8I), \
            8equal(8if(8is_0(8remainder(8X, 8S)), \
                       8S, \
                       8remainder(8X, 8S)), \
                   8J))))

// M: sequence/matrix of size SxS
// S: size of the matrix
// I: row to remove
// J: column to remove
#define ORDER_PP_DEF_8submx \
  ORDER_PP_FN( \
    8fn(8M, 8S, 8I, 8J, \
        8seq_unpair( \
          8seq_filter(8fn(8E, \
                          8not(8index_is_row_or_column(8seq_last(8E), \
                                                       8S, 8I, 8J))), \
                      8seq_pair_with_index(8M)))))

#define ORDER_PP_DEF_8even \
  ORDER_PP_FN(8fn(8N, 8is_0(8remainder(8N, 2))))

// M: sequence/matrix of size SxS
// S: size of the matrix
// N: name of the c array/matrix
#define ORDER_PP_DEF_8determinant_r \
  ORDER_PP_FN( \
    8fn(8M, 8S, 8N, \
        8if(8equal(8S, 1), \
            8print(8N 8lsquare 8seq_head(8M) 8rsquare), \
            8for_each_in_range( \
              8fn(8I, 8do( \
                8print(8lparen), \
                8if(8even(8I), \
                    8print(8quote((-1)*)), \
                    8nil), \
                8print(8N 8lsquare 8seq_at(8dec(8I), 8M) 8rsquare 8quote(*) \
                       8lparen), \
                8determinant_r(8submx(8M, 8S, 8I, 1), \
                               8dec(8S), 8N), \
                8print(8rparen 8rparen), \
                8if(8not(8equal(8I, 8S)), \
                    8print(8quote(+)), \
                    8nil))), \
              1, 8inc(8S)))))

#define ORDER_PP_DEF_8determinant \
  ORDER_PP_FN( \
    8fn(8S, 8N, \
        8determinant_r(8seq_iota(0, 8mul(8S,8S)), 8S, 8N)))

#define DEF_DETERMINANT_N(n) \
  int determinant_##n##_i(int * mx) \
    { return ORDER_PP(8determinant(n, 8quote(mx))); }

ORDER_PP(
  8for_each_in_range(
    8fn(8I,
        8emit(8quote(DEF_DETERMINANT_N),
              8tuple(8I))),
    1,
    6)) // 5 + 1

int main(int argc, char ** argv)
{
  int mx2[] = {1, 0,
               0, 1};
  int mx3[] = {1, 0, 0,
               0, 1, 0,
               0, 0, 1};
  int mx4[] = {1, 0, 0, 0,
               0, 1, 0, 0,
               0, 0, 1, 0,
               0, 0, 0, 1};
  int mx5[] = {1, 0, 0, 0, 0,
               0, 1, 0, 0, 0,
               0, 0, 1, 0, 0,
               0, 0, 0, 1, 0,
               0, 0, 0, 0, 1};
  printf("%i\n", determinant_2_i(mx2));
  printf("%i\n", determinant_3_i(mx3));
  printf("%i\n", determinant_4_i(mx4));
  printf("%i\n", determinant_5_i(mx5));
}