URL https://opencores.org/ocsvn/openrisc/openrisc/trunk

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [gnu-dev/] [or1k-gcc/] [gcc/] [testsuite/] [go.test/] [test/] [chan/] [powser2.go] - Blame information for rev 858

Go to most recent revision | Details | Compare with Previous | View Log


// $G $D/$F.go && $L $F.$A && ./$A.out
 
// Copyright 2009 The Go Authors. All rights reserved.
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.
 
// Power series package
// A power series is a channel, along which flow rational
// coefficients.  A denominator of zero signifies the end.
// Original code in Newsqueak by Doug McIlroy.
// See Squinting at Power Series by Doug McIlroy,
//   http://www.cs.bell-labs.com/who/rsc/thread/squint.pdf
// Like powser1.go but uses channels of interfaces.
// Has not been cleaned up as much as powser1.go, to keep
// it distinct and therefore a different test.
 
package main
 
import "os"
 
type rat struct  {
        num, den  int64 // numerator, denominator
}
 
type item interface {
        pr()
        eq(c item) bool
}
 
func (u *rat) pr(){
        if u.den==1 {
                print(u.num)
        } else {
                print(u.num, "/", u.den)
        }
        print(" ")
}
 
func (u *rat) eq(c item) bool {
        c1 := c.(*rat)
        return u.num == c1.num && u.den == c1.den
}
 
type dch struct {
        req chan  int
        dat chan  item
        nam int
}
 
type dch2 [2] *dch
 
var chnames string
var chnameserial int
var seqno int
 
func mkdch() *dch {
        c := chnameserial % len(chnames)
        chnameserial++
        d := new(dch)
        d.req = make(chan int)
        d.dat = make(chan item)
        d.nam = c
        return d
}
 
func mkdch2() *dch2 {
        d2 := new(dch2)
        d2[0] = mkdch()
        d2[1] = mkdch()
        return d2
}
 
// split reads a single demand channel and replicates its
// output onto two, which may be read at different rates.
// A process is created at first demand for an item and dies
// after the item has been sent to both outputs.
 
// When multiple generations of split exist, the newest
// will service requests on one channel, which is
// always renamed to be out[0]; the oldest will service
// requests on the other channel, out[1].  All generations but the
// newest hold queued data that has already been sent to
// out[0].  When data has finally been sent to out[1],
// a signal on the release-wait channel tells the next newer
// generation to begin servicing out[1].
 
func dosplit(in *dch, out *dch2, wait chan int ){
        both := false   // do not service both channels
 
        select {
        case <-out[0].req:
 
        case <-wait:
                both = true
                select {
                case <-out[0].req:
 
                case <-out[1].req:
                        out[0],out[1] = out[1], out[0]
                }
        }
 
        seqno++
        in.req <- seqno
        release := make(chan  int)
        go dosplit(in, out, release)
        dat := <-in.dat
        out[0].dat <- dat
        if !both {
                <-wait
        }
        <-out[1].req
        out[1].dat <- dat
        release <- 0
}
 
func split(in *dch, out *dch2){
        release := make(chan int)
        go dosplit(in, out, release)
        release <- 0
}
 
func put(dat item, out *dch){
        <-out.req
        out.dat <- dat
}
 
func get(in *dch) *rat {
        seqno++
        in.req <- seqno
        return (<-in.dat).(*rat)
}
 
// Get one item from each of n demand channels
 
func getn(in []*dch) []item {
        n:=len(in)
        if n != 2 { panic("bad n in getn") }
        req := make([] chan int, 2)
        dat := make([] chan item, 2)
        out := make([]item, 2)
        var i int
        var it item
        for i=0; i
                req[i] = in[i].req
                dat[i] = nil
        }
        for n=2*n; n>0; n-- {
                seqno++
 
                select{
                case req[0] <- seqno:
                        dat[0] = in[0].dat
                        req[0] = nil
                case req[1] <- seqno:
                        dat[1] = in[1].dat
                        req[1] = nil
                case it = <-dat[0]:
                        out[0] = it
                        dat[0] = nil
                case it = <-dat[1]:
                        out[1] = it
                        dat[1] = nil
                }
        }
        return out
}
 
// Get one item from each of 2 demand channels
 
func get2(in0 *dch, in1 *dch)  []item {
        return getn([]*dch{in0, in1})
}
 
func copy(in *dch, out *dch){
        for {
                <-out.req
                out.dat <- get(in)
        }
}
 
func repeat(dat item, out *dch){
        for {
                put(dat, out)
        }
}
 
type PS *dch    // power series
type PS2 *[2] PS // pair of power series
 
var Ones PS
var Twos PS
 
func mkPS() *dch {
        return mkdch()
}
 
func mkPS2() *dch2 {
        return mkdch2()
}
 
// Conventions
// Upper-case for power series.
// Lower-case for rationals.
// Input variables: U,V,...
// Output variables: ...,Y,Z
 
// Integer gcd; needed for rational arithmetic
 
func gcd (u, v int64) int64{
        if u < 0 { return gcd(-u, v) }
        if u == 0 { return v }
        return gcd(v%u, u)
}
 
// Make a rational from two ints and from one int
 
func i2tor(u, v int64) *rat{
        g := gcd(u,v)
        r := new(rat)
        if v > 0 {
                r.num = u/g
                r.den = v/g
        } else {
                r.num = -u/g
                r.den = -v/g
        }
        return r
}
 
