Wednesday, March 25, 2009

Higher Order Perl (Python Style) : Chapter 6 - Infinite Streams




TOC

### Chapter 6 - Infinite Streams

### 6.1 Linked Lists

# sub node {
# my ($h, $t) = @_;
# [$h, $t];
# }
# sub head {
# my ($ls) = @_;
# $ls->[0];
# }
# sub tail {
# my ($ls) = @_;
# $ls->[1];
# }
# sub set_head {
# my ($ls, $new_head) = @_;
# $ls->[0] = $new_head;
# }
# sub set_tail {
# my ($ls, $new_tail) = @_;
# $ls->[1] = $new_tail;
# }
def node(h,t):
return [h,t]
def head(ls):
return ls[0]
def tail(ls):
return ls[1]
def set_head(ls, new_head):
ls[0] = new_head
def set_tail(ls, new_tail):
ls[1] = new_tail


# $my_list = node($new_data, $my_list);
my_list = node(new_data, my_list)


# sub insert_after {
# my ($node, $new_data) = @_;
# my $new_node = node($new_data, tail($node));
# set_tail($node, $new_node);
# }
def insert_after(node, new_data):
new_node = node(new_data, tail(node))
set_tail(node, new_node)



### 6.2 Lazy Linked Lists

# package Stream;

# use base Exporter;
# @EXPORT_OK = qw(node head tail drop upto upfrom show promise
# filter transform merge list_to_stream cutsort
# iterate_function cut_loops);
# %EXPORT_TAGS = ('all' => \@EXPORT_OK);
# sub node {
# my ($h, $t) = @_;
# [$h, $t];
# }
# sub head {
# my ($s) = @_;
# $s->[0];
# }
# sub tail {
# my ($s) = @_;
# if (is_promise($s->[1])) {
# return $s->[1]->();
# }
# $s->[1];
# }
# sub is_promise {
# UNIVERSAL::isa($_[0], 'CODE');
# }
def node(h,t):
return [h,t]
def head(s):
return s[0]
def tail(s):
if is_promise(s[1]):
return s[1]()
return s[1]
def is_promise(s):
return callable(s)


# sub promise (&) { $_[0] }
# this is the silliest kind of syntactic sugar
# but this will avoid some scoping confusion and
# premature execution i was experiencing
# in later examples and makes it look more like the perl example
def promise(func):
return func


# sub upto_list {
# my ($m, $n) = @_;
# return if $m > $n;
# node($m, upto_list($m+1, $n));
# }
def upto_list(m, n):
if m > n:
return None
return node(m, upto_list(m+1, n))

# sub upto {
# my ($m, $n) = @_;
# return if $m > $n;
# node($m, promise { upto($m+1, $n) } );
# }
def upto(m, n):
if m > n:
return None
return node(m, promise(lambda: upto(m+1, n)))


# sub upfrom {
# my ($m) = @_;
# node($m, promise { upfrom($m+1) } );
# }
def upfrom(m):
return node(m, promise(lambda: upfrom(m+1)))


# sub upfrom_list {
# my ($m) = @_;
# node($m, upfrom_list($m+1) );
# }
def upfrom_list(m):
return node(m, upfrom_list(m+1))


# sub show {
# my $s = shift;
# while ($s) {
# print head($s), $";
# $s = tail($s);
# }
# print $/;
# }
def show(s):
while(s):
print head(s),
s = tail(s)
print


# sub show {
# my ($s, $n) = @_;
# while ($s && (! defined $n || $n-- > 0)) {
# print head($s), $";
# $s = tail($s);
# }
# print $/;
# }
# NOTE: i'm following the perlish code here
# but in real life i would use itertools
def show(s, n=None):
while s and (n == None or n > 0):
print head(s),
s = tail(s)
if n != None:
n -= 1
print


# sub drop {
# my $h = head($_[0]);
# $_[0] = tail($_[0]);
# return $h;
# }
# to do this in python we'd either
# need to make a stream class or
# return a pair (h,s)
# where s is the stream in it's new state
def drop(s):
h, tail = s
return h, tail # which really was a noop

# OR something like:
class Stream(object):
def __init__(self, contents):
self.contents = contents

def drop(self):
h = head(self.contents)
self.contents = tail(contents)
return h


