From 01abb351b58cb6b09c0033322f8ae2d14ed48b13 Mon Sep 17 00:00:00 2001
From: John Benediktsson <mrjbq7@gmail.com>
Date: Fri, 27 Jul 2012 12:44:51 -0700
Subject: [PATCH] math.extras: more maths functions.

---
 extra/math/extras/extras-docs.factor  | 32 ++++++++++++++++++++
 extra/math/extras/extras-tests.factor |  8 +++++
 extra/math/extras/extras.factor       | 43 +++++++++++++++++++++++++--
 3 files changed, 81 insertions(+), 2 deletions(-)

diff --git a/extra/math/extras/extras-docs.factor b/extra/math/extras/extras-docs.factor
index 930fb74857..e745caa3d1 100644
--- a/extra/math/extras/extras-docs.factor
+++ b/extra/math/extras/extras-docs.factor
@@ -21,3 +21,35 @@ HELP: chi2P
 { $values { "chi" real } { "df" real } { "p" real } }
 { $description "Returns the inverse chi-squared value according to " { $snippet "P(chi|df) = P(df/2,chi/2)" } "." } ;
 
+HELP: bartlett
+{ $values { "n" integer } { "seq" sequence } }
+{ $description "Return the Bartlett window." } ;
+
+HELP: hanning
+{ $values { "n" integer } { "seq" sequence } }
+{ $description "Return the Hanning window." } ;
+
+HELP: hamming
+{ $values { "n" integer } { "seq" sequence } }
+{ $description "Return the Hamming window." } ;
+
+HELP: blackman
+{ $values { "n" integer } { "seq" sequence } }
+{ $description "Return the Blackman window." } ;
+
+HELP: nan-sum
+{ $values { "seq" sequence } { "n" number } }
+{ $description "Return the " { $link sum } " of " { $snippet "seq" } " treating any NaNs as zero." } ;
+
+HELP: nan-min
+{ $values { "seq" sequence } { "n" number } }
+{ $description "Return the " { $link infimum } " of " { $snippet "seq" } " ignoring any NaNs." } ;
+
+HELP: nan-max
+{ $values { "seq" sequence } { "n" number } }
+{ $description "Return the " { $link supremum } " of " { $snippet "seq" } " ignoring any NaNs." } ;
+
+HELP: sinc
+{ $values { "x" number } { "y" number } }
+{ $description "Returns the " { $link sinc } " function, calculated according to " { $snippet "sin(pi * x) / (pi * x)" } ". The name " { $link sinc } " is short for \"sine cardinal\" or \"sinus cardinalis\"." }
+{ $notes { $snippet "0 sinc" } " is the limit value of 1." } ;
diff --git a/extra/math/extras/extras-tests.factor b/extra/math/extras/extras-tests.factor
index aa33351d59..fe797d916e 100644
--- a/extra/math/extras/extras-tests.factor
+++ b/extra/math/extras/extras-tests.factor
@@ -23,3 +23,11 @@ IN: math.extras.test
 
 { { } } [ { 0 0 } nonzero ] unit-test
 { { 1 2 3 } } [ { 0 1 0 2 0 3 0 } nonzero ] unit-test
+
+{ { } } [ 0 bartlett ] unit-test
+{ { 1 } } [ 1 bartlett ] unit-test
+{ { 0 0 } } [ 2 bartlett ] unit-test
+{ { 0 1 0 } } [ 3 bartlett ] unit-test
+{ { 0 2/3 2/3 0 } } [ 4 bartlett ] unit-test
+{ { 0 1/2 1 1/2 0 } } [ 5 bartlett ] unit-test
+{ { 0 2/5 4/5 4/5 2/5 0 } } [ 6 bartlett ] unit-test
diff --git a/extra/math/extras/extras.factor b/extra/math/extras/extras.factor
index 77c301e6f5..fa02ea8e15 100644
--- a/extra/math/extras/extras.factor
+++ b/extra/math/extras/extras.factor
@@ -2,8 +2,9 @@
 ! See http://factorcode.org/license.txt for BSD license
 
 USING: combinators.short-circuit grouping kernel math
-math.combinatorics math.functions math.order math.primes
-math.ranges math.statistics math.vectors memoize sequences ;
+math.combinatorics math.constants math.functions math.order
+math.primes math.ranges math.statistics math.vectors memoize
+sequences ;
 
 IN: math.extras
 
@@ -97,3 +98,41 @@ PRIVATE>
 
 : nonzero ( seq -- seq' )
     [ zero? not ] filter ;
+
+: bartlett ( n -- seq )
+    dup 1 <= [ 1 = { 1 } { } ? ] [
+        [ iota ] [ 1 - 2 / ] bi [
+            [ recip * ] [ >= ] 2bi [ 2 swap - ] when
+        ] curry map
+    ] if ;
+
+: hanning ( n -- seq )
+    dup 1 <= [ 1 = { 1 } { } ? ] [
+        [ iota ] [ 1 - 2pi swap / ] bi v*n
+        [ cos -0.5 * 0.5 + ] map!
+    ] if ;
+
+: hamming ( n -- seq )
+    dup 1 <= [ 1 = { 1 } { } ? ] [
+        [ iota ] [ 1 - 2pi swap / ] bi v*n
+        [ cos -0.46 * 0.54 + ] map!
+    ] if ;
+
+: blackman ( n -- seq )
+    dup 1 <= [ 1 = { 1 } { } ? ] [
+        [ iota ] [ 1 - 2pi swap / ] bi v*n
+        [ [ cos -0.5 * ] map ] [ [ 2 * cos 0.08 * ] map ] bi
+        v+ 0.42 v+n
+    ] if ;
+
+: nan-sum ( seq -- n )
+    0 [ dup fp-nan? [ drop ] [ + ] if ] binary-reduce ;
+
+: nan-min ( seq -- n )
+    [ fp-nan? not ] filter infimum ;
+
+: nan-max ( seq -- n )
+    [ fp-nan? not ] filter supremum ;
+
+: sinc ( x -- y )
+    [ 1 ] [ pi * [ sin ] [ / ] bi ] if-zero ;