Merge lp:~gang65/ubuntu-calculator-app/reboot-upgrade-math-js-to-1.5 into lp:ubuntu-calculator-app

[Bartosz Kosiorek]

Upgrade math.js to 1.5.2

Upgrade math.js to 1.5.2

[Bartosz Kosiorek (gang65)]

Unfortunately with math.js1.5.0 I noticed regression with complex numbers:
With math.js 1.4.0:
   66i*i = -66

With math.js 1.5.0:
   66i*i = -6.6

Bug already submitted:
https://github.com/josdejong/mathjs/issues/319

review: Needs Fixing

[Ubuntu Phone Apps Jenkins Bot (ubuntu-phone-apps-jenkins-bot)]

FAILED: Continuous integration, rev:152
http://91.189.93.70:8080/job/ubuntu-calculator-app-reboot-ci/139/
Executed test runs:
    UNSTABLE: http://91.189.93.70:8080/job/generic-mediumtests-utopic/2523
        deb: http://91.189.93.70:8080/job/generic-mediumtests-utopic/2523/artifact/work/output/*zip*/output.zip
    SUCCESS: http://91.189.93.70:8080/job/ubuntu-calculator-app-reboot-utopic-amd64-ci/32
    SUCCESS: http://91.189.93.70:8080/job/ubuntu-calculator-app-reboot-vivid-amd64-ci/139

Click here to trigger a rebuild:
http://91.189.93.70:8080/job/ubuntu-calculator-app-reboot-ci/139/rebuild

review: Needs Fixing (continuous-integration)

[lp:~gang65/ubuntu-calculator-app/reboot-upgrade-math-js-to-1.5]

153. By Bartosz Kosiorek on 2015-04-08

Upgrade math.js to 1.5.1 to fix regression

[Bartosz Kosiorek (gang65)]

I upgraded to fixed math.js version 1.5.1.
Now everything should be fine.

review: Approve

[lp:~gang65/ubuntu-calculator-app/reboot-upgrade-math-js-to-1.5]

154. By Bartosz Kosiorek on 2015-04-08

Comment fix

[Bartosz Kosiorek (gang65)]

review: Needs Information

[Ubuntu Phone Apps Jenkins Bot (ubuntu-phone-apps-jenkins-bot)]

PASSED: Continuous integration, rev:154
http://91.189.93.70:8080/job/ubuntu-calculator-app-reboot-ci/141/
Executed test runs:
    SUCCESS: http://91.189.93.70:8080/job/generic-mediumtests-utopic/2536
        deb: http://91.189.93.70:8080/job/generic-mediumtests-utopic/2536/artifact/work/output/*zip*/output.zip
    SUCCESS: http://91.189.93.70:8080/job/ubuntu-calculator-app-reboot-utopic-amd64-ci/34
    SUCCESS: http://91.189.93.70:8080/job/ubuntu-calculator-app-reboot-vivid-amd64-ci/141

Click here to trigger a rebuild:
http://91.189.93.70:8080/job/ubuntu-calculator-app-reboot-ci/141/rebuild

review: Approve (continuous-integration)

[lp:~gang65/ubuntu-calculator-app/reboot-upgrade-math-js-to-1.5]

155. By Bartosz Kosiorek on 2015-04-10

Upgrade math.js to 1.5.2

[Ubuntu Phone Apps Jenkins Bot (ubuntu-phone-apps-jenkins-bot)]

FAILED: Continuous integration, rev:155
http://91.189.93.70:8080/job/ubuntu-calculator-app-reboot-ci/143/
Executed test runs:
    UNSTABLE: http://91.189.93.70:8080/job/generic-mediumtests-utopic/2550
        deb: http://91.189.93.70:8080/job/generic-mediumtests-utopic/2550/artifact/work/output/*zip*/output.zip
    SUCCESS: http://91.189.93.70:8080/job/ubuntu-calculator-app-reboot-utopic-amd64-ci/36
    SUCCESS: http://91.189.93.70:8080/job/ubuntu-calculator-app-reboot-vivid-amd64-ci/143

Click here to trigger a rebuild:
http://91.189.93.70:8080/job/ubuntu-calculator-app-reboot-ci/143/rebuild

review: Needs Fixing (continuous-integration)

[Riccardo Padovani (rpadovani)]

Looks good to me, thanks :-)

review: Approve

[Ubuntu Phone Apps Jenkins Bot (ubuntu-phone-apps-jenkins-bot)]

review: Approve (continuous-integration)

=== modified file 'app/engine/math.js'
--- app/engine/math.js	2015-03-02 09:15:57 +0000
+++ app/engine/math.js	2015-04-10 14:51:16 +0000
@@ -14,11 +14,11 @@
  * It features real and complex numbers, units, matrices, a large set of
  * mathematical functions, and a flexible expression parser.
  *
- * @version 1.2.0
- * @date    2014-12-25
+ * @version 1.5.2
+ * @date    2015-04-09
  *
  * @license
- * Copyright (C) 2013-2014 Jos de Jong <[e-mail]>
+ * Copyright (C) 2013-2015 Jos de Jong <[e-mail]>
  *
  * Licensed under the Apache License, Version 2.0 (the "License"); you may not
  * use this file except in compliance with the License. You may obtain a copy