func itor(u int64) *rat{
        return i2tor(u, 1)
}
 
var zero *rat
var one *rat
 
 
// End mark and end test
 
var finis *rat
 
func end(u *rat) int64 {
        if u.den==0 { return 1 }
        return 0
}
 
// Operations on rationals
 
func add(u, v *rat) *rat {
        g := gcd(u.den,v.den)
        return  i2tor(u.num*(v.den/g)+v.num*(u.den/g),u.den*(v.den/g))
}
 
func mul(u, v *rat) *rat{
        g1 := gcd(u.num,v.den)
        g2 := gcd(u.den,v.num)
        r := new(rat)
        r.num =(u.num/g1)*(v.num/g2)
        r.den = (u.den/g2)*(v.den/g1)
        return r
}
 
func neg(u *rat) *rat{
        return i2tor(-u.num, u.den)
}
 
func sub(u, v *rat) *rat{
        return add(u, neg(v))
}
 
func inv(u *rat) *rat{  // invert a rat
        if u.num == 0 { panic("zero divide in inv") }
        return i2tor(u.den, u.num)
}
 
// print eval in floating point of PS at x=c to n terms
func Evaln(c *rat, U PS, n int) {
        xn := float64(1)
        x := float64(c.num)/float64(c.den)
        val := float64(0)
        for i:=0; i
                u := get(U)
                if end(u) != 0 {
                        break
                }
                val = val + x * float64(u.num)/float64(u.den)
                xn = xn*x
        }
        print(val, "\n")
}
 
// Print n terms of a power series
func Printn(U PS, n int){
        done := false
        for ; !done && n>0; n-- {
                u := get(U)
                if end(u) != 0 {
                        done = true
                } else {
                        u.pr()
                }
        }
        print(("\n"))
}
 
func Print(U PS){
        Printn(U,1000000000)
}
 
// Evaluate n terms of power series U at x=c
func eval(c *rat, U PS, n int) *rat{
        if n==0 { return zero }
        y := get(U)
        if end(y) != 0 { return zero }
        return add(y,mul(c,eval(c,U,n-1)))
}
 
// Power-series constructors return channels on which power
// series flow.  They start an encapsulated generator that
// puts the terms of the series on the channel.
 
// Make a pair of power series identical to a given power series
 
func Split(U PS) *dch2{
        UU := mkdch2()
        go split(U,UU)
        return UU
}
 
// Add two power series
func Add(U, V PS) PS{
        Z := mkPS()
        go func(U, V, Z PS){
                var uv [] item
                for {
                        <-Z.req
                        uv = get2(U,V)
                        switch end(uv[0].(*rat))+2*end(uv[1].(*rat)) {
                        case 0:
                                Z.dat <- add(uv[0].(*rat), uv[1].(*rat))
                        case 1:
                                Z.dat <- uv[1]
                                copy(V,Z)
                        case 2:
                                Z.dat <- uv[0]
                                copy(U,Z)
                        case 3:
                                Z.dat <- finis
                        }
                }
        }(U, V, Z)
        return Z
}
 
// Multiply a power series by a constant
func Cmul(c *rat,U PS) PS{
        Z := mkPS()
        go func(c *rat, U, Z PS){
                done := false
                for !done {
                        <-Z.req
                        u := get(U)
                        if end(u) != 0 {
                                done = true
                        } else {
                                Z.dat <- mul(c,u)
                        }
                }
                Z.dat <- finis
        }(c, U, Z)
        return Z
}
 
// Subtract
 
func Sub(U, V PS) PS{
        return Add(U, Cmul(neg(one), V))
}
 
// Multiply a power series by the monomial x^n
 
func Monmul(U PS, n int) PS{
        Z := mkPS()
        go func(n int, U PS, Z PS){
                for ; n>0; n-- { put(zero,Z) }
                copy(U,Z)
        }(n, U, Z)
        return Z
}
 
// Multiply by x
 
func Xmul(U PS) PS{
        return Monmul(U,1)
}
 
func Rep(c *rat) PS{
        Z := mkPS()
        go repeat(c,Z)
        return Z
}
 
// Monomial c*x^n
 
func Mon(c *rat, n int) PS{
        Z:=mkPS()
        go func(c *rat, n int, Z PS){
                if(c.num!=0) {
                        for ; n>0; n=n-1 { put(zero,Z) }
                        put(c,Z)
                }
                put(finis,Z)
        }(c, n, Z)
        return Z
}
 
func Shift(c *rat, U PS) PS{
        Z := mkPS()
        go func(c *rat, U, Z PS){
                put(c,Z)
                copy(U,Z)
        }(c, U, Z)
        return Z
}
 
// simple pole at 1: 1/(1-x) = 1 1 1 1 1 ...
 
// Convert array of coefficients, constant term first
// to a (finite) power series
 
/*
func Poly(a [] *rat) PS{
        Z:=mkPS()
        begin func(a [] *rat, Z PS){
                j:=0
                done:=0
                for j=len(a); !done&&j>0; j=j-1)
                        if(a[j-1].num!=0) done=1
                i:=0
                for(; i
                put(finis,Z)
        }()
        return Z
}
*/
 
// Multiply. The algorithm is
//      let U = u + x*UU
//      let V = v + x*VV
//      then UV = u*v + x*(u*VV+v*UU) + x*x*UU*VV
 
