1 // $G $D/$F.go && $L $F.$A && ./$A.out
3 // Copyright 2009 The Go Authors. All rights reserved.
4 // Use of this source code is governed by a BSD-style
5 // license that can be found in the LICENSE file.
7 // Power series package
8 // A power series is a channel, along which flow rational
9 // coefficients. A denominator of zero signifies the end.
10 // Original code in Newsqueak by Doug McIlroy.
11 // See Squinting at Power Series by Doug McIlroy,
12 // http://www.cs.bell-labs.com/who/rsc/thread/squint.pdf
19 num
, den
int64 // numerator, denominator
26 print(u
.num
, "/", u
.den
)
31 func (u rat
) eq(c rat
) bool {
32 return u
.num
== c
.num
&& u
.den
== c
.den
48 c
:= chnameserial
% len(chnames
)
51 d
.req
= make(chan int)
52 d
.dat
= make(chan rat
)
64 // split reads a single demand channel and replicates its
65 // output onto two, which may be read at different rates.
66 // A process is created at first demand for a rat and dies
67 // after the rat has been sent to both outputs.
69 // When multiple generations of split exist, the newest
70 // will service requests on one channel, which is
71 // always renamed to be out[0]; the oldest will service
72 // requests on the other channel, out[1]. All generations but the
73 // newest hold queued data that has already been sent to
74 // out[0]. When data has finally been sent to out[1],
75 // a signal on the release-wait channel tells the next newer
76 // generation to begin servicing out[1].
78 func dosplit(in
*dch
, out
*dch2
, wait
chan int ) {
79 both
:= false // do not service both channels
90 out
[0], out
[1] = out
[1], out
[0]
96 release
:= make(chan int)
97 go dosplit(in
, out
, release
)
108 func split(in
*dch
, out
*dch2
) {
109 release
:= make(chan int)
110 go dosplit(in
, out
, release
)
114 func put(dat rat
, out
*dch
) {
119 func get(in
*dch
) rat
{
125 // Get one rat from each of n demand channels
127 func getn(in
[]*dch
) []rat
{
129 if n
!= 2 { panic("bad n in getn") }
130 req
:= new([2] chan int)
131 dat
:= new([2] chan rat
)
132 out
:= make([]rat
, 2)
139 for n
=2*n
; n
>0; n
-- {
143 case req
[0] <- seqno
:
146 case req
[1] <- seqno
:
160 // Get one rat from each of 2 demand channels
162 func get2(in0
*dch
, in1
*dch
) []rat
{
163 return getn([]*dch
{in0
, in1
})
166 func copy(in
*dch
, out
*dch
) {
173 func repeat(dat rat
, out
*dch
) {
179 type PS
*dch
// power series
180 type PS2
*[2] PS
// pair of power series
194 // Upper-case for power series.
195 // Lower-case for rationals.
196 // Input variables: U,V,...
197 // Output variables: ...,Y,Z
199 // Integer gcd; needed for rational arithmetic
201 func gcd (u
, v
int64) int64 {
202 if u
< 0 { return gcd(-u
, v
) }
203 if u
== 0 { return v
}
207 // Make a rational from two ints and from one int
209 func i2tor(u
, v
int64) rat
{
222 func itor(u
int64) rat
{
230 // End mark and end test
234 func end(u rat
) int64 {
235 if u
.den
==0 { return 1 }
239 // Operations on rationals
241 func add(u
, v rat
) rat
{
242 g
:= gcd(u
.den
,v
.den
)
243 return i2tor(u
.num
*(v
.den
/g
)+v
.num
*(u
.den
/g
),u
.den
*(v
.den
/g
))
246 func mul(u
, v rat
) rat
{
247 g1
:= gcd(u
.num
,v
.den
)
248 g2
:= gcd(u
.den
,v
.num
)
250 r
.num
= (u
.num
/g1
)*(v
.num
/g2
)
251 r
.den
= (u
.den
/g2
)*(v
.den
/g1
)
255 func neg(u rat
) rat
{
256 return i2tor(-u
.num
, u
.den
)
259 func sub(u
, v rat
) rat
{
260 return add(u
, neg(v
))
263 func inv(u rat
) rat
{ // invert a rat
264 if u
.num
== 0 { panic("zero divide in inv") }
265 return i2tor(u
.den
, u
.num
)
268 // print eval in floating point of PS at x=c to n terms
269 func evaln(c rat
, U PS
, n
int) {
271 x
:= float64(c
.num
)/float64(c
.den
)
278 val
= val
+ x
* float64(u
.num
)/float64(u
.den
)
284 // Print n terms of a power series
285 func printn(U PS
, n
int) {
287 for ; !done
&& n
>0; n
-- {
298 // Evaluate n terms of power series U at x=c
299 func eval(c rat
, U PS
, n
int) rat
{
300 if n
==0 { return zero
}
302 if end(y
) != 0 { return zero
}
303 return add(y
,mul(c
,eval(c
,U
,n
-1)))
306 // Power-series constructors return channels on which power
307 // series flow. They start an encapsulated generator that
308 // puts the terms of the series on the channel.
310 // Make a pair of power series identical to a given power series
312 func Split(U PS
) *dch2
{
318 // Add two power series
319 func Add(U
, V PS
) PS
{
326 switch end(uv
[0])+2*end(uv
[1]) {
328 Z
.dat
<- add(uv
[0], uv
[1])
343 // Multiply a power series by a constant
344 func Cmul(c rat
,U PS
) PS
{
364 func Sub(U
, V PS
) PS
{
365 return Add(U
, Cmul(neg(one
), V
))
368 // Multiply a power series by the monomial x^n
370 func Monmul(U PS
, n
int) PS
{
373 for ; n
>0; n
-- { put(zero
,Z
) }
393 func Mon(c rat
, n
int) PS
{
397 for ; n
>0; n
=n
-1 { put(zero
,Z
) }
405 func Shift(c rat
, U PS
) PS
{
414 // simple pole at 1: 1/(1-x) = 1 1 1 1 1 ...
416 // Convert array of coefficients, constant term first
417 // to a (finite) power series
420 func Poly(a []rat) PS {
422 begin func(a []rat, Z PS) {
425 for j=len(a); !done&&j>0; j=j-1)
426 if(a[j-1].num!=0) done=1
428 for(; i<j; i=i+1) put(a[i],Z)
435 // Multiply. The algorithm is
438 // then UV = u*v + x*(u*VV+v*UU) + x*x*UU*VV
440 func Mul(U
, V PS
) PS
{
445 if end(uv
[0])!=0 ||
end(uv
[1]) != 0 {
448 Z
.dat
<- mul(uv
[0],uv
[1])
451 W
:= Add(Cmul(uv
[0],VV
[0]),Cmul(uv
[1],UU
[0]))
454 copy(Add(W
,Mul(UU
[1],VV
[1])),Z
)
469 for i
:=1; !done
; i
++ {
474 Z
.dat
<- mul(itor(int64(i
)),u
)
484 // Integrate, with const of integration
485 func Integ(c rat
,U PS
) PS
{
490 for i
:=1; !done
; i
++ {
493 if end(u
) != 0 { done
= true }
494 Z
.dat
<- mul(i2tor(1,int64(i
)),u
)
501 // Binomial theorem (1+x)^c
503 func Binom(c rat
) PS
{
510 t
= mul(mul(t
,c
),i2tor(1,int64(n
)))
519 // Reciprocal of a power series
522 // (u+x*UU)*(z+x*ZZ) = 1
524 // u*ZZ + z*UU +x*UU*ZZ = 0
525 // ZZ = -UU*(z+x*ZZ)/u
527 func Recip(U PS
) PS
{
534 split(Mul(Cmul(neg(z
),U
),Shift(z
,ZZ
[0])),ZZ
)
540 // Exponential of a power series with constant term 0
541 // (nonzero constant term would make nonrational coefficients)
542 // bug: the constant term is simply ignored
545 // integrate to get Z
549 split(Integ(one
,Mul(ZZ
[0],Diff(U
))),ZZ
)
553 // Substitute V for x in U, where the leading term of V is zero
556 // then S(U,V) = u + VV*S(V,UU)
557 // bug: a nonzero constant term is ignored
559 func Subst(U
, V PS
) PS
{
567 if end(get(VV
[0])) != 0 {
570 copy(Mul(VV
[0],Subst(U
,VV
[1])),Z
)
577 // Monomial Substition: U(c x^n)
578 // Each Ui is multiplied by c^i and followed by n-1 zeros
580 func MonSubst(U PS
, c0 rat
, n
int) PS
{
593 for i
:= 1; i
< n
; i
++ {
606 chnames
= "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz"
614 func check(U PS
, c rat
, count
int, str
string) {
615 for i
:= 0; i
< count
; i
++ {
629 func checka(U PS
, a
[]rat
, str
string) {
630 for i
:= 0; i
< N
; i
++ {
631 check(U
, a
[i
], 1, str
)
637 if len(os
.Args
) > 1 { // print
638 print("Ones: "); printn(Ones
, 10)
639 print("Twos: "); printn(Twos
, 10)
640 print("Add: "); printn(Add(Ones
, Twos
), 10)
641 print("Diff: "); printn(Diff(Ones
), 10)
642 print("Integ: "); printn(Integ(zero
, Ones
), 10)
643 print("CMul: "); printn(Cmul(neg(one
), Ones
), 10)
644 print("Sub: "); printn(Sub(Ones
, Twos
), 10)
645 print("Mul: "); printn(Mul(Ones
, Ones
), 10)
646 print("Exp: "); printn(Exp(Ones
), 15)
647 print("MonSubst: "); printn(MonSubst(Ones
, neg(one
), 2), 10)
648 print("ATan: "); printn(Integ(zero
, MonSubst(Ones
, neg(one
), 2)), 10)
650 check(Ones
, one
, 5, "Ones")
651 check(Add(Ones
, Ones
), itor(2), 0, "Add Ones Ones") // 1 1 1 1 1
652 check(Add(Ones
, Twos
), itor(3), 0, "Add Ones Twos") // 3 3 3 3 3
655 for i
:=0; i
< N
; i
++ {
656 a
[i
] = itor(int64(i
+1))
658 checka(d
, a
, "Diff") // 1 2 3 4 5
659 in
:= Integ(zero
, Ones
)
660 a
[0] = zero
// integration constant
661 for i
:=1; i
< N
; i
++ {
662 a
[i
] = i2tor(1, int64(i
))
664 checka(in
, a
, "Integ") // 0 1 1/2 1/3 1/4 1/5
665 check(Cmul(neg(one
), Twos
), itor(-2), 10, "CMul") // -1 -1 -1 -1 -1
666 check(Sub(Ones
, Twos
), itor(-1), 0, "Sub Ones Twos") // -1 -1 -1 -1 -1
668 for i
:=0; i
< N
; i
++ {
669 a
[i
] = itor(int64(i
+1))
671 checka(m
, a
, "Mul") // 1 2 3 4 5
679 a
[6] = i2tor(4051,720)
680 a
[7] = i2tor(37633,5040)
681 a
[8] = i2tor(43817,4480)
682 a
[9] = i2tor(4596553,362880)
683 checka(e
, a
, "Exp") // 1 1 3/2 13/6 73/24
684 at
:= Integ(zero
, MonSubst(Ones
, neg(one
), 2))
685 for c
, i
:= 1, 0; i
< N
; i
++ {
689 a
[i
] = i2tor(int64(c
), int64(i
))
693 checka(at
, a
, "ATan") // 0 -1 0 -1/3 0 -1/5
695 t := Revert(Integ(zero, MonSubst(Ones, neg(one), 2)))
705 a[9] = i2tor(62,2835)
706 checka(t, a, "Tan") // 0 1 0 1/3 0 2/15