# sub show {
# my ($s, $n) = @_;
# while ($s && (! defined $n || $n-- > 0)) {
# print drop($s), $";
# }
# print $/;
# }
### we don't really get any savings here so i'll skip this


# sub transform (&$) {
# my $f = shift;
# my $s = shift;
# return unless $s;
# node($f->(head($s)),
# promise { transform($f, tail($s)) });
# }
def transform(f, s):
if not s:
return None
return node(f(head(s)),
promise(lambda: transform(f, tail(s))))



# my $evens = transform { $_[0] * 2 } upfrom(1);
evens = transform(lambda x: x * 2, upfrom(1))


# sub filter (&$) {
# my $f = shift;
# my $s = shift;
# until (! $s || $f->(head($s))) {
# drop($s);
# }
# return if ! $s;
# node(head($s),
# promise { filter($f, tail($s)) });
# }
def filter(f, s):
while not (not s or f(head(s))):
s = tail(s)

if not s:
return None

return node(head(s),
promise(lambda: filter(f, tail(s))))


# sub iterate_function {
# my ($f, $x) = @_;
# node($x, promise { iterate_function($f, $f->($x)) });
# }
def iterate_function(f, x):
return node(x, promise(lambda: iterate_function(f, f(x))))



### 6.3 Recursive Streams

# sub carrots {
# node('carrot', promise { carrots() });
# }
# my $carrots = carrots();
def carrots():
return node("carrot", promise(carrots))


# my $carrots = node('carrot', promise { carrots() });
carrots = node("carrot", promise(carrots()))


# my $carrots = node('carrot', promise { $carrots });
carrots = node("carrot", promise(lambda: carrots))



# sub pow2_from {
# my $n = shift;
# node($n, promise {pow2_from($n*2)})
# }
# my $powers_of_2 = pow2_from(1);
def pow2_from(n):
return node(n, promise(lambda: pow2_from(n*2)))


# my $powers_of_2;
# $powers_of_2 =
# node(1, promise { transform {$_[0]*2} $powers_of_2 });
# def powers_of_2():
# return node(1, powers_of_2)
powers_of_2 = node(1, promise(lambda: transform(lambda x: x*2, powers_of_2)))


# sub tail {
# my ($s) = @_;
# if (is_promise($s->[1])) {
# $s->[1] = $s->[1]->();
# }
# $s->[1];
# }
def tail(s):
if is_promise(s[1]):
s[1] = s[1]()
return s[1]


### 6.4 The Hamming Problem



# sub is_hamming {
# my $n = shift;
# $n/=2 while $n%2 == 0;
# $n/=3 while $n%3 == 0;
# $n/=5 while $n%5 == 0;
# return $n == 1;
# }
# # Return the first $N hamming numbers
# sub hamming {
# my $N = shift;
# my @hamming;
# my $t = 1;
# until (@hamming == $N) {
# push @hamming, $t if is_hamming($t);
# $t++;
# }
# @hamming;
# }
def is_hamming(n):
while n % 2 == 0:
n /= 2
while n % 3 == 0:
n /= 3
while n % 5 == 0:
n /= 5

return n == 1
def hamming(N):
result = []
t = 1
while len(result) < N:
if is_hamming(t):
result.append(t)
t += 1
return result



# sub merge {
# my ($S, $T) = @_;
# return $T unless $S;
# return $S unless $T;
# my ($s, $t) = (head($S), head($T));
# if ($s > $t) {
# node($t, promise {merge( $S, tail($T))});
# } elsif ($s < $t) {
# node($s, promise {merge(tail($S), $T)});
# } else {
# node($s, promise {merge(tail($S), tail($T))});
# }
# }
def merge(S,T):
if not S:
return T
if not T:
return S
s, t = head(S), head(T)
if s > t:
return node(t, promise(lambda: merge(S, tail(T))))
elif s < t:
return node(s, promise(lambda: merge(tail(S), T)))
else:
return node(s, promise(lambda: merge(tail(S), tail(T))))


# sub scale {
# my ($s, $c) = @_;
# transform { $_[0]*$c } $s;
# }
def scale(s, c):
return transform(lambda x: x * c, s)


