integration using simpson's method

TODO.FACTOR.txt

@@ -22,7 +22,6 @@
 - vector-each/map examples
 - string construction examples
 - string construction ackward
-- read#
 
 + tests:
 
@@ -47,6 +46,8 @@
 
 + native:
 
+- broken pipe errors with httpd and telnetd in cfactor
+- read#
 - write buffer
 - to_fixnum and friends: error on float
 ERROR: I/O error: [ "primitive_read_line_fd_8" "Resource temporarily unavailable" ]
@@ -88,7 +89,6 @@ ERROR: I/O error: [ "primitive_read_line_fd_8" "Resource temporarily unavailable
 
 + httpd:
 
-- broken pipe errors with httpd and telnetd in cfactor
 - multitasking
 - inspect: always use inspect/ URL prefix, not responder name var
 - httpd: don't flush so much

library/lists.factor

@@ -267,6 +267,36 @@ DEFER: tree-contains?
         transp over >r >r call r> cons r>
     ] each drop nreverse ; inline interpret-only
 
+: 2uncons ( list1 list2 -- car1 car2 cdr1 cdr2 )
+    uncons >r >r uncons r> swap r> ;
+
+: 2each-step ( list list quot -- cdr cdr )
+    >r 2uncons r> -rot >r >r call r> r> ; inline interpret-only
+
+: 2each ( list list quot -- )
+    #! Apply the quotation to each pair of elements from the
+    #! two lists in turn. The quotation must have stack effect
+    #! ( x y -- ).
+    >r 2dup and [
+        r> dup >r 2each-step r> 2each
+    ] [
+        r> 3drop
+    ] ifte ; inline interpret-only
+
+: 2map-step ( accum quot elt elt -- accum )
+    2swap swap >r call r> cons ;
+
+: <2map ( list list quot -- accum quot list list )
+    >r f -rot r> -rot ;
+
+: 2map ( list list quot -- list )
+    #! Apply the quotation to each pair of elements from the
+    #! two lists in turn, collecting the return value into a
+    #! new list. The quotation must have stack effect
+    #! ( x y -- z ).
+    <2map [ pick >r 2map-step r> ] 2each drop nreverse ;
+    inline interpret-only
+
 : substitute ( new old list -- list )
     [ 2dup = [ drop over ] when ] inject ;
 

library/math/arithmetic.factor

@@ -35,6 +35,9 @@ USE: stack
 : rational? dup integer? swap ratio? or ;
 : real? dup number? swap complex? not and ;
 
+: odd? 2 mod 1 = ;
+: even? 2 mod 0 = ;
+
 : i #{ 0 1 } ; inline
 : -i #{ 0 -1 } ; inline
 : inf 1.0 0.0 / ; inline

library/math/list-math.factor

@@ -25,15 +25,28 @@
 ! OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
 ! ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
-IN: math USE: arithmetic USE: lists USE: stack
+IN: math
+USE: arithmetic
+USE: lists
+USE: stack
 
-: [+] ( list -- sum )
+: |+ ( list -- sum )
     #! sum all elements in a list.
     0 swap [ + ] each ;
 
-: [*] ( list -- sum )
+: +| ( list list -- list )
+    [ + ] 2map ;
+
+: |* ( list -- sum )
     #! multiply all elements in a list.
     1 swap [ * ] each ;
 
+: *| ( list list -- list )
+    [ * ] 2map ;
+
+: *|+ ( list list -- dot )
+    #! Dot product
+    *| |+ ;
+
 : average ( list -- avg )
-    dup [+] swap length / ;
+    dup |+ swap length / ;

library/math/simpson.factor