func Mul(U, V PS) PS{
        Z:=mkPS()
        go func(U, V, Z PS){
                <-Z.req
                uv := get2(U,V)
                if end(uv[0].(*rat))!=0 || end(uv[1].(*rat)) != 0 {
                        Z.dat <- finis
                } else {
                        Z.dat <- mul(uv[0].(*rat),uv[1].(*rat))
                        UU := Split(U)
                        VV := Split(V)
                        W := Add(Cmul(uv[0].(*rat),VV[0]),Cmul(uv[1].(*rat),UU[0]))
                        <-Z.req
                        Z.dat <- get(W)
                        copy(Add(W,Mul(UU[1],VV[1])),Z)
                }
        }(U, V, Z)
        return Z
}
 
// Differentiate
 
func Diff(U PS) PS{
        Z:=mkPS()
        go func(U, Z PS){
                <-Z.req
                u := get(U)
                if end(u) == 0 {
                        done:=false
                        for i:=1; !done; i++ {
                                u = get(U)
                                if end(u) != 0 {
                                        done=true
                                } else {
                                        Z.dat <- mul(itor(int64(i)),u)
                                        <-Z.req
                                }
                        }
                }
                Z.dat <- finis
        }(U, Z)
        return Z
}
 
// Integrate, with const of integration
func Integ(c *rat,U PS) PS{
        Z:=mkPS()
        go func(c *rat, U, Z PS){
                put(c,Z)
                done:=false
                for i:=1; !done; i++ {
                        <-Z.req
                        u := get(U)
                        if end(u) != 0 { done= true }
                        Z.dat <- mul(i2tor(1,int64(i)),u)
                }
                Z.dat <- finis
        }(c, U, Z)
        return Z
}
 
// Binomial theorem (1+x)^c
 
func Binom(c *rat) PS{
        Z:=mkPS()
        go func(c *rat, Z PS){
                n := 1
                t := itor(1)
                for c.num!=0 {
                        put(t,Z)
                        t = mul(mul(t,c),i2tor(1,int64(n)))
                        c = sub(c,one)
                        n++
                }
                put(finis,Z)
        }(c, Z)
        return Z
}
 
// Reciprocal of a power series
//      let U = u + x*UU
//      let Z = z + x*ZZ
//      (u+x*UU)*(z+x*ZZ) = 1
//      z = 1/u
//      u*ZZ + z*UU +x*UU*ZZ = 0
//      ZZ = -UU*(z+x*ZZ)/u
 
func Recip(U PS) PS{
        Z:=mkPS()
        go func(U, Z PS){
                ZZ:=mkPS2()
                <-Z.req
                z := inv(get(U))
                Z.dat <- z
                split(Mul(Cmul(neg(z),U),Shift(z,ZZ[0])),ZZ)
                copy(ZZ[1],Z)
        }(U, Z)
        return Z
}
 
// Exponential of a power series with constant term 0
// (nonzero constant term would make nonrational coefficients)
// bug: the constant term is simply ignored
//      Z = exp(U)
//      DZ = Z*DU
//      integrate to get Z
 
func Exp(U PS) PS{
        ZZ := mkPS2()
        split(Integ(one,Mul(ZZ[0],Diff(U))),ZZ)
        return ZZ[1]
}
 
// Substitute V for x in U, where the leading term of V is zero
//      let U = u + x*UU
//      let V = v + x*VV
//      then S(U,V) = u + VV*S(V,UU)
// bug: a nonzero constant term is ignored
 
func Subst(U, V PS) PS {
        Z:= mkPS()
        go func(U, V, Z PS) {
                VV := Split(V)
                <-Z.req
                u := get(U)
                Z.dat <- u
                if end(u) == 0 {
                        if end(get(VV[0])) != 0 {
                                put(finis,Z)
                        } else {
                                copy(Mul(VV[0],Subst(U,VV[1])),Z)
                        }
                }
        }(U, V, Z)
        return Z
}
 
// Monomial Substition: U(c x^n)
// Each Ui is multiplied by c^i and followed by n-1 zeros
 
func MonSubst(U PS, c0 *rat, n int) PS {
        Z:= mkPS()
        go func(U, Z PS, c0 *rat, n int) {
                c := one
                for {
                        <-Z.req
                        u := get(U)
                        Z.dat <- mul(u, c)
                        c = mul(c, c0)
                        if end(u) != 0 {
                                Z.dat <- finis
                                break
                        }
                        for i := 1; i < n; i++ {
                                <-Z.req
                                Z.dat <- zero
                        }
                }
        }(U, Z, c0, n)
        return Z
}
 
 
func Init() {
        chnameserial = -1
        seqno = 0
        chnames = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz"
        zero = itor(0)
        one = itor(1)
        finis = i2tor(1,0)
        Ones = Rep(one)
        Twos = Rep(itor(2))
}
 
func check(U PS, c *rat, count int, str string) {
        for i := 0; i < count; i++ {
                r := get(U)
                if !r.eq(c) {
                        print("got: ")
                        r.pr()
                        print("should get ")
                        c.pr()
                        print("\n")
                        panic(str)
                }
        }
}
 
const N=10
func checka(U PS, a []*rat, str string) {
        for i := 0; i < N; i++ {
                check(U, a[i], 1, str)
        }
}
 
