test/chan/powser2.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 // 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
  20 package main
  21
  22 import "os"
  23
  24 type rat struct {
  25         num, den int64 // numerator, denominator
  26 }
  27
  28 type item interface {
  29         pr()
  30         eq(c item) bool
  31 }
  32
  33 func (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
  42 func (u *rat) eq(c item) bool {
  43         c1 := c.(*rat)
  44         return u.num == c1.num && u.den == c1.den
  45 }
  46
  47 type dch struct {
  48         req chan int
  49         dat chan item
  50         nam int
  51 }
  52
  53 type dch2 [2]*dch
  54
  55 var chnames string
  56 var chnameserial int
  57 var seqno int
  58
  59 func 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
  69 func 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
  90 func 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
 120 func split(in *dch, out *dch2) {
 121         release := make(chan int)
 122         go dosplit(in, out, release)
 123         release <- 0
 124 }
 125
 126 func put(dat item, out *dch) {
 127         <-out.req
 128         out.dat <- dat
 129 }
 130
 131 func 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
 139 func 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
 176 func get2(in0 *dch, in1 *dch) []item {
 177         return getn([]*dch{in0, in1})
 178 }
 179
 180 func copy(in *dch, out *dch) {
 181         for {
 182                 <-out.req
 183                 out.dat <- get(in)
 184         }
 185 }
 186
 187 func repeat(dat item, out *dch) {
 188         for {
 189                 put(dat, out)
 190         }
 191 }
 192
 193 type PS *dch    // power series
 194 type PS2 *[2]PS // pair of power series
 195
 196 var Ones PS
 197 var Twos PS
 198
 199 func mkPS() *dch {
 200         return mkdch()
 201 }
 202
 203 func 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
 215 func 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
 227 func 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
 240 func itor(u int64) *rat {
 241         return i2tor(u, 1)
 242 }
 243
 244 var zero *rat
 245 var one *rat
 246
 247 // End mark and end test
 248
 249 var finis *rat
 250
 251 func end(u *rat) int64 {
 252         if u.den == 0 {
 253                 return 1
 254         }
 255         return 0
 256 }
 257
 258 // Operations on rationals
 259
 260 func 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
 265 func 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
 274 func neg(u *rat) *rat {
 275         return i2tor(-u.num, u.den)
 276 }
 277
 278 func sub(u, v *rat) *rat {
 279         return add(u, neg(v))
 280 }
 281
 282 func 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
 290 func 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
 306 func 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
 319 func Print(U PS) {
 320         Printn(U, 1000000000)
 321 }
 322
 323 // Evaluate n terms of power series U at x=c
 324 func 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
 341 func Split(U PS) *dch2 {
 342         UU := mkdch2()
 343         go split(U, UU)
 344         return UU
 345 }
 346
 347 // Add two power series
 348 func 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
 373 func 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
 393 func 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
 399 func 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
 412 func Xmul(U PS) PS {
 413         return Monmul(U, 1)
 414 }
 415
 416 func Rep(c *rat) PS {
 417         Z := mkPS()
 418         go repeat(c, Z)
 419         return Z
 420 }
 421
 422 // Monomial c*x^n
 423
 424 func 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
 438 func 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 /*
 453 func 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
 473 func 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
 495 func 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
 518 func 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
 538 func 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
 562 func 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
 582 func 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
 594 func 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
 615 func 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
 637 func 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
 648 func 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
 662 const N = 10
 663
 664 func checka(U PS, a []*rat, str string) {
 665         for i := 0; i < N; i++ {
 666                 check(U, a[i], 1, str)
 667         }
 668 }
 669
 670 func 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 }