@@ -0,0 +1,70 @@
+! :folding=indent:collapseFolds=0:
+
+! $Id$
+!
+! Copyright (C) 2004 Slava Pestov.
+! 
+! Redistribution and use in source and binary forms, with or without
+! modification, are permitted provided that the following conditions are met:
+! 
+! 1. Redistributions of source code must retain the above copyright notice,
+!    this list of conditions and the following disclaimer.
+! 
+! 2. Redistributions in binary form must reproduce the above copyright notice,
+!    this list of conditions and the following disclaimer in the documentation
+!    and/or other materials provided with the distribution.
+! 
+! THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES,
+! INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
+! FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
+! DEVELOPERS AND CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+! PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
+! OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
+! WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
+! OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
+! ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+IN: math
+USE: arithmetic
+USE: combinators
+USE: kernel
+USE: lists
+USE: logic
+USE: stack
+
+: multiplier ( n -- 2|4 )
+    odd? 4 2 ? ;
+
+: (multipliers) ( list n -- list )
+    dup 2 <= [
+        drop
+    ] [
+        dup >r multiplier swons r> pred (multipliers)
+    ] ifte ;
+
+: multipliers ( n -- list )
+    #! The value n must be odd. Makes a list like [ 1 4 2 4 1 ]
+    [ 1 ] swap (multipliers) 1 swons ;
+
+: x-values ( lower upper n -- list )
+    #! The value n must be odd.
+    pred >r over - r> dup succ count [
+        >r 3dup r> swap / * +
+    ] inject >r 3drop r> ;
+
+: y-values ( lower upper n quot -- values )
+    >r x-values r> inject ;
+
+: (simpson) ( lower upper n quot -- value )
+    over multipliers >r y-values r> *|+ ;
+
+: h ( lower upper n -- h )
+    transp - swap pred / 3 / ;
+
+: simpson ( lower upper n quot -- value )
+    #! Compute the integral between the lower and upper bound,
+    #! using Simpson's method with n steps. The value of n must
+    #! be odd. The quotation must have stack effect
+    #! ( x -- f(x) ).
+    >r 3dup r> (simpson) >r h r> * ;

library/platform/jvm/boot-sumo.factor

@@ -82,6 +82,7 @@ USE: parser
 "/library/math/arc-trig-hyp.factor"      run-resource ! math
 "/library/math/quadratic.factor"         run-resource ! math
 "/library/math/list-math.factor"         run-resource ! math
+"/library/math/simpson.factor"           run-resource ! math
 
 !!! Development tools.
 "/library/vocabulary-style.factor"         run-resource ! style

library/platform/native/boot.factor

@@ -96,6 +96,7 @@ primitives,
     "/library/math/pow.factor"
     "/library/math/quadratic.factor"
     "/library/math/trig-hyp.factor"
+    "/library/math/simpson.factor"
     "/library/test/test.factor"
     "/library/platform/native/errors.factor"
     "/library/platform/native/io-internals.factor"

library/test/math/simpson.factor

@@ -0,0 +1,20 @@
+IN: scratchpad
+USE: arithmetic
+USE: math
+USE: test
+
+[ 2 ] [ 0 multiplier ] unit-test
+[ 4 ] [ 1 multiplier ] unit-test
+
+[ [ 1 4 1 ] ] [ 3 multipliers ] unit-test
+[ [ 1 4 2 4 1 ] ] [ 5 multipliers ] unit-test
+[ [ 1 4 2 4 2 4 1 ] ] [ 7 multipliers ] unit-test
+
+[ [ 0 5 10 ] ] [ 0 10 3 x-values ] unit-test
+[ [ 10 15 20 ] ] [ 10 20 3 x-values ] unit-test
+
+[ [ 0 1 4 9 16 ] ] [ 0 4 5 [ sq ] y-values ] unit-test
+
+[ 5/3 ] [ 10 20 3 h ] unit-test
+
+[ 1000/3 ] [ 0 10 3 [ sq ] simpson-try ] unit-test

library/test/test.factor

@@ -85,6 +85,7 @@ USE: vocabularies
         "math/float"
         "math/complex"
         "math/irrational"
+        "math/simpson"
         "httpd/url-encoding"
         "httpd/html"
         "httpd/httpd"