func main() {
        Init()
        if len(os.Args) > 1 {  // print
                print("Ones: "); Printn(Ones, 10)
                print("Twos: "); Printn(Twos, 10)
                print("Add: "); Printn(Add(Ones, Twos), 10)
                print("Diff: "); Printn(Diff(Ones), 10)
                print("Integ: "); Printn(Integ(zero, Ones), 10)
                print("CMul: "); Printn(Cmul(neg(one), Ones), 10)
                print("Sub: "); Printn(Sub(Ones, Twos), 10)
                print("Mul: "); Printn(Mul(Ones, Ones), 10)
                print("Exp: "); Printn(Exp(Ones), 15)
                print("MonSubst: "); Printn(MonSubst(Ones, neg(one), 2), 10)
                print("ATan: "); Printn(Integ(zero, MonSubst(Ones, neg(one), 2)), 10)
        } else {  // test
                check(Ones, one, 5, "Ones")
                check(Add(Ones, Ones), itor(2), 0, "Add Ones Ones")  // 1 1 1 1 1
                check(Add(Ones, Twos), itor(3), 0, "Add Ones Twos") // 3 3 3 3 3
                a := make([]*rat, N)
                d := Diff(Ones)
                for i:=0; i < N; i++ {
                        a[i] = itor(int64(i+1))
                }
                checka(d, a, "Diff")  // 1 2 3 4 5
                in := Integ(zero, Ones)
                a[0] = zero  // integration constant
                for i:=1; i < N; i++ {
                        a[i] = i2tor(1, int64(i))
                }
                checka(in, a, "Integ")  // 0 1 1/2 1/3 1/4 1/5
                check(Cmul(neg(one), Twos), itor(-2), 10, "CMul")  // -1 -1 -1 -1 -1
                check(Sub(Ones, Twos), itor(-1), 0, "Sub Ones Twos")  // -1 -1 -1 -1 -1
                m := Mul(Ones, Ones)
                for i:=0; i < N; i++ {
                        a[i] = itor(int64(i+1))
                }
                checka(m, a, "Mul")  // 1 2 3 4 5
                e := Exp(Ones)
                a[0] = itor(1)
                a[1] = itor(1)
                a[2] = i2tor(3,2)
                a[3] = i2tor(13,6)
                a[4] = i2tor(73,24)
                a[5] = i2tor(167,40)
                a[6] = i2tor(4051,720)
                a[7] = i2tor(37633,5040)
                a[8] = i2tor(43817,4480)
                a[9] = i2tor(4596553,362880)
                checka(e, a, "Exp")  // 1 1 3/2 13/6 73/24
                at := Integ(zero, MonSubst(Ones, neg(one), 2))
                for c, i := 1, 0; i < N; i++ {
                        if i%2 == 0 {
                                a[i] = zero
                        } else {
                                a[i] = i2tor(int64(c), int64(i))
                                c *= -1
                        }
                }
                checka(at, a, "ATan");  // 0 -1 0 -1/3 0 -1/5
/*
                t := Revert(Integ(zero, MonSubst(Ones, neg(one), 2)))
                a[0] = zero
                a[1] = itor(1)
                a[2] = zero
                a[3] = i2tor(1,3)
                a[4] = zero
                a[5] = i2tor(2,15)
                a[6] = zero
                a[7] = i2tor(17,315)
                a[8] = zero
                a[9] = i2tor(62,2835)
                checka(t, a, "Tan")  // 0 1 0 1/3 0 2/15
*/
        }
}

Browse

Tools

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [gnu-dev/] [or1k-gcc/] [gcc/] [testsuite/] [go.test/] [test/] [chan/] [powser2.go] - Blame information for rev 858