# my $hamming;
# $hamming = node(1,
# promise {
# merge(scale($hamming, 2),
# merge(scale($hamming, 3),
# scale($hamming, 5),
# ))
# }
# );
# show($hamming, 3000);
hamming = node(1,
promise(lambda: merge(scale(hamming, 2),
merge(scale(hamming, 3),
scale(hamming, 5),
))
))
### just for kicks here are two solutions from the tubes:
### - OO version : http://aspn.activestate.com/ASPN/Mail/Message/python-list/905315
### - iterator verion (my preference) : http://mail.python.org/pipermail/python-list/2005-January/303480.html
### There is also a couple variations of this in python's test suite: test_generators.py

### but i'd wager that the haskell solution of this is still the king


### 6.5 Regex String Generation


# package Regex;
# use Stream ':all';
# use base 'Exporter';
# @EXPORT_OK = qw(literal union concat star plus charclass show
# matches);

# sub literal {
# my $string = shift;
# node($string, undef);
# }
# show(literal("foo"));
# foo
def literal(s):
return node(s, None)


# sub mingle2 {
# my ($s, $t) = @_;
# return $t unless $s;
# return $s unless $t;
# node(head($s),
# node(head($t),
# promise { mingle2(tail($s),
# tail($t))
# }
# ));
# }
def mingle2(s, t):
if not s:
return t
if not t:
return s
return node(head(s),
node(head(t),
promise(lambda: mingle2(tail(s), tail(t)))))


# sub union {
# my ($h, @s) = grep $_, @_;
# return unless $h;
# return $h unless @s;
# node(head($h),
# promise {
# union(@s, tail($h));
# });
# }
def union(*streams):
streams = [_s for _s in streams if _s != None]
if len(streams) == 0:
return None

if len(streams) == 1:
return streams[0]

return node(head(streams[0]),
promise(lambda: union(*(streams[1:]+[tail(streams[0])]))))



# # generate infinite stream ($k:1, $k:2, $k:3, ...)
# sub constant {
# my $k = shift;
# my $i = shift || 1;
# my $s = node("$k:$i", promise { constant($k, $i+1) });
# }
# my $fish = constant('fish');
# show($fish, 3);
# fish:1 fish:2 fish:3
# my $soup = union($fish, constant('dog'), constant('carrot'));
# show($soup, 10);
# fish:1 dog:1 carrot:1 fish:2 dog:2 carrot:2 fish:3 dog:3 carrot:3 fish:4
def constant(k, i=1):
return node("%s:%s" % (k,i),
promise(lambda: constant(k, i+1)))
fish = constant("fish")
soup = union(fish, constant("dog"), constant("carrot"))
show(soup, 10)



# sub concat {
# my ($S, $T) = @_;
# return unless $S && $T;
# my ($s, $t) = (head($S), head($T));
# node("$s$t", promise {
# union(postcat(tail($S), $t),
# precat(tail($T), $s),
# concat(tail($S), tail($T)),
# )
# });
# }
# sub precat {
# my ($s, $c) = @_;
# transform {"$c$_[0]"} $s;
# }
# sub postcat {
# my ($s, $c) = @_;
# transform {"$_[0]$c"} $s;
# }

def concat(S,T):
if None in (S, T):
return None

s, t = head(S), head(T)
return node("%s%s" % (s,t),
promise(lambda: (union(postcat(tail(S), t),
precat(tail(T),s),
concat(tail(S),tail(T))))))

def precat(s,c):
return transform(lambda x: "%s%s" % (c,x), s)

def postcat(s,c):
return transform(lambda x: "%s%s" % (x,c), s)


# # Im /(a|b)(c|d)$/
# my $z = concat(union(literal("a"), literal("b")),
# union(literal("c"), literal("d")),
# );
# show($z);
z = (concat(union(literal("a"), literal("b")),
union(literal("c"), literal("d")),
))
show(z)


# sub star {
# my $s = shift;
# my $r;
# $r = node("", promise { concat($s, $r) });
# }
# def star(s):
# _s = s[0]
# return iterate_function(lambda x: x + _s, "")
def star(s):
r = node("", promise(lambda: concat(s, r)))
return r



# sub show {
# my ($s, $n) = @_;
# while ($s && (! defined $n || $n-- > 0)) {
# print qq{"}, drop($s), qq{"\n};
# }
# print "\n";
# }
def show(s, n=None):
while s and (n == None or n > 0):
print repr(head(s))
s = tail(s)
if n != None:
n -= 1
print


