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