@@ -35,46 +35,45 @@
 
 (function webpackUniversalModuleDefinition(root, factory) {
     // UCA: we delete all exports, we don't need them, and we keep only our var
-    mathJs = factory();
-})(this, function() {
+    mathJs = factory();})(this, function() {
 return /******/ (function(modules) { // webpackBootstrap
 /******/ 	// The module cache
 /******/ 	var installedModules = {};
-/******/
+
 /******/ 	// The require function
 /******/ 	function __webpack_require__(moduleId) {
-/******/
+
 /******/ 		// Check if module is in cache
 /******/ 		if(installedModules[moduleId])
 /******/ 			return installedModules[moduleId].exports;
-/******/
+
 /******/ 		// Create a new module (and put it into the cache)
 /******/ 		var module = installedModules[moduleId] = {
 /******/ 			exports: {},
 /******/ 			id: moduleId,
 /******/ 			loaded: false
 /******/ 		};
-/******/
+
 /******/ 		// Execute the module function
 /******/ 		modules[moduleId].call(module.exports, module, module.exports, __webpack_require__);
-/******/
+
 /******/ 		// Flag the module as loaded
 /******/ 		module.loaded = true;
-/******/
+
 /******/ 		// Return the exports of the module
 /******/ 		return module.exports;
 /******/ 	}
-/******/
-/******/
+
+
 /******/ 	// expose the modules object (__webpack_modules__)
 /******/ 	__webpack_require__.m = modules;
-/******/
+
 /******/ 	// expose the module cache
 /******/ 	__webpack_require__.c = installedModules;
-/******/
+
 /******/ 	// __webpack_public_path__
 /******/ 	__webpack_require__.p = "";
-/******/
+
 /******/ 	// Load entry module and return exports
 /******/ 	return __webpack_require__(0);
 /******/ })
@@ -162,27 +161,27 @@
 
 	      // TODO: remove deprecated setting some day (deprecated since version 0.17.0)
 	      if (options.number && options.number.defaultType) {
-	        throw new Error('setting `number.defaultType` is deprecated. Use `number` instead.')
+	        throw new Error('setting `number.defaultType` is deprecated. Use `number` instead.');
 	      }
 
 	      // TODO: remove deprecated setting some day (deprecated since version 0.17.0)
 	      if (options.number && options.number.precision) {
-	        throw new Error('setting `number.precision` is deprecated. Use `precision` instead.')
+	        throw new Error('setting `number.precision` is deprecated. Use `precision` instead.');
 	      }
 
 	      // TODO: remove deprecated setting some day (deprecated since version 0.17.0)
 	      if (options.matrix && options.matrix.defaultType) {
-	        throw new Error('setting `matrix.defaultType` is deprecated. Use `matrix` instead.')
+	        throw new Error('setting `matrix.defaultType` is deprecated. Use `matrix` instead.');
 	      }
 
 	      // TODO: remove deprecated setting some day (deprecated since version 0.15.0)
 	      if (options.matrix && options.matrix['default']) {
-	        throw new Error('setting `matrix.default` is deprecated. Use `matrix` instead.')
+	        throw new Error('setting `matrix.default` is deprecated. Use `matrix` instead.');
 	      }
 
 	      // TODO: remove deprecated setting some day (deprecated since version 0.20.0)
 	      if (options.decimals) {
-	        throw new Error('setting `decimals` is deprecated. Use `precision` instead.')
+	        throw new Error('setting `decimals` is deprecated. Use `precision` instead.');
 	      }
 	    }
 
@@ -207,6 +206,29 @@
 	  // create a new BigNumber factory for this instance of math.js
 	  var BigNumber = __webpack_require__(5).constructor();
 
+	  /**
+	   * Get a JSON representation of a BigNumber containing
+	   * type information
+	   * @returns {Object} Returns a JSON object structured as:
+	   *                   `{"mathjs": "BigNumber", "value": "0.2"}`
+	   */
+	  BigNumber.prototype.toJSON = function () {
+	    return {
+	      mathjs: 'BigNumber',
+	      value: this.toString()
+	    };
+	  };
+
+	  /**
+	   * Instantiate a BigNumber from a JSON object
+	   * @param {Object} json  a JSON object structured as:
+	   *                       `{"mathjs": "BigNumber", "value": "0.2"}`
+	   * @return {BigNumber}
+	   */
+	  BigNumber.fromJSON = function (json) {
+	    return new BigNumber(json.value);
+	  };
+
 	  // extend BigNumber with a function clone
 	  if (typeof BigNumber.prototype.clone !== 'function') {
 	    /**
@@ -248,33 +270,39 @@
 	  math.type.Complex = __webpack_require__(7);
 	  math.type.Range = __webpack_require__(8);
 	  math.type.Index = __webpack_require__(9);
-	  math.type.Matrix = __webpack_require__(10);
+	  math.type.Matrix = __webpack_require__(10)(_config);
 	  math.type.Unit = __webpack_require__(11);
 	  math.type.Help = __webpack_require__(12);
 	  math.type.ResultSet = __webpack_require__(13);
 	  math.type.BigNumber = BigNumber;
 
-	  math.collection = __webpack_require__(14);
+	  math.collection = __webpack_require__(14)(math, _config);
+
+	  // matrix storage formats
+	  math.type.CcsMatrix = __webpack_require__(15)(math, _config);
+	  math.type.CrsMatrix = __webpack_require__(16)(math, _config);
+	  math.type.DenseMatrix = __webpack_require__(17)(math, _config);
+
+	  // matrix storage format registry
+	  math.type.Matrix._storage.ccs = math.type.CcsMatrix;
+	  math.type.Matrix._storage.crs = math.type.CrsMatrix;
+	  math.type.Matrix._storage.dense = math.type.DenseMatrix;
+	  math.type.Matrix._storage['default'] = math.type.DenseMatrix;
 
 	  // expression (parse, Parser, nodes, docs)
 	  math.expression = {};
-	  math.expression.node = __webpack_require__(15);
-	  math.expression.parse = __webpack_require__(16);
-	  math.expression.Parser = __webpack_require__(17);
-	  math.expression.docs = __webpack_require__(18);
+	  math.expression.node = __webpack_require__(18);
+	  math.expression.parse = __webpack_require__(19)(math, _config);
+	  math.expression.Parser = __webpack_require__(20)(math, _config);
+	  math.expression.docs = __webpack_require__(21);
 
 	  // serialization utilities
 	  math.json = {
-	    reviver: __webpack_require__(19)
+	    reviver: __webpack_require__(22)(math, _config)
 	  };
-
-	  // expression parser
-	  __webpack_require__(31)(math, _config);
-	  __webpack_require__(32)(math, _config);
-	  __webpack_require__(33)(math, _config);
+	  
+	  // functions - construction (must be defined before the rest of functions)
 	  __webpack_require__(34)(math, _config);
-
-	  // functions - arithmetic
 	  __webpack_require__(35)(math, _config);
 	  __webpack_require__(36)(math, _config);
 	  __webpack_require__(37)(math, _config);
@@ -284,10 +312,14 @@
 	  __webpack_require__(41)(math, _config);
 	  __webpack_require__(42)(math, _config);
 	  __webpack_require__(43)(math, _config);
+
+	  // expression parser
 	  __webpack_require__(44)(math, _config);
 	  __webpack_require__(45)(math, _config);
 	  __webpack_require__(46)(math, _config);
 	  __webpack_require__(47)(math, _config);
+
+	  // functions - arithmetic
 	  __webpack_require__(48)(math, _config);
 	  __webpack_require__(49)(math, _config);
 	  __webpack_require__(50)(math, _config);
@@ -304,8 +336,6 @@
 	  __webpack_require__(61)(math, _config);
 	  __webpack_require__(62)(math, _config);
 	  __webpack_require__(63)(math, _config);
-
-	  // functions - bitwise
 	  __webpack_require__(64)(math, _config);
 	  __webpack_require__(65)(math, _config);
 	  __webpack_require__(66)(math, _config);
@@ -313,16 +343,14 @@
 	  __webpack_require__(68)(math, _config);
 	  __webpack_require__(69)(math, _config);
 	  __webpack_require__(70)(math, _config);
-
-	  // functions - complex
 	  __webpack_require__(71)(math, _config);
 	  __webpack_require__(72)(math, _config);
 	  __webpack_require__(73)(math, _config);
 	  __webpack_require__(74)(math, _config);
-
-	  // functions - construction
 	  __webpack_require__(75)(math, _config);
 	  __webpack_require__(76)(math, _config);
+
+	  // functions - bitwise
 	  __webpack_require__(77)(math, _config);
 	  __webpack_require__(78)(math, _config);
 	  __webpack_require__(79)(math, _config);
@@ -330,18 +358,20 @@
 	  __webpack_require__(81)(math, _config);
 	  __webpack_require__(82)(math, _config);
 	  __webpack_require__(83)(math, _config);
+
+	  // functions - complex
 	  __webpack_require__(84)(math, _config);
-
-	  // functions - logical
 	  __webpack_require__(85)(math, _config);
 	  __webpack_require__(86)(math, _config);
 	  __webpack_require__(87)(math, _config);
+
+	  // functions - logical
 	  __webpack_require__(88)(math, _config);
-
-	  // functions - matrix
 	  __webpack_require__(89)(math, _config);
 	  __webpack_require__(90)(math, _config);
 	  __webpack_require__(91)(math, _config);
+
+	  // functions - matrix
 	  __webpack_require__(92)(math, _config);
 	  __webpack_require__(93)(math, _config);
 	  __webpack_require__(94)(math, _config);
@@ -356,41 +386,41 @@
 	  __webpack_require__(103)(math, _config);
 	  __webpack_require__(104)(math, _config);
 	  __webpack_require__(105)(math, _config);
-
-	  // functions - probability
-	  //require('./function/probability/distribution')(math, _config); // TODO: rethink math.distribution
 	  __webpack_require__(106)(math, _config);
 	  __webpack_require__(107)(math, _config);
 	  __webpack_require__(108)(math, _config);
+
+	  // functions - probability
+	  //require('./function/probability/distribution')(math, _config); // TODO: rethink math.distribution
 	  __webpack_require__(109)(math, _config);
 	  __webpack_require__(110)(math, _config);
 	  __webpack_require__(111)(math, _config);
 	  __webpack_require__(112)(math, _config);
-
-	  // functions - relational
 	  __webpack_require__(113)(math, _config);
 	  __webpack_require__(114)(math, _config);
 	  __webpack_require__(115)(math, _config);
+
+	  // functions - relational
 	  __webpack_require__(116)(math, _config);
 	  __webpack_require__(117)(math, _config);
 	  __webpack_require__(118)(math, _config);
 	  __webpack_require__(119)(math, _config);
 	  __webpack_require__(120)(math, _config);
-
-	  // functions - statistics
 	  __webpack_require__(121)(math, _config);
 	  __webpack_require__(122)(math, _config);
 	  __webpack_require__(123)(math, _config);
+
+	  // functions - statistics
 	  __webpack_require__(124)(math, _config);
 	  __webpack_require__(125)(math, _config);
 	  __webpack_require__(126)(math, _config);
 	  __webpack_require__(127)(math, _config);
 	  __webpack_require__(128)(math, _config);
-
-	  // functions - trigonometry
 	  __webpack_require__(129)(math, _config);
 	  __webpack_require__(130)(math, _config);
 	  __webpack_require__(131)(math, _config);
+
+	  // functions - trigonometry
 	  __webpack_require__(132)(math, _config);
 	  __webpack_require__(133)(math, _config);
 	  __webpack_require__(134)(math, _config);
@@ -413,20 +443,23 @@
 	  __webpack_require__(151)(math, _config);
 	  __webpack_require__(152)(math, _config);
 	  __webpack_require__(153)(math, _config);
-
-	  // functions - units
 	  __webpack_require__(154)(math, _config);
-
-	  // functions - utils
 	  __webpack_require__(155)(math, _config);
 	  __webpack_require__(156)(math, _config);
+
+	  // functions - units
 	  __webpack_require__(157)(math, _config);
+
+	  // functions - utils
 	  __webpack_require__(158)(math, _config);
 	  __webpack_require__(159)(math, _config);
 	  __webpack_require__(160)(math, _config);
 	  __webpack_require__(161)(math, _config);
 	  __webpack_require__(162)(math, _config);
 	  __webpack_require__(163)(math, _config);
+	  __webpack_require__(164)(math, _config);
+	  __webpack_require__(165)(math, _config);
+	  __webpack_require__(166)(math, _config);
 
 	  // TODO: deprecated since version 0.25.0, remove some day.
 	  math.ifElse = function () {
@@ -438,21 +471,21 @@
 
 	  // attach transform functions (for converting one-based indices to zero-based)
 	  math.expression.transform = {
-	    concat: __webpack_require__(20)(math, _config),
-	    filter: __webpack_require__(21)(math, _config),
-	    forEach:__webpack_require__(22)(math, _config),
-	    index:  __webpack_require__(23)(math, _config),
-	    map:    __webpack_require__(24)(math, _config),
-	    max:    __webpack_require__(25)(math, _config),
-	    mean:   __webpack_require__(26)(math, _config),
-	    min:    __webpack_require__(27)(math, _config),
-	    range:  __webpack_require__(28)(math, _config),
-	    subset: __webpack_require__(29)(math, _config)
+	    concat: __webpack_require__(23)(math, _config),
+	    filter: __webpack_require__(24)(math, _config),
+	    forEach:__webpack_require__(25)(math, _config),
+	    index:  __webpack_require__(26)(math, _config),
+	    map:    __webpack_require__(27)(math, _config),
+	    max:    __webpack_require__(28)(math, _config),
+	    mean:   __webpack_require__(29)(math, _config),
+	    min:    __webpack_require__(30)(math, _config),
+	    range:  __webpack_require__(31)(math, _config),
+	    subset: __webpack_require__(32)(math, _config)
 	  };
 
 	  // selector (we initialize after all functions are loaded)
 	  math.chaining = {};
-	  math.chaining.Chain = __webpack_require__(30)(math, _config);
+	  math.chaining.Chain = __webpack_require__(33)(math, _config);
 	  math.chaining.Selector = math.chaining.Chain; // TODO: deprecate in v2.0
 
 	  // apply provided configuration options
@@ -684,6 +717,8 @@
 
 	'use strict';
 
+	var NumberFormatter = __webpack_require__(167);
+
 	/**
 	 * Test whether value is a Number
 	 * @param {*} value
@@ -829,44 +864,15 @@
 	      return exports.toExponential(value, precision);
 
 	    case 'auto':
-	      // determine lower and upper bound for exponential notation.
-	        // TODO: implement support for upper and lower to be BigNumbers themselves
-	      var lower = 1e-3;
-	      var upper = 1e5;
-	      if (options && options.exponential) {
-	        if (options.exponential.lower !== undefined) {
-	          lower = options.exponential.lower;
-	        }
-	        if (options.exponential.upper !== undefined) {
-	          upper = options.exponential.upper;
-	        }
-	      }
-
-	      // handle special case zero
-	      if (value === 0) return '0';
-
-	      // determine whether or not to output exponential notation
-	      var str;
-	      var abs = Math.abs(value);
-	      if (abs >= lower && abs < upper) {
-	        // normal number notation
-	        // Note: IE7 does not allow value.toPrecision(undefined)
-	        var valueStr = precision ?
-	            value.toPrecision(Math.min(precision, 21)) :
-	            value.toPrecision();
-	        str = parseFloat(valueStr) + '';
-	      }
-	      else {
-	        // exponential notation
-	        str = exports.toExponential(value, precision);
-	      }
-
-	      // remove trailing zeros after the decimal point
-	      return str.replace(/((\.\d*?)(0+))($|e)/, function () {
-	        var digits = arguments[2];
-	        var e = arguments[4];
-	        return (digits !== '.') ? digits + e : e;
-	      });
+	      return exports
+	          .toPrecision(value, precision, options && options.exponential)
+
+	          // remove trailing zeros after the decimal point
+	          .replace(/((\.\d*?)(0+))($|e)/, function () {
+	            var digits = arguments[2];
+	            var e = arguments[4];
+	            return (digits !== '.') ? digits + e : e;
+	          });
 
 	    default:
 	      throw new Error('Unknown notation "' + notation + '". ' +
@@ -883,12 +889,7 @@
 	 * @returns {string} str
 	 */
 	exports.toExponential = function(value, precision) {
-	  if (precision !== undefined) {
-	    return value.toExponential(Math.min(precision - 1, 20));
-	  }
-	  else {
-	    return value.toExponential();
-	  }
+	  return new NumberFormatter(value).toExponential(precision);
 	};
 
 	/**
@@ -898,7 +899,20 @@
 	 *                                      decimal point. Zero by default.
 	 */
 	exports.toFixed = function(value, precision) {
-	  return value.toFixed(Math.min(precision, 20));
+	  return new NumberFormatter(value).toFixed(precision);
+	};
+
+	/**
+	 * Format a number with a certain precision
+	 * @param {Number} value
+	 * @param {Number} [precision=undefined] Optional number of digits.
+	 * @param {{lower: number, upper: number}} [options]  By default:
+	 *                                                    lower = 1e-3 (excl)
+	 *                                                    upper = 1e+5 (incl)
+	 * @return {string}
+	 */
+	exports.toPrecision = function(value, precision, options) {
+	  return new NumberFormatter(value).toPrecision(precision, options);
 	};
 
 	/**
@@ -970,14 +984,14 @@
 
 	module.exports = function (math, config) {
 	  var object = __webpack_require__(2);
-	  var bignumber = __webpack_require__(164);
+	  var bignumber = __webpack_require__(168);
 	  var Complex = __webpack_require__(7);
 	  var BigNumber = math.type.BigNumber;
 
 	  math['true']     = true;
 	  math['false']    = false;
 	  math['null']     = null;
-	  math['uninitialized'] = __webpack_require__(165).UNINITIALIZED;
+	  math['uninitialized'] = __webpack_require__(169).UNINITIALIZED;
 
 	  if (config.number === 'bignumber') {
 	    math['Infinity'] = new BigNumber(Infinity);
@@ -1022,7 +1036,7 @@
 	  math.i = new Complex(0, 1);
 
 	  // meta information
-	  math.version = __webpack_require__(166);
+	  math.version = __webpack_require__(170);
 	};
 
 
@@ -1030,33 +1044,10 @@
 /* 5 */
 /***/ function(module, exports, __webpack_require__) {
 
-	var BigNumber = __webpack_require__(336);
+	var BigNumber = __webpack_require__(340);
 
 	// FIXME: replace all require('decimal.js') with require('./BigNumber').
 
-	/**
-	 * Get a JSON representation of a BigNumber containing
-	 * type information
-	 * @returns {Object} Returns a JSON object structured as:
-	 *                   `{"mathjs": "BigNumber", "value": "0.2"}`
-	 */
-	BigNumber.prototype.toJSON = function () {
-	  return {
-	    mathjs: 'BigNumber',
-	    value: this.toString()
-	  };
-	};
-
-	/**
-	 * Instantiate a BigNumber from a JSON object
-	 * @param {Object} json  a JSON object structured as:
-	 *                       `{"mathjs": "BigNumber", "value": "0.2"}`
-	 * @return {BigNumber}
-	 */
-	BigNumber.fromJSON = function (json) {
-	  return new BigNumber(json.value);
-	};
-
 	module.exports = BigNumber;
 
 
@@ -1066,10 +1057,10 @@
 
 	'use strict';
 
-	exports.ArgumentsError = __webpack_require__(167);
-	exports.DimensionError = __webpack_require__(168);
-	exports.IndexError = __webpack_require__(169);
-	exports.UnsupportedTypeError = __webpack_require__(170);
+	exports.ArgumentsError = __webpack_require__(171);
+	exports.DimensionError = __webpack_require__(172);
+	exports.IndexError = __webpack_require__(173);
+	exports.UnsupportedTypeError = __webpack_require__(174);
 
 	// TODO: implement an InvalidValueError?
 
@@ -1080,7 +1071,7 @@
 
 	'use strict';
 
-	var util = __webpack_require__(171),
+	var util = __webpack_require__(175),
 	    Unit = __webpack_require__(11),
 	    number = util.number,
 
@@ -1559,7 +1550,7 @@
 
 	'use strict';
 
-	var util = __webpack_require__(171);
+	var util = __webpack_require__(175);
 
 	var number = util.number;
 	var string = util.string;
@@ -1858,7 +1849,7 @@
 
 	'use strict';
 
-	var util = __webpack_require__(171),
+	var util = __webpack_require__(175),
 
 	    Range = __webpack_require__(8),
 
@@ -2155,542 +2146,268 @@
 
 	'use strict';
 
-	var util = __webpack_require__(171);
-	var DimensionError = __webpack_require__(168);
-
-	var Index = __webpack_require__(9);
-
-	var number = util.number;
-	var string = util.string;
-	var array = util.array;
-	var object = util.object;
-
-	var isArray = Array.isArray;
-	var validateIndex = array.validateIndex;
-
-	/**
-	 * @constructor Matrix
-	 *
-	 * A Matrix is a wrapper around an Array. A matrix can hold a multi dimensional
-	 * array. A matrix can be constructed as:
-	 *     var matrix = new Matrix(data)
-	 *
-	 * Matrix contains the functions to resize, get and set values, get the size,
-	 * clone the matrix and to convert the matrix to a vector, array, or scalar.
-	 * Furthermore, one can iterate over the matrix using map and forEach.
-	 * The internal Array of the Matrix can be accessed using the function valueOf.
-	 *
-	 * Example usage:
-	 *     var matrix = new Matrix([[1, 2], [3, 4]);
-	 *     matix.size();              // [2, 2]
-	 *     matrix.resize([3, 2], 5);
-	 *     matrix.valueOf();          // [[1, 2], [3, 4], [5, 5]]
-	 *     matrix.subset([1,2])       // 3 (indexes are zero-based)
-	 *
-	 * @param {Array | Matrix} [data]    A multi dimensional array
-	 */
-	function Matrix(data) {
-	  if (!(this instanceof Matrix)) {
-	    throw new SyntaxError('Constructor must be called with the new operator');
-	  }
-
-	  if (data instanceof Matrix) {
-	    // clone data from a Matrix
-	    this._data = data.clone()._data;
-	  }
-	  else if (isArray(data)) {
-	    // use array
-	    // replace nested Matrices with Arrays
-	    this._data = preprocess(data);
-	  }
-	  else if (data != null) {
-	    // unsupported type
-	    throw new TypeError('Unsupported type of data (' + util.types.type(data) + ')');
-	  }
-	  else {
-	    // nothing provided
-	    this._data = [];
-	  }
-
-	  // verify the size of the array
-	  this._size = array.size(this._data);
-	}
-
-	/**
-	 * Test whether an object is a Matrix
-	 * @param {*} object
-	 * @return {Boolean} isMatrix
-	 */
-	Matrix.isMatrix = function (object) {
-	  return (object instanceof Matrix);
-	};
-
-	/**
-	 * Get a subset of the matrix, or replace a subset of the matrix.
-	 *
-	 * Usage:
-	 *     var subset = matrix.subset(index)               // retrieve subset
-	 *     var value = matrix.subset(index, replacement)   // replace subset
-	 *
-	 * @param {Index} index
-	 * @param {Array | Matrix | *} [replacement]
-	 * @param {*} [defaultValue=0]      Default value, filled in on new entries when
-	 *                                  the matrix is resized. If not provided,
-	 *                                  new matrix elements will be filled with zeros.
-	 */
-	Matrix.prototype.subset = function (index, replacement, defaultValue) {
-	  switch (arguments.length) {
-	    case 1:
-	      return _get(this, index);
-
-	    // intentional fall through
-	    case 2:
-	    case 3:
-	      return _set(this, index, replacement, defaultValue);
-
-	    default:
-	      throw new SyntaxError('Wrong number of arguments');
-	  }
-	};
-
-	/**
-	 * Get a single element from the matrix.
-	 * @param {Number[]} index   Zero-based index
-	 * @return {*} value
-	 */
-	Matrix.prototype.get = function (index) {
-	  if (!isArray(index)) {
-	    throw new TypeError('Array expected');
-	  }
-	  if (index.length != this._size.length) {
-	    throw new DimensionError(index.length, this._size.length);
-	  }
-
-	  var data = this._data;
-	  for (var i = 0, ii = index.length; i < ii; i++) {
-	    var index_i = index[i];
-	    validateIndex(index_i, data.length);
-	    data = data[index_i];
-	  }
-
-	  return object.clone(data);
-	};
-
-	/**
-	 * Replace a single element in the matrix.
-	 * @param {Number[]} index   Zero-based index
-	 * @param {*} value
-	 * @param {*} [defaultValue]        Default value, filled in on new entries when
-	 *                                  the matrix is resized. If not provided,
-	 *                                  new matrix elements will be left undefined.
-	 * @return {Matrix} self
-	 */
-	Matrix.prototype.set = function (index, value, defaultValue) {
-	  var i, ii;
-
-	  // validate input type and dimensions
-	  if (!isArray(index)) {
-	    throw new Error('Array expected');
-	  }
-	  if (index.length < this._size.length) {
-	    throw new DimensionError(index.length, this._size.length, '<');
-	  }
-
-	  // enlarge matrix when needed
-	  var size = index.map(function (i) {
-	    return i + 1;
-	  });
-	  _fit(this, size, defaultValue);
-
-	  // traverse over the dimensions
-	  var data = this._data;
-	  for (i = 0, ii = index.length - 1; i < ii; i++) {
-	    var index_i = index[i];
-	    validateIndex(index_i, data.length);
-	    data = data[index_i];
-	  }
-
-	  // set new value
-	  index_i = index[index.length - 1];
-	  validateIndex(index_i, data.length);
-	  data[index_i] = value;
-
-	  return this;
-	};
-
-	/**
-	 * Get a submatrix of this matrix
-	 * @param {Matrix} matrix
-	 * @param {Index} index   Zero-based index
-	 * @private
-	 */
-	function _get (matrix, index) {
-	  if (!(index instanceof Index)) {
-	    throw new TypeError('Invalid index');
-	  }
-
-	  var isScalar = index.isScalar();
-	  if (isScalar) {
-	    // return a scalar
-	    return matrix.get(index.min());
-	  }
-	  else {
-	    // validate dimensions
-	    var size = index.size();
-	    if (size.length != matrix._size.length) {
-	      throw new DimensionError(size.length, matrix._size.length);
-	    }
-
-	    // validate if any of the ranges in the index is out of range
-	    var min = index.min();
-	    var max = index.max();
-	    for (var i = 0, ii = matrix._size.length; i < ii; i++) {
-	      validateIndex(min[i], matrix._size[i]);
-	      validateIndex(max[i], matrix._size[i]);
-	    }
-
-	    // retrieve submatrix
-	    // TODO: more efficient when creating an empty matrix and setting _data and _size manually
-	    return new Matrix(_getSubmatrix(matrix._data, index, size.length, 0));
-	  }
-	}
-
-	/**
-	 * Recursively get a submatrix of a multi dimensional matrix.
-	 * Index is not checked for correct number or length of dimensions.
-	 * @param {Array} data
-	 * @param {Index} index
-	 * @param {number} dims   Total number of dimensions
-	 * @param {number} dim    Current dimension
-	 * @return {Array} submatrix
-	 * @private
-	 */
-	function _getSubmatrix (data, index, dims, dim) {
-	  var last = (dim == dims - 1);
-	  var range = index.range(dim);
-
-	  if (last) {
-	    return range.map(function (i) {
-	      return data[i];
-	    });
-	  }
-	  else {
-	    return range.map(function (i) {
-	      var child = data[i];
-	      return _getSubmatrix(child, index, dims, dim + 1);
-	    });
-	  }
-	}
-
-	/**
-	 * Replace a submatrix in this matrix
-	 * Indexes are zero-based.
-	 * @param {Matrix} matrix
-	 * @param {Index} index
-	 * @param {Matrix | Array | *} submatrix
-	 * @param {*} defaultValue          Default value, filled in on new entries when
-	 *                                  the matrix is resized.
-	 * @return {Matrix} matrix
-	 * @private
-	 */
-	function _set (matrix, index, submatrix, defaultValue) {
-	  if (!(index instanceof Index)) {
-	    throw new TypeError('Invalid index');
-	  }
-
-	  // get index size and check whether the index contains a single value
-	  var iSize = index.size(),
-	      isScalar = index.isScalar();
-
-	  // calculate the size of the submatrix, and convert it into an Array if needed
-	  var sSize;
-	  if (submatrix instanceof Matrix) {
-	    sSize = submatrix.size();
-	    submatrix = submatrix.valueOf();
-	  }
-	  else {
-	    sSize = array.size(submatrix);
-	  }
-
-	  if (isScalar) {
-	    // set a scalar
-
-	    // check whether submatrix is a scalar
-	    if (sSize.length != 0) {
-	      throw new TypeError('Scalar expected');
-	    }
-
-	    matrix.set(index.min(), submatrix, defaultValue);
-	  }
-	  else {
-	    // set a submatrix
-
-	    // validate dimensions
-	    if (iSize.length < matrix._size.length) {
-	      throw new DimensionError(iSize.length, matrix._size.length, '<');
-	    }
-
-	    if (sSize.length < iSize.length) {
-	      // calculate number of missing outer dimensions
-	      var i = 0;
-	      var outer = 0;
-	      while (iSize[i] === 1 && sSize[i] === 1) {
-	        i++;
-	      }
-	      while (iSize[i] === 1) {
-	        outer++;
-	        i++;
-	      }
-
-	      // unsqueeze both outer and inner dimensions
-	      submatrix = array.unsqueeze(submatrix, iSize.length, outer, sSize);
-	    }
-
-	    // check whether the size of the submatrix matches the index size
-	    if (!object.deepEqual(iSize, sSize)) {
-	      throw new DimensionError(iSize, sSize, '>');
-	    }
-
-	    // enlarge matrix when needed
-	    var size = index.max().map(function (i) {
-	      return i + 1;
-	    });
-	    _fit(matrix, size, defaultValue);
-
-	    // insert the sub matrix
-	    var dims = iSize.length,
-	        dim = 0;
-	    _setSubmatrix (matrix._data, index, submatrix, dims, dim);
-	  }
-
-	  return matrix;
-	}
-
-	/**
-	 * Replace a submatrix of a multi dimensional matrix.
-	 * @param {Array} data
-	 * @param {Index} index
-	 * @param {Array} submatrix
-	 * @param {number} dims   Total number of dimensions
-	 * @param {number} dim
-	 * @private
-	 */
-	function _setSubmatrix (data, index, submatrix, dims, dim) {
-	  var last = (dim == dims - 1),
-	      range = index.range(dim);
-
-	  if (last) {
-	    range.forEach(function (dataIndex, subIndex) {
-	      validateIndex(dataIndex);
-	      data[dataIndex] = submatrix[subIndex];
-	    });
-	  }
-	  else {
-	    range.forEach(function (dataIndex, subIndex) {
-	      validateIndex(dataIndex);
-	      _setSubmatrix(data[dataIndex], index, submatrix[subIndex], dims, dim + 1);
-	    });
-	  }
-	}
-
-	/**
-	 * Resize the matrix
-	 * @param {Number[]} size
-	 * @param {*} [defaultValue=0]      Default value, filled in on new entries.
-	 *                                  If not provided, the matrix elements will
-	 *                                  be filled with zeros.
-	 * @return {Matrix} self            The matrix itself is returned
-	 */
-	Matrix.prototype.resize = function (size, defaultValue) {
-	  this._size = object.clone(size);
-	  this._data = array.resize(this._data, this._size, defaultValue);
-
-	  // return the matrix itself
-	  return this;
-	};
-
-	/**
-	 * Enlarge the matrix when it is smaller than given size.
-	 * If the matrix is larger or equal sized, nothing is done.
-	 * @param {Matrix} matrix           The matrix to be resized
-	 * @param {Number[]} size
-	 * @param {*} defaultValue          Default value, filled in on new entries.
-	 * @private
-	 */
-	function _fit(matrix, size, defaultValue) {
-	  var newSize = object.clone(matrix._size),
-	      changed = false;
-
-	  // add dimensions when needed
-	  while (newSize.length < size.length) {
-	    newSize.push(0);
-	    changed = true;
-	  }
-
-	  // enlarge size when needed
-	  for (var i = 0, ii = size.length; i < ii; i++) {
-	    if (size[i] > newSize[i]) {
-	      newSize[i] = size[i];
-	      changed = true;
-	    }
-	  }
-
-	  if (changed) {
-	    // resize only when size is changed
-	    matrix.resize(newSize, defaultValue);
-	  }
-	}
-
-	/**
-	 * Create a clone of the matrix
-	 * @return {Matrix} clone
-	 */
-	Matrix.prototype.clone = function () {
-	  var matrix = new Matrix();
-	  matrix._data = object.clone(this._data);
-	  matrix._size = object.clone(this._size);
-	  return matrix;
-	};
-
-	/**
-	 * Retrieve the size of the matrix.
-	 * @returns {Number[]} size
-	 */
-	Matrix.prototype.size = function() {
-	  return this._size;
-	};
-
-	/**
-	 * Create a new matrix with the results of the callback function executed on
-	 * each entry of the matrix.
-	 * @param {function} callback   The callback function is invoked with three
-	 *                              parameters: the value of the element, the index
-	 *                              of the element, and the Matrix being traversed.
-	 * @return {Matrix} matrix
-	 */
-	Matrix.prototype.map = function (callback) {
-	  var me = this;
-	  var matrix = new Matrix();
-
-	  var recurse = function (value, index) {
-	    if (isArray(value)) {
-	      return value.map(function (child, i) {
-	        return recurse(child, index.concat(i));
-	      });
-	    }
-	    else {
-	      return callback(value, index, me);
-	    }
-	  };
-	  matrix._data = recurse(this._data, []);
-	  matrix._size = object.clone(this._size);
-
-	  return matrix;
-	};
-
-	/**
-	 * Execute a callback function on each entry of the matrix.
-	 * @param {function} callback   The callback function is invoked with three
-	 *                              parameters: the value of the element, the index
-	 *                              of the element, and the Matrix being traversed.
-	 */
-	Matrix.prototype.forEach = function (callback) {
-	  var me = this;
-
-	  var recurse = function (value, index) {
-	    if (isArray(value)) {
-	      value.forEach(function (child, i) {
-	        recurse(child, index.concat(i));
-	      });
-	    }
-	    else {
-	      callback(value, index, me);
-	    }
-	  };
-	  recurse(this._data, []);
-	};
-
-	/**
-	 * Create an Array with a copy of the data of the Matrix
-	 * @returns {Array} array
-	 */
-	Matrix.prototype.toArray = function () {
-	  return object.clone(this._data);
-	};
-
-	/**
-	 * Get the primitive value of the Matrix: a multidimensional array
-	 * @returns {Array} array
-	 */
-	Matrix.prototype.valueOf = function () {
-	  return this._data;
-	};
-
-	/**
-	 * Get a string representation of the matrix, with optional formatting options.
-	 * @param {Object | Number | Function} [options]  Formatting options. See
-	 *                                                lib/util/number:format for a
-	 *                                                description of the available
-	 *                                                options.
-	 * @returns {String} str
-	 */
-	Matrix.prototype.format = function (options) {
-	  return string.format(this._data, options);
-	};
-
-	/**
-	 * Get a string representation of the matrix
-	 * @returns {String} str
-	 */
-	Matrix.prototype.toString = function () {
-	  return string.format(this._data);
-	};
-
-	/**
-	 * Get a JSON representation of the matrix
-	 * @returns {Object}
-	 */
-	Matrix.prototype.toJSON = function () {
-	  return {
-	    mathjs: 'Matrix',
-	    data: this._data
-	  }
-	};
-
-	/**
-	 * Generate a matrix from a JSON object
-	 * @param {Object} json  An object structured like
-	 *                       `{"mathjs": "Matrix", data: []}`,
-	 *                       where mathjs is optional
-	 * @returns {Matrix}
-	 */
-	Matrix.fromJSON = function (json) {
-	  return new Matrix(json.data);
-	};
-
-	/**
-	 * Preprocess data, which can be an Array or Matrix with nested Arrays and
-	 * Matrices. Replaces all nested Matrices with Arrays
-	 * @param {Array} data
-	 * @return {Array} data
-	 */
-	function preprocess(data) {
-	  for (var i = 0, ii = data.length; i < ii; i++) {
-	    var elem = data[i];
-	    if (isArray(elem)) {
-	      data[i] = preprocess(elem);
-	    }
-	    else if (elem instanceof Matrix) {
-	      data[i] = preprocess(elem._data);
-	    }
-	  }
-
-	  return data;
-	}
-
-	// exports
-	module.exports = Matrix;
-
+	var string = __webpack_require__(176),
+
+	    isString = string.isString;
+
+	module.exports = function (config) {
+
+	  /**
+	   * @constructor Matrix
+	   *
+	   * A Matrix is a wrapper around an Array. A matrix can hold a multi dimensional
+	   * array. A matrix can be constructed as:
+	   *     var matrix = math.matrix(data)
+	   *
+	   * Matrix contains the functions to resize, get and set values, get the size,
+	   * clone the matrix and to convert the matrix to a vector, array, or scalar.
+	   * Furthermore, one can iterate over the matrix using map and forEach.
+	   * The internal Array of the Matrix can be accessed using the function valueOf.
+	   *
+	   * Example usage:
+	   *     var matrix = math.matrix([[1, 2], [3, 4]]);
+	   *     matix.size();              // [2, 2]
+	   *     matrix.resize([3, 2], 5);
+	   *     matrix.valueOf();          // [[1, 2], [3, 4], [5, 5]]
+	   *     matrix.subset([1,2])       // 3 (indexes are zero-based)
+	   *
+	   */
+	  function Matrix() {
+	    if (!(this instanceof Matrix)) {
+	      throw new SyntaxError('Constructor must be called with the new operator');
+	    }
+	  }
+
+	  /**
+	   * Test whether an object is a Matrix
+	   * @param {*} object
+	   * @return {Boolean} isMatrix
+	   */
+	  Matrix.isMatrix = function (object) {
+	    return (object instanceof Matrix);
+	  };
+
+	  /**
+	   * Get the Matrix storage constructor for the given format.
+	   *
+	   * @param {string} format       The Matrix storage format.
+	   *
+	   * @return {Function}           The Matrix storage constructor.
+	   */
+	  Matrix.storage = function (format) {
+	    // check storage format is a string
+	    if (!isString(format)) {
+	      throw new TypeError('format must be a string value');
+	    }
+
+	    // get storage format constructor
+	    var constructor = Matrix._storage[format];
+	    if (!constructor) {
+	      throw new SyntaxError('Unsupported matrix storage format: ' + format);
+	    }
+
+	    // return storage constructor
+	    return constructor;
+	  };
+
+	  // a map with all constructors for all storage types
+	  Matrix._storage = {};
+
+	  /**
+	   * Get the storage format used by the matrix.
+	   *
+	   * Usage:
+	   *     var format = matrix.storage()                   // retrieve storage format
+	   *
+	   * @return {string}           The storage format.
+	   */
+	  Matrix.prototype.storage = function () {
+	    // must be implemented by each of the Matrix implementations
+	    throw new Error('Cannot invoke storage on a Matrix interface');
+	  };
+	  
+	  /**
+	   * Get a subset of the matrix, or replace a subset of the matrix.
+	   *
+	   * Usage:
+	   *     var subset = matrix.subset(index)               // retrieve subset
+	   *     var value = matrix.subset(index, replacement)   // replace subset
+	   *
+	   * @param {Index} index
+	   * @param {Array | Matrix | *} [replacement]
+	   * @param {*} [defaultValue=0]      Default value, filled in on new entries when
+	   *                                  the matrix is resized. If not provided,
+	   *                                  new matrix elements will be filled with zeros.
+	   */
+	  Matrix.prototype.subset = function (index, replacement, defaultValue) {
+	    // must be implemented by each of the Matrix implementations
+	    throw new Error('Cannot invoke subset on a Matrix interface');
+	  };
+
+	  /**
+	   * Get a single element from the matrix.
+	   * @param {Number[]} index   Zero-based index
+	   * @return {*} value
+	   */
+	  Matrix.prototype.get = function (index) {
+	    // must be implemented by each of the Matrix implementations
+	    throw new Error('Cannot invoke get on a Matrix interface');
+	  };
+
+	  /**
+	   * Replace a single element in the matrix.
+	   * @param {Number[]} index   Zero-based index
+	   * @param {*} value
+	   * @param {*} [defaultValue]        Default value, filled in on new entries when
+	   *                                  the matrix is resized. If not provided,
+	   *                                  new matrix elements will be left undefined.
+	   * @return {Matrix} self
+	   */
+	  Matrix.prototype.set = function (index, value, defaultValue) {
+	    // must be implemented by each of the Matrix implementations
+	    throw new Error('Cannot invoke set on a Matrix interface');
+	  };
+
+	  /**
+	   * Resize the matrix to the given size. Returns a copy of the matrix when 
+	   * `copy=true`, otherwise return the matrix itself (resize in place).
+	   *
+	   * @param {Number[]} size           The new size the matrix should have.
+	   * @param {*} [defaultValue=0]      Default value, filled in on new entries.
+	   *                                  If not provided, the matrix elements will
+	   *                                  be filled with zeros.
+	   * @param {boolean} [copy]          Return a resized copy of the matrix
+	   *
+	   * @return {Matrix}                 The resized matrix
+	   */
+	  Matrix.prototype.resize = function (size, defaultValue) {
+	    // must be implemented by each of the Matrix implementations
+	    throw new Error('Cannot invoke resize on a Matrix interface');
+	  };
+
+	  /**
+	   * Create a clone of the matrix
+	   * @return {Matrix} clone
+	   */
+	  Matrix.prototype.clone = function () {
+	    // must be implemented by each of the Matrix implementations
+	    throw new Error('Cannot invoke clone on a Matrix interface');
+	  };
+
+	  /**
+	   * Retrieve the size of the matrix.
+	   * @returns {Number[]} size
+	   */
+	  Matrix.prototype.size = function() {
+	    // must be implemented by each of the Matrix implementations
+	    throw new Error('Cannot invoke size on a Matrix interface');
+	  };
+
+	  /**
+	   * Create a new matrix with the results of the callback function executed on
+	   * each entry of the matrix.
+	   * @param {function} callback   The callback function is invoked with three
+	   *                              parameters: the value of the element, the index
+	   *                              of the element, and the Matrix being traversed.
+	   * @param {boolean} [skipZeros] Invoke callback function for non-zero values only.
+	   *
+	   * @return {Matrix} matrix
+	   */
+	  Matrix.prototype.map = function (callback, skipZeros) {
+	    // must be implemented by each of the Matrix implementations
+	    throw new Error('Cannot invoke map on a Matrix interface');
+	  };
+
+	  /**
+	   * Execute a callback function on each entry of the matrix.
+	   * @param {function} callback   The callback function is invoked with three
+	   *                              parameters: the value of the element, the index
+	   *                              of the element, and the Matrix being traversed.
+	   */
+	  Matrix.prototype.forEach = function (callback) {
+	    // must be implemented by each of the Matrix implementations
+	    throw new Error('Cannot invoke forEach on a Matrix interface');
+	  };
+
+	  /**
+	   * Create an Array with a copy of the data of the Matrix
+	   * @returns {Array} array
+	   */
+	  Matrix.prototype.toArray = function () {
+	    // must be implemented by each of the Matrix implementations
+	    throw new Error('Cannot invoke toArray on a Matrix interface');
+	  };
+
+	  /**
+	   * Get the primitive value of the Matrix: a multidimensional array
+	   * @returns {Array} array
+	   */
+	  Matrix.prototype.valueOf = function () {
+	    // must be implemented by each of the Matrix implementations
+	    throw new Error('Cannot invoke valueOf on a Matrix interface');
+	  };
+
+	  /**
+	   * Get a string representation of the matrix, with optional formatting options.
+	   * @param {Object | Number | Function} [options]  Formatting options. See
+	   *                                                lib/util/number:format for a
+	   *                                                description of the available
+	   *                                                options.
+	   * @returns {String} str
+	   */
+	  Matrix.prototype.format = function (options) {
+	    // must be implemented by each of the Matrix implementations
+	    throw new Error('Cannot invoke format on a Matrix interface');
+	  };
+
+	  /**
+	   * Get a string representation of the matrix
+	   * @returns {String} str
+	   */
+	  Matrix.prototype.toString = function () {
+	    // must be implemented by each of the Matrix implementations
+	    throw new Error('Cannot invoke toString on a Matrix interface');
+	  };
+
+	  /**
+	   * Calculates the transpose of the matrix
+	   * @returns {Matrix}
+	   */
+	  Matrix.prototype.transpose = function () {
+	    // must be implemented by each of the Matrix implementations
+	    throw new Error('Cannot invoke transpose on a Matrix interface');
+	  };
+
+	  /**
+	   * Calculate the trace of a matrix: the sum of the elements on the main
+	   * diagonal of a square matrix.
+	   *
+	   * See also:
+	   *
+	   *    diagonal
+	   *
+	   * @returns {Number}       The matrix trace
+	   */
+	  Matrix.prototype.trace = function () {
+	    // must be implemented by each of the Matrix implementations
+	    throw new Error('Cannot invoke transpose on a Matrix interface');
+	  };
+
+	  /**
+	   * Multiply the matrix values times the argument.
+	   *
+	   * @param  {Number | BigNumber | Boolean | Complex | Unit | Array | Matrix | null} Value to multiply.
+	   *
+	   * @return {Number | BigNumber | Complex | Unit | Matrix}
+	   */
+	  Matrix.prototype.multiply = function (value) {
+	    // must be implemented by each of the Matrix implementations
+	    throw new Error('Cannot invoke multiply on a Matrix interface');
+	  };
+	  
+	  // exports
+	  return Matrix;
+	};
 
 /***/ },
 /* 11 */
@@ -2698,7 +2415,7 @@
 
 	'use strict';
 
-	var util = __webpack_require__(171),
+	var util = __webpack_require__(175),
 
 	    number = util.number,
 	    string = util.string,
@@ -3293,14 +3010,14 @@
 	  },
 	  BINARY_SHORT: {
 	    '': {name: '', value: 1, scientific: true},
-	    'k': {name: 'k', value: 1024, scientific: true},
-	    'M': {name: 'M', value: Math.pow(1024, 2), scientific: true},
-	    'G': {name: 'G', value: Math.pow(1024, 3), scientific: true},
-	    'T': {name: 'T', value: Math.pow(1024, 4), scientific: true},
-	    'P': {name: 'P', value: Math.pow(1024, 5), scientific: true},
-	    'E': {name: 'E', value: Math.pow(1024, 6), scientific: true},
-	    'Z': {name: 'Z', value: Math.pow(1024, 7), scientific: true},
-	    'Y': {name: 'Y', value: Math.pow(1024, 8), scientific: true},
+	    'k': {name: 'k', value: 1e3, scientific: true},
+	    'M': {name: 'M', value: 1e6, scientific: true},
+	    'G': {name: 'G', value: 1e9, scientific: true},
+	    'T': {name: 'T', value: 1e12, scientific: true},
+	    'P': {name: 'P', value: 1e15, scientific: true},
+	    'E': {name: 'E', value: 1e18, scientific: true},
+	    'Z': {name: 'Z', value: 1e21, scientific: true},
+	    'Y': {name: 'Y', value: 1e24, scientific: true},
 
 	    'Ki': {name: 'Ki', value: 1024, scientific: true},
 	    'Mi': {name: 'Mi', value: Math.pow(1024, 2), scientific: true},
@@ -3313,21 +3030,21 @@
 	  },
 	  BINARY_LONG: {
 	    '': {name: '', value: 1, scientific: true},
-	    'kilo': {name: 'kilo', value: 1024, scientific: true},
-	    'mega': {name: 'mega', value: Math.pow(1024, 2), scientific: true},
-	    'giga': {name: 'giga', value: Math.pow(1024, 3), scientific: true},
-	    'tera': {name: 'tera', value: Math.pow(1024, 4), scientific: true},
-	    'peta': {name: 'peta', value: Math.pow(1024, 5), scientific: true},
-	    'exa': {name: 'exa', value: Math.pow(1024, 6), scientific: true},
-	    'zetta': {name: 'zetta', value: Math.pow(1024, 7), scientific: true},
-	    'yotta': {name: 'yotta', value: Math.pow(1024, 8), scientific: true},
+	    'kilo': {name: 'kilo', value: 1e3, scientific: true},
+	    'mega': {name: 'mega', value: 1e6, scientific: true},
+	    'giga': {name: 'giga', value: 1e9, scientific: true},
+	    'tera': {name: 'tera', value: 1e12, scientific: true},
+	    'peta': {name: 'peta', value: 1e15, scientific: true},
+	    'exa':  {name: 'exa', value: 1e18, scientific: true},
+	    'zetta': {name: 'zetta', value: 1e21, scientific: true},
+	    'yotta': {name: 'yotta', value: 1e24, scientific: true},
 
 	    'kibi': {name: 'kibi', value: 1024, scientific: true},
 	    'mebi': {name: 'mebi', value: Math.pow(1024, 2), scientific: true},
 	    'gibi': {name: 'gibi', value: Math.pow(1024, 3), scientific: true},
 	    'tebi': {name: 'tebi', value: Math.pow(1024, 4), scientific: true},
 	    'pebi': {name: 'pebi', value: Math.pow(1024, 5), scientific: true},
-	    'exi': {name: 'exi', value: Math.pow(1024, 6), scientific: true},
+	    'exi':  {name: 'exi', value: Math.pow(1024, 6), scientific: true},
 	    'zebi': {name: 'zebi', value: Math.pow(1024, 7), scientific: true},
 	    'yobi': {name: 'yobi', value: Math.pow(1024, 8), scientific: true}
 	  }
@@ -3442,6 +3159,7 @@
 	  poundmass: {name: 'poundmass', base: BASE_UNITS.MASS, prefixes: PREFIXES.NONE, value: 453.59237e-3, offset: 0},
 	  hundredweight: {name: 'hundredweight', base: BASE_UNITS.MASS, prefixes: PREFIXES.NONE, value: 45.359237, offset: 0},
 	  stick: {name: 'stick', base: BASE_UNITS.MASS, prefixes: PREFIXES.NONE, value: 115e-3, offset: 0},
+	  stone: {name: 'stone', base: BASE_UNITS.MASS, prefixes: PREFIXES.NONE, value: 6350, offset: 0},
 
 	  gr: {name: 'gr', base: BASE_UNITS.MASS, prefixes: PREFIXES.NONE, value: 64.79891e-6, offset: 0},
 	  dr: {name: 'dr', base: BASE_UNITS.MASS, prefixes: PREFIXES.NONE, value: 1.7718451953125e-3, offset: 0},
@@ -3473,9 +3191,9 @@
 	  ampere: {name: 'ampere', base: BASE_UNITS.CURRENT, prefixes: PREFIXES.LONG, value: 1, offset: 0},
 
 	  // Temperature
-	  // K(C) = °C + 273.15
-	  // K(F) = (°F + 459.67) / 1.8
-	  // K(R) = °R / 1.8
+	  // K(C) = °C + 273.15
+	  // K(F) = (°F + 459.67) / 1.8
+	  // K(R) = °R / 1.8
 	  K: {name: 'K', base: BASE_UNITS.TEMPERATURE, prefixes: PREFIXES.NONE, value: 1, offset: 0},
 	  degC: {name: 'degC', base: BASE_UNITS.TEMPERATURE, prefixes: PREFIXES.NONE, value: 1, offset: 273.15},
 	  degF: {name: 'degF', base: BASE_UNITS.TEMPERATURE, prefixes: PREFIXES.NONE, value: 1/1.8, offset: 459.67},
@@ -3596,7 +3314,7 @@
 
 	'use strict';
 
-	var util = __webpack_require__(171);
+	var util = __webpack_require__(175);
 	var object = util.object;
 	var string = util.string;
 
@@ -3790,265 +3508,272 @@
 	// utility methods for arrays and matrices
 	'use strict';
 
-	var util = __webpack_require__(171),
-
-	    IndexError = __webpack_require__(169),
-	    DimensionError = __webpack_require__(168),
-
-	    Matrix = __webpack_require__(10),
+	var util = __webpack_require__(175),
+
+	    IndexError = __webpack_require__(173),
+	    DimensionError = __webpack_require__(172),
 
 	    array = util.array,
-	    isArray = util.array.isArray,
-	    isString = util.string.isString;
-
-	/**
-	 * Convert function arguments to an array. Arguments can have the following
-	 * signature:
-	 *     fn()
-	 *     fn(n)
-	 *     fn(m, n, p, ...)
-	 *     fn([m, n, p, ...])
-	 * @param {...Number | Array | Matrix} args
-	 * @returns {Array} array
-	 */
-	exports.argsToArray = function(args) {
-	  if (args.length == 0) {
-	    // fn()
-	    return [];
-	  }
-	  else if (args.length == 1) {
-	    // fn(n)
-	    // fn([m, n, p, ...])
-	    var array = args[0];
+	    isArray = util.array.isArray;
+
+	module.exports = function (math) {
+
+	  var Matrix = math.type.Matrix;
+
+	  var collection = {};
+
+	  /**
+	   * Convert function arguments to an array. Arguments can have the following
+	   * signature:
+	   *     fn()
+	   *     fn(n)
+	   *     fn(m, n, p, ...)
+	   *     fn([m, n, p, ...])
+	   * @param {...Number | Array | Matrix} args
+	   * @returns {Array} array
+	   */
+	  collection.argsToArray = function(args) {
+	    if (args.length === 0) {
+	      // fn()
+	      return [];
+	    }
+	    else if (args.length == 1) {
+	      // fn(n)
+	      // fn([m, n, p, ...])
+	      var array = args[0];
+	      if (array instanceof Matrix) {
+	        array = array.valueOf();
+	      }
+	      if (!isArray(array)) {
+	        array = [array];
+	      }
+	      return array;
+	    }
+	    else {
+	      // fn(m, n, p, ...)
+	      return util.array.argsToArray(args);
+	    }
+	  };
+
+
+	  /**
+	   * Test whether a value is a collection: an Array or Matrix
+	   * @param {*} x
+	   * @returns {boolean} isCollection
+	   */
+	  collection.isCollection = function(x) {
+	    return (isArray(x) || (x instanceof Matrix));
+	  };
+
+	  /**
+	   * Execute the callback function element wise for each element in array and any
+	   * nested array
+	   * Returns an array with the results
+	   * @param {Array | Matrix} array
+	   * @param {function} callback   The callback is called with two parameters:
+	   *                              value1 and value2, which contain the current
+	   *                              element of both arrays.
+	   * @param {boolean} [skipZeros] Invoke callback function for non-zero values only.
+	   *
+	   * @return {Array | Matrix} res
+	   */
+	  collection.deepMap = function deepMap(array, callback, skipZeros) {
+	    if (array && (typeof array.map === 'function')) {
+	      return array.map(function (x) {
+	        return deepMap(x, callback, skipZeros);
+	      });
+	    }
+	    else {
+	      return callback(array);
+	    }
+	  };
+
+	  /**
+	   * Execute the callback function element wise for each entry in two given arrays,
+	   * and for any nested array. Objects can also be scalar objects.
+	   * Returns an array with the results.
+	   * @param {Array | Matrix | Object} array1
+	   * @param {Array | Matrix | Object} array2
+	   * @param {function} callback   The callback is called with two parameters:
+	   *                              value1 and value2, which contain the current
+	   *                              element of both arrays.
+	   * @return {Array | Matrix} res
+	   */
+	  collection.deepMap2 = function deepMap2(array1, array2, callback) {
+	    var res, len, i;
+
+	    if (isArray(array1)) {
+	      if (isArray(array2)) {
+	        // callback(array, array)
+	        if (array1.length != array2.length) {
+	          throw new DimensionError(array1.length, array2.length);
+	        }
+
+	        res = [];
+	        len = array1.length;
+	        for (i = 0; i < len; i++) {
+	          res[i] = deepMap2(array1[i], array2[i], callback);
+	        }
+	      }
+	      else if (array2 instanceof Matrix) {
+	        // callback(array, matrix)
+	        res = deepMap2(array1, array2.valueOf(), callback);
+	        return math.matrix(res);
+	      }
+	      else {
+	        // callback(array, object)
+	        res = [];
+	        len = array1.length;
+	        for (i = 0; i < len; i++) {
+	          res[i] = deepMap2(array1[i], array2, callback);
+	        }
+	      }
+	    }
+	    else if (array1 instanceof Matrix) {
+	      if (array2 instanceof Matrix) {
+	        // callback(matrix, matrix)
+	        res = deepMap2(array1.valueOf(), array2.valueOf(), callback);
+	        return math.matrix(res);
+	      }
+	      else {
+	        // callback(matrix, array)
+	        // callback(matrix, object)
+	        res = deepMap2(array1.valueOf(), array2, callback);
+	        return math.matrix(res);
+	      }
+	    }
+	    else {
+	      if (isArray(array2)) {
+	        // callback(object, array)
+	        res = [];
+	        len = array2.length;
+	        for (i = 0; i < len; i++) {
+	          res[i] = deepMap2(array1, array2[i], callback);
+	        }
+	      }
+	      else if (array2 instanceof Matrix) {
+	        // callback(object, matrix)
+	        res = deepMap2(array1, array2.valueOf(), callback);
+	        return math.matrix(res);
+	      }
+	      else {
+	        // callback(object, object)
+	        res = callback(array1, array2);
+	      }
+	    }
+
+	    return res;
+	  };
+
+	  /**
+	   * Reduce a given matrix or array to a new matrix or
+	   * array with one less dimension, applying the given
+	   * callback in the selected dimension.
+	   * @param {Array | Matrix} mat
+	   * @param {Number} dim
+	   * @param {function} callback
+	   * @return {Array | Matrix} res
+	   */
+	  collection.reduce = function(mat, dim, callback) {
+	    var size = isArray(mat) ? array.size(mat) : mat.size();
+	    if (dim < 0) {
+	      // TODO: would be more clear when throwing a DimensionError here
+	      throw new IndexError(dim);
+	    }
+	    if (dim >= size.length) {
+	      // TODO: would be more clear when throwing a DimensionError here
+	      throw new IndexError(dim, size.length);
+	    }
+
+	      if (mat instanceof Matrix) {
+	          return math.matrix(_reduce(mat.valueOf(), dim, callback));
+	      }else {
+	          return _reduce(mat, dim, callback);
+	      }
+	  };
+
+	  /**
+	   * Recursively reduce a matrix
+	   * @param {Array} mat
+	   * @param {Number} dim
+	   * @param {Function} callback
+	   * @returns {Array} ret
+	   * @private
+	   */
+	  function _reduce(mat, dim, callback){
+	    var i, ret, val, tran;
+
+	      if(dim<=0){
+	          if( !isArray(mat[0]) ){
+	              val = mat[0];
+	              for(i=1; i<mat.length; i++){
+	                  val = callback(val, mat[i]);
+	              }
+	              return val;
+	          }else{
+	              tran = _switch(mat);
+	              ret = [];
+	              for(i=0; i<tran.length; i++){
+	                  ret[i] = _reduce(tran[i], dim-1, callback);
+	              }
+	              return ret;
+	          }
+	      }else{
+	          ret = [];
+	          for(i=0; i<mat.length; i++){
+	              ret[i] = _reduce(mat[i], dim-1, callback);
+	          }
+	          return ret;
+	      }
+	  }
+
+	  /**
+	   * Transpose a matrix
+	   * @param {Array} mat
+	   * @returns {Array} ret
+	   * @private
+	   */
+	  function _switch(mat){
+	    var I = mat.length;
+	    var J = mat[0].length;
+	    var i, j;
+	    var ret = [];
+	    for( j=0; j<J; j++) {
+	      var tmp = [];
+	      for( i=0; i<I; i++) {
+	        tmp.push(mat[i][j]);
+	      }
+	      ret.push(tmp);
+	    }
+	    return ret;
+	  }
+
+	  /**
+	   * Recursively loop over all elements in a given multi dimensional array
+	   * and invoke the callback on each of the elements.
+	   * @param {Array | Matrix} array
+	   * @param {function} callback     The callback method is invoked with one
+	   *                                parameter: the current element in the array
+	   */
+	  collection.deepForEach = function deepForEach (array, callback) {
 	    if (array instanceof Matrix) {
 	      array = array.valueOf();
 	    }
-	    if (!isArray(array)) {
-	      array = [array];
-	    }
-	    return array;
-	  }
-	  else {
-	    // fn(m, n, p, ...)
-	    return util.array.argsToArray(args);
-	  }
-	};
-
-
-	/**
-	 * Test whether a value is a collection: an Array or Matrix
-	 * @param {*} x
-	 * @returns {boolean} isCollection
-	 */
-	exports.isCollection = function(x) {
-	  return (isArray(x) || (x instanceof Matrix));
-	};
-
-	/**
-	 * Execute the callback function element wise for each element in array and any
-	 * nested array
-	 * Returns an array with the results
-	 * @param {Array | Matrix} array
-	 * @param {function} callback   The callback is called with two parameters:
-	 *                              value1 and value2, which contain the current
-	 *                              element of both arrays.
-	 * @return {Array | Matrix} res
-	 */
-	exports.deepMap = function deepMap(array, callback) {
-	  if (array && (typeof array.map === 'function')) {
-	    return array.map(function (x) {
-	      return deepMap(x, callback);
-	    });
-	  }
-	  else {
-	    return callback(array);
-	  }
-	};
-
-	/**
-	 * Execute the callback function element wise for each entry in two given arrays,
-	 * and for any nested array. Objects can also be scalar objects.
-	 * Returns an array with the results.
-	 * @param {Array | Matrix | Object} array1
-	 * @param {Array | Matrix | Object} array2
-	 * @param {function} callback   The callback is called with two parameters:
-	 *                              value1 and value2, which contain the current
-	 *                              element of both arrays.
-	 * @return {Array | Matrix} res
-	 */
-	exports.deepMap2 = function deepMap2(array1, array2, callback) {
-	  var res, len, i;
-
-	  if (isArray(array1)) {
-	    if (isArray(array2)) {
-	      // callback(array, array)
-	      if (array1.length != array2.length) {
-	        throw new DimensionError(array1.length, array2.length);
-	      }
-
-	      res = [];
-	      len = array1.length;
-	      for (i = 0; i < len; i++) {
-	        res[i] = deepMap2(array1[i], array2[i], callback);
-	      }
-	    }
-	    else if (array2 instanceof Matrix) {
-	      // callback(array, matrix)
-	      res = deepMap2(array1, array2.valueOf(), callback);
-	      return new Matrix(res);
-	    }
-	    else {
-	      // callback(array, object)
-	      res = [];
-	      len = array1.length;
-	      for (i = 0; i < len; i++) {
-	        res[i] = deepMap2(array1[i], array2, callback);
-	      }
-	    }
-	  }
-	  else if (array1 instanceof Matrix) {
-	    if (array2 instanceof Matrix) {
-	      // callback(matrix, matrix)
-	      res = deepMap2(array1.valueOf(), array2.valueOf(), callback);
-	      return new Matrix(res);
-	    }
-	    else {
-	      // callback(matrix, array)
-	      // callback(matrix, object)
-	      res = deepMap2(array1.valueOf(), array2, callback);
-	      return new Matrix(res);
-	    }
-	  }
-	  else {
-	    if (isArray(array2)) {
-	      // callback(object, array)
-	      res = [];
-	      len = array2.length;
-	      for (i = 0; i < len; i++) {
-	        res[i] = deepMap2(array1, array2[i], callback);
-	      }
-	    }
-	    else if (array2 instanceof Matrix) {
-	      // callback(object, matrix)
-	      res = deepMap2(array1, array2.valueOf(), callback);
-	      return new Matrix(res);
-	    }
-	    else {
-	      // callback(object, object)
-	      res = callback(array1, array2);
-	    }
-	  }
-
-	  return res;
-	};
-
-	/**
-	 * Reduce a given matrix or array to a new matrix or
-	 * array with one less dimension, applying the given
-	 * callback in the selected dimension.
-	 * @param {Array | Matrix} mat
-	 * @param {Number} dim
-	 * @param {function} callback
-	 * @return {Array | Matrix} res
-	 */
-	exports.reduce = function(mat, dim, callback) {
-	  var size = isArray(mat) ? array.size(mat) : mat.size();
-	  if (dim < 0) {
-	    // TODO: would be more clear when throwing a DimensionError here
-	    throw new IndexError(dim);
-	  }
-	  if (dim >= size.length) {
-	    // TODO: would be more clear when throwing a DimensionError here
-	    throw new IndexError(dim, size.length);
-	  }
-
-		if (mat instanceof Matrix) {
-			return new Matrix(_reduce(mat.valueOf(), dim, callback));
-		}else {
-			return _reduce(mat, dim, callback);
-		}
-	};
-
-	/**
-	 * Recursively reduce a matrix
-	 * @param {Array} mat
-	 * @param {Number} dim
-	 * @param {Function} callback
-	 * @returns {Array} ret
-	 * @private
-	 */
-	function _reduce(mat, dim, callback){
-	  var i, ret, val, tran;
-
-		if(dim<=0){
-			if( !isArray(mat[0]) ){
-				val = mat[0];
-				for(i=1; i<mat.length; i++){
-					val = callback(val, mat[i]);
-				}
-				return val;
-			}else{
-				tran = _switch(mat);
-				ret = [];
-				for(i=0; i<tran.length; i++){
-					ret[i] = _reduce(tran[i], dim-1, callback);
-				}
-				return ret
-			}
-		}else{
-			ret = [];
-			for(i=0; i<mat.length; i++){
-				ret[i] = _reduce(mat[i], dim-1, callback);
-			}
-			return ret;
-		}
-	}
-
-	/**
-	 * Transpose a matrix
-	 * @param {Array} mat
-	 * @returns {Array} ret
-	 * @private
-	 */
-	function _switch(mat){
-	  var I = mat.length;
-	  var J = mat[0].length;
-	  var i, j;
-	  var ret = [];
-	  for( j=0; j<J; j++) {
-	    var tmp = [];
-	    for( i=0; i<I; i++) {
-	      tmp.push(mat[i][j]);
-	    }
-	    ret.push(tmp);
-	  }
-	  return ret;
-	}
-
-	/**
-	 * Recursively loop over all elements in a given multi dimensional array
-	 * and invoke the callback on each of the elements.
-	 * @param {Array | Matrix} array
-	 * @param {function} callback     The callback method is invoked with one
-	 *                                parameter: the current element in the array
-	 */
-	exports.deepForEach = function deepForEach (array, callback) {
-	  if (array instanceof Matrix) {
-	    array = array.valueOf();
-	  }
-
-	  for (var i = 0, ii = array.length; i < ii; i++) {
-	    var value = array[i];
-
-	    if (isArray(value)) {
-	      deepForEach(value, callback);
-	    }
-	    else {
-	      callback(value);
-	    }
-	  }
-	};
-
+
+	    for (var i = 0, ii = array.length; i < ii; i++) {
+	      var value = array[i];
+
+	      if (isArray(value)) {
+	        deepForEach(value, callback);
+	      }
+	      else {
+	        callback(value);
+	      }
+	    }
+	  };
+	  
+	  return collection;
+	};
 
 /***/ },
 /* 15 */