# # charclass('abc') = /[abc]$/
# sub charclass {
# my $class = shift;
# union(map literal($_), split(//, $class));
# }
def charclass(_class):
return union(*[literal(c) for c in _class])


# # plus($s) = /s+$/
# sub plus {
# my $s = shift;
# concat($s, star($s));
# }
def plus(s):
return concat(s, star(s))

# use Regex qw(concat star literal show);
# # I represent /ab*$/
# my $regex1 = concat( literal("a"),
# star(literal("b"))
# );
# show($regex1, 10);
regex1 = concat( literal("a"),
star(literal("b"))
)
show(regex1, 10)


# # I represent /(aa|b)*$/
# my $regex2 = star(union(literal("aa"),
# literal("b"),
# ));
# show($regex2, 16);
regex2 = star(union(literal("aa"),
literal("b"),
))
show(regex2, 16)


# # I represent /(ab+|c)*$/
# my $regex3 = star(union(concat( literal("a"),
# plus(literal("b"))),
# literal("c")
# ));
# show($regex3, 20);
regex3 = star(union(concat( literal("a"),
plus(literal("b"))),
literal("c")
))
show(regex3, 20)


# sub union {
# my (@s) = grep $_, @_;
# return unless @s;
# return $s[0] if @s == 1;
# my $si = index_of_shortest(@s);
# node(head($s[$si]),
# promise {
# union(map $_ == $si ? tail($s[$_]) : $s[$_],
# 0 .. $#s);
# });
# }
# curse python's lame lambdas!
def union(*streams):
streams = [_s for _s in streams if _s != None]
if len(streams) == 0:
return None

if len(streams) == 1:
return streams[0]

si = index_of_shortest(streams)

return node(head(streams[si]),
promise(lambda: union(*[_s if _i != si else tail(_s) for (_i, _s) in enumerate(streams)])))


# sub index_of_shortest {
# my @s = @_;
# my $minlen = length(head($s[0]));
# my $si = 0;
# for (1 .. $#s) {
# my $h = head($s[$_]);
# if (length($h) < $minlen) {
# $minlen = length($h);
# $si = $_;
# }
# }
# $si;
# }
def index_of_shortest(*streams):
minlin = len(head(streams[0]))
si = 0
for i in range(1,len(streams)):
h = head(streams[i])
if len(h) < minlin:
minlen = len(h)
si = i
return si
### ugg. i can't seem to get correct ordering to work. :(



# sub matches {
# my ($string, $regex) = @_;
# while ($regex) {
# my $s = drop($regex);
# return 1 if $s eq $string;
# return 0 if length($s) > length($string);
# }
# return 0;
# }
def matches(string, regex):
while regex:
s = head(regex)
regex = tail(regex)
if s == string:
return True
if len(s) > len(string):
return False
return False
# ugg. this fails when star is involved


# sub bal {
# my $contents = shift;
# my $bal;
# $bal = node("", promise {
# concat($bal,
# union($contents,
# transform {"($_[0])"} $bal,
# )
# )
# });
# }
def bal(contents):
_bal = node("",
promise(lambda: union(contents,
transform(lambda x: "(%s)" % x, _bal))))


return _bal



# sub cut_bylen {
# my ($a, $b) = @_;
# # Its OK to emit item $a if the next item in the stream is $b
# length($a) < length($b);
# }
def cut_bylen(a,b):
return len(a) < len(b)

# sub list_to_stream {
# my $node = pop;
# while (@_) {
# $node = node(pop, $node);
# }
# $node;
# }
def list_to_stream(nodes):
_nodes = nodes[:]
_node = _nodes.pop()
while _nodes:
_node = node(_nodes.pop(), _node)
return _node

# sub insert (\@$$);
# sub cutsort {
# my ($s, $cmp, $cut, @pending) = @_;
# my @emit;
# while ($s) {
# while (@pending && $cut->($pending[0], head($s))) {
# push @emit, shift @pending;
# }
# if (@emit) {
# return list_to_stream(@emit,
# promise { cutsort($s, $cmp, $cut, @pending) });
# } else {
# insert(@pending, head($s), $cmp);
# $s = tail($s);
# }
# }
# return list_to_stream(@pending, undef);
# }
def cutsort(s, _cmp, cut, _pending=[]):
pending = _pending[:]
emit = []
while s:
while pending and cut(pending[0],head(s)):
emit.append(pending.pop(0))
if emit:
return list_to_stream(emit,
promise(lambda: cutsort(s, _cmp, cut, pending)))
else:
insert(pending, head(s), _cmp)
s = tail(s)

