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
16package main
17
18import "os"
19
20type rat struct  {
21	num, den  int64	// numerator, denominator
22}
23
24func (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
33func (u rat) eq(c rat) bool {
34	return u.num == c.num && u.den == c.den
35}
36
37type dch struct {
38	req chan  int
39	dat chan  rat
40	nam int
41}
42
43type dch2 [2] *dch
44
45var chnames string
46var chnameserial int
47var seqno int
48
49func 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
59func 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
80func 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
110func split(in *dch, out *dch2) {
111	release := make(chan int)
112	go dosplit(in, out, release)
113	release <- 0
114}
115
116func put(dat rat, out *dch) {
117	<-out.req
118	out.dat <- dat
119}
120
121func 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
129func 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
164func get2(in0 *dch, in1 *dch) []rat {
165	return getn([]*dch{in0, in1})
166}
167
168func copy(in *dch, out *dch) {
169	for {
170		<-out.req
171		out.dat <- get(in)
172	}
173}
174
175func repeat(dat rat, out *dch) {
176	for {
177		put(dat, out)
178	}
179}
180
181type PS *dch	// power series
182type PS2 *[2] PS // pair of power series
183
184var Ones PS
185var Twos PS
186
187func mkPS() *dch {
188	return mkdch()
189}
190
191func 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
203func 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
211func 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
224func itor(u int64) rat {
225	return i2tor(u, 1)
226}
227
228var zero rat
229var one rat
230
231
232// End mark and end test
233
234var finis rat
235
236func end(u rat) int64 {
237	if u.den==0 { return 1 }
238	return 0
239}
240
241// Operations on rationals
242
243func 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
248func 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
257func neg(u rat) rat {
258	return i2tor(-u.num, u.den)
259}
260
261func sub(u, v rat) rat {
262	return add(u, neg(v))
263}
264
265func 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
271func 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
287func 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
301func 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
314func Split(U PS) *dch2 {
315	UU := mkdch2()
316	go split(U,UU)
317	return UU
318}
319
320// Add two power series
321func 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
346func 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
366func 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
372func 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
383func Xmul(U PS) PS {
384	return Monmul(U,1)
385}
386
387func Rep(c rat) PS {
388	Z := mkPS()
389	go repeat(c,Z)
390	return Z
391}
392
393// Monomial c*x^n
394
395func 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
407func 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/*
422func 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
442func 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
464func 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
487func 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
505func 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
529func 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
549func 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
561func 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
582func 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
605func 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
616func 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
630const N=10
631func checka(U PS, a []rat, str string) {
632	for i := 0; i < N; i++ {
633		check(U, a[i], 1, str)
634	}
635}
636
637func 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