add (^) as an evaluable factor, re: #39

Cargo.toml

@@ -4,9 +4,9 @@ version = "0.7.9"
 authors = ["Mark Thom"]
 
 [dependencies]
+downcast = "0.9.1"
 num = "0.2"
 ordered-float = "0.5.0"
-downcast = "0.9.1"
 
 [dependencies.termion]
 version = "1.4.0"

README.md

@@ -28,7 +28,7 @@ Extend rusty-wam to include the following, among other features:
 * Built-in control operators (`,`, `;`, `->`, etc.) (_done_).
 * A revised, not-terrible module system (_done, I think_).
 * Built-in predicates for list processing and top-level declarative
-  control (`setup_call_control/3`, `call_with_inference_limit/3`,
+  control (`setup_call_cleanup/3`, `call_with_inference_limit/3`,
   etc.) (_done_)
 * Default representation of strings as list of chars, using a packed
   internal representation (_done_).
@@ -110,7 +110,7 @@ information.
 The following predicates are built-in to rusty-wam.
 
 * Arithmetic support:
-    * `is/2` works for `(+)/2`, `(-)/{1,2}`, `(*)/2`, `(//)/2`, `(div)/2`, `(/)/2`, `(rdiv)/2`,
+    * `is/2` works for `(+)/2`, `(-)/{1,2}`, `(*)/2`, `(//)/2`, `(^)/2`, `(div)/2`, `(/)/2`, `(rdiv)/2`,
       `(xor)/2`, `(rem)/2`, `(mod)/2`, `(/\)/2`, `(\/)/2`, `(>>)/2`, `(<<)/2`.
     * Comparison operators: `>`, `<`, `=<`, `>=`, `=:=`, `=\=`.
 * `(:)/2`

src/prolog/arithmetic.rs

@@ -105,6 +105,7 @@ impl<'a> ArithmeticEvaluator<'a>
                        -> Result<ArithmeticInstruction, ArithmeticError>
     {
         match name.as_str() {
+            "abs" => Ok(ArithmeticInstruction::Abs(a1, t)),
             "-" => Ok(ArithmeticInstruction::Neg(a1, t)),
              _  => Err(ArithmeticError::InvalidOp)
         }
@@ -121,6 +122,7 @@ impl<'a> ArithmeticEvaluator<'a>
             "div"  => Ok(ArithmeticInstruction::FIDiv(a1, a2, t)),
             "rdiv" => Ok(ArithmeticInstruction::RDiv(a1, a2, t)),
             "*"    => Ok(ArithmeticInstruction::Mul(a1, a2, t)),
+            "^"    => Ok(ArithmeticInstruction::Pow(a1, a2, t)),
             ">>"   => Ok(ArithmeticInstruction::Shr(a1, a2, t)),
             "<<"   => Ok(ArithmeticInstruction::Shl(a1, a2, t)),
             "/\\"  => Ok(ArithmeticInstruction::And(a1, a2, t)),

src/prolog/ast.rs

@@ -1,6 +1,6 @@
-use prolog::num::bigint::BigInt;
+use prolog::num::bigint::{BigInt, BigUint};
-use prolog::num::{Float, ToPrimitive, Zero};
+use prolog::num::{Float, Integer, One, Signed, ToPrimitive, Zero};
-use prolog::num::rational::Ratio;
+use prolog::num::rational::{BigRational, Ratio};
 use prolog::ordered_float::*;
 use prolog::string_list::*;
 use prolog::tabled_rc::*;
@@ -11,7 +11,7 @@ use std::collections::{BTreeSet, HashMap, VecDeque};
 use std::fmt;
 use std::hash::{Hash, Hasher};
 use std::io::Error as IOError;
-use std::ops::{Add, AddAssign, Div, Index, IndexMut, Sub, Mul, Neg};
+use std::ops::{Add, AddAssign, Div, Index, IndexMut, Sub, Mul, MulAssign, Neg};
 use std::rc::Rc;
 use std::str::Utf8Error;
 use std::vec::Vec;
@@ -476,6 +476,7 @@ pub enum ArithmeticError {
     InvalidAtom,
     InvalidOp,
     InvalidTerm,
+    NoRoots,
     UninstantiatedVar
 }
 
@@ -1173,7 +1174,136 @@ impl Default for Number {
     }
 }
 
+fn binary_pow<T>(mut n: T, mut power: BigUint) -> T
+    where T: Clone + Mul + One,
+    for<'a> T: MulAssign<&'a T>
+{
+    if power.is_zero() {
+        return T::one();
+    }
+
+    let mut oddand = T::one();
+    let one = BigUint::one();
+
+    while power > one {
+        if power.is_odd() {
+            oddand *= &n;
+        }
+
+        n = n.clone() * n;
+        power >>= 1;
+    }
+
+    n * oddand
+}
+
+fn rational_pow(r1: BigRational, r2: BigRational) -> Result<BigRational, ArithmeticError>
+{
+    #[inline]
+    fn to_unsigned(n: &BigInt) -> Result<BigUint, ArithmeticError> {
+        n.abs().to_biguint().ok_or(ArithmeticError::NoRoots)
+    };
+
+    #[inline]
+    fn to_big_rational(n: BigUint) -> BigRational {
+        BigRational::from_integer(BigInt::from(n))
+    };
+
+    let r2 = r2.reduced(); // so that gcd(numer, denom) = 1
+    let n = to_unsigned(r2.denom())?;
+
+    if n == BigUint::one() {
+        return if r2.is_positive() {
+            Ok(binary_pow(r1, to_unsigned(&r2.numer())?))
+        } else if r2.is_negative() {
+            Ok(binary_pow(r1, to_unsigned(&r2.numer())?).recip())
+        } else {
+            Ok(BigRational::one())
+        };
+    }
+
+    if (n.is_even() && r1.is_negative()) || (r2.is_negative() && r1.is_zero()) {
+        return Err(ArithmeticError::NoRoots);
+    }
+
+    let sgn = r1.signum();
+    let r1 = r1 * &sgn; // set r1 to its absolute value.
+
+    let epsilon = BigRational::new_raw(BigInt::one(), BigInt::from(10000));
+    let n1      = n.clone() - BigUint::one(); // n -1
+
+    // 1 + r1 / (n-1) is a good initial point.
+    let mut x_i = BigRational::one() + r1.clone() / to_big_rational(n1.clone());
+    let mut x_i_n1 = binary_pow(x_i.clone(), n1.clone()); // x_i^{n-1}
+    let mut delta_x_i = BigRational::one();
+
+    while delta_x_i.abs() > epsilon {
+        x_i = x_i.reduced();
+        x_i_n1 = x_i_n1.reduced();
+
+        let r_quot = r1.clone() / &x_i_n1; // r1 / x_i^{n-1}
+        let r_n    = to_big_rational(n.clone());
+
+        delta_x_i = ( r_quot - &x_i ) / &r_n;
+
+        x_i += &delta_x_i;
+        x_i_n1 = binary_pow(x_i.clone(), n1.clone());
+    }
+
+    if r2.is_positive() {
+        Ok(binary_pow(sgn * x_i, to_unsigned(r2.numer())?))
+    } else {
+        Ok(binary_pow(sgn * x_i, to_unsigned(r2.numer())?).recip())
+    }
+}
+
+fn pow_float(f1: f64, f2: f64) -> Result<Number, ArithmeticError> {
+    let result = OrderedFloat(f1.powf(f2));
+
+    if result.is_finite() {
+        Ok(Number::Float(result))
+    } else {
+        Err(ArithmeticError::NoRoots)
+    }
+}
+
+fn rational_to_f64(r: &BigRational) -> Option<f64> {
+    match (r.numer().to_f64(), r.denom().to_f64()) {
+        (Some(ref f1), Some(ref f2)) if f2.is_normal() => Some(*f1 / *f2),
+        _ => None
+    }
+}
+
 impl Number {
+    pub fn pow(self, other: Number) -> Result<Self, ArithmeticError> {
+        match NumberPair::from(self, other) {
+            NumberPair::Integer(n1, n2) =>
+                if let Some(n2) = n2.to_biguint() {
+                    Ok(Number::Integer(Rc::new(binary_pow((*n1).clone(), n2))))
+                } else if n1.is_zero() {
+                    Err(ArithmeticError::NoRoots)
+                } else {
+                    let r1 = Ratio::new(BigInt::one(), (*n1).clone());
+                    let n2 = n2.abs().to_biguint().unwrap();
+
+                    Ok(Number::Rational(Rc::new(binary_pow(r1, n2))))
+                },
+            NumberPair::Float(n1, n2) =>
+                pow_float(n1.into_inner(), n2.into_inner()),
+            NumberPair::Rational(r1, r2) => {
+                if let (Some(f1), Some(f2)) = (rational_to_f64(&r1), rational_to_f64(&r2)) {                    
+                    if let Ok(result) = pow_float(f1, f2) {
+                        return Ok(result);
+                    }
+                }
+                
+                let root = rational_pow((*r1).clone(), (*r2).clone())?;
+                Ok(Number::Rational(Rc::new(root)))
+            }                    
+        }
+    }
+
+    #[inline]
     pub fn is_zero(&self) -> bool {
         match self {
             &Number::Float(fl)       => fl.into_inner().is_zero(),
@@ -1181,6 +1311,15 @@ impl Number {
             &Number::Rational(ref r) => r.is_zero()
         }
     }
+
+    #[inline]
+    pub fn abs(&self) -> Self {
+        match self {
+            &Number::Float(ref fl)   => Number::Float(OrderedFloat(fl.into_inner().abs())),
+            &Number::Integer(ref n)  => Number::Integer(Rc::new((*n).clone().abs())),
+            &Number::Rational(ref r) => Number::Rational(Rc::new((*r).clone().abs()))
+        }
+    }
 }
 
 pub enum NumberPair {
@@ -1364,6 +1503,7 @@ pub enum ArithmeticInstruction {
     Add(ArithmeticTerm, ArithmeticTerm, usize),
     Sub(ArithmeticTerm, ArithmeticTerm, usize),
     Mul(ArithmeticTerm, ArithmeticTerm, usize),
+    Pow(ArithmeticTerm, ArithmeticTerm, usize),
     IDiv(ArithmeticTerm, ArithmeticTerm, usize),
     FIDiv(ArithmeticTerm, ArithmeticTerm, usize),
     RDiv(ArithmeticTerm, ArithmeticTerm, usize),
@@ -1375,7 +1515,8 @@ pub enum ArithmeticInstruction {
     Or(ArithmeticTerm, ArithmeticTerm, usize),
     Mod(ArithmeticTerm, ArithmeticTerm, usize),
     Rem(ArithmeticTerm, ArithmeticTerm, usize),
-    Neg(ArithmeticTerm, usize)
+    Abs(ArithmeticTerm, usize),
+    Neg(ArithmeticTerm, usize),
 }
 
 #[derive(Clone)]

src/prolog/io.rs

@@ -244,12 +244,16 @@ impl fmt::Display for ArithmeticTerm {
 impl fmt::Display for ArithmeticInstruction {
     fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
         match self {
+            &ArithmeticInstruction::Abs(ref a1, ref t) =>
+                write!(f, "abs {}, @{}", a1, t),
             &ArithmeticInstruction::Add(ref a1, ref a2, ref t) =>
                 write!(f, "add {}, {}, @{}", a1, a2, t),
             &ArithmeticInstruction::Sub(ref a1, ref a2, ref t) =>
                 write!(f, "sub {}, {}, @{}", a1, a2, t),
             &ArithmeticInstruction::Mul(ref a1, ref a2, ref t) =>
                 write!(f, "mul {}, {}, @{}", a1, a2, t),
+            &ArithmeticInstruction::Pow(ref a1, ref a2, ref t) =>
+                write!(f, "pow {}, {}, @{}", a1, a2, t),
             &ArithmeticInstruction::Div(ref a1, ref a2, ref t) =>
                 write!(f, "div {}, {}, @{}", a1, a2, t),
             &ArithmeticInstruction::IDiv(ref a1, ref a2, ref t) =>

src/prolog/lib/builtins.pl

@@ -1,6 +1,6 @@
 :- op(400, yfx, /).
 
-:- module(builtins, [(=)/2, (+)/2, (*)/2, (-)/2, (/)/2, (/\)/2,
+:- module(builtins, [(=)/2, (+)/2, (^)/2, (*)/2, (-)/2, (/)/2, (/\)/2,
 	(\/)/2, (is)/2, (xor)/2, (div)/2, (//)/2, (rdiv)/2, (<<)/2,
 	(>>)/2, (mod)/2, (rem)/2, (>)/2, (<)/2, (=\=)/2, (=:=)/2,
 	(-)/1, (>=)/2, (=<)/2, (,)/2, (->)/2, (;)/2, (=..)/2, (==)/2,
@@ -19,6 +19,7 @@
 :- op(500, yfx, +).
 :- op(500, yfx, -).
 :- op(400, yfx, *).
+:- op(200, xfy, ^).
 :- op(500, yfx, /\).
 :- op(500, yfx, \/).
 :- op(500, yfx, xor).

src/prolog/machine/machine_errors.rs

@@ -233,7 +233,8 @@ pub enum EvalError {
 //    IntOverflow,
 //    Undefined,
 //    Underflow,
-    ZeroDivisor
+    ZeroDivisor,
+    NoRoots
 }
 
 impl EvalError {
@@ -243,7 +244,8 @@ impl EvalError {
 //            EvalError::IntOverflow => "int_overflow",
 //            EvalError::Undefined => "undefined",
 //            EvalError::Underflow => "underflow",
-            EvalError::ZeroDivisor => "zero_divisor"
+            EvalError::ZeroDivisor => "zero_divisor",
+            EvalError::NoRoots => "no_roots"
         }
     }
 }

src/prolog/machine/machine_state_impl.rs

@@ -429,6 +429,7 @@ impl MachineState {
                         "-" => interms.push(a1 - a2),
                         "*" => interms.push(a1 * a2),
                         "/" => interms.push(self.div(a1, a2)?),
+                        "^" => interms.push(self.pow(a1, a2)?),
                         "rdiv" => {
                             let r1 = self.get_rational(&ArithmeticTerm::Number(a1), &caller)?;
                             let r2 = self.get_rational(&ArithmeticTerm::Number(a2), &caller)?;
@@ -535,6 +536,18 @@ impl MachineState {
         }
     }
 
+    fn pow(&self, n1: Number, n2: Number) -> Result<Number, MachineStub>
+    {
+        match n1.pow(n2) {
+            Ok(result) => Ok(result),
+            Err(_) => {
+                let stub = MachineError::functor_stub(clause_name!("^"), 2);
+                Err(self.error_form(MachineError::evaluation_error(EvalError::NoRoots),
+                                    stub))
+            }
+        }
+    }
+
     fn shr(&self, n1: Number, n2: Number) -> Result<Rc<BigInt>, MachineStub>
     {
         let stub = MachineError::functor_stub(clause_name!("(>>)"), 2);
@@ -700,6 +713,13 @@ impl MachineState {
                 self.interms[t - 1] = n1 * n2;
                 self.p += 1;
             },
+            &ArithmeticInstruction::Pow(ref a1, ref a2, t) => {
+                let n1 = try_or_fail!(self, self.get_number(a1));
+                let n2 = try_or_fail!(self, self.get_number(a2));
+
+                self.interms[t - 1] = try_or_fail!(self, self.pow(n1, n2));
+                self.p += 1;
+            },
             &ArithmeticInstruction::RDiv(ref a1, ref a2, t) => {
                 let stub = MachineError::functor_stub(clause_name!("(rdiv)"), 2);
 
@@ -723,6 +743,12 @@ impl MachineState {
                 self.interms[t - 1] = Number::Integer(try_or_fail!(self, self.idiv(n1, n2)));
                 self.p += 1;
             },
+            &ArithmeticInstruction::Abs(ref a1, t) => {
+                let n1 = try_or_fail!(self, self.get_number(a1));
+
+                self.interms[t - 1] = n1.abs();
+                self.p += 1;
+            },
             &ArithmeticInstruction::Neg(ref a1, t) => {
                 let n1 = try_or_fail!(self, self.get_number(a1));
 

src/tests.rs

@@ -1019,6 +1019,75 @@ fn test_queries_on_arithmetic()
     assert_prolog_success!(&mut wam, "?- f(5, Sum).", [["Sum = 33"]]);
     assert_prolog_success!(&mut wam, "?- f(5, 33).");
     assert_prolog_failure!(&mut wam, "?- f(5, 32).");
+
+    // exponentiation.    
+
+    // the ~ operators tests whether |X - Y| <= 1/10000...
+    // or whatever degree of approximation used by Newton's method in rational_pow.
+    submit(&mut wam, ":- op(900, xfx, ~).");
+    submit(&mut wam, "X ~ Y :- abs(X - Y) =< 1 rdiv 10000.");
+
+    assert_prolog_success!(&mut wam, "?- X is 3 ^ 3.",
+                           [["X = 27"]]);
+    assert_prolog_success!(&mut wam, "?- X is 3 ^ 0.",
+                           [["X = 1"]]);
+    assert_prolog_success!(&mut wam, "?- X is 3 ^ -0.",
+                           [["X = 1"]]);
+    assert_prolog_success!(&mut wam, "?- X is 3 ^ 1.",
+                           [["X = 3"]]);
+    assert_prolog_success!(&mut wam, "?- X is 3 ^ -3.",
+                           [["X = 1/27"]]);
+    assert_prolog_success!(&mut wam, "?- X is (-3) ^ 3.",
+                           [["X = -27"]]);
+    assert_prolog_success!(&mut wam, "?- X is (-3) ^ 3.",
+                           [["X = -27"]]);
+    assert_prolog_success!(&mut wam, "?- X is (-3) ^ 0.",
+                           [["X = 1"]]);
+    assert_prolog_success!(&mut wam, "?- X is (-3) ^ -0.",
+                           [["X = 1"]]);
+    assert_prolog_success!(&mut wam, "?- X is (-3) ^ 1.",
+                           [["X = -3"]]);
+    assert_prolog_success!(&mut wam, "?- X is (-3) ^ -3.",
+                           [["X = -1/27"]]);
+    assert_prolog_success!(&mut wam, "?- X is (1 rdiv 27) ^ -3, X ~ 19683.");
+    assert_prolog_success!(&mut wam, "?- X is (-1 rdiv 27) ^ -3, X ~ -19683.");
+
+    assert_prolog_success!(&mut wam, "?- X is 0.0 ^ 0.",
+                           [["X = 1"]]);
+    assert_prolog_success!(&mut wam, "?- catch(_ is 0.0 ^ -2342, error(E, _), true).",
+                           [["E = evaluation_error(no_roots)"]]);
+    assert_prolog_success!(&mut wam, "?- X is 0.0 ^ 2342.",
+                           [["X = 0"]]);
+
+    assert_prolog_success!(&mut wam, "?- catch(_ is (-3) ^ (1 rdiv 2), error(E, _), true).",
+                           [["E = evaluation_error(no_roots)"]]);
+    assert_prolog_success!(&mut wam, "?- catch(_ is (-3/2) ^ (1 rdiv 2), error(E, _), true).",
+                           [["E = evaluation_error(no_roots)"]]);
+    assert_prolog_success!(&mut wam, "?- catch(_ is (-3 rdiv 2) ^ (1 rdiv 4), error(E, _), true).",
+                           [["E = evaluation_error(no_roots)"]]);
+    assert_prolog_success!(&mut wam, "?- catch(_ is (-3 rdiv 2) ^ (-1 rdiv 4), error(E, _), true).",
+                           [["E = evaluation_error(no_roots)"]]);
+    assert_prolog_success!(&mut wam, "?- catch(_ is 0 ^ (-5 rdiv 4), error(E, _), true).",
+                           [["E = evaluation_error(no_roots)"]]);
+
+    assert_prolog_success!(&mut wam, "?- X is 3 ^ (1 rdiv 3), Y is X ^ 3, Y ~ 3.");
+    assert_prolog_success!(&mut wam, "?- X is (-3) ^ (1 rdiv 3), Y is X ^ 3, Y ~ -3.");
+    assert_prolog_failure!(&mut wam, "?- X is (-5) ^ (1 rdiv 3), Y is X ^ 3, Y ~ -3.");
+    assert_prolog_failure!(&mut wam, "?- X is 5 ^ (1 rdiv 3), Y is X ^ 3, Y ~ 3.");
+    assert_prolog_failure!(&mut wam, "?- X is (1 rdiv 3) ^ 0.5, Y is X ^ 2, X ~ Y.");
+    assert_prolog_success!(&mut wam, "?- X is (1 rdiv 3) ^ 0.5, Y is X ^ 2, 1 rdiv 3 ~ Y.");
+
+    assert_prolog_success!(&mut wam, "?- X is (-5) ^ (-1 rdiv 3), Y is X ^ 3, Y ~ -1 rdiv 5.");
+    assert_prolog_failure!(&mut wam, "?- X is (-5) ^ (-1 rdiv 3), Y is X ^ 3, Y ~ 1 rdiv 5.");    
+    
+    assert_prolog_success!(&mut wam, "?- X is (0 rdiv 5) ^ 5.",
+                           [["X = 0"]]);
+    assert_prolog_success!(&mut wam, "?- X is (-0 rdiv 5) ^ 5.",
+                           [["X = 0"]]);
+    assert_prolog_success!(&mut wam, "?- X is (0 rdiv 5) ^ 0.",
+                           [["X = 1"]]);
+    assert_prolog_success!(&mut wam, "?- catch(_ is (0 rdiv 0) ^ 5, error(E, _), true).",
+                           [["E = evaluation_error(zero_divisor)"]]);
 }
 
 #[test]