return list_to_stream(pending+[None])



# my $sorted =
# cutsort($regex3,
# sub { $_[0] cmp $_[1] }, # comparator
# \&cut_bylen # cutting function
# );
sorted = cutsort(regex3,
lambda x,y: cmp(x,y),
cut_bylen,
)

# sub insert (\@$$) {
# my ($a, $e, $cmp) = @_;
# my ($lo, $hi) = (0, scalar(@$a));
# while ($lo < $hi) {
# my $med = int(($lo + $hi) / 2);
# my $d = $cmp->($a->[$med], $e);
# if ($d <= 0) {
# $lo = $med+1;
# } else {
# $hi = $med;
# }
# }
# splice(@$a, $lo, 0, $e);
# }
def insert(a, e, _cmp):
lo, hi = 0, len(a)
while lo < hi:
med = int((lo+hi)/2)
d = _cmp(a[med],e)
if d <= 0:
lo = med+1
else:
hi = med

splice(a, lo, 0, e)
##
## above not tested
##



# sub _devino {
# my $f = shift;
# my ($dev, $ino) = stat($f);
# return unless defined $dev;
# "$dev;$ino";
# }
def _devino(f):
stat_tuple = os.state(f)
dev, ino = stat_tuple.st_dev, stat_tuple.st_ino
if not dev:
return None
return "%s;%s" % (dev, ino)


##
## ugg... sorry got lazy again and skipped the rest of this section
##


### 6.6 The Newton-Raphson Method

# sub sqrt2 {
# my $g = 2; # Initial guess
# until (close_enough($g*$g, 2)) {
# $g = ($g*$g + 2) / (2*$g);
# }
# $g;
# }
def sqrt2():
g = 2
while not close_enough(g*g, 2):
g = float(g*g + 2) / (2*g)

return g

# sub close_enough {
# my ($a, $b) = @_;
# return abs($a - $b) < 1e-12;
# }
def close_enough(a,b):
return abs(a-b) < 1e-12


# sub sqrtn {
# my $n = shift;
# my $g = $n; # Initial guess
# until (close_enough($g*$g, $n)) {
# $g = ($g*$g + $n) / (2*$g);
# }
# $g;
# }
def sqrtn(n):
g = n
while not close_enough(g*g, n):
g = float(g*g + n) / (2*g)

return g


# use Stream 'iterate_function';
# sub sqrt_stream {
# my $n = shift;
# iterate_function (sub { my $g = shift;
# ($g*$g + $n) / (2*$g);
# },
# $n);
# }
# 1;
def sqrt_stream(n):
return iterate_function(lambda g: float(g*g + n)/(2*g),
n)

# sub iterate_function {
# my ($f, $x) = @_;
# my $s;
# $s = node($x, promise { &transform($f, $s) });
# }
def iterate_function(f, x):
s = node(x, promise(lambda: transform(f, s)))
return s


# sub slope {
# my ($f, $x) = @_;
# my $e = 0.00000095367431640625;
# ($f->($x+$e) - $f->($x-$e)) / (2*$e);
# }
def slope(f, x):
e = 0.00000095367431640625
return (f(x+e) - f(x-e)) / (2*e)


# # Return a stream of numbers $x that make $f->($x) close to 0
# sub solve {
# my $f = shift;
# my $guess = shift || 1;
# iterate_function(sub { my $g = shift;
# $g - $f->($g)/slope($f, $g);
# },
# $guess);
# }
def solve(f, guess=1):
return iterate_function(lambda g: g - f(g)/slope(f,g),
guess)


