...
Run Format

Source file test/chan/powser1.go

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

View as plain text