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