# use Math::BigFloat;
# my $sqrt2 = solve(sub { $_[0] * $_[0] - 2 },
# Math::BigFloat->new(2));
# if want to use Decimal then need to retrofit a few functions:
import decimal
decimal.getcontext().prec = 64
def slope(f, x):
x = decimal.Decimal(x)
e = decimal.Decimal("0.00000095367431640625")
return (f(x+e) - f(x-e)) / (2*e)
def solve(f, guess=1):
return iterate_function(lambda g: g - f(g)/slope(f,g),
decimal.Decimal(guess))
sqrt2 = solve(lambda x: x*x - decimal.Decimal(2),
decimal.Decimal(2))


# sub cut_loops {
# my $s = shift;
# return unless $s;
# my @previous_values = @_;
# for (@previous_values) {
# if (head($s) == $_ ) {
# return;
# }
# }
# node(head($s),
# promise { cut_loops(tail($s), head($s), @previous_values) });
# }
def cut_loops(s, *args):
if s == None:
return None
previous_values = list(args)
for v in previous_values:
if head(s) == v:
return None
return node(head(s),
promise(lambda: cut_loops(tail(s), *([head(s)] + previous_values))))



# sub cut_loops {
# my ($tortoise, $hare) = @_;
# return unless $tortoise;
# # The hare and tortoise start at the same place
# $hare = $tortoise unless defined $hare;
# # The hare moves two steps every time the tortoise moves one
# $hare = tail(tail($hare));
# # If the hare and the tortoise are in the same place, cut the loop
# return if head($tortoise) == head($hare);
# return node(head($tortoise),
# promise { cut_loops(tail($tortoise), $hare) });
# }
def cut_loops(tortoise, hare=None):
if tortoise == None:
return None
# The hare and tortoise start at the same place
if hare == None:
hare = tortoise
# The hare moves two steps every time the tortoise moves one
hare = tail(tail(hare))
# If the hare and the tortoise are in the same place, cut the loop
if head(tortoise) == head(hare):
return None
return node(head(tortoise),
promise(lambda: cut_loops(tail(tortoise), hare)))



# sub cut_loops2 {
# my ($tortoise, $hare, $n) = @_;
# return unless $tortoise;
# $hare = $tortoise unless defined $hare;
# $hare = tail(tail($hare));
# return if head($tortoise) == head($hare)
# && $n++;
# return node(head($tortoise),
# promise { cut_loops(tail($tortoise), $hare, $n) });
# }
def cut_loops2(tortoise, hare=None, n=0):
if tortoise == None:
return None
if hare == None:
hare = tortoise
hare = tail(tail(hare))
if head(tortoise) == head(hare):
if n > 0:
return None
else:
n += 1
return node(head(tortoise),
promise(lambda: cut_loops2(tail(tortoise), hare, n)))


# sub owed {
# my ($P, $N, $pmt, $i) = @_;
# my $payment_factor = 0;
# for (0 .. $N-1) {
# $payment_factor += (1+$i) ** $_;
# }
# return $P * (1+$i)**$N - $pmt * $payment_factor;
# }
def owed(P,N,pmt,i):
payment_factor = 0
for x in range(0,N):
payment_factor += (1+i) ** x
return P * (1+i)**N - pmt * payment_factor



# sub owed {
# my ($P, $N, $pmt, $i) = @_;
# return $P * (1+$i)**$N - $pmt * ((1+$i)**$N - 1) / $i;
# }
def owed(P,N,pmt,i):
return P * (1+i)**N - pmt * ((1+i)**N - 1) / i


# sub owed_after_n_months {
# my $N = shift;
# owed(100_000, $N, 1_000, 0.005);
# }
# my $stream = cut_loops(solve(\&owed_after_n_months));
# my $n;
# $n = drop($stream) while $stream;
# print "You will be paid off in only $n months!\n";
def owed_after_n_months(N):
return owed(100000, N, 1000, 0.005)
stream = cut_loops(solve(owed_after_n_months))
n = head(stream)
while stream:
n = head(stream)
stream = tail(stream)
print "You will be paid off in only %s months!" % n


# sub affordable_mortgage {
# my $mortgage = shift;
# owed($mortgage, 30, 15_600, 0.0675);
# }
# my $stream = cut_loops(solve(\&affordable_mortgage));
# my $n;
# $n = drop($stream) while $stream;
# print "You can afford a \$$n mortgage.\n";
def affordable_mortgage(mortgage):
return owed(mortgage, 30, 15600, 0.0675)
stream = cut_loops(solve(affordable_mortgage))
n = head(stream)
while stream:
n = head(stream)
stream = tail(stream)
print "You can afford a $%s mortgage." % n