@@ -4056,20 +3781,1214 @@
 
 	'use strict';
 
-	exports.ArrayNode = __webpack_require__(172);
-	exports.AssignmentNode = __webpack_require__(173);
-	exports.BlockNode = __webpack_require__(174);
-	exports.ConditionalNode = __webpack_require__(175);
-	exports.ConstantNode = __webpack_require__(176);
-	exports.IndexNode = __webpack_require__(177);
-	exports.FunctionAssignmentNode = __webpack_require__(178);
-	exports.FunctionNode = __webpack_require__(179);
-	exports.Node = __webpack_require__(180);
-	exports.OperatorNode = __webpack_require__(181);
-	exports.RangeNode = __webpack_require__(182);
-	exports.SymbolNode = __webpack_require__(183);
-	exports.UpdateNode = __webpack_require__(184);
-
+	var util = __webpack_require__(175);
+	var DimensionError = __webpack_require__(172);
+
+	var array = util.array;
+	var object = util.object;
+	var string = util.string;
+	var number = util.number;
+
+	var isArray = Array.isArray;
+	var isNumber = util.number.isNumber;
+	var isInteger = util.number.isInteger;
+
+	var validateIndex = array.validateIndex;
+
+	module.exports = function (math) {
+	  
+	  var Index = math.type.Index,
+	      BigNumber = math.type.BigNumber,
+	      Matrix = math.type.Matrix;
+	  
+	  function CcsMatrix(data) {
+	    if (!(this instanceof CcsMatrix))
+	      throw new SyntaxError('Constructor must be called with the new operator');
+
+	    if (data instanceof Matrix) {
+	      // check data is a CcsMatrix
+	      if (data.type === 'CcsMatrix') {
+	        // clone arrays
+	        this._values = object.clone(data._values);
+	        this._index = object.clone(data._index);
+	        this._ptr = object.clone(data._ptr);
+	        this._size = object.clone(data._size);
+	      }
+	      else {
+	        // build from matrix data
+	        _createFromArray(this, data.valueOf());
+	      }
+	    }
+	    else if (data && isArray(data.values) && isArray(data.index) && isArray(data.ptr) && isArray(data.size)) {
+	      // initialize fields
+	      this._values = data.values;
+	      this._index = data.index;
+	      this._ptr = data.ptr;
+	      this._size = data.size;
+	    }
+	    else if (isArray(data)) {
+	      // create from array
+	      _createFromArray(this, data);
+	    }
+	    else if (data) {
+	      // unsupported type
+	      throw new TypeError('Unsupported type of data (' + util.types.type(data) + ')');
+	    }
+	    else {
+	      // nothing provided
+	      this._values = [];
+	      this._index = [];
+	      this._ptr = [0];
+	      this._size = [0];
+	    }
+	  }
+	  
+	  var _createFromArray = function (matrix, data) {
+	    // initialize fields
+	    matrix._values = [];
+	    matrix._index = [];
+	    matrix._ptr = [];
+	    // discover rows & columns, do not use math.size() to avoid looping array twice
+	    var rows = data.length;
+	    var columns = 0;
+
+	    // check we have rows (empty array)
+	    if (rows > 0) {
+	      // column index
+	      var j = 0;
+	      do {
+	        // store pointer to values index
+	        matrix._ptr.push(matrix._values.length);
+	        // loop rows
+	        for (var i = 0; i < rows; i++) {
+	          // current row
+	          var row = data[i];
+	          // check row is an array
+	          if (isArray(row)) {
+	            // update columns if needed (only on first column)
+	            if (j ===0 && columns < row.length)
+	              columns = row.length;
+	            // check row has column
+	            if (j < row.length) {
+	              // value
+	              var v = row[j];
+	              // check value != 0
+	              if (!math.equal(v, 0)) {
+	                // store value
+	                matrix._values.push(v);
+	                // index
+	                matrix._index.push(i);
+	              }
+	            }
+	          }
+	          else {
+	            // update columns if needed (only on first column)
+	            if (j === 0 && columns < 1)
+	              columns = 1;
+	            // check value != 0 (row is a scalar)
+	            if (!math.equal(row, 0)) {
+	              // store value
+	              matrix._values.push(row);
+	              // index
+	              matrix._index.push(i);
+	            }
+	          }
+	        }
+	        // increment index
+	        j++;      
+	      }
+	      while (j < columns);
+	    }
+	    // store number of values in ptr
+	    matrix._ptr.push(matrix._values.length);
+	    // size
+	    matrix._size = [rows, columns];
+	  };
+	  
+	  CcsMatrix.prototype = new math.type.Matrix();
+
+	  CcsMatrix.prototype.type = 'CcsMatrix';
+	  
+	  /**
+	   * Get the storage format used by the matrix.
+	   *
+	   * Usage:
+	   *     var format = matrix.storage()                   // retrieve storage format
+	   *
+	   * @return {string}           The storage format.
+	   */
+	  CcsMatrix.prototype.storage = function () {
+	    return 'ccs';
+	  };
+	  
+	  /**
+	   * Get a subset of the matrix, or replace a subset of the matrix.
+	   *
+	   * Usage:
+	   *     var subset = matrix.subset(index)               // retrieve subset
+	   *     var value = matrix.subset(index, replacement)   // replace subset
+	   *
+	   * @param {Index} index
+	   * @param {Array | Maytrix | *} [replacement]
+	   * @param {*} [defaultValue=0]      Default value, filled in on new entries when
+	   *                                  the matrix is resized. If not provided,
+	   *                                  new matrix elements will be filled with zeros.
+	   */
+	  CcsMatrix.prototype.subset = function (index, replacement, defaultValue) {
+	    // check arguments
+	    switch (arguments.length) {
+	      case 1:
+	        return _getsubset(this, index);
+
+	        // intentional fall through
+	      case 2:
+	      case 3:
+	        return _setsubset(this, index, replacement, defaultValue);
+
+	      default:
+	        throw new SyntaxError('Wrong number of arguments');
+	    }
+	  };
+	  
+	  var _getsubset = function (matrix, index) {
+	    // check index
+	    if (!(index instanceof Index)) {
+	      throw new TypeError('Invalid index');
+	    }
+
+	    var isScalar = index.isScalar();
+	    if (isScalar) {
+	      // return a scalar
+	      return matrix.get(index.min());
+	    }
+	    // validate dimensions
+	    var size = index.size();
+	    if (size.length != matrix._size.length) {
+	      throw new DimensionError(size.length, matrix._size.length);
+	    }
+
+	    // validate if any of the ranges in the index is out of range
+	    var min = index.min();
+	    var max = index.max();
+	    for (var i = 0, ii = matrix._size.length; i < ii; i++) {
+	      validateIndex(min[i], matrix._size[i]);
+	      validateIndex(max[i], matrix._size[i]);
+	    }
+
+	    // map callback
+	    var callback = function (v) {
+	      // return value
+	      return v;
+	    };
+	    // get sub-matrix
+	    return _map(matrix, min[0], max[0], min[1], max[1], callback, false);
+	  };
+	  
+	  var _setsubset = function (matrix, index, submatrix, defaultValue) {
+	    // check index
+	    if (!(index instanceof Index)) {
+	      throw new TypeError('Invalid index');
+	    }
+	    
+	    // get index size and check whether the index contains a single value
+	    var iSize = index.size(),
+	        isScalar = index.isScalar();
+	    
+	    // calculate the size of the submatrix, and convert it into an Array if needed
+	    var sSize;
+	    if (submatrix instanceof Matrix) {
+	      // submatrix size
+	      sSize = submatrix.size();
+	      // use array representation
+	      submatrix = submatrix.toArray();
+	    }
+	    else {
+	      // get submatrix size (array, scalar)
+	      sSize = array.size(submatrix);
+	    }
+	    
+	    // check index is a scalar
+	    if (isScalar) {
+	      // verify submatrix is a scalar
+	      if (sSize.length !== 0) {
+	        throw new TypeError('Scalar expected');
+	      }
+	      // set value
+	      matrix.set(index.min(), submatrix, defaultValue);
+	    }
+	    else {
+	      // validate dimensions, index size must be one or two dimensions
+	      if (iSize.length !== 1 && iSize.length !== 2) {
+	        throw new DimensionError(iSize.length, matrix._size.length, '<');
+	      }
+	      
+	      // check submatrix and index have the same dimensions
+	      if (sSize.length < iSize.length) {
+	        // calculate number of missing outer dimensions
+	        var i = 0;
+	        var outer = 0;
+	        while (iSize[i] === 1 && sSize[i] === 1) {
+	          i++;
+	        }
+	        while (iSize[i] === 1) {
+	          outer++;
+	          i++;
+	        }
+	        // unsqueeze both outer and inner dimensions
+	        submatrix = array.unsqueeze(submatrix, iSize.length, outer, sSize);
+	      }
+	      
+	      // check whether the size of the submatrix matches the index size
+	      if (!object.deepEqual(iSize, sSize)) {
+	        throw new DimensionError(iSize, sSize, '>');
+	      }
+	      
+	      // offsets
+	      var x0 = index.min()[0];
+	      var y0 = index.min()[1];      
+	      
+	      // submatrix rows and columns
+	      var m = sSize[0];
+	      var n = sSize[1];
+
+	      // loop submatrix
+	      for (var x = 0; x < m; x++) {
+	        // loop columns
+	        for (var y = 0; y < n; y++) {
+	          // value at i, j
+	          var v = submatrix[x][y];
+	          // invoke set (zero value will remove entry from matrix)
+	          matrix.set([x + x0, y + y0], v, defaultValue);
+	        }
+	      }
+	    }
+	    return matrix;
+	  };
+
+	  /**
+	   * Get a single element from the matrix.
+	   * @param {Number[]} index   Zero-based index
+	   * @return {*} value
+	   */
+	  CcsMatrix.prototype.get = function (index) {
+	    if (!isArray(index))
+	      throw new TypeError('Array expected');
+	    if (index.length != this._size.length)
+	      throw new DimensionError(index.length, this._size.length);
+
+	    // row and column
+	    var i = index[0];
+	    var j = index[1];
+
+	    // check i, j are valid
+	    validateIndex(i, this._size[0]);
+	    validateIndex(j, this._size[1]);
+
+	    // find value index
+	    var k = _getValueIndex(i, this._ptr[j], this._ptr[j + 1], this._index);
+	    // check k is prior to next column k and it is in the correct row
+	    if (k < this._ptr[j + 1] && this._index[k] === i)
+	      return object.clone(this._values[k]);
+
+	    return 0;
+	  };
+	  
+	  /**
+	   * Replace a single element in the matrix.
+	   * @param {Number[]} index   Zero-based index
+	   * @param {*} value
+	   * @param {*} [defaultValue]        Default value, filled in on new entries when
+	   *                                  the matrix is resized. If not provided,
+	   *                                  new matrix elements will be set to zero.
+	   * @return {CcsMatrix} self
+	   */
+	  CcsMatrix.prototype.set = function (index, v, defaultValue) {
+	    if (!isArray(index))
+	      throw new TypeError('Array expected');
+	    if (index.length != this._size.length)
+	      throw new DimensionError(index.length, this._size.length);
+
+	    // row and column
+	    var i = index[0];
+	    var j = index[1];
+
+	    // rows & columns
+	    var rows = this._size[0];
+	    var columns = this._size[1];
+
+	    // check we need to resize matrix
+	    if (i > rows - 1 || j > columns - 1) {
+	      // resize matrix
+	      _resize(this, Math.max(i + 1, rows), Math.max(j + 1, columns), defaultValue);
+	      // update rows & columns
+	      rows = this._size[0];
+	      columns = this._size[1];
+	    }
+
+	    // check i, j are valid
+	    validateIndex(i, rows);
+	    validateIndex(j, columns);
+
+	    // find value index
+	    var k = _getValueIndex(i, this._ptr[j], this._ptr[j + 1], this._index);
+	    // check k is prior to next column k and it is in the correct row
+	    if (k < this._ptr[j + 1] && this._index[k] === i) {
+	      // check value != 0
+	      if (!math.equal(v, 0)) {
+	        // update value
+	        this._values[k] = v;
+	      }
+	      else {
+	        // remove value from matrix
+	        _remove(k, j, this._values, this._index, this._ptr);
+	      }
+	    }
+	    else {
+	      // insert value @ (i, j)
+	      _insert(k, i, j, v, this._values, this._index, this._ptr);
+	    }
+
+	    return this;
+	  };
+	  
+	  var _getValueIndex = function(i, top, bottom, index) {
+	    // check row is on the bottom side
+	    if (bottom - top === 0 || i > index[bottom - 1])
+	      return bottom;
+	    // loop until we find row index
+	    while (top < bottom) {
+	      // point in the middle (fast integer division)
+	      var p = ~~((top + bottom) / 2);
+	      // row @ p
+	      var r = index[p];
+	      // check we have to look on the top side, bottom side or we found the row
+	      if (i < r)
+	        bottom = p;
+	      else if (i > r)
+	        top = p + 1;
+	      else
+	        return p;
+	    }
+	    return top;
+	  };
+
+	  var _remove = function (k, j, values, index, ptr) {
+	    // remove value @ k
+	    values.splice(k, 1);
+	    index.splice(k, 1);
+	    // update pointers
+	    for (var x = j + 1; x < ptr.length; x++)
+	      ptr[x]--;
+	  };
+
+	  var _insert = function (k, i, j, v, values, index, ptr) {
+	    // insert value
+	    values.splice(k, 0, v);
+	    // update row for k
+	    index.splice(k, 0, i);
+	    // update column pointers
+	    for (var x = j + 1; x < ptr.length; x++)
+	      ptr[x]++;
+	  };
+	  
+	  /**
+	   * Resize the matrix to the given size. Returns a copy of the matrix when 
+	   * `copy=true`, otherwise return the matrix itself (resize in place).
+	   *
+	   * @param {Number[]} size           The new size the matrix should have.
+	   * @param {*} [defaultValue=0]      Default value, filled in on new entries.
+	   *                                  If not provided, the matrix elements will
+	   *                                  be filled with zeros.
+	   * @param {boolean} [copy]          Return a resized copy of the matrix
+	   *
+	   * @return {Matrix}                 The resized matrix
+	   */
+	  CcsMatrix.prototype.resize = function (size, defaultValue, copy) {    
+	    // validate arguments
+	    if (!isArray(size))
+	      throw new TypeError('Array expected');
+	    if (size.length !== 2)
+	      throw new Error('Only two dimensions matrix are supported');
+
+	    // check sizes
+	    size.forEach(function (value) {
+	      if (!number.isNumber(value) || !number.isInteger(value) || value < 0) {
+	        throw new TypeError('Invalid size, must contain positive integers ' +
+	                            '(size: ' + string.format(size) + ')');
+	      }
+	    });
+	    
+	    // matrix to resize
+	    var m = copy ? this.clone() : this;
+	    // resize matrix
+	    return _resize(m, size[0], size[1], defaultValue);
+	  };
+	  
+	  var _resize = function (matrix, rows, columns, defaultValue) {
+	    // value to insert at the time of growing matrix
+	    var value = defaultValue || 0;
+	    // should we insert the value?
+	    var ins = !math.equal(value, 0);
+
+	    // old columns and rows
+	    var r = matrix._size[0];
+	    var c = matrix._size[1];
+
+	    var i, j, k;
+
+	    // check we need to increase columns
+	    if (columns > c) {
+	      // loop new columns
+	      for (j = c; j < columns; j++) {
+	        // update matrix._ptr for current column
+	        matrix._ptr[j] = matrix._values.length;
+	        // check we need to insert matrix._values
+	        if (ins) {
+	          // loop rows
+	          for (i = 0; i < r; i++) {
+	            // add new matrix._values
+	            matrix._values.push(value);
+	            // update matrix._index
+	            matrix._index.push(i);
+	          }
+	        }        
+	      }
+	      // store number of matrix._values in matrix._ptr
+	      matrix._ptr[columns] = matrix._values.length;
+	    }
+	    else if (columns < c) {
+	      // truncate matrix._ptr
+	      matrix._ptr.splice(columns + 1, c - columns);
+	      // truncate matrix._values and matrix._index
+	      matrix._values.splice(matrix._ptr[columns], matrix._values.length);
+	      matrix._index.splice(matrix._ptr[columns], matrix._index.length);
+	    }
+	    // update columns
+	    c = columns;
+
+	    // check we need to increase rows
+	    if (rows > r) {
+	      // check we have to insert values
+	      if (ins) {
+	        // inserts
+	        var n = 0;
+	        // loop columns
+	        for (j = 0; j < c; j++) {
+	          // update matrix._ptr for current column
+	          matrix._ptr[j] = matrix._ptr[j] + n;
+	          // where to insert matrix._values
+	          k = matrix._ptr[j + 1] + n;
+	          // pointer
+	          var p = 0;
+	          // loop new rows, initialize pointer
+	          for (i = r; i < rows; i++, p++) {
+	            // add value
+	            matrix._values.splice(k + p, 0, value);
+	            // update matrix._index
+	            matrix._index.splice(k + p, 0, i);
+	            // increment inserts
+	            n++;
+	          }
+	        }
+	        // store number of matrix._values in matrix._ptr
+	        matrix._ptr[c] = matrix._values.length;
+	      }
+	    }
+	    else if (rows < r) {
+	      // deletes
+	      var d = 0;
+	      // loop columns
+	      for (j = 0; j < c; j++) {
+	        // update matrix._ptr for current column
+	        matrix._ptr[j] = matrix._ptr[j] - d;
+	        // where matrix._values start for next column
+	        var k0 = matrix._ptr[j];
+	        var k1 = matrix._ptr[j + 1] - d;
+	        // loop matrix._index
+	        for (k = k0; k < k1; k++) {
+	          // row
+	          i = matrix._index[k];
+	          // check we need to delete value and matrix._index
+	          if (i > rows - 1) {
+	            // remove value
+	            matrix._values.splice(k, 1);
+	            // remove item from matrix._index
+	            matrix._index.splice(k, 1);
+	            // increase deletes
+	            d++;
+	          }
+	        }
+	      }
+	      // update matrix._ptr for current column
+	      matrix._ptr[j] = matrix._values.length;
+	    }
+	    // update matrix._size
+	    matrix._size[0] = rows;
+	    matrix._size[1] = columns;
+	    // return matrix
+	    return matrix;
+	  };
+	  
+	  /**
+	   * Create a clone of the matrix
+	   * @return {CcsMatrix} clone
+	   */
+	  CcsMatrix.prototype.clone = function () {
+	    var m = new CcsMatrix({
+	      values: object.clone(this._values),
+	      index: object.clone(this._index),
+	      ptr: object.clone(this._ptr),
+	      size: object.clone(this._size)
+	    });
+	    return m;
+	  };
+	  
+	  /**
+	   * Retrieve the size of the matrix.
+	   * @returns {Number[]} size
+	   */
+	  CcsMatrix.prototype.size = function() {
+	    return object.clone(this._size);
+	  };
+	  
+	  /**
+	   * Create a new matrix with the results of the callback function executed on
+	   * each entry of the matrix.
+	   * @param {function} callback   The callback function is invoked with three
+	   *                              parameters: the value of the element, the index
+	   *                              of the element, and the Matrix being traversed.
+	   * @param {boolean} [skipZeros] Invoke callback function for non-zero values only.
+	   *
+	   * @return {CcsMatrix} matrix
+	   */
+	  CcsMatrix.prototype.map = function (callback, skipZeros) {
+	    // matrix instance
+	    var me = this;
+	    // rows and columns
+	    var rows = this._size[0];
+	    var columns = this._size[1];
+	    // invoke callback
+	    var invoke = function (v, i, j) {
+	      // invoke callback
+	      return callback(v, [i, j], me);
+	    };
+	    // invoke _map
+	    return _map(this, 0, rows - 1, 0, columns - 1, invoke, skipZeros);
+	  };
+
+	  /**
+	   * Create a new matrix with the results of the callback function executed on the interval
+	   * [minRow..maxRow, minColumn..maxColumn].
+	   */
+	  var _map = function (matrix, minRow, maxRow, minColumn, maxColumn, callback, skipZeros) {
+	    // result arrays
+	    var values = [];
+	    var index = [];
+	    var ptr = [];
+	    // invoke callback
+	    var invoke = function (v, x, y) {
+	      // invoke callback
+	      v = callback(v, x, y);
+	      // check value != 0
+	      if (!math.equal(v, 0)) {
+	        // store value
+	        values.push(v);
+	        // index
+	        index.push(x);
+	      }
+	    };
+	    // loop columns
+	    for (var j = minColumn; j <= maxColumn; j++) {
+	      // store pointer to values index
+	      ptr.push(values.length);
+	      // k0 <= k < k1 where k0 = _ptr[j] && k1 = _ptr[j+1]
+	      var k0 = matrix._ptr[j];
+	      var k1 = matrix._ptr[j + 1];
+	      // row pointer
+	      var p = minRow;
+	      // loop k within [k0, k1[
+	      for (var k = k0; k < k1; k++) {
+	        // row index
+	        var i = matrix._index[k];
+	        // check i is in range
+	        if (i >= minRow && i <= maxRow) {
+	          // zero values
+	          if (!skipZeros) {
+	           for (var x = p; x < i; x++)
+	             invoke(0, x - minRow, j - minColumn);
+	          }
+	          // value @ k
+	          invoke(matrix._values[k], i - minRow, j - minColumn);
+	        }
+	        // update pointer
+	        p = i + 1;
+	      }
+	      // zero values
+	      if (!skipZeros) {
+	        for (var y = p; y <= maxRow; y++)
+	          invoke(0, y - minRow, j - minColumn);
+	      }
+	    }
+	    // store number of values in ptr
+	    ptr.push(values.length);
+	    // return ccs
+	    return new CcsMatrix({
+	      values: values,
+	      index: index,
+	      ptr: ptr,
+	      size: [maxRow - minRow + 1, maxColumn - minColumn + 1]
+	    });
+	  };
+	  
+	  /**
+	   * Execute a callback function on each entry of the matrix.
+	   * @param {function} callback   The callback function is invoked with three
+	   *                              parameters: the value of the element, the index
+	   *                              of the element, and the Matrix being traversed.
+	   * @param {boolean} [skipZeros] Invoke callback function for non-zero values only.
+	   */
+	  CcsMatrix.prototype.forEach = function (callback, skipZeros) {
+	    // matrix instance
+	    var me = this;
+	    // rows and columns
+	    var rows = this._size[0];
+	    var columns = this._size[1];
+	    // loop columns
+	    for (var j = 0; j < columns; j++) {
+	      // k0 <= k < k1 where k0 = _ptr[j] && k1 = _ptr[j+1]
+	      var k0 = this._ptr[j];
+	      var k1 = this._ptr[j + 1];
+	      // column pointer
+	      var p = 0;
+	      // loop k within [k0, k1[
+	      for (var k = k0; k < k1; k++) {
+	        // row index
+	        var i = this._index[k];
+	        // check we need to process zeros
+	        if (!skipZeros) {
+	          // zero values
+	          for (var x = p; x < i; x++)
+	            callback(0, [x, j], me);
+	        }
+	        // value @ k
+	        callback(this._values[k], [i, j], me);
+	        // update pointer
+	        p = i + 1;
+	      }
+	      // check we need to process zeros
+	      if (!skipZeros) {
+	        // zero values
+	        for (var y = p; y < rows; y++)
+	          callback(0, [y, j], me);
+	      }
+	    }
+	  };
+	  
+	  /**
+	   * Create an Array with a copy of the data of the CcsMatrix
+	   * @returns {Array} array
+	   */
+	  CcsMatrix.prototype.toArray = function () {
+	    return _toArray(this, true);
+	  };
+
+	  /**
+	   * Get the primitive value of the CcsMatrix: a two dimensions array
+	   * @returns {Array} array
+	   */
+	  CcsMatrix.prototype.valueOf = function () {
+	    return _toArray(this, false);
+	  };
+	  
+	  var _toArray = function (matrix, copy) {
+	    // result
+	    var a = [];
+	    // rows and columns
+	    var rows = matrix._size[0];
+	    var columns = matrix._size[1];
+	    // loop columns
+	    for (var j = 0; j < columns; j++) {
+	      // k0 <= k < k1 where k0 = _ptr[j] && k1 = _ptr[j+1]
+	      var k0 = matrix._ptr[j];
+	      var k1 = matrix._ptr[j + 1];
+	      // row pointer
+	      var p = 0;
+	      // loop k within [k0, k1[
+	      for (var k = k0; k < k1; k++) {
+	        // row index
+	        var i = matrix._index[k];
+	        // zeros
+	        for (var x = p; x < i; x++)
+	          (a[x] = (a[x] || []))[j] = 0;
+	        // set value
+	        (a[i] = (a[i] || []))[j] = copy ? object.clone(matrix._values[k]) : matrix._values[k];
+	        // update pointer
+	        p = i + 1;
+	      }
+	      // zero values
+	      for (var y = p; y < rows; y++)
+	        (a[y] = (a[y] || []))[j] = 0;
+	    }
+	    return a;
+	  };
+	  
+	  /**
+	   * Get a string representation of the matrix, with optional formatting options.
+	   * @param {Object | Number | Function} [options]  Formatting options. See
+	   *                                                lib/util/number:format for a
+	   *                                                description of the available
+	   *                                                options.
+	   * @returns {String} str
+	   */
+	  CcsMatrix.prototype.format = function (options) {
+	    // rows and columns
+	    var rows = this._size[0];
+	    var columns = this._size[1];
+	    // rows & columns
+	    var str = 'CCS [' + string.format(rows, options) + ' x ' + string.format(columns, options) + '] density: ' + string.format(this._values.length / (rows * columns), options) + '\n';
+	    // loop columns
+	    for (var j = 0; j < columns; j++) {
+	      // k0 <= k < k1 where k0 = _ptr[j] && k1 = _ptr[j+1]
+	      var k0 = this._ptr[j];
+	      var k1 = this._ptr[j + 1];
+	      // loop k within [k0, k1[
+	      for (var k = k0; k < k1; k++) {
+	        // row index
+	        var i = this._index[k];
+	        // append value
+	        str += '\n    (' + string.format(i, options) + ', ' + string.format(j, options) + ') ==> ' + string.format(this._values[k], options);
+	      }
+	    }
+	    return str;
+	  };
+	  
+	  /**
+	   * Get a string representation of the matrix
+	   * @returns {String} str
+	   */
+	  CcsMatrix.prototype.toString = function () {
+	    return string.format(this.toArray());
+	  };
+	  
+	  /**
+	   * Get a JSON representation of the matrix
+	   * @returns {Object}
+	   */
+	  CcsMatrix.prototype.toJSON = function () {
+	    return {
+	      mathjs: 'CcsMatrix',
+	      values: this._values,
+	      index: this._index,
+	      ptr: this._ptr,
+	      size: this._size
+	    };
+	  };
+
+	  /**
+	   * Calculates the transpose of the matrix
+	   * @returns {Matrix}
+	   */
+	  CcsMatrix.prototype.transpose = function () {
+	    // rows and columns
+	    var rows = this._size[0];
+	    var columns = this._size[1];
+	    // check columns
+	    if (columns === 0) {
+	      // throw exception
+	      throw new RangeError('Cannot transpose a 2D matrix with no columns (size: ' + string.format(this._size) + ')');
+	    }
+	    // ccs transpose is a crs matrix with the same structure
+	    return new math.type.CrsMatrix({
+	      values: object.clone(this._values),
+	      index: object.clone(this._index),
+	      ptr: object.clone(this._ptr),
+	      size: [columns, rows]
+	    });
+	  };
+
+	  /**
+	   * Get the kth Matrix diagonal.
+	   *
+	   * @param {Number | BigNumber} [k=0]     The kth diagonal where the vector will retrieved.
+	   *
+	   * @returns {Array}                      The array vector with the diagonal values.
+	   */
+	  CcsMatrix.prototype.diagonal = function(k) {
+	    // validate k if any
+	    if (k) {
+	      // convert BigNumber to a number
+	      if (k instanceof BigNumber) 
+	        k = k.toNumber();
+	      // is must be an integer
+	      if (!isNumber(k) || !isInteger(k)) {
+	        throw new TypeError ('The parameter k must be an integer number');
+	      }
+	    }
+	    else {
+	      // default value
+	      k = 0;
+	    }
+
+	    var kSuper = k > 0 ? k : 0;
+	    var kSub = k < 0 ? -k : 0;
+	    
+	    // rows & columns
+	    var rows = this._size[0];
+	    var columns = this._size[1];
+	    
+	    // number diagonal values
+	    var n = Math.min(rows - kSub, columns -  kSuper);
+	    
+	    // diagonal
+	    var values = [];
+	    // loop columns
+	    for (var j = kSuper; j < columns && values.length < n; j++) {
+	      // k0 <= k < k1 where k0 = _ptr[j] && k1 = _ptr[j+1]
+	      var k0 = this._ptr[j];
+	      var k1 = this._ptr[j + 1];
+	      // column value flag
+	      var cv = false;
+	      // loop x within [k0, k1[
+	      for (var x = k0; x < k1; x++) {
+	        // row index
+	        var i = this._index[x];
+	        // check row
+	        if (i === j - kSuper + kSub) {
+	          // set flag
+	          cv = true;
+	          // value on this column
+	          values.push(object.clone(this._values[x]));
+	          // exit loop
+	          break;
+	        }
+	        else if (i > j - kSuper + kSub) {
+	          // exit loop, no value on the diagonal for column j
+	          break;
+	        }
+	      }
+	      // check this column has a value set
+	      if (!cv && values.length < n) {
+	        // zero on this column
+	        values.push(0);
+	      }
+	    }
+	    return values;
+	  };
+	  
+	  /**
+	   * Generate a matrix from a JSON object
+	   * @param {Object} json  An object structured like
+	   *                       `{"mathjs": "CcsMatrix", "values": [], "index": [], "ptr": [], "size": []}`,
+	   *                       where mathjs is optional
+	   * @returns {CcsMatrix}
+	   */
+	  CcsMatrix.fromJSON = function (json) {
+	    return new CcsMatrix(json);
+	  };
+
+	  /**
+	   * Create a diagonal matrix.
+	   *
+	   * @param {Array} size                   The matrix size.
+	   * @param {Number, Array} value          The values for the diagonal.
+	   * @param {Number | BigNumber} [k=0]     The kth diagonal where the vector will be filled in.
+	   *
+	   * @returns {CcsMatrix}
+	   */
+	  CcsMatrix.diagonal = function (size, value, k) {
+	    if (!isArray(size))
+	      throw new TypeError('Array expected, size parameter');
+	    if (size.length !== 2)
+	      throw new Error('Only two dimensions matrix are supported');
+	    
+	    // map size & validate
+	    size = size.map(function (s) {
+	      // check it is a big number
+	      if (s instanceof BigNumber) {
+	        // convert it
+	        s = s.toNumber();
+	      }
+	      // validate arguments
+	      if (!isNumber(s) || !isInteger(s) || s < 1) {
+	        throw new Error('Size values must be positive integers');
+	      } 
+	      return s;
+	    });
+	    
+	    // validate k if any
+	    if (k) {
+	      // convert BigNumber to a number
+	      if (k instanceof BigNumber) 
+	        k = k.toNumber();
+	      // is must be an integer
+	      if (!isNumber(k) || !isInteger(k)) {
+	        throw new TypeError ('The parameter k must be an integer number');
+	      }
+	    }
+	    else {
+	      // default value
+	      k = 0;
+	    }
+
+	    var kSuper = k > 0 ? k : 0;
+	    var kSub = k < 0 ? -k : 0;
+	    
+	    // rows and columns
+	    var rows = size[0];
+	    var columns = size[1];
+	    
+	    // number of non-zero items
+	    var n = Math.min(rows - kSub, columns -  kSuper);
+	    
+	    // value extraction function
+	    var _value;
+	      
+	    // check value
+	    if (isArray(value)) {
+	      // validate array
+	      if (value.length !== n) {
+	        // number of values in array must be n
+	        throw new Error('Invalid value array length');
+	      }
+	      // define function
+	      _value = function (i) {
+	        // return value @ i
+	        return value[i];
+	      };
+	    }
+	    else {
+	      // define function
+	      _value = function () {
+	        // return value
+	        return value;
+	      };
+	    }
+	    
+	    // create arrays
+	    var values = [];
+	    var index = [];
+	    var ptr = [];
+	    
+	    // loop items
+	    for (var j = 0; j < columns; j++) {
+	      // number of rows with value
+	      ptr.push(values.length);
+	      // diagonal index
+	      var i = j - kSuper;      
+	      // check we need to set diagonal value
+	      if (i >= 0 && i < n) {
+	        // get value @ i
+	        var v = _value(i);
+	        // check for zero
+	        if (!math.equal(v, 0)) {
+	          // column
+	          index.push(i + kSub);
+	          // add value
+	          values.push(v);
+	        }
+	      }
+	    }
+	    // last value should be number of values
+	    ptr.push(values.length);
+	    // create CcsMatrix
+	    return new CcsMatrix({
+	      values: values,
+	      index: index,
+	      ptr: ptr,
+	      size: [rows, columns]
+	    });
+	  };
+	  
+	  /**
+	   * Calculate the trace of a matrix: the sum of the elements on the main
+	   * diagonal of a square matrix.
+	   *
+	   * See also:
+	   *
+	   *    diagonal
+	   *
+	   * @returns {Number}       The matrix trace
+	   */
+	  CcsMatrix.prototype.trace = function () {
+	    // size
+	    var size = this._size;
+	    // check dimensions
+	    var rows = size[0];
+	    var columns = size[1];
+	    // matrix must be square
+	    if (rows === columns) {
+	      // calulate sum
+	      var sum = 0;
+	      // check we have data (avoid looping columns)
+	      if (this._values.length > 0) {
+	        // loop columns
+	        for (var j = 0; j < columns; j++) {
+	          // k0 <= k < k1 where k0 = _ptr[j] && k1 = _ptr[j+1]
+	          var k0 = this._ptr[j];
+	          var k1 = this._ptr[j + 1];
+	          // loop k within [k0, k1[
+	          for (var k = k0; k < k1; k++) {
+	            // row index
+	            var i = this._index[k];
+	            // check row
+	            if (i === j) {
+	              // accumulate value
+	              sum = math.add(sum, this._values[k]);
+	              // exit loop
+	              break;
+	            }
+	            if (i > j) {
+	              // exit loop, no value on the diagonal for column j
+	              break;
+	            }
+	          }
+	        }
+	      }
+	      // return trace
+	      return sum;
+	    }
+	    throw new RangeError('Matrix must be square (size: ' + string.format(size) + ')');        
+	  };
+	  
+	  /**
+	   * Multiply the matrix values times the argument.
+	   *
+	   * @param  {Number | BigNumber | Boolean | Complex | Unit | Array | Matrix | null} Value to multiply.
+	   *
+	   * @return {Number | BigNumber | Complex | Unit | Matrix}
+	   */
+	  CcsMatrix.prototype.multiply = function (value) {
+	    // check dimensions
+	    var rows = this._size[0];
+	    var columns = this._size[1];
+	    
+	    // check value is a matrix
+	    if (value instanceof Matrix) {
+	      // matrix size
+	      var z = value.size();
+	      // check value is a vector
+	      if (z.length === 1) {
+	        // mutiply matrix x vector array
+	        return _multiply(this, z[0], 1, function (i) {
+	          // value[i]
+	          return value.get([i]);
+	        });
+	      }
+	      // check two dimensions matrix
+	      if (z.length === 2) {        
+	        // mutiply matrix x matrix
+	        return _multiply(this, z[0], z[1], function (i, j) {
+	          // value[i, j]
+	          return value.get([i, j]);
+	        });
+	      }
+	      throw new Error('Can only multiply a 1 or 2 dimensional matrix ' +
+	                      '(value has ' + z.length + ' dimensions)');
+	    }
+
+	    // check value is an array
+	    if (isArray(value)) {
+	      // array size
+	      var s = array.size(value);
+	      // check value is a vector
+	      if (s.length === 1) {
+	        // mutiply matrix x vector array
+	        return _multiply(this, s[0], 1, function (i) {
+	          // value[i]
+	          return value[i];
+	        });
+	      }
+	      if (s.length === 2) {
+	        // mutiply matrix x array
+	        return _multiply(this, s[0], s[1], function (i, j) {
+	          // value[i, j]
+	          return value[i][j];
+	        });
+	      }
+	      throw new Error('Can only multiply a 1 or 2 dimensional matrix ' +
+	                      '(value has ' + s.length + ' dimensions)');
+	    }
+	    
+	    var callback = function (v) {
+	      return math.multiply(v, value);
+	    };
+
+	    // map non zero elements
+	    return _map(this, 0, rows - 1, 0, columns - 1, callback, false);
+	  };
+
+	  var _multiply = function (matrix, r, c, get) {
+	  
+	    // matrix dimensions
+	    var rows = matrix._size[0];
+	    var columns = matrix._size[1];
+	    
+	    // check dimensions match
+	    if (columns !== r) {
+	      // throw error
+	      throw new RangeError('Dimension mismatch in multiplication. ' +
+	                           'Columns of A must match length of B ' +
+	                           '(A is ' + rows + 'x' + columns +
+	                           ', B is ' + r + ', ' +
+	                           columns + ' != ' + r + ')');
+	    }
+	    
+	    // result arrays
+	    var values = [];
+	    var index = [];
+	    var ptr = [];
+
+	    // create array with rows entries
+	    var data = [];
+	    for (var x = 0; x < rows; x++)
+	      data[x] = 0;
+	    // loop value columns
+	    for (var z = 0; z < c; z++) {
+	      // update ptr
+	      ptr.push(values.length);
+	      // do not traverse rows in matrix, it is not efficient in CCS
+	      for (var j = 0; j < columns; j++) {          
+	        // k0 <= k < k1 where k0 = _ptr[j] && k1 = _ptr[j+1]
+	        var k0 = matrix._ptr[j];
+	        var k1 = matrix._ptr[j + 1];
+	        // loop k within [k0, k1[
+	        for (var k = k0; k < k1; k++) {
+	          // row
+	          var i = matrix._index[k];
+	          // multiply & aggregate
+	          data[i] = math.add(data[i], math.multiply(matrix._values[k], get(j, z)));
+	        }          
+	      }
+	      // finished processing column z, compress results
+	      for (var y = 0; y < rows; y++) {
+	        // check value is different than zero
+	        if (!math.equal(data[y], 0)) {
+	          // push value
+	          values.push(data[y]);
+	          index.push(y);          
+	        }
+	        // reset value
+	        data[y] = 0;
+	      }
+	    }
+	    // update ptr
+	    ptr.push(values.length);
+	    
+	    // check we need to squeeze the result into a scalar
+	    if (rows === 1 && c === 1)
+	      return values.length === 1 ? values[0] : 0;
+	    
+	    // return CCS matrix
+	    return new CcsMatrix({
+	      values: values,
+	      index: index,
+	      ptr: ptr,
+	      size: [rows, c]
+	    });
+	  };
+	  
+	  return CcsMatrix;
+	};
 
 /***/ },
 /* 16 */