Line No.	Rev	Author	Line
1	700	jeremybenn	`// $G $D/$F.go && $L $F.$A && ./$A.out`
2
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.`
6
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`
13			`// Like powser1.go but uses channels of interfaces.`
14			`// Has not been cleaned up as much as powser1.go, to keep`
15			`// it distinct and therefore a different test.`
16
17			`package main`
18
19			`import "os"`
20
21			`type rat struct {`
22			`num, den int64 // numerator, denominator`
23			`}`
24
25			`type item interface {`
26			`pr()`
27			`eq(c item) bool`
28			`}`
29
30			`func (u *rat) pr(){`
31			`if u.den==1 {`
32			`print(u.num)`
33			`} else {`
34			`print(u.num, "/", u.den)`
35			`}`
36			`print(" ")`
37			`}`
38
39			`func (u *rat) eq(c item) bool {`
40			`c1 := c.(*rat)`
41			`return u.num == c1.num && u.den == c1.den`
42			`}`
43
44			`type dch struct {`
45			`req chan int`
46			`dat chan item`
47			`nam int`
48			`}`
49
50			`type dch2 [2] *dch`
51
52			`var chnames string`
53			`var chnameserial int`
54			`var seqno int`
55
56			`func mkdch() *dch {`
57			`c := chnameserial % len(chnames)`
58			`chnameserial++`
59			`d := new(dch)`
60			`d.req = make(chan int)`
61			`d.dat = make(chan item)`
62			`d.nam = c`
63			`return d`
64			`}`
65
66			`func mkdch2() *dch2 {`
67			`d2 := new(dch2)`
68			`d2[0] = mkdch()`
69			`d2[1] = mkdch()`
70			`return d2`
71			`}`
72
73			`// split reads a single demand channel and replicates its`
74			`// output onto two, which may be read at different rates.`
75			`// A process is created at first demand for an item and dies`
76			`// after the item has been sent to both outputs.`
77
78			`// When multiple generations of split exist, the newest`
79			`// will service requests on one channel, which is`
80			`// always renamed to be out[0]; the oldest will service`
81			`// requests on the other channel, out[1]. All generations but the`
82			`// newest hold queued data that has already been sent to`
83			`// out[0]. When data has finally been sent to out[1],`
84			`// a signal on the release-wait channel tells the next newer`
85			`// generation to begin servicing out[1].`
86
87			`func dosplit(in dch, out dch2, wait chan int ){`
88			`both := false // do not service both channels`
89
90			`select {`
91			`case <-out[0].req:`
92
93			`case <-wait:`
94			`both = true`
95			`select {`
96			`case <-out[0].req:`
97
98			`case <-out[1].req:`
99			`out[0],out[1] = out[1], out[0]`
100			`}`
101			`}`
102
103			`seqno++`
104			`in.req <- seqno`
105			`release := make(chan int)`
106			`go dosplit(in, out, release)`
107			`dat := <-in.dat`
108			`out[0].dat <- dat`
109			`if !both {`
110			`<-wait`
111			`}`
112			`<-out[1].req`
113			`out[1].dat <- dat`
114			`release <- 0`
115			`}`
116
117			`func split(in dch, out dch2){`
118			`release := make(chan int)`
119			`go dosplit(in, out, release)`
120			`release <- 0`
121			`}`
122
123			`func put(dat item, out *dch){`
124			`<-out.req`
125			`out.dat <- dat`
126			`}`
127
128			`func get(in dch) rat {`
129			`seqno++`
130			`in.req <- seqno`
131			`return (<-in.dat).(*rat)`
132			`}`
133
134			`// Get one item from each of n demand channels`
135
136			`func getn(in []*dch) []item {`
137			`n:=len(in)`
138			`if n != 2 { panic("bad n in getn") }`
139			`req := make([] chan int, 2)`
140			`dat := make([] chan item, 2)`
141			`out := make([]item, 2)`
142			`var i int`
143			`var it item`
144			`for i=0; i`
145			`req[i] = in[i].req`
146			`dat[i] = nil`
147			`}`
148			`for n=2*n; n>0; n-- {`
149			`seqno++`
150
151			`select{`
152			`case req[0] <- seqno:`
153			`dat[0] = in[0].dat`
154			`req[0] = nil`
155			`case req[1] <- seqno:`
156			`dat[1] = in[1].dat`
157			`req[1] = nil`
158			`case it = <-dat[0]:`
159			`out[0] = it`
160			`dat[0] = nil`
161			`case it = <-dat[1]:`
162			`out[1] = it`
163			`dat[1] = nil`
164			`}`
165			`}`
166			`return out`
167			`}`
168
169			`// Get one item from each of 2 demand channels`
170
171			`func get2(in0 dch, in1 dch) []item {`
172			`return getn([]*dch{in0, in1})`
173			`}`
174
175			`func copy(in dch, out dch){`
176			`for {`
177			`<-out.req`
178			`out.dat <- get(in)`
179			`}`
180			`}`
181
182			`func repeat(dat item, out *dch){`
183			`for {`
184			`put(dat, out)`
185			`}`
186			`}`
187
188			`type PS *dch // power series`
189			`type PS2 *[2] PS // pair of power series`
190
191			`var Ones PS`
192			`var Twos PS`
193
194			`func mkPS() *dch {`
195			`return mkdch()`
196			`}`
197
198			`func mkPS2() *dch2 {`
199			`return mkdch2()`
200			`}`
201
202			`// Conventions`
203			`// Upper-case for power series.`
204			`// Lower-case for rationals.`
205			`// Input variables: U,V,...`
206			`// Output variables: ...,Y,Z`
207
208			`// Integer gcd; needed for rational arithmetic`
209
210			`func gcd (u, v int64) int64{`
211			`if u < 0 { return gcd(-u, v) }`
212			`if u == 0 { return v }`
213			`return gcd(v%u, u)`
214			`}`
215
216			`// Make a rational from two ints and from one int`
217
218			`func i2tor(u, v int64) *rat{`
219			`g := gcd(u,v)`
220			`r := new(rat)`
221			`if v > 0 {`
222			`r.num = u/g`
223			`r.den = v/g`
224			`} else {`
225			`r.num = -u/g`
226			`r.den = -v/g`
227			`}`
228			`return r`
229			`}`
230
231			`func itor(u int64) *rat{`
232			`return i2tor(u, 1)`
233			`}`
234
235			`var zero *rat`
236			`var one *rat`
237
238
239			`// End mark and end test`
240
241			`var finis *rat`
242
243			`func end(u *rat) int64 {`
244			`if u.den==0 { return 1 }`
245			`return 0`
246			`}`
247
248			`// Operations on rationals`
249
250			`func add(u, v rat) rat {`
251			`g := gcd(u.den,v.den)`
252			`return i2tor(u.num(v.den/g)+v.num(u.den/g),u.den*(v.den/g))`
253			`}`
254
255			`func mul(u, v rat) rat{`
256			`g1 := gcd(u.num,v.den)`
257			`g2 := gcd(u.den,v.num)`
258			`r := new(rat)`
259			`r.num =(u.num/g1)*(v.num/g2)`
260			`r.den = (u.den/g2)*(v.den/g1)`
261			`return r`
262			`}`
263
264			`func neg(u rat) rat{`
265			`return i2tor(-u.num, u.den)`
266			`}`
267
268			`func sub(u, v rat) rat{`
269			`return add(u, neg(v))`
270			`}`
271
272			`func inv(u rat) rat{ // invert a rat`
273			`if u.num == 0 { panic("zero divide in inv") }`
274			`return i2tor(u.den, u.num)`
275			`}`
276
277			`// print eval in floating point of PS at x=c to n terms`
278			`func Evaln(c *rat, U PS, n int) {`
279			`xn := float64(1)`
280			`x := float64(c.num)/float64(c.den)`
281			`val := float64(0)`
282			`for i:=0; i`
283			`u := get(U)`
284			`if end(u) != 0 {`
285			`break`
286			`}`
287			`val = val + x * float64(u.num)/float64(u.den)`
288			`xn = xn*x`
289			`}`
290			`print(val, "\n")`
291			`}`
292
293			`// Print n terms of a power series`
294			`func Printn(U PS, n int){`
295			`done := false`
296			`for ; !done && n>0; n-- {`
297			`u := get(U)`
298			`if end(u) != 0 {`
299			`done = true`
300			`} else {`
301			`u.pr()`
302			`}`
303			`}`
304			`print(("\n"))`
305			`}`
306
307			`func Print(U PS){`
308			`Printn(U,1000000000)`
309			`}`
310
311			`// Evaluate n terms of power series U at x=c`
312			`func eval(c rat, U PS, n int) rat{`
313			`if n==0 { return zero }`
314			`y := get(U)`
315			`if end(y) != 0 { return zero }`
316			`return add(y,mul(c,eval(c,U,n-1)))`
317			`}`
318
319			`// Power-series constructors return channels on which power`
320			`// series flow. They start an encapsulated generator that`
321			`// puts the terms of the series on the channel.`
322
323			`// Make a pair of power series identical to a given power series`
324
325			`func Split(U PS) *dch2{`
326			`UU := mkdch2()`
327			`go split(U,UU)`
328			`return UU`
329			`}`
330
331			`// Add two power series`
332			`func Add(U, V PS) PS{`
333			`Z := mkPS()`
334			`go func(U, V, Z PS){`
335			`var uv [] item`
336			`for {`
337			`<-Z.req`
338			`uv = get2(U,V)`
339			`switch end(uv[0].(rat))+2end(uv[1].(*rat)) {`
340			`case 0:`
341			`Z.dat <- add(uv[0].(rat), uv[1].(rat))`
342			`case 1:`
343			`Z.dat <- uv[1]`
344			`copy(V,Z)`
345			`case 2:`
346			`Z.dat <- uv[0]`
347			`copy(U,Z)`
348			`case 3:`
349			`Z.dat <- finis`
350			`}`
351			`}`
352			`}(U, V, Z)`
353			`return Z`
354			`}`
355
356			`// Multiply a power series by a constant`
357			`func Cmul(c *rat,U PS) PS{`
358			`Z := mkPS()`
359			`go func(c *rat, U, Z PS){`
360			`done := false`
361			`for !done {`
362			`<-Z.req`
363			`u := get(U)`
364			`if end(u) != 0 {`
365			`done = true`
366			`} else {`
367			`Z.dat <- mul(c,u)`
368			`}`
369			`}`
370			`Z.dat <- finis`
371			`}(c, U, Z)`
372			`return Z`
373			`}`
374
375			`// Subtract`
376
377			`func Sub(U, V PS) PS{`
378			`return Add(U, Cmul(neg(one), V))`
379			`}`
380
381			`// Multiply a power series by the monomial x^n`
382
383			`func Monmul(U PS, n int) PS{`
384			`Z := mkPS()`
385			`go func(n int, U PS, Z PS){`
386			`for ; n>0; n-- { put(zero,Z) }`
387			`copy(U,Z)`
388			`}(n, U, Z)`
389			`return Z`
390			`}`
391
392			`// Multiply by x`
393
394			`func Xmul(U PS) PS{`
395			`return Monmul(U,1)`
396			`}`
397
398			`func Rep(c *rat) PS{`
399			`Z := mkPS()`
400			`go repeat(c,Z)`
401			`return Z`
402			`}`
403
404			`// Monomial c*x^n`
405
406			`func Mon(c *rat, n int) PS{`
407			`Z:=mkPS()`
408			`go func(c *rat, n int, Z PS){`
409			`if(c.num!=0) {`
410			`for ; n>0; n=n-1 { put(zero,Z) }`
411			`put(c,Z)`
412			`}`
413			`put(finis,Z)`
414			`}(c, n, Z)`
415			`return Z`
416			`}`
417
418			`func Shift(c *rat, U PS) PS{`
419			`Z := mkPS()`
420			`go func(c *rat, U, Z PS){`
421			`put(c,Z)`
422			`copy(U,Z)`
423			`}(c, U, Z)`
424			`return Z`
425			`}`
426
427			`// simple pole at 1: 1/(1-x) = 1 1 1 1 1 ...`
428
429			`// Convert array of coefficients, constant term first`
430			`// to a (finite) power series`
431
432			`/*`
433			`func Poly(a [] *rat) PS{`
434			`Z:=mkPS()`
435			`begin func(a [] *rat, Z PS){`
436			`j:=0`
437			`done:=0`
438			`for j=len(a); !done&&j>0; j=j-1)`
439			`if(a[j-1].num!=0) done=1`
440			`i:=0`
441			`for(; i`
442			`put(finis,Z)`
443			`}()`
444			`return Z`
445			`}`
446			`*/`
447
448			`// Multiply. The algorithm is`
449			`// let U = u + x*UU`
450			`// let V = v + x*VV`
451			`// then UV = uv + x(uVV+vUU) + xxUU*VV`
452
453			`func Mul(U, V PS) PS{`
454			`Z:=mkPS()`
455			`go func(U, V, Z PS){`
456			`<-Z.req`
457			`uv := get2(U,V)`
458			`if end(uv[0].(rat))!=0 \|\| end(uv[1].(rat)) != 0 {`
459			`Z.dat <- finis`
460			`} else {`
461			`Z.dat <- mul(uv[0].(rat),uv[1].(rat))`
462			`UU := Split(U)`
463			`VV := Split(V)`
464			`W := Add(Cmul(uv[0].(rat),VV[0]),Cmul(uv[1].(rat),UU[0]))`
465			`<-Z.req`
466			`Z.dat <- get(W)`
467			`copy(Add(W,Mul(UU[1],VV[1])),Z)`
468			`}`
469			`}(U, V, Z)`
470			`return Z`
471			`}`
472
473			`// Differentiate`
474
475			`func Diff(U PS) PS{`
476			`Z:=mkPS()`
477			`go func(U, Z PS){`
478			`<-Z.req`
479			`u := get(U)`
480			`if end(u) == 0 {`
481			`done:=false`
482			`for i:=1; !done; i++ {`
483			`u = get(U)`
484			`if end(u) != 0 {`
485			`done=true`
486			`} else {`
487			`Z.dat <- mul(itor(int64(i)),u)`
488			`<-Z.req`
489			`}`
490			`}`
491			`}`
492			`Z.dat <- finis`
493			`}(U, Z)`
494			`return Z`
495			`}`
496
497			`// Integrate, with const of integration`
498			`func Integ(c *rat,U PS) PS{`
499			`Z:=mkPS()`
500			`go func(c *rat, U, Z PS){`
501			`put(c,Z)`
502			`done:=false`
503			`for i:=1; !done; i++ {`
504			`<-Z.req`
505			`u := get(U)`
506			`if end(u) != 0 { done= true }`
507			`Z.dat <- mul(i2tor(1,int64(i)),u)`
508			`}`
509			`Z.dat <- finis`
510			`}(c, U, Z)`
511			`return Z`
512			`}`
513
514			`// Binomial theorem (1+x)^c`
515
516			`func Binom(c *rat) PS{`
517			`Z:=mkPS()`
518			`go func(c *rat, Z PS){`
519			`n := 1`
520			`t := itor(1)`
521			`for c.num!=0 {`
522			`put(t,Z)`
523			`t = mul(mul(t,c),i2tor(1,int64(n)))`
524			`c = sub(c,one)`
525			`n++`
526			`}`
527			`put(finis,Z)`
528			`}(c, Z)`
529			`return Z`
530			`}`
531
532			`// Reciprocal of a power series`
533			`// let U = u + x*UU`
534			`// let Z = z + x*ZZ`
535			`// (u+xUU)(z+x*ZZ) = 1`
536			`// z = 1/u`
537			`// uZZ + zUU +xUUZZ = 0`
538			`// ZZ = -UU(z+xZZ)/u`
539
540			`func Recip(U PS) PS{`
541			`Z:=mkPS()`
542			`go func(U, Z PS){`
543			`ZZ:=mkPS2()`
544			`<-Z.req`
545			`z := inv(get(U))`
546			`Z.dat <- z`
547			`split(Mul(Cmul(neg(z),U),Shift(z,ZZ[0])),ZZ)`
548			`copy(ZZ[1],Z)`
549			`}(U, Z)`
550			`return Z`
551			`}`
552
553			`// Exponential of a power series with constant term 0`
554			`// (nonzero constant term would make nonrational coefficients)`
555			`// bug: the constant term is simply ignored`
556			`// Z = exp(U)`
557			`// DZ = Z*DU`
558			`// integrate to get Z`
559
560			`func Exp(U PS) PS{`
561			`ZZ := mkPS2()`
562			`split(Integ(one,Mul(ZZ[0],Diff(U))),ZZ)`
563			`return ZZ[1]`
564			`}`
565
566			`// Substitute V for x in U, where the leading term of V is zero`
567			`// let U = u + x*UU`
568			`// let V = v + x*VV`
569			`// then S(U,V) = u + VV*S(V,UU)`
570			`// bug: a nonzero constant term is ignored`
571
572			`func Subst(U, V PS) PS {`
573			`Z:= mkPS()`
574			`go func(U, V, Z PS) {`
575			`VV := Split(V)`
576			`<-Z.req`
577			`u := get(U)`
578			`Z.dat <- u`
579			`if end(u) == 0 {`
580			`if end(get(VV[0])) != 0 {`
581			`put(finis,Z)`
582			`} else {`
583			`copy(Mul(VV[0],Subst(U,VV[1])),Z)`
584			`}`
585			`}`
586			`}(U, V, Z)`
587			`return Z`
588			`}`
589
590			`// Monomial Substition: U(c x^n)`
591			`// Each Ui is multiplied by c^i and followed by n-1 zeros`
592
593			`func MonSubst(U PS, c0 *rat, n int) PS {`
594			`Z:= mkPS()`
595			`go func(U, Z PS, c0 *rat, n int) {`
596			`c := one`
597			`for {`
598			`<-Z.req`
599			`u := get(U)`
600			`Z.dat <- mul(u, c)`
601			`c = mul(c, c0)`
602			`if end(u) != 0 {`
603			`Z.dat <- finis`
604			`break`
605			`}`
606			`for i := 1; i < n; i++ {`
607			`<-Z.req`
608			`Z.dat <- zero`
609			`}`
610			`}`
611			`}(U, Z, c0, n)`
612			`return Z`
613			`}`
614
615
616			`func Init() {`
617			`chnameserial = -1`
618			`seqno = 0`
619			`chnames = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz"`
620			`zero = itor(0)`
621			`one = itor(1)`
622			`finis = i2tor(1,0)`
623			`Ones = Rep(one)`
624			`Twos = Rep(itor(2))`
625			`}`
626
627			`func check(U PS, c *rat, count int, str string) {`
628			`for i := 0; i < count; i++ {`
629			`r := get(U)`
630			`if !r.eq(c) {`
631			`print("got: ")`
632			`r.pr()`
633			`print("should get ")`
634			`c.pr()`
635			`print("\n")`
636			`panic(str)`
637			`}`
638			`}`
639			`}`
640
641			`const N=10`
642			`func checka(U PS, a []*rat, str string) {`
643			`for i := 0; i < N; i++ {`
644			`check(U, a[i], 1, str)`
645			`}`
646			`}`
647
648			`func main() {`
649			`Init()`
650			`if len(os.Args) > 1 { // print`
651			`print("Ones: "); Printn(Ones, 10)`
652			`print("Twos: "); Printn(Twos, 10)`
653			`print("Add: "); Printn(Add(Ones, Twos), 10)`
654			`print("Diff: "); Printn(Diff(Ones), 10)`
655			`print("Integ: "); Printn(Integ(zero, Ones), 10)`
656			`print("CMul: "); Printn(Cmul(neg(one), Ones), 10)`
657			`print("Sub: "); Printn(Sub(Ones, Twos), 10)`
658			`print("Mul: "); Printn(Mul(Ones, Ones), 10)`
659			`print("Exp: "); Printn(Exp(Ones), 15)`
660			`print("MonSubst: "); Printn(MonSubst(Ones, neg(one), 2), 10)`
661			`print("ATan: "); Printn(Integ(zero, MonSubst(Ones, neg(one), 2)), 10)`
662			`} else { // test`
663			`check(Ones, one, 5, "Ones")`
664			`check(Add(Ones, Ones), itor(2), 0, "Add Ones Ones") // 1 1 1 1 1`
665			`check(Add(Ones, Twos), itor(3), 0, "Add Ones Twos") // 3 3 3 3 3`
666			`a := make([]*rat, N)`
667			`d := Diff(Ones)`
668			`for i:=0; i < N; i++ {`
669			`a[i] = itor(int64(i+1))`
670			`}`
671			`checka(d, a, "Diff") // 1 2 3 4 5`
672			`in := Integ(zero, Ones)`
673			`a[0] = zero // integration constant`
674			`for i:=1; i < N; i++ {`
675			`a[i] = i2tor(1, int64(i))`
676			`}`
677			`checka(in, a, "Integ") // 0 1 1/2 1/3 1/4 1/5`
678			`check(Cmul(neg(one), Twos), itor(-2), 10, "CMul") // -1 -1 -1 -1 -1`
679			`check(Sub(Ones, Twos), itor(-1), 0, "Sub Ones Twos") // -1 -1 -1 -1 -1`
680			`m := Mul(Ones, Ones)`
681			`for i:=0; i < N; i++ {`
682			`a[i] = itor(int64(i+1))`
683			`}`
684			`checka(m, a, "Mul") // 1 2 3 4 5`
685			`e := Exp(Ones)`
686			`a[0] = itor(1)`
687			`a[1] = itor(1)`
688			`a[2] = i2tor(3,2)`
689			`a[3] = i2tor(13,6)`
690			`a[4] = i2tor(73,24)`
691			`a[5] = i2tor(167,40)`
692			`a[6] = i2tor(4051,720)`
693			`a[7] = i2tor(37633,5040)`
694			`a[8] = i2tor(43817,4480)`
695			`a[9] = i2tor(4596553,362880)`
696			`checka(e, a, "Exp") // 1 1 3/2 13/6 73/24`
697			`at := Integ(zero, MonSubst(Ones, neg(one), 2))`
698			`for c, i := 1, 0; i < N; i++ {`
699			`if i%2 == 0 {`
700			`a[i] = zero`
701			`} else {`
702			`a[i] = i2tor(int64(c), int64(i))`
703			`c *= -1`
704			`}`
705			`}`
706			`checka(at, a, "ATan"); // 0 -1 0 -1/3 0 -1/5`
707			`/*`
708			`t := Revert(Integ(zero, MonSubst(Ones, neg(one), 2)))`
709			`a[0] = zero`
710			`a[1] = itor(1)`
711			`a[2] = zero`
712			`a[3] = i2tor(1,3)`
713			`a[4] = zero`
714			`a[5] = i2tor(2,15)`
715			`a[6] = zero`
716			`a[7] = i2tor(17,315)`
717			`a[8] = zero`
718			`a[9] = i2tor(62,2835)`
719			`checka(t, a, "Tan") // 0 1 0 1/3 0 2/15`
720			`*/`
721			`}`
722			`}`