### 6.7 Power Series

# # Approximate sin(x) using the first n terms of the power series
# sub approx_sin {
# my $n = shift;
# my $x = shift;
# my ($denom, $c, $num, $total) = (1, 1, $x, 0);
# while ($n--) {
# $total += $num / $denom;
# $num *= $x*$x * -1;
# $denom *= ($c+1) * ($c+2);
# $c += 2;
# }
# $total;
# }
# 1;
def approx_sin(n, x):
denom, c, num, total = 1,1,x,0
while n:
n -= 1
total += num / float(denom)
num *= x*x * - 1
denom *= (c+1) * (c+2)
c += 2
return total



# package PowSeries;
# use base 'Exporter';
# @EXPORT_OK = qw(add2 mul2 partial_sums powers_of term_values
# evaluate derivative multiply recip divide
# $sin $cos $exp $log_ $tan);
# use Stream ':all';
# sub tabulate {
# my $f = shift;
# &transform($f, upfrom(0));
# }
def tabulate(f):
return transform(f, upfrom(0))



# my @fact = (1);
# sub factorial {
# my $n = shift;
# return $fact[$n] if defined $fact[$n];
# $fact[$n] = $n * factorial($n-1);
# }

# $sin = tabulate(sub { my $N = shift;
# return 0 if $N % 2 == 0;
# my $sign = int($N/2) % 2 ? -1 : 1;
# $sign/factorial($N)
# });

# $cos = tabulate(sub { my $N = shift;
# return 0 if $N % 2 != 0;
# my $sign = int($N/2) % 2 ? -1 : 1;
# $sign/factorial($N)
# });
fact = {0:1}
def factorial(n):
if fact.get(n):
return fact[n]
fact[n] = n * factorial(n-1)
return fact[n]

def sin_lambda(N):
if N % 2 == 0:
return 0
sign = 1
if (N/2) % 2 != 0:
-1
return sign / float(factorial(N))
sin = tabulate(sin_lambda)

def cos_lambda(N):
if N % 2 != 0:
return 0
sign = 1
if (N/2) % 2 != 0:
-1
return sign / float(factorial(N))
cos = tabulate(cos_lambda)




# sub add2 {
# my ($s, $t) = @_;
# return $s unless $t;
# return $t unless $s;
# node(head($s) + head($t),
# promise { add2(tail($s), tail($t)) });
# }
def add2(s,t):
if not t:
return s
if not s:
return t
return node(head(s) + head(t),
promise(lambda: add2(tail(s),tail(t))))

# sub mul2 {
# my ($s, $t) = @_;
# return unless $s && $t;
# node(head($s) * head($t),
# promise { mul2(tail($s), tail($t)) });
# }
def mul2(s,t):
if not (s and t):
return None
return node(head(s)*head(t),
promise(lambda: mul2(tail(s),tail(t))))


# sub partial_sums {
# my $s = shift;
# my $r;
# $r = node(head($s), promise { add2($r, tail($s)) });
# }
def partial_sums(s):
r = node(head(s), promise(lambda: add2(r, tail(s))))
return r


# sub powers_of {
# my $x = shift;
# iterate_function(sub {$_[0] * $x}, 1);
# }
def powers_of(x):
return iterate_function(lambda i: i*x, 1)


# sub term_values {
# my ($s, $x) = @_;
# mul2($s, powers_of($x));
# }
def term_values(s,x):
return mul2(s, powers_of(x))


# sub evaluate {
# my ($s, $x) = @_;
# partial_sums(term_values($s, $x));
# }
def evaluate(s,x):
return partial_sums(term_values(s,x))


# my $pi = 3.1415926535897932;
# show(evaluate($cos, $pi/6), 20);
pi = 3.1415926535897932
show(evaluate(cos, pi/6.0), 20)


### grrr..... each of the pieces i test
### above seems to be behaving correctly
### but for some reason the cos(pi/6) example
### is not coming out correctly. no idea :(