@@ -4077,1769 +4996,4033 @@
 
 	'use strict';
 
-	var util = __webpack_require__(171),
-
-	    ArgumentsError = __webpack_require__(167),
+	var util = __webpack_require__(175);
+	var DimensionError = __webpack_require__(172);
+
+	var array = util.array;
+	var object = util.object;
+	var string = util.string;
+	var number = util.number;
+
+	var isArray = Array.isArray;
+	var isNumber = util.number.isNumber;
+	var isInteger = util.number.isInteger;
+
+	var validateIndex = array.validateIndex;
+
+	module.exports = function (math) {
+
+	  var Index = math.type.Index,
+	      BigNumber = math.type.BigNumber,
+	      Matrix = math.type.Matrix;
+
+	  function CrsMatrix(data) {
+	    if (!(this instanceof CrsMatrix))
+	      throw new SyntaxError('Constructor must be called with the new operator');
+
+	    if (data instanceof Matrix) {
+	      // check data is a CrsMatrix
+	      if (data.type === 'CrsMatrix') {
+	        // clone arrays
+	        this._values = object.clone(data._values);
+	        this._index = object.clone(data._index);
+	        this._ptr = object.clone(data._ptr);
+	        this._size = object.clone(data._size);
+	      }
+	      else {
+	        // build from matrix data
+	        _createFromArray(this, data.valueOf());
+	      }
+	    }
+	    else if (data && isArray(data.values) && isArray(data.index) && isArray(data.ptr) && isArray(data.size)) {
+	      // initialize fields
+	      this._values = data.values;
+	      this._index = data.index;
+	      this._ptr = data.ptr;
+	      this._size = data.size;
+	    }
+	    else if (isArray(data)) {
+	      // create from array
+	      _createFromArray(this, data);
+	    }
+	    else if (data) {
+	      // unsupported type
+	      throw new TypeError('Unsupported type of data (' + util.types.type(data) + ')');
+	    }
+	    else {
+	      // nothing provided
+	      this._values = [];
+	      this._index = [];
+	      this._ptr = [0];
+	      this._size = [0];
+	    }
+	  }
+
+	  var _createFromArray = function (matrix, data) {
+	    // initialize fields
+	    matrix._values = [];
+	    matrix._index = [];
+	    matrix._ptr = [];
+	    // discover rows & columns, do not use math.size() to avoid looping array twice
+	    var rows = data.length;
+	    var columns = 0;
+
+	    // loop rows
+	    for (var i = 0; i < rows; i++) {
+	      // store value index in ptr
+	      matrix._ptr.push(matrix._values.length);
+	      // current row
+	      var row = data[i];      
+	      // check row is an array
+	      if (isArray(row)) {
+	        // update columns if needed
+	        if (row.length > columns)
+	          columns = row.length;
+	        // loop columns
+	        for (var j = 0; j < row.length; j++) {
+	          // value at data[i][j]
+	          var v = row[j];
+	          // check value != 0
+	          if (!math.equal(v, 0)) {
+	            // store value
+	            matrix._values.push(v);
+	            // add column index
+	            matrix._index.push(j);
+	          }
+	        }
+	      }
+	      else {
+	        // update columns if needed (only on first row)
+	        if (i === 0 && columns < 1)
+	          columns = 1;
+	        // check value != 0 (row is a scalar)
+	        if (!math.equal(row, 0)) {
+	          // store value
+	          matrix._values.push(row);
+	          // index
+	          matrix._index.push(0);
+	        }
+	      }
+	    }
+	    // store number of values in ptr
+	    matrix._ptr.push(matrix._values.length);
+	    // size
+	    matrix._size = [rows, columns];
+	  };
+
+	  CrsMatrix.prototype = new math.type.Matrix();
+
+	  CrsMatrix.prototype.type = 'CrsMatrix';
+
+	  /**
+	   * Get the storage format used by the matrix.
+	   *
+	   * Usage:
+	   *     var format = matrix.storage()                   // retrieve storage format
+	   *
+	   * @return {string}           The storage format.
+	   */
+	  CrsMatrix.prototype.storage = function () {
+	    return 'crs';
+	  };
+
+	  /**
+	   * Get a subset of the matrix, or replace a subset of the matrix.
+	   *
+	   * Usage:
+	   *     var subset = matrix.subset(index)               // retrieve subset
+	   *     var value = matrix.subset(index, replacement)   // replace subset
+	   *
+	   * @param {Index} index
+	   * @param {Array | Maytrix | *} [replacement]
+	   * @param {*} [defaultValue=0]      Default value, filled in on new entries when
+	   *                                  the matrix is resized. If not provided,
+	   *                                  new matrix elements will be filled with zeros.
+	   */
+	  CrsMatrix.prototype.subset = function (index, replacement, defaultValue) {
+	    // check arguments
+	    switch (arguments.length) {
+	      case 1:
+	        return _getsubset(this, index);
+
+	        // intentional fall through
+	      case 2:
+	      case 3:
+	        return _setsubset(this, index, replacement, defaultValue);
+
+	      default:
+	        throw new SyntaxError('Wrong number of arguments');
+	    }
+	  };
+
+	  var _getsubset = function (matrix, index) {
+	    // check index
+	    if (!(index instanceof Index)) {
+	      throw new TypeError('Invalid index');
+	    }
+
+	    var isScalar = index.isScalar();
+	    if (isScalar) {
+	      // return a scalar
+	      return matrix.get(index.min());
+	    }
+	    // validate dimensions
+	    var size = index.size();
+	    if (size.length != matrix._size.length) {
+	      throw new DimensionError(size.length, matrix._size.length);
+	    }
+
+	    // validate if any of the ranges in the index is out of range
+	    var min = index.min();
+	    var max = index.max();
+	    for (var i = 0, ii = matrix._size.length; i < ii; i++) {
+	      validateIndex(min[i], matrix._size[i]);
+	      validateIndex(max[i], matrix._size[i]);
+	    }
+
+	    // map callback
+	    var callback = function (v) {
+	      // return value
+	      return v;
+	    };
+	    // get sub-matrix
+	    return _map(matrix, min[0], max[0], min[1], max[1], callback, false);
+	  };
+
+	  var _setsubset = function (matrix, index, submatrix, defaultValue) {
+	    // check index
+	    if (!(index instanceof Index)) {
+	      throw new TypeError('Invalid index');
+	    }
+
+	    // get index size and check whether the index contains a single value
+	    var iSize = index.size(),
+	        isScalar = index.isScalar();
+
+	    // calculate the size of the submatrix, and convert it into an Array if needed
+	    var sSize;
+	    if (submatrix instanceof Matrix) {
+	      // submatrix size
+	      sSize = submatrix.size();
+	      // use array representation
+	      submatrix = submatrix.toArray();
+	    }
+	    else {
+	      // get submatrix size (array, scalar)
+	      sSize = array.size(submatrix);
+	    }
+
+	    // check index is a scalar
+	    if (isScalar) {
+	      // verify submatrix is a scalar
+	      if (sSize.length !== 0) {
+	        throw new TypeError('Scalar expected');
+	      }
+	      // set value
+	      matrix.set(index.min(), submatrix, defaultValue);
+	    }
+	    else {
+	      // validate dimensions, index size must be one or two dimensions
+	      if (iSize.length !== 1 && iSize.length !== 2) {
+	        throw new DimensionError(iSize.length, matrix._size.length, '<');
+	      }
+
+	      // check submatrix and index have the same dimensions
+	      if (sSize.length < iSize.length) {
+	        // calculate number of missing outer dimensions
+	        var i = 0;
+	        var outer = 0;
+	        while (iSize[i] === 1 && sSize[i] === 1) {
+	          i++;
+	        }
+	        while (iSize[i] === 1) {
+	          outer++;
+	          i++;
+	        }
+	        // unsqueeze both outer and inner dimensions
+	        submatrix = array.unsqueeze(submatrix, iSize.length, outer, sSize);
+	      }
+
+	      // check whether the size of the submatrix matches the index size
+	      if (!object.deepEqual(iSize, sSize)) {
+	        throw new DimensionError(iSize, sSize, '>');
+	      }
+
+	      // offsets
+	      var x0 = index.min()[0];
+	      var y0 = index.min()[1];      
+
+	      // submatrix rows and columns
+	      var m = sSize[0];
+	      var n = sSize[1];
+
+	      // loop submatrix
+	      for (var x = 0; x < m; x++) {
+	        // loop columns
+	        for (var y = 0; y < n; y++) {
+	          // value at i, j
+	          var v = submatrix[x][y];
+	          // invoke set (zero value will remove entry from matrix)
+	          matrix.set([x + x0, y + y0], v, defaultValue);
+	        }
+	      }
+	    }
+	    return matrix;
+	  };
+
+	  /**
+	   * Get a single element from the matrix.
+	   * @param {Number[]} index   Zero-based index
+	   * @return {*} value
+	   */
+	  CrsMatrix.prototype.get = function (index) {
+	    if (!isArray(index))
+	      throw new TypeError('Array expected');
+	    if (index.length != this._size.length)
+	      throw new DimensionError(index.length, this._size.length);
+
+	    // row and column
+	    var i = index[0];
+	    var j = index[1];
+
+	    // check i, j are valid
+	    validateIndex(i, this._size[0]);
+	    validateIndex(j, this._size[1]);
+
+	    // find value index
+	    var k = _getValueIndex(j, this._ptr[i], this._ptr[i + 1], this._index);
+	    // check k is prior to next row k and it is in the correct row
+	    if (k < this._ptr[i + 1] && this._index[k] === j)
+	      return object.clone(this._values[k]);
+
+	    return 0;
+	  };
+
+	  /**
+	   * Replace a single element in the matrix.
+	   * @param {Number[]} index   Zero-based index
+	   * @param {*} value
+	   * @param {*} [defaultValue]        Default value, filled in on new entries when
+	   *                                  the matrix is resized. If not provided,
+	   *                                  new matrix elements will be set to zero.
+	   * @return {CrsMatrix} self
+	   */
+	  CrsMatrix.prototype.set = function (index, v, defaultValue) {
+	    if (!isArray(index))
+	      throw new TypeError('Array expected');
+	    if (index.length != this._size.length)
+	      throw new DimensionError(index.length, this._size.length);
+
+	    // row and column
+	    var i = index[0];
+	    var j = index[1];
+
+	    // rows & columns
+	    var rows = this._size[0];
+	    var columns = this._size[1];
+
+	    // check we need to resize matrix
+	    if (i > rows - 1 || j > columns - 1) {
+	      // resize matrix
+	      _resize(this, Math.max(i + 1, rows), Math.max(j + 1, columns), defaultValue);
+	      // update rows & columns
+	      rows = this._size[0];
+	      columns = this._size[1];
+	    }
+
+	    // check i, j are valid
+	    validateIndex(i, rows);
+	    validateIndex(j, columns);
+
+	    // find value index
+	    var k = _getValueIndex(j, this._ptr[i], this._ptr[i + 1], this._index);
+	    // check k is prior to next row k and it is in the correct column
+	    if (k < this._ptr[i + 1] && this._index[k] === j) {
+	      // check value != 0
+	      if (!math.equal(v, 0)) {
+	        // update value
+	        this._values[k] = v;
+	      }
+	      else {
+	        // remove value from matrix
+	        _remove(k, i, this._values, this._index, this._ptr);
+	      }
+	    }
+	    else {
+	      // insert value @ (i, j)
+	      _insert(k, i, j, v, this._values, this._index, this._ptr);
+	    }
+
+	    return this;
+	  };
+
+	  var _getValueIndex = function(j, left, right, index) {
+	    // check column is on the right side
+	    if (right - left === 0 || j > index[right - 1])
+	      return right;
+	    // loop until we find row index
+	    while (left < right) {
+	      // point in the middle (fast integer division)
+	      var p = ~~((left + right) / 2);
+	      // column @ p
+	      var c = index[p];
+	      // check we have to look on the left side, right side or we found the column
+	      if (j < c)
+	        right = p;
+	      else if (j > c)
+	        left = p + 1;
+	      else
+	        return p;
+	    }
+	    return left;
+	  };
+
+	  var _remove = function (k, i, values, index, ptr) {
+	    // remove value @ k
+	    values.splice(k, 1);
+	    index.splice(k, 1);
+	    // update pointers
+	    for (var x = i + 1; x < ptr.length; x++)
+	      ptr[x]--;
+	  };
+
+	  var _insert = function (k, i, j, v, values, index, ptr) {
+	    // insert value
+	    values.splice(k, 0, v);
+	    // update column for k
+	    index.splice(k, 0, j);
+	    // update row pointers
+	    for (var x = i + 1; x < ptr.length; x++)
+	      ptr[x]++;
+	  };
+
+	  /**
+	   * Resize the matrix to the given size. Returns a copy of the matrix when 
+	   * `copy=true`, otherwise return the matrix itself (resize in place).
+	   *
+	   * @param {Number[]} size           The new size the matrix should have.
+	   * @param {*} [defaultValue=0]      Default value, filled in on new entries.
+	   *                                  If not provided, the matrix elements will
+	   *                                  be filled with zeros.
+	   * @param {boolean} [copy]          Return a resized copy of the matrix
+	   *
+	   * @return {Matrix}                 The resized matrix
+	   */
+	  CrsMatrix.prototype.resize = function (size, defaultValue, copy) {    
+	    // validate arguments
+	    if (!isArray(size))
+	      throw new TypeError('Array expected');
+	    if (size.length !== 2)
+	      throw new Error('Only two dimensions matrix are supported');
+
+	    // check sizes
+	    size.forEach(function (value) {
+	      if (!number.isNumber(value) || !number.isInteger(value) || value < 0) {
+	        throw new TypeError('Invalid size, must contain positive integers ' +
+	                            '(size: ' + string.format(size) + ')');
+	      }
+	    });
+
+	    // matrix to resize
+	    var m = copy ? this.clone() : this;
+	    // resize matrix
+	    return _resize(m, size[0], size[1], defaultValue);
+	  };
+
+	  var _resize = function (matrix, rows, columns, defaultValue) {
+	    // value to insert at the time of growing matrix
+	    var value = defaultValue || 0;
+	    // should we insert the value?
+	    var ins = !math.equal(value, 0);
+
+	    // old columns and rows
+	    var r = matrix._size[0];
+	    var c = matrix._size[1];
+
+	    var i, j, k;
+
+	    // check we need to increase rows
+	    if (rows > r) {
+	      // loop new rows
+	      for (i = r; i < rows; i++) {
+	        // update matrix._ptr for current column
+	        matrix._ptr[i] = matrix._values.length;
+	        // check we need to insert matrix._values
+	        if (ins) {
+	          // loop columns
+	          for (j = 0; j < c; j++) {
+	            // add new matrix._values
+	            matrix._values.push(value);
+	            // update matrix._index
+	            matrix._index.push(j);
+	          }
+	        }        
+	      }
+	      // store number of matrix._values in matrix._ptr
+	      matrix._ptr[rows] = matrix._values.length;
+	    }
+	    else if (rows < r) {
+	      // truncate matrix._ptr
+	      matrix._ptr.splice(rows + 1, r - rows);
+	      // truncate matrix._values and matrix._index
+	      matrix._values.splice(matrix._ptr[rows], matrix._values.length);
+	      matrix._index.splice(matrix._ptr[rows], matrix._index.length);
+	    }
+	    // update rows
+	    r = rows;
+
+	    // check we need to increase columns
+	    if (columns > c) {
+	      // check we have to insert values
+	      if (ins) {
+	        // inserts
+	        var n = 0;
+	        // loop rows
+	        for (i = 0; i < r; i++) {
+	          // update matrix._ptr for current row
+	          matrix._ptr[i] = matrix._ptr[i] + n;
+	          // where to insert matrix._values
+	          k = matrix._ptr[i + 1] + n;
+	          // pointer
+	          var p = 0;
+	          // loop new columns, initialize pointer
+	          for (j = c; j < columns; j++, p++) {
+	            // add value
+	            matrix._values.splice(k + p, 0, value);
+	            // update matrix._index
+	            matrix._index.splice(k + p, 0, j);
+	            // increment inserts
+	            n++;
+	          }
+	        }
+	        // store number of matrix._values in matrix._ptr
+	        matrix._ptr[r] = matrix._values.length;
+	      }
+	    }
+	    else if (columns < c) {
+	      // deletes
+	      var d = 0;
+	      // loop rows
+	      for (i = 0; i < r; i++) {
+	        // update matrix._ptr for current row
+	        matrix._ptr[i] = matrix._ptr[i] - d;
+	        // where matrix._values start for next column
+	        var k0 = matrix._ptr[i];
+	        var k1 = matrix._ptr[i + 1] - d;
+	        // loop matrix._index
+	        for (k = k0; k < k1; k++) {
+	          // column
+	          j = matrix._index[k];
+	          // check we need to delete value and matrix._index
+	          if (j > columns - 1) {
+	            // remove value
+	            matrix._values.splice(k, 1);
+	            // remove item from matrix._index
+	            matrix._index.splice(k, 1);
+	            // increase deletes
+	            d++;
+	          }
+	        }
+	      }
+	      // update matrix._ptr for current column
+	      matrix._ptr[i] = matrix._values.length;
+	    }
+	    // update matrix._size
+	    matrix._size[0] = rows;
+	    matrix._size[1] = columns;
+	    // return matrix
+	    return matrix;
+	  };
+
+	  /**
+	   * Create a clone of the matrix
+	   * @return {CrsMatrix} clone
+	   */
+	  CrsMatrix.prototype.clone = function () {
+	    var m = new CrsMatrix({
+	      values: object.clone(this._values),
+	      index: object.clone(this._index),
+	      ptr: object.clone(this._ptr),
+	      size: object.clone(this._size)
+	    });
+	    return m;
+	  };
+
+	  /**
+	   * Retrieve the size of the matrix.
+	   * @returns {Number[]} size
+	   */
+	  CrsMatrix.prototype.size = function() {
+	    return object.clone(this._size);
+	  };
+
+	  /**
+	   * Create a new matrix with the results of the callback function executed on
+	   * each entry of the matrix.
+	   * @param {function} callback   The callback function is invoked with three
+	   *                              parameters: the value of the element, the index
+	   *                              of the element, and the Matrix being traversed.
+	   * @param {boolean} [skipZeros] Invoke callback function for non-zero values only.
+	   *
+	   * @return {Matrix} matrix
+	   */
+	  CrsMatrix.prototype.map = function (callback, skipZeros) {
+	    // matrix instance
+	    var me = this;
+	    // rows and columns
+	    var rows = this._size[0];
+	    var columns = this._size[1];
+	    // invoke callback
+	    var invoke = function (v, i, j) {
+	      // invoke callback
+	      return callback(v, [i, j], me);
+	    };
+	    // invoke _map
+	    return _map(this, 0, rows - 1, 0, columns - 1, invoke, skipZeros);
+	  };
+
+	  /**
+	   * Create a new matrix with the results of the callback function executed on the interval
+	   * [minRow..maxRow, minColumn..maxColumn].
+	   */
+	  var _map = function (matrix, minRow, maxRow, minColumn, maxColumn, callback, skipZeros) {
+	    // result arrays
+	    var values = [];
+	    var index = [];
+	    var ptr = [];
+	    // invoke callback
+	    var invoke = function (v, x, y) {
+	      // invoke callback
+	      v = callback(v, x, y);
+	      // check value != 0
+	      if (!math.equal(v, 0)) {
+	        // store value
+	        values.push(v);
+	        // index
+	        index.push(y);
+	      }
+	    };
+	    // loop rows
+	    for (var i = minRow; i <= maxRow; i++) {
+	      // store pointer to values index
+	      ptr.push(values.length);
+	      // k0 <= k < k1 where k0 = _ptr[j] && k1 = _ptr[j+1]
+	      var k0 = matrix._ptr[i];
+	      var k1 = matrix._ptr[i + 1];
+	      // column pointer
+	      var p = minColumn;
+	      // loop k within [k0, k1[
+	      for (var k = k0; k < k1; k++) {
+	        // column index
+	        var j = matrix._index[k];
+	        // check j is in range
+	        if (j >= minColumn && j <= maxColumn) {
+	          // zero values
+	          if (!skipZeros) {
+	            // write zeros from column p to j
+	            for (var x = p; x < j; x++)
+	              invoke(0, i - minRow, x - minColumn);
+	          }
+	          // value @ k
+	          invoke(matrix._values[k], i - minRow, j - minColumn);
+	        }
+	        // update pointer
+	        p = j + 1;
+	      }
+	      // zero values
+	      if (!skipZeros) {
+	        // write zeros from column p to maxColumn
+	        for (var y = p; y <= maxColumn; y++)
+	          invoke(0, i - minRow, y - minColumn);
+	      }
+	    }
+	    // store number of values in ptr
+	    ptr.push(values.length);
+	    // return ccs
+	    return new CrsMatrix({
+	      values: values,
+	      index: index,
+	      ptr: ptr,
+	      size: [maxRow - minRow + 1, maxColumn - minColumn + 1]
+	    });
+	  };
+
+	  /**
+	   * Execute a callback function on each entry of the matrix.
+	   * @param {function} callback   The callback function is invoked with three
+	   *                              parameters: the value of the element, the index
+	   *                              of the element, and the Matrix being traversed.
+	   * @param {boolean} [skipZeros] Invoke callback function for non-zero values only.
+	   */
+	  CrsMatrix.prototype.forEach = function (callback, skipZeros) {
+	    // matrix instance
+	    var me = this;
+	    // rows and columns
+	    var rows = this._size[0];
+	    var columns = this._size[1];
+	    // loop rows
+	    for (var i = 0; i < rows; i++) {
+	      // k0 <= k < k1 where k0 = _ptr[i] && k1 = _ptr[i+1]
+	      var k0 = this._ptr[i];
+	      var k1 = this._ptr[i + 1];
+	      // column pointer
+	      var p = 0;
+	      // loop k within [k0, k1[
+	      for (var k = k0; k < k1; k++) {
+	        // column index
+	        var j = this._index[k];
+	        // check we need to process zeros
+	        if (!skipZeros) {
+	          // zero values
+	          for (var x = p; x < j; x++)
+	            callback(0, [i, x], me);
+	        }
+	        // value @ k
+	        callback(this._values[k], [i, j], me);
+	        // update pointer
+	        p = j + 1;
+	      }
+	      // check we need to process zeros
+	      if (!skipZeros) {
+	        // zero values
+	        for (var y = p; y < columns; y++)
+	          callback(0, [i, y], me);
+	      }
+	    }
+	  };
+
+	  /**
+	   * Create an Array with a copy of the data of the CrsMatrix
+	   * @returns {Array} array
+	   */
+	  CrsMatrix.prototype.toArray = function () {
+	    return _toArray(this, true);
+	  };
+
+	  /**
+	   * Get the primitive value of the CrsMatrix: a two dimensions array
+	   * @returns {Array} array
+	   */
+	  CrsMatrix.prototype.valueOf = function () {
+	    return _toArray(this, false);
+	  };
+
+	  var _toArray = function (matrix, copy) {
+	    // result
+	    var a = [];
+	    // rows and columns
+	    var rows = matrix._size[0];
+	    var columns = matrix._size[1];
+	    // loop rows
+	    for (var i = 0; i < rows; i++) {
+	      // push row
+	      var r = a[i] = [];
+	      // k0 <= k < k1 where k0 = _ptr[i] && k1 = _ptr[i+1]
+	      var k0 = matrix._ptr[i];
+	      var k1 = matrix._ptr[i + 1];
+	      // column pointer
+	      var p = 0;
+	      // loop k is within [k0, k1[
+	      for (var k = k0; k < k1; k++) {
+	        // column index
+	        var j = matrix._index[k];
+	        // zero values
+	        for (var x = p; x < j; x++)
+	          r[x] = 0;
+	        // set value
+	        r[j] = copy ? object.clone(matrix._values[k]) : matrix._values[k];
+	        // update pointer
+	        p = j + 1;
+	      }
+	      // zero values
+	      for (var y = p; y < columns; y++)
+	        r[y] = 0;
+	    }
+	    return a;
+	  };
+
+	  /**
+	   * Get a string representation of the matrix, with optional formatting options.
+	   * @param {Object | Number | Function} [options]  Formatting options. See
+	   *                                                lib/util/number:format for a
+	   *                                                description of the available
+	   *                                                options.
+	   * @returns {String} str
+	   */
+	  CrsMatrix.prototype.format = function (options) {
+	    // rows and columns
+	    var rows = this._size[0];
+	    var columns = this._size[1];
+	    // rows & columns
+	    var str = 'CRS [' + string.format(rows, options) + ' x ' + string.format(columns, options) + '] density: ' + string.format(this._values.length / (rows * columns), options) + '\n';
+	    // loop rows
+	    for (var i = 0; i < rows; i++) {
+	      // k0 <= k < k1 where k0 = _ptr[i] && k1 = _ptr[i+1]
+	      var k0 = this._ptr[i];
+	      var k1 = this._ptr[i + 1];
+	      // loop k within [k0, k1[
+	      for (var k = k0; k < k1; k++) {
+	        // column index
+	        var j = this._index[k];
+	        // append value
+	        str += '\n    (' + string.format(i, options) + ', ' + string.format(j, options) + ') ==> ' + string.format(this._values[k], options);
+	      }
+	    }
+	    return str;
+	  };
+
+	  /**
+	   * Get a string representation of the matrix
+	   * @returns {String} str
+	   */
+	  CrsMatrix.prototype.toString = function () {
+	    return string.format(this.toArray());
+	  };
+
+	  /**
+	   * Get a JSON representation of the matrix
+	   * @returns {Object}
+	   */
+	  CrsMatrix.prototype.toJSON = function () {
+	    return {
+	      mathjs: 'CrsMatrix',
+	      values: this._values,
+	      index: this._index,
+	      ptr: this._ptr,
+	      size: this._size
+	    };
+	  };
+
+	  /**
+	   * Calculates the transpose of the matrix
+	   * @returns {Matrix}
+	   */
+	  CrsMatrix.prototype.transpose = function () {
+	    // rows and columns
+	    var rows = this._size[0];
+	    var columns = this._size[1];
+	    // check columns
+	    if (columns === 0) {
+	      // throw exception
+	      throw new RangeError('Cannot transpose a 2D matrix with no columns (size: ' + string.format(this._size) + ')');
+	    }
+	    // crs transpose is a ccs matrix with the same structure
+	    return new math.type.CcsMatrix({
+	      values: object.clone(this._values),
+	      index: object.clone(this._index),
+	      ptr: object.clone(this._ptr),
+	      size: [columns, rows]
+	    });
+	  };
+
+	  /**
+	   * Get the kth Matrix diagonal.
+	   *
+	   * @param {Number | BigNumber} [k=0]     The kth diagonal where the vector will retrieved.
+	   *
+	   * @returns {Array}                      The array vector with the diagonal values.
+	   */
+	  CrsMatrix.prototype.diagonal = function(k) {
+	    // validate k if any
+	    if (k) {
+	      // convert BigNumber to a number
+	      if (k instanceof BigNumber) 
+	        k = k.toNumber();
+	      // is must be an integer
+	      if (!isNumber(k) || !isInteger(k)) {
+	        throw new TypeError ('The parameter k must be an integer number');
+	      }
+	    }
+	    else {
+	      // default value
+	      k = 0;
+	    }
+
+	    var kSuper = k > 0 ? k : 0;
+	    var kSub = k < 0 ? -k : 0;
+
+	    // rows & columns
+	    var rows = this._size[0];
+	    var columns = this._size[1];
+
+	    // number diagonal values
+	    var n = Math.min(rows - kSub, columns -  kSuper);
+
+	    // diagonal
+	    var values = [];
+	    // loop rows
+	    for (var i = kSub; i < rows && values.length < n; i++) {
+	      // k0 <= k < k1 where k0 = _ptr[i] && k1 = _ptr[i+1]
+	      var k0 = this._ptr[i];
+	      var k1 = this._ptr[i + 1];
+	      // row value flag
+	      var rv = false;
+	      // loop x within [k0, k1[
+	      for (var x = k0; x < k1; x++) {
+	        // column index
+	        var j = this._index[x];
+	        // check column
+	        if (j === i + kSuper - kSub) {
+	          // set flag
+	          rv = true;
+	          // value on this column
+	          values.push(object.clone(this._values[x]));
+	          // exit loop
+	          break;
+	        }
+	        else if (j > i + kSuper - kSub) {
+	          // exit loop, no value on the diagonal for row i
+	          break;
+	        }
+	      }
+	      // check this row has a value set
+	      if (!rv && values.length < n) {
+	        // zero on this column
+	        values.push(0);
+	      }
+	    }
+	    return values;
+	  };
+
+	  /**
+	   * Generate a matrix from a JSON object
+	   * @param {Object} json  An object structured like
+	   *                       `{"mathjs": "CrsMatrix", "values": [], "index": [], "ptr": [], "size": []}`,
+	   *                       where mathjs is optional
+	   * @returns {CrsMatrix}
+	   */
+	  CrsMatrix.fromJSON = function (json) {
+	    return new CrsMatrix(json);
+	  };
+
+	  /**
+	   * Create a diagonal matrix.
+	   *
+	   * @param {Array} size                   The matrix size.
+	   * @param {Number, Array} value          The values for the diagonal.
+	   * @param {Number | BigNumber} [k=0]     The kth diagonal where the vector will be filled in.
+	   *
+	   * @returns {CrsMatrix}
+	   */
+	  CrsMatrix.diagonal = function (size, value, k) {
+	    if (!isArray(size))
+	      throw new TypeError('Array expected, size parameter');
+	    if (size.length !== 2)
+	      throw new Error('Only two dimensions matrix are supported');
+
+	    // map size & validate
+	    size = size.map(function (s) {
+	      // check it is a big number
+	      if (s instanceof BigNumber) {
+	        // convert it
+	        s = s.toNumber();
+	      }
+	      // validate arguments
+	      if (!isNumber(s) || !isInteger(s) || s < 1) {
+	        throw new Error('Size values must be positive integers');
+	      } 
+	      return s;
+	    });
+
+	    // validate k if any
+	    if (k) {
+	      // convert BigNumber to a number
+	      if (k instanceof BigNumber) 
+	        k = k.toNumber();
+	      // is must be an integer
+	      if (!isNumber(k) || !isInteger(k)) {
+	        throw new TypeError ('The parameter k must be an integer number');
+	      }
+	    }
+	    else {
+	      // default value
+	      k = 0;
+	    }
+
+	    var kSuper = k > 0 ? k : 0;
+	    var kSub = k < 0 ? -k : 0;
+
+	    // rows and columns
+	    var rows = size[0];
+	    var columns = size[1];
+
+	    // number of non-zero items
+	    var n = Math.min(rows - kSub, columns -  kSuper);
+
+	    // value extraction function
+	    var _value;
+
+	    // check value
+	    if (isArray(value)) {
+	      // validate array
+	      if (value.length !== n) {
+	        // number of values in array must be n
+	        throw new Error('Invalid value array length');
+	      }
+	      // define function
+	      _value = function (i) {
+	        // return value @ i
+	        return value[i];
+	      };
+	    }
+	    else {
+	      // define function
+	      _value = function () {
+	        // return value
+	        return value;
+	      };
+	    }
+
+	    // create arrays
+	    var values = [];
+	    var index = [];
+	    var ptr = [];
+
+	    // loop items
+	    for (var i = 0; i < rows; i++) {
+	      // number of rows with value
+	      ptr.push(values.length);
+	      // diagonal index
+	      var j = i - kSub;
+	      // check we need to set diagonal value
+	      if (j >= 0 && j < n) {
+	        // get value @ j
+	        var v = _value(j);
+	        // check for zero
+	        if (!math.equal(v, 0)) {
+	          // column
+	          index.push(j + kSuper);
+	          // add value
+	          values.push(v);
+	        }
+	      }
+	    }
+	    // last value should be number of values
+	    ptr.push(values.length);
+	    // create CrsMatrix
+	    return new CrsMatrix({
+	      values: values,
+	      index: index,
+	      ptr: ptr,
+	      size: [rows, columns]
+	    });
+	  };
+
+	  /**
+	   * Calculate the trace of a matrix: the sum of the elements on the main
+	   * diagonal of a square matrix.
+	   *
+	   * See also:
+	   *
+	   *    diagonal
+	   *
+	   * @returns {Number}       The matrix trace
+	   */
+	  CrsMatrix.prototype.trace = function () {
+	    // size
+	    var size = this._size;
+	    // check dimensions
+	    var rows = size[0];
+	    var columns = size[1];
+	    // matrix must be square
+	    if (rows === columns) {
+	      // calulate sum
+	      var sum = 0;
+	      // check we have data (avoid looping rows)
+	      if (this._values.length > 0) {
+	        // loop rows
+	        for (var i = 0; i < rows; i++) {
+	          // k0 <= k < k1 where k0 = _ptr[i] && k1 = _ptr[i+1]
+	          var k0 = this._ptr[i];
+	          var k1 = this._ptr[i + 1];
+	          // loop k within [k0, k1[
+	          for (var k = k0; k < k1; k++) {
+	            // column index
+	            var j = this._index[k];
+	            // check row
+	            if (i === j) {
+	              // accumulate value
+	              sum = math.add(sum, this._values[k]);
+	              // exit loop
+	              break;
+	            }
+	            if (j > i) {
+	              // exit loop, no value on the diagonal for column j
+	              break;
+	            }
+	          }
+	        }
+	      }
+	      // return trace
+	      return sum;
+	    }
+	    throw new RangeError('Matrix must be square (size: ' + string.format(size) + ')');        
+	  };
+	  
+	  /**
+	   * Multiply the matrix values times the argument.
+	   *
+	   * @param  {Number | BigNumber | Boolean | Complex | Unit | Array | Matrix | null} Value to multiply.
+	   *
+	   * @return {Number | BigNumber | Complex | Unit | Matrix}
+	   */
+	  CrsMatrix.prototype.multiply = function (value) {
+	    // check dimensions
+	    var rows = this._size[0];
+	    var columns = this._size[1];
+
+	    // check value is a matrix
+	    if (value instanceof Matrix) {
+	      // matrix size
+	      var z = value.size();
+	      // check value is a vector
+	      if (z.length === 1) {
+	        // mutiply matrix x vector array
+	        return _multiply(this, z[0], 1, function (i) {
+	          // value[i]
+	          return value.get([i]);
+	        });
+	      }
+	      // check two dimensions matrix
+	      if (z.length === 2) {        
+	        // mutiply matrix x matrix
+	        return _multiply(this, z[0], z[1], function (i, j) {
+	          // value[i, j]
+	          return value.get([i, j]);
+	        });
+	      }
+	      throw new Error('Can only multiply a 1 or 2 dimensional matrix ' +
+	                      '(value has ' + z.length + ' dimensions)');
+	    }
+
+	    // check value is an array
+	    if (isArray(value)) {
+	      // array size
+	      var s = array.size(value);
+	      // check value is a vector
+	      if (s.length === 1) {
+	        // mutiply matrix x vector array
+	        return _multiply(this, s[0], 1, function (i) {
+	          // value[i]
+	          return value[i];
+	        });
+	      }
+	      if (s.length === 2) {
+	        // mutiply matrix x array
+	        return _multiply(this, s[0], s[1], function (i, j) {
+	          // value[i, j]
+	          return value[i][j];
+	        });
+	      }
+	      throw new Error('Can only multiply a 1 or 2 dimensional matrix ' +
+	                      '(value has ' + s.length + ' dimensions)');
+	    }
+
+	    var callback = function (v) {
+	      return math.multiply(v, value);
+	    };
+
+	    // map non zero elements
+	    return _map(this, 0, rows - 1, 0, columns - 1, callback, false);
+	  };
+	  
+	  var _multiply = function (matrix, r, c, get) {
+
+	    // matrix dimensions
+	    var rows = matrix._size[0];
+	    var columns = matrix._size[1];
+
+	    // check dimensions match
+	    if (columns !== r) {
+	      // throw error
+	      throw new RangeError('Dimension mismatch in multiplication. ' +
+	                           'Columns of A must match length of B ' +
+	                           '(A is ' + rows + 'x' + columns +
+	                           ', B is ' + r + ', ' +
+	                           columns + ' != ' + r + ')');
+	    }
+
+	    // result arrays
+	    var values = [];
+	    var index = [];
+	    var ptr = [];
+
+	    // loop rows
+	    for (var i = 0; i < rows; i++) {
+	      // update ptr
+	      ptr.push(values.length);
+	      // k0 <= k < k1 where k0 = _ptr[i] && k1 = _ptr[i+1]
+	      var k0 = matrix._ptr[i];
+	      var k1 = matrix._ptr[i + 1];
+	      // loop value columns
+	      for (var z = 0; z < c; z++) {
+	        // value @ (i, x)
+	        var value = 0;
+	        // loop k within [k0, k1[
+	        for (var k = k0; k < k1; k++) {
+	          // column
+	          var j = matrix._index[k];
+	          // multiply & aggregate
+	          value = math.add(value, math.multiply(matrix._values[k], get(j, z)));
+	        }
+	        // check value is different than zero
+	        if (!math.equal(value, 0)) {
+	          // push value & column
+	          values.push(value);
+	          index.push(z);
+	        }        
+	      }
+	    }
+	    // update ptr
+	    ptr.push(values.length);
+
+	    // check we need to squeeze the result into a scalar
+	    if (rows === 1 && c === 1)
+	      return values.length === 1 ? values[0] : 0;
+
+	    // return CRS matrix
+	    return new CrsMatrix({
+	      values: values,
+	      index: index,
+	      ptr: ptr,
+	      size: [rows, c]
+	    });
+	  };
+
+	  return CrsMatrix;
+	};
+
+/***/ },
+/* 17 */
+/***/ function(module, exports, __webpack_require__) {
+
+	'use strict';
+
+	var util = __webpack_require__(175);
+	var DimensionError = __webpack_require__(172);
+
+	var string = util.string;
+	var array = util.array;
+	var object = util.object;
+
+	var isArray = Array.isArray;
+	var isNumber = util.number.isNumber;
+	var isInteger = util.number.isInteger;
+
+	var validateIndex = array.validateIndex;
+
+	module.exports = function (math) {
+
+	  var Index = math.type.Index,
+	      BigNumber = math.type.BigNumber,
+	      Matrix = math.type.Matrix;
+	  
+	  function DenseMatrix(data) {
+	    if (!(this instanceof DenseMatrix))
+	      throw new SyntaxError('Constructor must be called with the new operator');
+
+	    if (data instanceof Matrix) {
+	      // check data is a DenseMatrix
+	      if (data.type === 'DenseMatrix') {
+	        // clone data & size
+	        this._data = object.clone(data._data);
+	        this._size = object.clone(data._size);
+	      }
+	      else {
+	        // build data from existing matrix
+	        this._data = data.toArray();
+	        this._size = data.size();
+	      }
+	    }
+	    else if (data && isArray(data.data) && isArray(data.size)) {
+	      // initialize fields from JSON representation
+	      this._data = data.data;
+	      this._size = data.size;
+	    }
+	    else if (isArray(data)) {
+	      // replace nested Matrices with Arrays
+	      this._data = preprocess(data);
+	      // verify the size of the array, TODO: compute size while processing array
+	      this._size = array.size(this._data);
+	    }
+	    else if (data) {
+	      // unsupported type
+	      throw new TypeError('Unsupported type of data (' + util.types.type(data) + ')');
+	    }
+	    else {
+	      // nothing provided
+	      this._data = [];
+	      this._size = [0];
+	    }
+	  }
+	  
+	  DenseMatrix.prototype = new math.type.Matrix();
+
+	  DenseMatrix.prototype.type = 'DenseMatrix';
+
+	  /**
+	   * Get the storage format used by the matrix.
+	   *
+	   * Usage:
+	   *     var format = matrix.storage()                   // retrieve storage format
+	   *
+	   * @return {string}           The storage format.
+	   */
+	  DenseMatrix.prototype.storage = function () {
+	    return 'dense';
+	  };
+	  
+	  /**
+	   * Get a subset of the matrix, or replace a subset of the matrix.
+	   *
+	   * Usage:
+	   *     var subset = matrix.subset(index)               // retrieve subset
+	   *     var value = matrix.subset(index, replacement)   // replace subset
+	   *
+	   * @param {Index} index
+	   * @param {Array | DenseMatrix | *} [replacement]
+	   * @param {*} [defaultValue=0]      Default value, filled in on new entries when
+	   *                                  the matrix is resized. If not provided,
+	   *                                  new matrix elements will be filled with zeros.
+	   */
+	  DenseMatrix.prototype.subset = function (index, replacement, defaultValue) {
+	    switch (arguments.length) {
+	      case 1:
+	        return _get(this, index);
+
+	        // intentional fall through
+	      case 2:
+	      case 3:
+	        return _set(this, index, replacement, defaultValue);
+
+	      default:
+	        throw new SyntaxError('Wrong number of arguments');
+	    }
+	  };
+	  
+	  /**
+	   * Get a single element from the matrix.
+	   * @param {Number[]} index   Zero-based index
+	   * @return {*} value
+	   */
+	  DenseMatrix.prototype.get = function (index) {
+	    if (!isArray(index))
+	      throw new TypeError('Array expected');
+	    if (index.length != this._size.length)
+	      throw new DimensionError(index.length, this._size.length);
+
+	    // check index
+	    for (var x = 0; x < index.length; x++)
+	      validateIndex(index[x], this._size[x]);
+
+	    var data = this._data;
+	    for (var i = 0, ii = index.length; i < ii; i++) {
+	      var index_i = index[i];
+	      validateIndex(index_i, data.length);
+	      data = data[index_i];
+	    }
+
+	    return object.clone(data);
+	  };
+	  
+	  /**
+	   * Replace a single element in the matrix.
+	   * @param {Number[]} index   Zero-based index
+	   * @param {*} value
+	   * @param {*} [defaultValue]        Default value, filled in on new entries when
+	   *                                  the matrix is resized. If not provided,
+	   *                                  new matrix elements will be left undefined.
+	   * @return {DenseMatrix} self
+	   */
+	  DenseMatrix.prototype.set = function (index, value, defaultValue) {
+	    if (!isArray(index))
+	      throw new TypeError('Array expected');
+	    if (index.length < this._size.length)
+	      throw new DimensionError(index.length, this._size.length, '<');
+
+	    var i, ii, index_i;
+
+	    // enlarge matrix when needed
+	    var size = index.map(function (i) {
+	      return i + 1;
+	    });
+	    _fit(this, size, defaultValue);
+
+	    // traverse over the dimensions
+	    var data = this._data;
+	    for (i = 0, ii = index.length - 1; i < ii; i++) {
+	      index_i = index[i];
+	      validateIndex(index_i, data.length);
+	      data = data[index_i];
+	    }
+
+	    // set new value
+	    index_i = index[index.length - 1];
+	    validateIndex(index_i, data.length);
+	    data[index_i] = value;
+
+	    return this;
+	  };
+	  
+	  /**
+	   * Get a submatrix of this matrix
+	   * @param {DenseMatrix} matrix
+	   * @param {Index} index   Zero-based index
+	   * @private
+	   */
+	  function _get (matrix, index) {
+	    if (!(index instanceof Index)) {
+	      throw new TypeError('Invalid index');
+	    }
+
+	    var isScalar = index.isScalar();
+	    if (isScalar) {
+	      // return a scalar
+	      return matrix.get(index.min());
+	    }
+	    else {
+	      // validate dimensions
+	      var size = index.size();
+	      if (size.length != matrix._size.length) {
+	        throw new DimensionError(size.length, matrix._size.length);
+	      }
+
+	      // validate if any of the ranges in the index is out of range
+	      var min = index.min();
+	      var max = index.max();
+	      for (var i = 0, ii = matrix._size.length; i < ii; i++) {
+	        validateIndex(min[i], matrix._size[i]);
+	        validateIndex(max[i], matrix._size[i]);
+	      }
+
+	      // retrieve submatrix
+	      // TODO: more efficient when creating an empty matrix and setting _data and _size manually
+	      return new DenseMatrix(_getSubmatrix(matrix._data, index, size.length, 0));
+	    }
+	  }
+	  
+	  /**
+	   * Recursively get a submatrix of a multi dimensional matrix.
+	   * Index is not checked for correct number or length of dimensions.
+	   * @param {Array} data
+	   * @param {Index} index
+	   * @param {number} dims   Total number of dimensions
+	   * @param {number} dim    Current dimension
+	   * @return {Array} submatrix
+	   * @private
+	   */
+	  function _getSubmatrix (data, index, dims, dim) {
+	    var last = (dim == dims - 1);
+	    var range = index.range(dim);
+
+	    if (last) {
+	      return range.map(function (i) {
+	        return data[i];
+	      });
+	    }
+	    else {
+	      return range.map(function (i) {
+	        var child = data[i];
+	        return _getSubmatrix(child, index, dims, dim + 1);
+	      });
+	    }
+	  }
+	  
+	  /**
+	   * Replace a submatrix in this matrix
+	   * Indexes are zero-based.
+	   * @param {DenseMatrix} matrix
+	   * @param {Index} index
+	   * @param {DenseMatrix | Array | *} submatrix
+	   * @param {*} defaultValue          Default value, filled in on new entries when
+	   *                                  the matrix is resized.
+	   * @return {DenseMatrix} matrix
+	   * @private
+	   */
+	  function _set (matrix, index, submatrix, defaultValue) {
+	    if (!(index instanceof Index)) {
+	      throw new TypeError('Invalid index');
+	    }
+
+	    // get index size and check whether the index contains a single value
+	    var iSize = index.size(),
+	        isScalar = index.isScalar();
+
+	    // calculate the size of the submatrix, and convert it into an Array if needed
+	    var sSize;
+	    if (submatrix instanceof math.type.Matrix) {
+	      sSize = submatrix.size();
+	      submatrix = submatrix.valueOf();
+	    }
+	    else {
+	      sSize = array.size(submatrix);
+	    }
+
+	    if (isScalar) {
+	      // set a scalar
+
+	      // check whether submatrix is a scalar
+	      if (sSize.length !== 0) {
+	        throw new TypeError('Scalar expected');
+	      }
+
+	      matrix.set(index.min(), submatrix, defaultValue);
+	    }
+	    else {
+	      // set a submatrix
+
+	      // validate dimensions
+	      if (iSize.length < matrix._size.length) {
+	        throw new DimensionError(iSize.length, matrix._size.length, '<');
+	      }
+
+	      if (sSize.length < iSize.length) {
+	        // calculate number of missing outer dimensions
+	        var i = 0;
+	        var outer = 0;
+	        while (iSize[i] === 1 && sSize[i] === 1) {
+	          i++;
+	        }
+	        while (iSize[i] === 1) {
+	          outer++;
+	          i++;
+	        }
+
+	        // unsqueeze both outer and inner dimensions
+	        submatrix = array.unsqueeze(submatrix, iSize.length, outer, sSize);
+	      }
+
+	      // check whether the size of the submatrix matches the index size
+	      if (!object.deepEqual(iSize, sSize)) {
+	        throw new DimensionError(iSize, sSize, '>');
+	      }
+
+	      // enlarge matrix when needed
+	      var size = index.max().map(function (i) {
+	        return i + 1;
+	      });
+	      _fit(matrix, size, defaultValue);
+
+	      // insert the sub matrix
+	      var dims = iSize.length,
+	          dim = 0;
+	      _setSubmatrix (matrix._data, index, submatrix, dims, dim);
+	    }
+
+	    return matrix;
+	  }
+	  
+	  /**
+	   * Replace a submatrix of a multi dimensional matrix.
+	   * @param {Array} data
+	   * @param {Index} index
+	   * @param {Array} submatrix
+	   * @param {number} dims   Total number of dimensions
+	   * @param {number} dim
+	   * @private
+	   */
+	  function _setSubmatrix (data, index, submatrix, dims, dim) {
+	    var last = (dim == dims - 1),
+	        range = index.range(dim);
+
+	    if (last) {
+	      range.forEach(function (dataIndex, subIndex) {
+	        validateIndex(dataIndex);
+	        data[dataIndex] = submatrix[subIndex];
+	      });
+	    }
+	    else {
+	      range.forEach(function (dataIndex, subIndex) {
+	        validateIndex(dataIndex);
+	        _setSubmatrix(data[dataIndex], index, submatrix[subIndex], dims, dim + 1);
+	      });
+	    }
+	  }
+	  
+	  /**
+	   * Resize the matrix to the given size. Returns a copy of the matrix when
+	   * `copy=true`, otherwise return the matrix itself (resize in place).
+	   *
+	   * @param {Number[]} size           The new size the matrix should have.
+	   * @param {*} [defaultValue=0]      Default value, filled in on new entries.
+	   *                                  If not provided, the matrix elements will
+	   *                                  be filled with zeros.
+	   * @param {boolean} [copy]          Return a resized copy of the matrix
+	   *
+	   * @return {Matrix}                 The resized matrix
+	   */
+	  DenseMatrix.prototype.resize = function (size, defaultValue, copy) {
+	    // validate arguments
+	    if (!isArray(size))
+	      throw new TypeError('Array expected');
+
+	    // matrix to resize
+	    var m = copy ? this.clone() : this;
+	    // resize matrix
+	    return _resize(m, size, defaultValue);
+	  };
+	  
+	  var _resize = function (matrix, size, defaultValue) {
+	    // check size
+	    if (size.length === 0) {
+	      // first value in matrix
+	      var v = matrix._data;
+	      // go deep
+	      while (isArray(v)) {
+	        v = v[0];
+	      }
+	      return object.clone(v);
+	    }
+	    // resize matrix
+	    matrix._size = object.clone(size);
+	    matrix._data = array.resize(matrix._data, matrix._size, defaultValue);
+	    // return matrix
+	    return matrix;
+	  };
+	  
+	  /**
+	   * Enlarge the matrix when it is smaller than given size.
+	   * If the matrix is larger or equal sized, nothing is done.
+	   * @param {DenseMatrix} matrix           The matrix to be resized
+	   * @param {Number[]} size
+	   * @param {*} defaultValue          Default value, filled in on new entries.
+	   * @private
+	   */
+	  function _fit(matrix, size, defaultValue) {
+	    var newSize = object.clone(matrix._size),
+	        changed = false;
+
+	    // add dimensions when needed
+	    while (newSize.length < size.length) {
+	      newSize.push(0);
+	      changed = true;
+	    }
+
+	    // enlarge size when needed
+	    for (var i = 0, ii = size.length; i < ii; i++) {
+	      if (size[i] > newSize[i]) {
+	        newSize[i] = size[i];
+	        changed = true;
+	      }
+	    }
+
+	    if (changed) {
+	      // resize only when size is changed
+	      _resize(matrix, newSize, defaultValue);
+	    }
+	  }
+	  
+	  /**
+	   * Create a clone of the matrix
+	   * @return {DenseMatrix} clone
+	   */
+	  DenseMatrix.prototype.clone = function () {
+	    var m = new DenseMatrix({
+	      data: object.clone(this._data),
+	      size: object.clone(this._size)
+	    });
+	    return m;
+	  };
+	  
+	  /**
+	   * Retrieve the size of the matrix.
+	   * @returns {Number[]} size
+	   */
+	  DenseMatrix.prototype.size = function() {
+	    return this._size;
+	  };
+	  
+	  /**
+	   * Create a new matrix with the results of the callback function executed on
+	   * each entry of the matrix.
+	   * @param {function} callback   The callback function is invoked with three
+	   *                              parameters: the value of the element, the index
+	   *                              of the element, and the Matrix being traversed.
+	   *
+	   * @return {DenseMatrix} matrix
+	   */
+	  DenseMatrix.prototype.map = function (callback) {
+	    // matrix instance
+	    var me = this;
+	    var recurse = function (value, index) {
+	      if (isArray(value)) {
+	        return value.map(function (child, i) {
+	          return recurse(child, index.concat(i));
+	        });
+	      }
+	      else {
+	        return callback(value, index, me);
+	      }
+	    };
+	    // return dense format
+	    return new DenseMatrix({
+	      data: recurse(this._data, []),
+	      size: object.clone(this._size)
+	    });
+	  };
+	  
+	  /**
+	   * Execute a callback function on each entry of the matrix.
+	   * @param {function} callback   The callback function is invoked with three
+	   *                              parameters: the value of the element, the index
+	   *                              of the element, and the Matrix being traversed.
+	   */
+	  DenseMatrix.prototype.forEach = function (callback) {
+	    // matrix instance
+	    var me = this;
+	    var recurse = function (value, index) {
+	      if (isArray(value)) {
+	        value.forEach(function (child, i) {
+	          recurse(child, index.concat(i));
+	        });
+	      }
+	      else {
+	        callback(value, index, me);
+	      }
+	    };
+	    recurse(this._data, []);
+	  };
+	  
+	  /**
+	   * Create an Array with a copy of the data of the DenseMatrix
+	   * @returns {Array} array
+	   */
+	  DenseMatrix.prototype.toArray = function () {
+	    return object.clone(this._data);
+	  };
+	  
+	  /**
+	   * Get the primitive value of the DenseMatrix: a multidimensional array
+	   * @returns {Array} array
+	   */
+	  DenseMatrix.prototype.valueOf = function () {
+	    return this._data;
+	  };
+	  
+	  /**
+	   * Get a string representation of the matrix, with optional formatting options.
+	   * @param {Object | Number | Function} [options]  Formatting options. See
+	   *                                                lib/util/number:format for a
+	   *                                                description of the available
+	   *                                                options.
+	   * @returns {String} str
+	   */
+	  DenseMatrix.prototype.format = function (options) {
+	    return string.format(this._data, options);
+	  };
+	  
+	  /**
+	   * Get a string representation of the matrix
+	   * @returns {String} str
+	   */
+	  DenseMatrix.prototype.toString = function () {
+	    return string.format(this._data);
+	  };
+	  
+	  /**
+	   * Get a JSON representation of the matrix
+	   * @returns {Object}
+	   */
+	  DenseMatrix.prototype.toJSON = function () {
+	    return {
+	      mathjs: 'DenseMatrix',
+	      data: this._data,
+	      size: this._size
+	    };
+	  };
+	  
+	  /**
+	   * Calculates the transpose of the matrix
+	   * @returns {Matrix}
+	   */
+	  DenseMatrix.prototype.transpose = function () {
+	    // check dimensions
+	    switch (this._size.length) {
+	        case 1:
+	          // vector
+	          return this.clone();
+	        case 2:
+	          // rows and columns
+	          var rows = this._size[0];
+	          var columns = this._size[1];
+	          // check columns
+	          if (columns === 0) {
+	            // throw exception
+	            throw new RangeError('Cannot transpose a 2D matrix with no columns (size: ' + string.format(this._size) + ')');
+	          }
+	          // transposed matrix data
+	          var transposed = [];
+	          var transposedRow;
+	          // loop columns
+	          for (var j = 0; j < columns; j++) {
+	            // initialize row
+	            transposedRow = transposed[j] = [];
+	            // loop rows
+	            for (var i = 0; i < rows; i++) {