s operator

examples/dna_mass.txt

@@ -0,0 +1,3 @@
+303.7 10 * 79 + (Da)
+10 (ul) 100 (uM)
+3 (ng) sp

examples/mpg.txt

@@ -0,0 +1,6 @@
+@base(usd) # define usd as a base unit
+0 (mi) 1 (gal) / 1 (mi) 1 (gal) / @derived(mpg) # define mpg as 1 mile per gallon
+3.4 (usd) 1 (gal) / # gas price
+0.36 (usd) 1 (kWh) / # electricity price
+36 (kWh) 100 (mi) / # electricity consumption
+3 (mpg) sp # solve for mpg

src/interpreter/mod.rs

@@ -47,6 +47,8 @@ pub enum InterpreterError {
     UndefinedVariable(String),
     #[error("Incompatible units: {0}")]
     IncompatibleUnits(UnitCombo),
+    #[error("No solution: {0}")]
+    NoSolution(String),
     #[error("Already defined: {0}")]
     AlreadyDefined(String),
 }
@@ -83,6 +85,7 @@ impl<'a> Interpreter<'a> {
                 Token::Operator('c') => self.op_c()?,
                 Token::Operator('d') => self.op_d()?,
                 Token::Operator('r') => self.op_r()?,
+                Token::Operator('s') => self.op_s()?,
                 Token::VarRecall(name) => self.op_recall(&name)?,
                 Token::VarStore(name) => self.op_store(&name)?,
                 Token::MacroInvoke((name, args)) => match name.as_str() {

src/interpreter/ops.rs

@@ -1,8 +1,12 @@
 use num_rational::BigRational;
+use num_traits::{FromPrimitive, ToPrimitive};
 
-use crate::quantity::{
-    units::{Unit, UnitCombo},
-    Quantity,
+use crate::{
+    linear_system::{transpose, LinearSystem},
+    quantity::{
+        units::{Unit, UnitCombo},
+        Quantity,
+    },
 };
 
 use super::{Interpreter, InterpreterError, InterpreterResult, Output};
@@ -138,6 +142,82 @@ impl<'a> Interpreter<'a> {
         self.stack.push(a);
         self.stack.push(b);
 
+        Ok(())
+    }
+    pub fn op_s(&mut self) -> InterpreterResult<()> {
+        let target = self.stack.pop().ok_or(InterpreterError::StackUnderflow)?;
+
+        let n_src_quantities = target.number_in_derived_unit().to_integer();
+        if n_src_quantities.to_usize().unwrap_or(0) > self.stack.len() {
+            return Err(InterpreterError::StackUnderflow);
+        }
+        let src_quantities = self
+            .stack
+            .split_off(self.stack.len() - n_src_quantities.to_usize().unwrap_or(0));
+        let dst_unit = target.unit.reduce();
+        let mut units_involved = Vec::new();
+        for q in &src_quantities {
+            for u in &q.unit.0 {
+                if !units_involved.contains(&u.unit) {
+                    units_involved.push(u.unit.clone());
+                }
+            }
+        }
+        for u in &dst_unit.0 {
+            if !units_involved.contains(&u.unit) {
+                return Err(InterpreterError::IncompatibleUnits(dst_unit));
+            }
+        }
+        let mut result_coefs = Vec::with_capacity(units_involved.len());
+        for u in &units_involved {
+            let mut coef = 0;
+            for u2 in &dst_unit.0 {
+                if u2.unit == *u {
+                    coef = u2.exponent;
+                    break;
+                }
+            }
+            result_coefs.push(coef);
+        }
+        let mut unit_coefs = Vec::with_capacity(src_quantities.len());
+        for q in &src_quantities {
+            let mut unit_coef = Vec::with_capacity(units_involved.len());
+            for u in &units_involved {
+                let mut coef = 0;
+                for u2 in &q.unit.0 {
+                    if u2.unit == *u {
+                        coef = u2.exponent;
+                        break;
+                    }
+                }
+                unit_coef.push(coef);
+            }
+            unit_coefs.push(unit_coef);
+        }
+        let mut lin = LinearSystem::new_equation_system(transpose(&unit_coefs), result_coefs);
+        let soln = lin.solve();
+        if soln.is_none() {
+            return Err(InterpreterError::NoSolution(
+                "failed to solve unit conversion".to_string(),
+            ));
+        }
+        if lin.is_overdetermined() {
+            return Err(InterpreterError::NoSolution(
+                "Linear system is overdetermined".to_string(),
+            ));
+        }
+        if lin.is_underdetermined() {
+            return Err(InterpreterError::NoSolution(
+                "Linear system is underdetermined".to_string(),
+            ));
+        }
+        let mut result = Quantity::new(BigRational::from_usize(1).unwrap(), dst_unit);
+        for (q, coef) in src_quantities.into_iter().zip(soln.unwrap()) {
+            result.number *= q.number.pow(coef);
+        }
+        result.use_derived_unit = target.use_derived_unit;
+        self.stack.push(result);
+
         Ok(())
     }
 }

src/lib.rs

@@ -1,3 +1,4 @@
 pub mod interpreter;
+pub mod linear_system;
 pub mod quantity;
 pub mod tokenizer;

src/linear_system.rs

@@ -0,0 +1,149 @@
+use log::debug;
+use num_traits::Signed;
+use std::{
+    fmt::Debug,
+    ops::{Add, Div, Mul, Neg, Sub},
+};
+
+pub fn transpose<T: Clone>(matrix: &Vec<Vec<T>>) -> Vec<Vec<T>> {
+    let mut result = Vec::new();
+    for j in 0..matrix[0].len() {
+        let mut row = Vec::new();
+        for i in 0..matrix.len() {
+            row.push(matrix[i][j].clone());
+        }
+        result.push(row);
+    }
+    result
+}
+
+#[derive(Debug, Clone)]
+pub struct LinearSystem<T>
+where
+    T: Add<Output = T>
+        + Sub<Output = T>
+        + Mul<Output = T>
+        + Div<Output = T>
+        + Neg<Output = T>
+        + Clone
+        + Signed,
+{
+    matrix: Vec<Vec<T>>,
+}
+
+impl<T> LinearSystem<T>
+where
+    T: Add<Output = T>
+        + Sub<Output = T>
+        + Mul<Output = T>
+        + Div<Output = T>
+        + Neg<Output = T>
+        + PartialOrd
+        + Clone
+        + Signed
+        + Debug,
+{
+    pub fn n(&self) -> usize {
+        self.matrix.len()
+    }
+    pub fn m(&self) -> usize {
+        if self.matrix.len() == 0 {
+            return 0;
+        }
+        self.matrix[0].len()
+    }
+    pub fn new(matrix: Vec<Vec<T>>) -> Self {
+        Self { matrix }
+    }
+    pub fn new_equation_system(left: Vec<Vec<T>>, right: Vec<T>) -> Self {
+        let mut matrix = left;
+        for (i, row) in matrix.iter_mut().enumerate() {
+            row.push(right[i].clone());
+        }
+        Self { matrix }
+    }
+    pub fn is_pivoted(&self, row: usize) -> bool {
+        let mut cur_col = 0;
+        while cur_col < self.m() - 1 && self.matrix[row][cur_col] == T::zero() {
+            cur_col += 1;
+        }
+        if self.matrix[row][cur_col] != T::one() {
+            return false;
+        }
+        for i in 0..self.n() {
+            if i != row && self.matrix[i][cur_col] != T::zero() {
+                return false;
+            }
+        }
+        true
+    }
+    pub fn n_pivoted(&self) -> usize {
+        self.matrix
+            .iter()
+            .enumerate()
+            .filter(|(i, _)| self.is_pivoted(*i))
+            .count()
+    }
+    fn scale_row(&mut self, i: usize, factor: T) {
+        for j in 0..self.m() {
+            self.matrix[i][j] = self.matrix[i][j].clone() * factor.clone();
+        }
+    }
+    fn add_row(&mut self, src: usize, dst: usize, factor: T) {
+        for j in 0..self.m() {
+            self.matrix[dst][j] =
+                self.matrix[dst][j].clone() + self.matrix[src][j].clone() * factor.clone();
+        }
+    }
+    fn reduce_column(&mut self, column: usize) {
+        let mut pivot_row = 0;
+        while pivot_row < self.n()
+            && (self.is_pivoted(pivot_row) || self.matrix[pivot_row][column] == T::zero())
+        {
+            pivot_row += 1;
+        }
+        if pivot_row == self.n() || self.matrix[pivot_row][column] == T::zero() {
+            return;
+        }
+        let n_pivoted = self.n_pivoted();
+        if pivot_row != n_pivoted {
+            self.matrix.swap(pivot_row, n_pivoted);
+            pivot_row = n_pivoted;
+        }
+        let scale_factor = self.matrix[pivot_row][column].clone();
+        self.scale_row(n_pivoted, T::one() / scale_factor);
+        for row in 0..self.n() {
+            if row != n_pivoted && self.matrix[row][column] != T::zero() {
+                let factor = self.matrix[row][column].clone();
+                self.add_row(n_pivoted, row, -factor);
+            }
+        }
+    }
+    pub fn is_underdetermined(&self) -> bool {
+        self.n_pivoted() > self.m() - 1
+    }
+    pub fn is_overdetermined(&self) -> bool {
+        self.n_pivoted() < self.m() - 1
+    }
+    pub fn solve(&mut self) -> Option<Vec<T>> {
+        debug!("starting matrix: {:?}", self.matrix);
+        for col in 0..if self.n() < self.m() {
+            self.n()
+        } else {
+            self.m()
+        } {
+            debug!("col: {}", col);
+            self.reduce_column(col);
+            debug!("matrix: {:?}", self.matrix);
+        }
+        let mut result = Vec::new();
+        for row in 0..self.n() {
+            if self.is_pivoted(row) {
+                result.push(self.matrix[row][self.m() - 1].clone());
+            } else if self.matrix[row][self.m() - 1] != T::zero() {
+                return None;
+            }
+        }
+        Some(result)
+    }
+}

src/quantity/mod.rs

@@ -48,6 +48,19 @@ impl Quantity {
             use_derived_unit: Vec::new(),
         }
     }
+    pub fn number_in_derived_unit(&self) -> BigRational {
+        let mut number = self.number.clone();
+
+        for d in &self.use_derived_unit {
+            if d.exponents == self.unit {
+                number -= d.offset.clone();
+                number /= d.scale.clone();
+                return number;
+            }
+        }
+
+        number
+    }
 }
 
 impl Display for Quantity {

src/tokenizer/mod.rs

@@ -151,6 +151,7 @@ impl<R: std::io::Read> Tokenizer<R> {
             Some('c') => Ok(Some(Token::Operator('c'))),
             Some('d') => Ok(Some(Token::Operator('d'))),
             Some('r') => Ok(Some(Token::Operator('r'))),
+            Some('s') => Ok(Some(Token::Operator('s'))),
             Some('#') => {
                 while let Some(c) = self.next_char().map_err(TokenizerError::IOError)? {
                     match c {