# # Get the n'th term from a stream
# sub nth {
# my $s = shift;
# my $n = shift;
# return $n == 0 ? head($s) : nth(tail($s), $n-1);
# }
# # Calculate the approximate cosine of x
# sub cosine {
# my $x = shift;
# nth(evaluate($cos, $x), 20);
# }
def nth(s,n):
if n == 0:
return head(s)
else:
return nth(tail(s), n-1)

def cosine(x):
return nth(evaluate(cos, x), 20)


# sub is_zero_when_x_is_pi {
# my $x = shift;
# my $c = cosine($x/6);
# $c * $c - 3/4;
# }
# show(solve(\&is_zero_when_x_is_pi), 20);
def is_zero_when_x_is_pi(x):
c = cosine(x/6.0)
return c * c - 3/4.0

show(solve(is_zero_when_x_is_pi), 20)

# sub derivative {
# my $s = shift;
# mul2(upfrom(1), tail($s));
# }
def derivative(s):
return mul2(upfrom(1), tail(s))


# show(derivative($sin), 20);
show(derivative(sin), 20)


# exp = tabulate(sub { my $N = shift; 1/factorial($N) });
exp = tabulate(lambda N: 1.0/factorial(N))

# $log_ = tabulate(sub { my $N = shift;
# $N==0 ? 0 : (-1)**$N/-$N });
log_ = tabulate(lambda N: N!=0 and (-1)**N/-N or 0)



# sub multiply {
# my ($S, $T) = @_;
# my ($s, $t) = (head($S), head($T));
# node($s*$t,
# promise { add2(scale(tail($T), $s),
# add2(scale(tail($S), $t),
# node(0,
# promise {multiply(tail($S), tail($T))}),
# ))
# }
# );
# }
def multiply(S,T):
s, t = head(S), head(T)
return node(s*t,
promise(add2(scale(tail(T),s),
add2(scale(tail(S),t),
node(0,
promise(multiply(tail(S), tail(T))))))))


# sub scale {
# my ($s, $c) = @_;
# return if $c == 0;
# return $s if $c == 1;
# transform { $_[0]*$c } $s;
# }
def scale(s, c):
if c == 0:
return None
if c == 1:
return s
return transform(lambda i: i*c, s)


# my $one = add2(multiply($cos, $cos), multiply($sin, $sin));
# show($one, 20);
one = add2(multiply(cos, cos), multiply(sin, sin))
show(one, 20)
### unfortunately i get:
### OverflowError: long int too large to convert to float
### perhaps come back later and figure this out

# sub sum {
# my @s = grep $_, @_;
# my $total = 0;
# $total += head($_ ) for @s;
# node($total,
# promise { sum(map tail($_ ), @s) }
# );
# }
def sum_(_s):
s = [s_ for s_ in _s if s_]
total = 0
for x in s:
total += head(x)
return node(total, promise(lambda: sum_([tail(x) for x in s])))


# sub multiply {
# my ($S, $T) = @_;
# my ($s, $t) = (head($S), head($T));
# node($s*$t,
# promise { sum(scale(tail($T), $s),
# scale(tail($S), $t),
# node(0,
# promise {multiply(tail($S), tail($T))}),
# )
# }
# );
# }
def multiply(S,T):
s,t = head(S), head(T)
return node(s*t,
promise(lambda: sum_(scale(tail(T),s), scale(tail(S),t), node(0, promise(lambda: multiply(tail(S), tail(T)))))))


# # Works only if head($s) = 1
# sub recip {
# my ($s) = shift;
# my $r;
# $r = node(1,
# promise { scale(multiply($r, tail($s)), -1) });
# }
def recip(s):
r = node(1,
promise(lambda: scale(multiply(r,tail(s)), -1)))
return r


# # Works only if head($t) = 1
# sub divide {
# my ($s, $t) = @_;
# multiply($s, recip($t));
# }
# $tan = divide($sin, $cos);
# show($tan, 10);
def divide(s,t):
return multiply(s, recip(t))
tan = divide(sin, cos)
show(tan, 10)


# my @fact = (Math::BigRat->new(1));
# sub factorial {
# my $n = shift;
# return $fact[$n] if defined $fact[$n];
# $fact[$n] = $n * factorial($n-1);
# }
# python 2.6 has "fractions" library
fact = [fraction.Fraction(1,1)]
def factorial(n):
if fact.get(n):
return fact[n]
fact[n] = n * factorial(n-1)
return fact[n]

No comments: