changeset 445:38c681e9d9ec

Add pow, round, ceil, isPowerOfTwo, nextPowerOfTwo

--- a/src/may/mathmisc.yeti	Wed Oct 23 13:59:21 2013 +0100
+++ b/src/may/mathmisc.yeti	Wed Oct 23 17:28:36 2013 +0100
@@ -1,12 +1,66 @@
 
 module may.mathmisc;
 
+/**
+ * The greatest common divisor of natural numbers a and b.
+ */
 gcd a b = (c = a % b; if c == 0 then b else gcd b c fi);
 
 factorial x = if x < 0 then 0 else fold do x y: x*y done 1 [1..x] fi;
 
+pow x y = 
+    if y < 0 then 1 / (pow x (-y))
+    elif y == 0 then 1
+    elif y == 1 then x
+    elif y == 2 then x*x
+    elif y == int(y) and y > 0
+    then
+        ipow a b = fold do x _: x*a done 1 [1..b]; // a^b where b∈ℕ
+        ipow x y
+    else
+        Math#pow(x, y) 
+    fi;
+
+round x =
+    if x < 0 then
+        -(round (-x))
+    else
+        int (x + 0.5)
+    fi;
+
+ceil x =
+    Math#ceil(x);
+
+/** 
+ * True if x is 2^n for some integer n >= 0.
+ */
+isPowerOfTwo x =
+    x >= 1 and
+       (x == 1 or
+           (half = x/2;
+            half == int half and isPowerOfTwo half));
+
+/**
+ * The next higher integer power of two from x, e.g. 1300 -> 2048, 2048 -> 2048.
+ */
+nextPowerOfTwo x =
+    if x <= 1 then 1
+    elif isPowerOfTwo x then x
+    else
+        np i n =
+            if i == 0 then n
+            else np (int (i/2)) (n * 2)
+            fi;
+        np x 1;
+    fi;
+
 {
     gcd,
     factorial,
+    pow,
+    round,
+    ceil,
+    isPowerOfTwo,
+    nextPowerOfTwo,
 }
 

--- a/src/may/mathmisc/test/test_mathmisc.yeti	Wed Oct 23 13:59:21 2013 +0100
+++ b/src/may/mathmisc/test/test_mathmisc.yeti	Wed Oct 23 17:28:36 2013 +0100
@@ -29,5 +29,73 @@
     compare (mm.factorial 20) 2432902008176640000
 ),
 
+"pow": \(
+    compare (mm.pow 0 1) 0 and
+    compare (mm.pow 1 0) 1 and
+    compare (mm.pow 2 0) 1 and
+    compare (mm.pow 3 1) 3 and
+    compare (mm.pow 5 2) 25 and
+    compare (mm.pow 5.5 2) 30.25 and
+    compare (mm.pow 5.5 3) 166.375 and
+    compare (mm.pow 4 (-1)) (1/4) and
+    compare (mm.pow 4 (-2)) (1/16) and
+    compare (mm.pow 4.1 5) 1158.56201 and
+    compare (mm.pow 4.1 5.1) 1334.1285230128055
+),
+    
+"round": \(
+    compare (mm.round 0.5) 1.0 and
+    compare (mm.round 0.49) 0.0 and
+    compare (mm.round 0.99) 1.0 and
+    compare (mm.round 0.01) 0.0 and
+    compare (mm.round 0.0) 0.0 and
+    compare (mm.round 100.0) 100.0 and
+    compare (mm.round (-0.2)) 0.0 and
+    compare (mm.round (-0.5)) (-1.0) and
+    compare (mm.round (-0.99)) (-1.0) and
+    compare (mm.round (-1.0)) (-1.0) and
+    compare (mm.round (-1.1)) (-1.0) and
+    compare (mm.round (-1.5)) (-2.0)
+),
+
+"ceil": \(
+    compare (mm.ceil 0.5) 1.0 and
+    compare (mm.ceil 0.49) 1.0 and
+    compare (mm.ceil 0.99) 1.0 and
+    compare (mm.ceil 0.01) 1.0 and
+    compare (mm.ceil 0.0) 0.0 and
+    compare (mm.ceil 100.0) 100.0 and
+    compare (mm.ceil (-0.2)) 0.0 and
+    compare (mm.ceil (-0.5)) 0.0 and
+    compare (mm.ceil (-0.99)) 0.0 and
+    compare (mm.ceil (-1.0)) (-1.0) and
+    compare (mm.ceil (-1.1)) (-1.0) and
+    compare (mm.ceil (-1.5)) (-1.0)
+),
+
+"isPowerOfTwo": \(
+    compare (mm.isPowerOfTwo 0) false and
+    compare (mm.isPowerOfTwo 1) true and
+    compare (mm.isPowerOfTwo (-2)) false and
+    compare (mm.isPowerOfTwo 2) true and
+    compare (mm.isPowerOfTwo 3) false and
+    compare (mm.isPowerOfTwo 12) false and
+    compare (mm.isPowerOfTwo 12.3) false and
+    compare (mm.isPowerOfTwo 16) true and
+    compare (mm.isPowerOfTwo 15.9) false
+),
+
+"nextPowerOfTwo": \(
+    compare (mm.nextPowerOfTwo 0) 1 and
+    compare (mm.nextPowerOfTwo 1) 1 and
+    compare (mm.nextPowerOfTwo (-2)) 1 and
+    compare (mm.nextPowerOfTwo 2) 2 and
+    compare (mm.nextPowerOfTwo 3) 4 and
+    compare (mm.nextPowerOfTwo 12) 16 and
+    compare (mm.nextPowerOfTwo 12.3) 16 and
+    compare (mm.nextPowerOfTwo 16) 16 and
+    compare (mm.nextPowerOfTwo 15.9) 16
+),
+
 ] is hash<string, () -> boolean>;
 

--- a/src/may/signal/window.yeti	Wed Oct 23 13:59:21 2013 +0100
+++ b/src/may/signal/window.yeti	Wed Oct 23 17:28:36 2013 +0100
@@ -102,13 +102,12 @@
 //!!! if kaiser is going to take a structure, so I guess should sinc and cosineWindow
 calculateKaiser { length, beta } =
    (terms = 20;
-    ipow a b = fold do x _: x*a done 1 [1..b]; // a^b where b∈ℕ
     factorials = array (map maths.factorial [1..terms]);
     bes0 x = 
         1 + sum
            (map do i:
                 f = factorials[i-1]; // this is i!, not (i-1)!
-                (ipow (x/2) (i*2)) / (f * f);
+                (maths.pow (x/2) (i*2)) / (f * f);
                 done [1..terms]);
     denominator = bes0 beta;
     even = (length % 2 == 0);