Source file test/chan/powser1.go

```     1  // run
2
4  // Use of this source code is governed by a BSD-style
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,
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 {
133  	}
134  	req := new([2]chan int)
135  	dat := new([2]chan rat)
136  	out := make([]rat, 2)
137  	var i int
138  	var it rat
139  	for i = 0; i < n; i++ {
140  		req[i] = in[i].req
141  		dat[i] = nil
142  	}
143  	for n = 2 * n; n > 0; n-- {
144  		seqno++
145
146  		select {
147  		case req[0] <- seqno:
148  			dat[0] = in[0].dat
149  			req[0] = nil
150  		case req[1] <- seqno:
151  			dat[1] = in[1].dat
152  			req[1] = nil
153  		case it = <-dat[0]:
154  			out[0] = it
155  			dat[0] = nil
156  		case it = <-dat[1]:
157  			out[1] = it
158  			dat[1] = nil
159  		}
160  	}
161  	return out
162  }
163
164  // Get one rat from each of 2 demand channels
165
166  func get2(in0 *dch, in1 *dch) []rat {
167  	return getn([]*dch{in0, in1})
168  }
169
170  func copy(in *dch, out *dch) {
171  	for {
172  		<-out.req
173  		out.dat <- get(in)
174  	}
175  }
176
177  func repeat(dat rat, out *dch) {
178  	for {
179  		put(dat, out)
180  	}
181  }
182
183  type PS *dch    // power series
184  type PS2 *[2]PS // pair of power series
185
186  var Ones PS
187  var Twos PS
188
189  func mkPS() *dch {
190  	return mkdch()
191  }
192
193  func mkPS2() *dch2 {
194  	return mkdch2()
195  }
196
197  // Conventions
198  // Upper-case for power series.
199  // Lower-case for rationals.
200  // Input variables: U,V,...
201  // Output variables: ...,Y,Z
202
203  // Integer gcd; needed for rational arithmetic
204
205  func gcd(u, v int64) int64 {
206  	if u < 0 {
207  		return gcd(-u, v)
208  	}
209  	if u == 0 {
210  		return v
211  	}
212  	return gcd(v%u, u)
213  }
214
215  // Make a rational from two ints and from one int
216
217  func i2tor(u, v int64) rat {
218  	g := gcd(u, v)
219  	var r rat
220  	if v > 0 {
221  		r.num = u / g
222  		r.den = v / g
223  	} else {
224  		r.num = -u / g
225  		r.den = -v / g
226  	}
227  	return r
228  }
229
230  func itor(u int64) rat {
231  	return i2tor(u, 1)
232  }
233
234  var zero rat
235  var one rat
236
237  // End mark and end test
238
239  var finis rat
240
241  func end(u rat) int64 {
242  	if u.den == 0 {
243  		return 1
244  	}
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  	var r 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 {
270  }
271
272  func inv(u rat) rat { // invert a rat
273  	if u.num == 0 {
274  		panic("zero divide in inv")
275  	}
276  	return i2tor(u.den, u.num)
277  }
278
279  // print eval in floating point of PS at x=c to n terms
280  func evaln(c rat, U PS, n int) {
281  	xn := float64(1)
282  	x := float64(c.num) / float64(c.den)
283  	val := float64(0)
284  	for i := 0; i < n; i++ {
285  		u := get(U)
286  		if end(u) != 0 {
287  			break
288  		}
289  		val = val + x*float64(u.num)/float64(u.den)
290  		xn = xn * x
291  	}
292  	print(val, "\n")
293  }
294
295  // Print n terms of a power series
296  func printn(U PS, n int) {
297  	done := false
298  	for ; !done && n > 0; n-- {
299  		u := get(U)
300  		if end(u) != 0 {
301  			done = true
302  		} else {
303  			u.pr()
304  		}
305  	}
306  	print(("\n"))
307  }
308
309  // Evaluate n terms of power series U at x=c
310  func eval(c rat, U PS, n int) rat {
311  	if n == 0 {
312  		return zero
313  	}
314  	y := get(U)
315  	if end(y) != 0 {
316  		return zero
317  	}
318  	return add(y, mul(c, eval(c, U, n-1)))
319  }
320
321  // Power-series constructors return channels on which power
322  // series flow.  They start an encapsulated generator that
323  // puts the terms of the series on the channel.
324
325  // Make a pair of power series identical to a given power series
326
327  func Split(U PS) *dch2 {
328  	UU := mkdch2()
329  	go split(U, UU)
330  	return UU
331  }
332
333  // Add two power series
334  func Add(U, V PS) PS {
335  	Z := mkPS()
336  	go func() {
337  		var uv []rat
338  		for {
339  			<-Z.req
340  			uv = get2(U, V)
341  			switch end(uv[0]) + 2*end(uv[1]) {
342  			case 0:
344  			case 1:
345  				Z.dat <- uv[1]
346  				copy(V, Z)
347  			case 2:
348  				Z.dat <- uv[0]
349  				copy(U, Z)
350  			case 3:
351  				Z.dat <- finis
352  			}
353  		}
354  	}()
355  	return Z
356  }
357
358  // Multiply a power series by a constant
359  func Cmul(c rat, U PS) PS {
360  	Z := mkPS()
361  	go func() {
362  		done := false
363  		for !done {
364  			<-Z.req
365  			u := get(U)
366  			if end(u) != 0 {
367  				done = true
368  			} else {
369  				Z.dat <- mul(c, u)
370  			}
371  		}
372  		Z.dat <- finis
373  	}()
374  	return Z
375  }
376
377  // Subtract
378
379  func Sub(U, V PS) PS {
381  }
382
383  // Multiply a power series by the monomial x^n
384
385  func Monmul(U PS, n int) PS {
386  	Z := mkPS()
387  	go func() {
388  		for ; n > 0; n-- {
389  			put(zero, Z)
390  		}
391  		copy(U, Z)
392  	}()
393  	return Z
394  }
395
396  // Multiply by x
397
398  func Xmul(U PS) PS {
399  	return Monmul(U, 1)
400  }
401
402  func Rep(c rat) PS {
403  	Z := mkPS()
404  	go repeat(c, Z)
405  	return Z
406  }
407
408  // Monomial c*x^n
409
410  func Mon(c rat, n int) PS {
411  	Z := mkPS()
412  	go func() {
413  		if c.num != 0 {
414  			for ; n > 0; n = n - 1 {
415  				put(zero, Z)
416  			}
417  			put(c, Z)
418  		}
419  		put(finis, Z)
420  	}()
421  	return Z
422  }
423
424  func Shift(c rat, U PS) PS {
425  	Z := mkPS()
426  	go func() {
427  		put(c, Z)
428  		copy(U, Z)
429  	}()
430  	return Z
431  }
432
433  // simple pole at 1: 1/(1-x) = 1 1 1 1 1 ...
434
435  // Convert array of coefficients, constant term first
436  // to a (finite) power series
437
438  /*
439  func Poly(a []rat) PS {
440  	Z:=mkPS()
441  	begin func(a []rat, Z PS) {
442  		j:=0
443  		done:=0
444  		for j=len(a); !done&&j>0; j=j-1)
445  			if(a[j-1].num!=0) done=1
446  		i:=0
447  		for(; i<j; i=i+1) put(a[i],Z)
448  		put(finis,Z)
449  	}()
450  	return Z
451  }
452  */
453
454  // Multiply. The algorithm is
455  //	let U = u + x*UU
456  //	let V = v + x*VV
457  //	then UV = u*v + x*(u*VV+v*UU) + x*x*UU*VV
458
459  func Mul(U, V PS) PS {
460  	Z := mkPS()
461  	go func() {
462  		<-Z.req
463  		uv := get2(U, V)
464  		if end(uv[0]) != 0 || end(uv[1]) != 0 {
465  			Z.dat <- finis
466  		} else {
467  			Z.dat <- mul(uv[0], uv[1])
468  			UU := Split(U)
469  			VV := Split(V)
470  			W := Add(Cmul(uv[0], VV[0]), Cmul(uv[1], UU[0]))
471  			<-Z.req
472  			Z.dat <- get(W)
474  		}
475  	}()
476  	return Z
477  }
478
479  // Differentiate
480
481  func Diff(U PS) PS {
482  	Z := mkPS()
483  	go func() {
484  		<-Z.req
485  		u := get(U)
486  		if end(u) == 0 {
487  			done := false
488  			for i := 1; !done; i++ {
489  				u = get(U)
490  				if end(u) != 0 {
491  					done = true
492  				} else {
493  					Z.dat <- mul(itor(int64(i)), u)
494  					<-Z.req
495  				}
496  			}
497  		}
498  		Z.dat <- finis
499  	}()
500  	return Z
501  }
502
503  // Integrate, with const of integration
504  func Integ(c rat, U PS) PS {
505  	Z := mkPS()
506  	go func() {
507  		put(c, Z)
508  		done := false
509  		for i := 1; !done; i++ {
510  			<-Z.req
511  			u := get(U)
512  			if end(u) != 0 {
513  				done = true
514  			}
515  			Z.dat <- mul(i2tor(1, int64(i)), u)
516  		}
517  		Z.dat <- finis
518  	}()
519  	return Z
520  }
521
522  // Binomial theorem (1+x)^c
523
524  func Binom(c rat) PS {
525  	Z := mkPS()
526  	go func() {
527  		n := 1
528  		t := itor(1)
529  		for c.num != 0 {
530  			put(t, Z)
531  			t = mul(mul(t, c), i2tor(1, int64(n)))
532  			c = sub(c, one)
533  			n++
534  		}
535  		put(finis, Z)
536  	}()
537  	return Z
538  }
539
540  // Reciprocal of a power series
541  //	let U = u + x*UU
542  //	let Z = z + x*ZZ
543  //	(u+x*UU)*(z+x*ZZ) = 1
544  //	z = 1/u
545  //	u*ZZ + z*UU +x*UU*ZZ = 0
546  //	ZZ = -UU*(z+x*ZZ)/u
547
548  func Recip(U PS) PS {
549  	Z := mkPS()
550  	go func() {
551  		ZZ := mkPS2()
552  		<-Z.req
553  		z := inv(get(U))
554  		Z.dat <- z
555  		split(Mul(Cmul(neg(z), U), Shift(z, ZZ[0])), ZZ)
556  		copy(ZZ[1], Z)
557  	}()
558  	return Z
559  }
560
561  // Exponential of a power series with constant term 0
562  // (nonzero constant term would make nonrational coefficients)
563  // bug: the constant term is simply ignored
564  //	Z = exp(U)
565  //	DZ = Z*DU
566  //	integrate to get Z
567
568  func Exp(U PS) PS {
569  	ZZ := mkPS2()
570  	split(Integ(one, Mul(ZZ[0], Diff(U))), ZZ)
571  	return ZZ[1]
572  }
573
574  // Substitute V for x in U, where the leading term of V is zero
575  //	let U = u + x*UU
576  //	let V = v + x*VV
577  //	then S(U,V) = u + VV*S(V,UU)
578  // bug: a nonzero constant term is ignored
579
580  func Subst(U, V PS) PS {
581  	Z := mkPS()
582  	go func() {
583  		VV := Split(V)
584  		<-Z.req
585  		u := get(U)
586  		Z.dat <- u
587  		if end(u) == 0 {
588  			if end(get(VV[0])) != 0 {
589  				put(finis, Z)
590  			} else {
591  				copy(Mul(VV[0], Subst(U, VV[1])), Z)
592  			}
593  		}
594  	}()
595  	return Z
596  }
597
598  // Monomial Substitution: U(c x^n)
599  // Each Ui is multiplied by c^i and followed by n-1 zeros
600
601  func MonSubst(U PS, c0 rat, n int) PS {
602  	Z := mkPS()
603  	go func() {
604  		c := one
605  		for {
606  			<-Z.req
607  			u := get(U)
608  			Z.dat <- mul(u, c)
609  			c = mul(c, c0)
610  			if end(u) != 0 {
611  				Z.dat <- finis
612  				break
613  			}
614  			for i := 1; i < n; i++ {
615  				<-Z.req
616  				Z.dat <- zero
617  			}
618  		}
619  	}()
620  	return Z
621  }
622
623  func Init() {
624  	chnameserial = -1
625  	seqno = 0
626  	chnames = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz"
627  	zero = itor(0)
628  	one = itor(1)
629  	finis = i2tor(1, 0)
630  	Ones = Rep(one)
631  	Twos = Rep(itor(2))
632  }
633
634  func check(U PS, c rat, count int, str string) {
635  	for i := 0; i < count; i++ {
636  		r := get(U)
637  		if !r.eq(c) {
638  			print("got: ")
639  			r.pr()
640  			print("should get ")
641  			c.pr()
642  			print("\n")
643  			panic(str)
644  		}
645  	}
646  }
647
648  const N = 10
649
650  func checka(U PS, a []rat, str string) {
651  	for i := 0; i < N; i++ {
652  		check(U, a[i], 1, str)
653  	}
654  }
655
656  func main() {
657  	Init()
658  	if len(os.Args) > 1 { // print
659  		print("Ones: ")
660  		printn(Ones, 10)
661  		print("Twos: ")
662  		printn(Twos, 10)
665  		print("Diff: ")
666  		printn(Diff(Ones), 10)
667  		print("Integ: ")
668  		printn(Integ(zero, Ones), 10)
669  		print("CMul: ")
670  		printn(Cmul(neg(one), Ones), 10)
671  		print("Sub: ")
672  		printn(Sub(Ones, Twos), 10)
673  		print("Mul: ")
674  		printn(Mul(Ones, Ones), 10)
675  		print("Exp: ")
676  		printn(Exp(Ones), 15)
677  		print("MonSubst: ")
678  		printn(MonSubst(Ones, neg(one), 2), 10)
679  		print("ATan: ")
680  		printn(Integ(zero, MonSubst(Ones, neg(one), 2)), 10)
681  	} else { // test
682  		check(Ones, one, 5, "Ones")
683  		check(Add(Ones, Ones), itor(2), 0, "Add Ones Ones") // 1 1 1 1 1
684  		check(Add(Ones, Twos), itor(3), 0, "Add Ones Twos") // 3 3 3 3 3
685  		a := make([]rat, N)
686  		d := Diff(Ones)
687  		for i := 0; i < N; i++ {
688  			a[i] = itor(int64(i + 1))
689  		}
690  		checka(d, a, "Diff") // 1 2 3 4 5
691  		in := Integ(zero, Ones)
692  		a[0] = zero // integration constant
693  		for i := 1; i < N; i++ {
694  			a[i] = i2tor(1, int64(i))
695  		}
696  		checka(in, a, "Integ")                               // 0 1 1/2 1/3 1/4 1/5
697  		check(Cmul(neg(one), Twos), itor(-2), 10, "CMul")    // -1 -1 -1 -1 -1
698  		check(Sub(Ones, Twos), itor(-1), 0, "Sub Ones Twos") // -1 -1 -1 -1 -1
699  		m := Mul(Ones, Ones)
700  		for i := 0; i < N; i++ {
701  			a[i] = itor(int64(i + 1))
702  		}
703  		checka(m, a, "Mul") // 1 2 3 4 5
704  		e := Exp(Ones)
705  		a[0] = itor(1)
706  		a[1] = itor(1)
707  		a[2] = i2tor(3, 2)
708  		a[3] = i2tor(13, 6)
709  		a[4] = i2tor(73, 24)
710  		a[5] = i2tor(167, 40)
711  		a[6] = i2tor(4051, 720)
712  		a[7] = i2tor(37633, 5040)
713  		a[8] = i2tor(43817, 4480)
714  		a[9] = i2tor(4596553, 362880)
715  		checka(e, a, "Exp") // 1 1 3/2 13/6 73/24
716  		at := Integ(zero, MonSubst(Ones, neg(one), 2))
717  		for c, i := 1, 0; i < N; i++ {
718  			if i%2 == 0 {
719  				a[i] = zero
720  			} else {
721  				a[i] = i2tor(int64(c), int64(i))
722  				c *= -1
723  			}
724  		}
725  		checka(at, a, "ATan") // 0 -1 0 -1/3 0 -1/5
726  		/*
727  			t := Revert(Integ(zero, MonSubst(Ones, neg(one), 2)))
728  			a[0] = zero
729  			a[1] = itor(1)
730  			a[2] = zero
731  			a[3] = i2tor(1,3)
732  			a[4] = zero
733  			a[5] = i2tor(2,15)
734  			a[6] = zero
735  			a[7] = i2tor(17,315)
736  			a[8] = zero
737  			a[9] = i2tor(62,2835)
738  			checka(t, a, "Tan")  // 0 1 0 1/3 0 2/15
739  		*/
740  	}
741  }
742
```

View as plain text