File: waveout.go 1 /* 2 The MIT License (MIT) 3 4 Copyright © 2020-2025 pacman64 5 6 Permission is hereby granted, free of charge, to any person obtaining a copy of 7 this software and associated documentation files (the “Software”), to deal 8 in the Software without restriction, including without limitation the rights to 9 use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies 10 of the Software, and to permit persons to whom the Software is furnished to do 11 so, subject to the following conditions: 12 13 The above copyright notice and this permission notice shall be included in all 14 copies or substantial portions of the Software. 15 16 THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR 17 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 18 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE 19 AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 20 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, 21 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE 22 SOFTWARE. 23 */ 24 25 /* 26 Single-file source-code for waveout. 27 28 To compile a smaller-sized command-line app, you can use the `go` command as 29 follows: 30 31 go build -ldflags "-s -w" -trimpath waveout.go 32 */ 33 34 package main 35 36 import ( 37 "bufio" 38 "encoding/base64" 39 "encoding/binary" 40 "errors" 41 "fmt" 42 "io" 43 "math" 44 "math/rand" 45 "os" 46 "strconv" 47 "strings" 48 "time" 49 "unicode" 50 "unicode/utf8" 51 ) 52 53 // aiff header format 54 // 55 // http://paulbourke.net/dataformats/audio/ 56 // 57 // wav header format 58 // 59 // http://soundfile.sapp.org/doc/WaveFormat/ 60 // http://www-mmsp.ece.mcgill.ca/Documents/AudioFormats/WAVE/WAVE.html 61 // https://docs.fileformat.com/audio/wav/ 62 63 const ( 64 // maxInt helps convert float64 values into int16 ones 65 maxInt = 1<<15 - 1 66 67 // wavIntPCM declares integer PCM sound-data in a wav header 68 wavIntPCM = 1 69 70 // wavFloatPCM declares floating-point PCM sound-data in a wav header 71 wavFloatPCM = 3 72 ) 73 74 // emitInt16LE writes a 16-bit signed integer in little-endian byte order 75 func emitInt16LE(w io.Writer, f float64) { 76 // binary.Write(w, binary.LittleEndian, int16(maxInt*f)) 77 var buf [2]byte 78 binary.LittleEndian.PutUint16(buf[:2], uint16(int16(maxInt*f))) 79 w.Write(buf[:2]) 80 } 81 82 // emitFloat32LE writes a 32-bit float in little-endian byte order 83 func emitFloat32LE(w io.Writer, f float64) { 84 var buf [4]byte 85 binary.LittleEndian.PutUint32(buf[:4], math.Float32bits(float32(f))) 86 w.Write(buf[:4]) 87 } 88 89 // emitInt16BE writes a 16-bit signed integer in big-endian byte order 90 func emitInt16BE(w io.Writer, f float64) { 91 // binary.Write(w, binary.BigEndian, int16(maxInt*f)) 92 var buf [2]byte 93 binary.BigEndian.PutUint16(buf[:2], uint16(int16(maxInt*f))) 94 w.Write(buf[:2]) 95 } 96 97 // emitFloat32BE writes a 32-bit float in big-endian byte order 98 func emitFloat32BE(w io.Writer, f float64) { 99 var buf [4]byte 100 binary.BigEndian.PutUint32(buf[:4], math.Float32bits(float32(f))) 101 w.Write(buf[:4]) 102 } 103 104 // wavSettings is an item in the type2wavSettings table 105 type wavSettings struct { 106 Type byte 107 BitsPerSample byte 108 } 109 110 // type2wavSettings encodes values used when emitting wav headers 111 var type2wavSettings = map[sampleFormat]wavSettings{ 112 int16LE: {wavIntPCM, 16}, 113 float32LE: {wavFloatPCM, 32}, 114 } 115 116 // emitWaveHeader writes the start of a valid .wav file: since it also starts 117 // the wav data section and emits its size, you only need to write all samples 118 // after calling this func 119 func emitWaveHeader(w io.Writer, cfg outputConfig) error { 120 const fmtChunkSize = 16 121 duration := cfg.MaxTime 122 numchan := uint32(len(cfg.Scripts)) 123 sampleRate := cfg.SampleRate 124 125 ws, ok := type2wavSettings[cfg.Samples] 126 if !ok { 127 const fs = `internal error: invalid output-format code %d` 128 return fmt.Errorf(fs, cfg.Samples) 129 } 130 kind := uint16(ws.Type) 131 bps := uint32(ws.BitsPerSample) 132 133 // byte rate 134 br := sampleRate * bps * numchan / 8 135 // data size in bytes 136 dataSize := uint32(float64(br) * duration) 137 // total file size 138 totalSize := uint32(dataSize + 44) 139 140 // general descriptor 141 w.Write([]byte(`RIFF`)) 142 binary.Write(w, binary.LittleEndian, uint32(totalSize)) 143 w.Write([]byte(`WAVE`)) 144 145 // fmt chunk 146 w.Write([]byte(`fmt `)) 147 binary.Write(w, binary.LittleEndian, uint32(fmtChunkSize)) 148 binary.Write(w, binary.LittleEndian, uint16(kind)) 149 binary.Write(w, binary.LittleEndian, uint16(numchan)) 150 binary.Write(w, binary.LittleEndian, uint32(sampleRate)) 151 binary.Write(w, binary.LittleEndian, uint32(br)) 152 binary.Write(w, binary.LittleEndian, uint16(bps*numchan/8)) 153 binary.Write(w, binary.LittleEndian, uint16(bps)) 154 155 // start data chunk 156 w.Write([]byte(`data`)) 157 binary.Write(w, binary.LittleEndian, uint32(dataSize)) 158 return nil 159 } 160 161 // config has all the parsed cmd-line options 162 type config struct { 163 // Scripts has the source codes of all scripts for all channels 164 Scripts []string 165 166 // To is the output format 167 To string 168 169 // MaxTime is the play duration of the resulting sound 170 MaxTime float64 171 172 // SampleRate is the number of samples per second for all channels 173 SampleRate uint 174 } 175 176 // parseFlags is the constructor for type config 177 func parseFlags(usage string) (config, error) { 178 cfg := config{ 179 To: `wav`, 180 MaxTime: math.NaN(), 181 SampleRate: 48_000, 182 } 183 184 args := os.Args[1:] 185 if len(args) == 0 { 186 fmt.Fprint(os.Stderr, usage) 187 os.Exit(0) 188 } 189 190 for _, s := range args { 191 switch s { 192 case `help`, `-h`, `--h`, `-help`, `--help`: 193 fmt.Fprint(os.Stdout, usage) 194 os.Exit(0) 195 } 196 197 err := cfg.handleArg(s) 198 if err != nil { 199 return cfg, err 200 } 201 } 202 203 if math.IsNaN(cfg.MaxTime) { 204 cfg.MaxTime = 1 205 } 206 if cfg.MaxTime < 0 { 207 const fs = `error: given negative duration %f` 208 return cfg, fmt.Errorf(fs, cfg.MaxTime) 209 } 210 return cfg, nil 211 } 212 213 func (c *config) handleArg(s string) error { 214 switch s { 215 case `44.1k`, `44.1K`: 216 c.SampleRate = 44_100 217 return nil 218 219 case `48k`, `48K`: 220 c.SampleRate = 48_000 221 return nil 222 223 case `dat`, `DAT`: 224 c.SampleRate = 48_000 225 return nil 226 227 case `cd`, `cda`, `CD`, `CDA`: 228 c.SampleRate = 44_100 229 return nil 230 } 231 232 // handle output-format names and their aliases 233 if kind, ok := name2type[s]; ok { 234 c.To = kind 235 return nil 236 } 237 238 // handle time formats, except when they're pure numbers 239 if math.IsNaN(c.MaxTime) { 240 dur, derr := time.ParseDuration(s) 241 if derr == nil { 242 c.MaxTime = float64(dur) / float64(time.Second) 243 return nil 244 } 245 f, err := strconv.ParseFloat(s, 64) 246 if err == nil && !math.IsNaN(f) && !math.IsInf(f, 0) && f >= 0 { 247 c.MaxTime = f 248 return nil 249 } 250 } 251 252 // handle sample-rate, given either in hertz or kilohertz 253 lc := strings.ToLower(s) 254 if strings.HasSuffix(lc, `khz`) { 255 lc = strings.TrimSuffix(lc, `khz`) 256 khz, err := strconv.ParseFloat(lc, 64) 257 if err != nil || isBadNumber(khz) || khz <= 0 { 258 const fs = `invalid sample-rate frequency %q` 259 return fmt.Errorf(fs, s) 260 } 261 c.SampleRate = uint(1_000 * khz) 262 return nil 263 } else if strings.HasSuffix(lc, `hz`) { 264 lc = strings.TrimSuffix(lc, `hz`) 265 hz, err := strconv.ParseUint(lc, 10, 64) 266 if err != nil { 267 const fs = `invalid sample-rate frequency %q` 268 return fmt.Errorf(fs, s) 269 } 270 c.SampleRate = uint(hz) 271 return nil 272 } 273 274 c.Scripts = append(c.Scripts, s) 275 return nil 276 } 277 278 type encoding byte 279 type headerType byte 280 type sampleFormat byte 281 282 const ( 283 directEncoding encoding = 1 284 uriEncoding encoding = 2 285 286 noHeader headerType = 1 287 wavHeader headerType = 2 288 289 int16BE sampleFormat = 1 290 int16LE sampleFormat = 2 291 float32BE sampleFormat = 3 292 float32LE sampleFormat = 4 293 ) 294 295 // name2type normalizes keys used for type2settings 296 var name2type = map[string]string{ 297 `datauri`: `data-uri`, 298 `dataurl`: `data-uri`, 299 `data-uri`: `data-uri`, 300 `data-url`: `data-uri`, 301 `uri`: `data-uri`, 302 `url`: `data-uri`, 303 304 `raw`: `raw`, 305 `raw16be`: `raw16be`, 306 `raw16le`: `raw16le`, 307 `raw32be`: `raw32be`, 308 `raw32le`: `raw32le`, 309 310 `audio/x-wav`: `wave-16`, 311 `audio/x-wave`: `wave-16`, 312 `wav`: `wave-16`, 313 `wave`: `wave-16`, 314 `wav16`: `wave-16`, 315 `wave16`: `wave-16`, 316 `wav-16`: `wave-16`, 317 `wave-16`: `wave-16`, 318 `x-wav`: `wave-16`, 319 `x-wave`: `wave-16`, 320 321 `wav16uri`: `wave-16-uri`, 322 `wave-16-uri`: `wave-16-uri`, 323 324 `wav32uri`: `wave-32-uri`, 325 `wave-32-uri`: `wave-32-uri`, 326 327 `wav32`: `wave-32`, 328 `wave32`: `wave-32`, 329 `wav-32`: `wave-32`, 330 `wave-32`: `wave-32`, 331 } 332 333 // outputSettings are format-specific settings which are controlled by the 334 // output-format option on the cmd-line 335 type outputSettings struct { 336 Encoding encoding 337 Header headerType 338 Samples sampleFormat 339 } 340 341 // type2settings translates output-format names into the specific settings 342 // these imply 343 var type2settings = map[string]outputSettings{ 344 ``: {directEncoding, wavHeader, int16LE}, 345 346 `data-uri`: {uriEncoding, wavHeader, int16LE}, 347 `raw`: {directEncoding, noHeader, int16LE}, 348 `raw16be`: {directEncoding, noHeader, int16BE}, 349 `raw16le`: {directEncoding, noHeader, int16LE}, 350 `wave-16`: {directEncoding, wavHeader, int16LE}, 351 `wave-16-uri`: {uriEncoding, wavHeader, int16LE}, 352 353 `raw32be`: {directEncoding, noHeader, float32BE}, 354 `raw32le`: {directEncoding, noHeader, float32LE}, 355 `wave-32`: {directEncoding, wavHeader, float32LE}, 356 `wave-32-uri`: {uriEncoding, wavHeader, float32LE}, 357 } 358 359 // outputConfig has all the info the core of this app needs to make sound 360 type outputConfig struct { 361 // Scripts has the source codes of all scripts for all channels 362 Scripts []string 363 364 // MaxTime is the play duration of the resulting sound 365 MaxTime float64 366 367 // SampleRate is the number of samples per second for all channels 368 SampleRate uint32 369 370 // all the configuration details needed to emit output 371 outputSettings 372 } 373 374 // newOutputConfig is the constructor for type outputConfig, translating the 375 // cmd-line info from type config 376 func newOutputConfig(cfg config) (outputConfig, error) { 377 oc := outputConfig{ 378 Scripts: cfg.Scripts, 379 MaxTime: cfg.MaxTime, 380 SampleRate: uint32(cfg.SampleRate), 381 } 382 383 if len(oc.Scripts) == 0 { 384 return oc, errors.New(`no formulas given`) 385 } 386 387 outFmt := strings.ToLower(strings.TrimSpace(cfg.To)) 388 if alias, ok := name2type[outFmt]; ok { 389 outFmt = alias 390 } 391 392 set, ok := type2settings[outFmt] 393 if !ok { 394 const fs = `unsupported output format %q` 395 return oc, fmt.Errorf(fs, cfg.To) 396 } 397 398 oc.outputSettings = set 399 return oc, nil 400 } 401 402 // mimeType gives the format's corresponding MIME type, or an empty string 403 // if the type isn't URI-encodable 404 func (oc outputConfig) mimeType() string { 405 if oc.Header == wavHeader { 406 return `audio/x-wav` 407 } 408 return `` 409 } 410 411 func isBadNumber(f float64) bool { 412 return math.IsNaN(f) || math.IsInf(f, 0) 413 } 414 415 const usage = ` 416 waveout [options...] [duration...] [formulas...] 417 418 419 This app emits wave-sound binary data using the script(s) given. Scripts 420 give you the float64-related functionality you may expect, from numeric 421 operations to several math functions. When given 1 formula, the result is 422 mono; when given 2 formulas (left and right), the result is stereo, and so 423 on. 424 425 Output is always uncompressed audio: "waveout" can emit that as is, or as a 426 base64-encoded data-URI, which you can use as a "src" attribute value in an 427 HTML audio tag. Output duration is 1 second by default, but you can change 428 that too by using a recognized time format. 429 430 The first recognized time format is the familiar hh:mm:ss, where the hours 431 are optional, and where seconds can have a decimal part after it. 432 433 The second recognized time format uses 1-letter shortcuts instead of colons 434 for each time component, each of which is optional: "h" stands for hour, "m" 435 for minutes, and "s" for seconds. 436 437 438 Output Formats 439 440 encoding header samples endian more info 441 442 wav direct wave int16 little default format 443 444 wav16 direct wave int16 little alias for "wav" 445 wav32 direct wave float32 little 446 uri data-URI wave int16 little MIME type is audio/x-wav 447 448 raw direct none int16 little 449 raw16le direct none int16 little alias for "raw" 450 raw32le direct none float32 little 451 raw16be direct none int16 big 452 raw32be direct none float32 big 453 454 455 Concrete Examples 456 457 # low-tones commonly used in club music as beats 458 waveout 2s 'sin(10 * tau * exp(-20 * u)) * exp(-2 * u)' > club-beats.wav 459 460 # 1 minute and 5 seconds of static-like random noise 461 waveout 1m5s 'rand()' > random-noise.wav 462 463 # many bell-like clicks in quick succession; can be a cellphone's ringtone 464 waveout 'sin(2048 * tau * t) * exp(-50 * (t%0.1))' > ringtone.wav 465 466 # similar to the door-opening sound from a dsc powerseries home alarm 467 waveout 'sin(4096 * tau * t) * exp(-10 * (t%0.1))' > home-alarm.wav 468 469 # watch your ears: quickly increases frequency up to 2khz 470 waveout 'sin(2_000 * t * tau * t)' > frequency-sweep.wav 471 472 # 1-second 400hz test tone 473 waveout 'sin(400 * tau * t)' > test-tone-400.wav 474 475 # 2s of a 440hz test tone, also called an A440 sound 476 waveout 2s 'sin(440 * tau * t)' > a440.wav 477 478 # 1s 400hz test tone with sudden volume drop at the end, to avoid clip 479 waveout 'sin(400 * tau * t) * min(1, exp(-100*(t-0.9)))' > nice-tone.wav 480 481 # old ringtone used in north america 482 waveout '0.5*sin(350 * tau * t) + 0.5*sin(450 * tau * t)' > na-ringtone.wav 483 484 # 20 seconds of periodic pings 485 waveout 20s 'sin(800 * tau * u) * exp(-20 * u)' > pings.wav 486 487 # 2 seconds of a european-style dial-tone 488 waveout 2s '(sin(350 * tau * t) + sin(450 * tau * t)) / 2' > dial-tone.wav 489 490 # 4 seconds of a north-american-style busy-phone signal 491 waveout 4s '(u < 0.5) * (sin(480*tau * t) + sin(620*tau * t)) / 2' > na-busy.wav 492 493 # hit the 51st key on a synthetic piano-like instrument 494 waveout 'sin(tau * 440 * 2**((51 - 49)/12) * t) * exp(-10*u)' > piano-key.wav 495 496 # hit of a synthetic snare-like sound 497 waveout 'random() * exp(-10 * t)' > synth-snare.wav 498 499 # a stereotypical "laser" sound 500 waveout 'sin(100 * tau * exp(-40 * t))' > laser.wav 501 ` 502 503 func main() { 504 cfg, err := parseFlags(usage[1:]) 505 if err != nil { 506 fmt.Fprintln(os.Stderr, err.Error()) 507 os.Exit(1) 508 } 509 510 oc, err := newOutputConfig(cfg) 511 if err != nil { 512 fmt.Fprintln(os.Stderr, err.Error()) 513 os.Exit(1) 514 } 515 516 addDetermFuncs() 517 518 if err := run(oc); err != nil { 519 fmt.Fprintln(os.Stderr, err.Error()) 520 os.Exit(1) 521 } 522 } 523 524 func run(cfg outputConfig) error { 525 // f, err := os.Create(`waveout.prof`) 526 // if err != nil { 527 // return err 528 // } 529 // defer f.Close() 530 531 // pprof.StartCPUProfile(f) 532 // defer pprof.StopCPUProfile() 533 534 w := bufio.NewWriterSize(os.Stdout, 64*1024) 535 defer w.Flush() 536 537 switch cfg.Encoding { 538 case directEncoding: 539 return runDirect(w, cfg) 540 541 case uriEncoding: 542 mtype := cfg.mimeType() 543 if mtype == `` { 544 return errors.New(`internal error: no MIME type`) 545 } 546 547 fmt.Fprintf(w, `data:%s;base64,`, mtype) 548 enc := base64.NewEncoder(base64.StdEncoding, w) 549 defer enc.Close() 550 return runDirect(enc, cfg) 551 552 default: 553 const fs = `internal error: wrong output-encoding code %d` 554 return fmt.Errorf(fs, cfg.Encoding) 555 } 556 } 557 558 // type2emitter chooses sample-emitter funcs from the format given 559 var type2emitter = map[sampleFormat]func(io.Writer, float64){ 560 int16LE: emitInt16LE, 561 int16BE: emitInt16BE, 562 float32LE: emitFloat32LE, 563 float32BE: emitFloat32BE, 564 } 565 566 // runDirect emits sound-data bytes: this func can be called with writers 567 // which keep bytes as given, or with re-encoders, such as base64 writers 568 func runDirect(w io.Writer, cfg outputConfig) error { 569 switch cfg.Header { 570 case noHeader: 571 // do nothing, while avoiding error 572 573 case wavHeader: 574 emitWaveHeader(w, cfg) 575 576 default: 577 const fs = `internal error: wrong header code %d` 578 return fmt.Errorf(fs, cfg.Header) 579 } 580 581 emitter, ok := type2emitter[cfg.Samples] 582 if !ok { 583 const fs = `internal error: wrong output-format code %d` 584 return fmt.Errorf(fs, cfg.Samples) 585 } 586 587 if len(cfg.Scripts) == 1 { 588 return emitMono(w, cfg, emitter) 589 } 590 return emit(w, cfg, emitter) 591 } 592 593 // makeDefs makes extra funcs and values available to scripts 594 func makeDefs(cfg outputConfig) map[string]any { 595 // copy extra built-in funcs 596 defs := make(map[string]any, len(extras)+6+5) 597 for k, v := range extras { 598 defs[k] = v 599 } 600 601 // add extra variables 602 defs[`t`] = 0.0 603 defs[`u`] = 0.0 604 defs[`d`] = cfg.MaxTime 605 defs[`dur`] = cfg.MaxTime 606 defs[`duration`] = cfg.MaxTime 607 defs[`end`] = cfg.MaxTime 608 609 // add pseudo-random funcs 610 611 seed := time.Now().UnixNano() 612 r := rand.New(rand.NewSource(seed)) 613 614 rand := func() float64 { 615 return random01(r) 616 } 617 randomf := func() float64 { 618 return random(r) 619 } 620 rexpf := func(scale float64) float64 { 621 return rexp(r, scale) 622 } 623 rnormf := func(mu, sigma float64) float64 { 624 return rnorm(r, mu, sigma) 625 } 626 627 defs[`rand`] = rand 628 defs[`rand01`] = rand 629 defs[`random`] = randomf 630 defs[`rexp`] = rexpf 631 defs[`rnorm`] = rnormf 632 633 return defs 634 } 635 636 type emitFunc = func(io.Writer, float64) 637 638 // emit runs the formulas given to emit all wave samples 639 func emit(w io.Writer, cfg outputConfig, emit emitFunc) error { 640 var c Compiler 641 defs := makeDefs(cfg) 642 643 programs := make([]Program, 0, len(cfg.Scripts)) 644 tvars := make([]*float64, 0, len(cfg.Scripts)) 645 uvars := make([]*float64, 0, len(cfg.Scripts)) 646 647 for _, s := range cfg.Scripts { 648 p, err := c.Compile(s, defs) 649 if err != nil { 650 return err 651 } 652 programs = append(programs, p) 653 t, _ := p.Get(`t`) 654 u, _ := p.Get(`u`) 655 tvars = append(tvars, t) 656 uvars = append(uvars, u) 657 } 658 659 dt := 1.0 / float64(cfg.SampleRate) 660 end := cfg.MaxTime 661 662 for i := 0.0; true; i++ { 663 now := dt * i 664 if now >= end { 665 return nil 666 } 667 668 _, u := math.Modf(now) 669 670 for j, p := range programs { 671 *tvars[j] = now 672 *uvars[j] = u 673 emit(w, p.Run()) 674 } 675 } 676 return nil 677 } 678 679 // emitMono runs the formula given to emit all single-channel wave samples 680 func emitMono(w io.Writer, cfg outputConfig, emit emitFunc) error { 681 var c Compiler 682 mono, err := c.Compile(cfg.Scripts[0], makeDefs(cfg)) 683 if err != nil { 684 return err 685 } 686 687 t, _ := mono.Get(`t`) 688 u, needsu := mono.Get(`u`) 689 690 dt := 1.0 / float64(cfg.SampleRate) 691 end := cfg.MaxTime 692 693 // update variable `u` only if script uses it: this can speed things 694 // up considerably when that variable isn't used 695 if needsu { 696 for i := 0.0; true; i++ { 697 now := dt * i 698 if now >= end { 699 return nil 700 } 701 702 *t = now 703 _, *u = math.Modf(now) 704 emit(w, mono.Run()) 705 } 706 return nil 707 } 708 709 for i := 0.0; true; i++ { 710 now := dt * i 711 if now >= end { 712 return nil 713 } 714 715 *t = now 716 emit(w, mono.Run()) 717 } 718 return nil 719 } 720 721 // tau is exactly 1 loop around a circle, which is handy to turn frequencies 722 // into trigonometric angles, since they're measured in radians 723 const tau = 2 * math.Pi 724 725 // extras has funcs beyond what the script built-ins offer: those built-ins 726 // are for general math calculations, while these are specific for sound 727 // effects, other sound-related calculations, or to make pseudo-random values 728 var extras = map[string]any{ 729 `hihat`: hihat, 730 } 731 732 // addDetermFuncs does what it says, ensuring these funcs are optimizable when 733 // they're given all-constant expressions as inputs 734 func addDetermFuncs() { 735 DefineDetFuncs(map[string]any{ 736 `ascale`: AnchoredScale, 737 `awrap`: AnchoredWrap, 738 `clamp`: Clamp, 739 `epa`: Epanechnikov, 740 `epanechnikov`: Epanechnikov, 741 `fract`: Fract, 742 `gauss`: Gauss, 743 `horner`: Polyval, 744 `logistic`: Logistic, 745 `mix`: Mix, 746 `polyval`: Polyval, 747 `scale`: Scale, 748 `sign`: Sign, 749 `sinc`: Sinc, 750 `smoothstep`: SmoothStep, 751 `step`: Step, 752 `tricube`: Tricube, 753 `unwrap`: Unwrap, 754 `wrap`: Wrap, 755 756 `drop`: dropsince, 757 `dropfrom`: dropsince, 758 `dropoff`: dropsince, 759 `dropsince`: dropsince, 760 `kick`: kick, 761 `kicklow`: kicklow, 762 `piano`: piano, 763 `pianokey`: piano, 764 `pickval`: pickval, 765 `pickvalue`: pickval, 766 `sched`: schedule, 767 `schedule`: schedule, 768 `timeval`: timeval, 769 `timevalues`: timeval, 770 }) 771 } 772 773 // random01 returns a random value in 0 .. 1 774 func random01(r *rand.Rand) float64 { 775 return r.Float64() 776 } 777 778 // random returns a random value in -1 .. +1 779 func random(r *rand.Rand) float64 { 780 return (2 * r.Float64()) - 1 781 } 782 783 // rexp returns an exponentially-distributed random value using the scale 784 // (expected value) given 785 func rexp(r *rand.Rand, scale float64) float64 { 786 return scale * r.ExpFloat64() 787 } 788 789 // rnorm returns a normally-distributed random value using the mean and 790 // standard deviation given 791 func rnorm(r *rand.Rand, mu, sigma float64) float64 { 792 return r.NormFloat64()*sigma + mu 793 } 794 795 // make sample for a synthetic-drum kick 796 func kick(t float64, f, k float64) float64 { 797 const p = 0.085 798 return math.Sin(tau*f*math.Pow(p, t)) * math.Exp(-k*t) 799 } 800 801 // make sample for a heavier-sounding synthetic-drum kick 802 func kicklow(t float64, f, k float64) float64 { 803 const p = 0.08 804 return math.Sin(tau*f*math.Pow(p, t)) * math.Exp(-k*t) 805 } 806 807 // make sample for a synthetic hi-hat hit 808 func hihat(t float64, k float64) float64 { 809 return rand.Float64() * math.Exp(-k*t) 810 } 811 812 // schedule rearranges time, without being a time machine 813 func schedule(t float64, period, delay float64) float64 { 814 v := t + (1 - delay) 815 if v < 0 { 816 return 0 817 } 818 return math.Mod(v*period, period) 819 } 820 821 // make sample for a synthetic piano key being hit 822 func piano(t float64, n float64) float64 { 823 p := (math.Floor(n) - 49) / 12 824 f := 440 * math.Pow(2, p) 825 return math.Sin(tau * f * t) 826 } 827 828 // multiply rest of a formula with this for a quick volume drop at the end: 829 // this is handy to avoid clips when sounds end playing 830 func dropsince(t float64, start float64) float64 { 831 // return math.Min(1, math.Exp(-100*(t-start))) 832 if t <= start { 833 return 1 834 } 835 return math.Exp(-100 * (t - start)) 836 } 837 838 // pickval requires at least 3 args, the first 2 being the current time and 839 // each slot's duration, respectively: these 2 are followed by all the values 840 // to pick for all time slots 841 func pickval(args ...float64) float64 { 842 if len(args) < 3 { 843 return 0 844 } 845 846 t := args[0] 847 slotdur := args[1] 848 values := args[2:] 849 850 u, _ := math.Modf(t / slotdur) 851 n := len(values) 852 i := int(u) % n 853 if 0 <= i && i < n { 854 return values[i] 855 } 856 return 0 857 } 858 859 // timeval requires at least 2 args, the first 2 being the current time and 860 // the total looping-period, respectively: these 2 are followed by pairs of 861 // numbers, each consisting of a timestamp and a matching value, in order 862 func timeval(args ...float64) float64 { 863 if len(args) < 2 { 864 return 0 865 } 866 867 t := args[0] 868 period := args[1] 869 u, _ := math.Modf(t / period) 870 871 // find the first value whose periodic timestamp is due 872 for rest := args[2:]; len(rest) >= 2; rest = rest[2:] { 873 if u >= rest[0]/period { 874 return rest[1] 875 } 876 } 877 return 0 878 } 879 880 // unary2op turns unary operators into their corresponding basic operations; 881 // some entries are only for the optimizer, and aren't accessible directly 882 // from valid source code 883 var unary2op = map[string]opcode{ 884 "-": neg, 885 "!": not, 886 "&": abs, 887 "*": square, 888 "^": square, 889 "/": rec, 890 "%": mod1, 891 } 892 893 // binary2op turns binary operators into their corresponding basic operations 894 var binary2op = map[string]opcode{ 895 "+": add, 896 "-": sub, 897 "*": mul, 898 "/": div, 899 "%": mod, 900 "&&": and, 901 "&": and, 902 "||": or, 903 "|": or, 904 "==": equal, 905 "!=": notequal, 906 "<>": notequal, 907 "<": less, 908 "<=": lessoreq, 909 ">": more, 910 ">=": moreoreq, 911 "**": pow, 912 "^": pow, 913 } 914 915 // Compiler lets you create Program objects, which you can then run. The whole 916 // point of this is to create quicker-to-run numeric scripts. You can even add 917 // variables with initial values, as well as functions for the script to use. 918 // 919 // Common math funcs and constants are automatically detected, and constant 920 // results are optimized away, unless builtins are redefined. In other words, 921 // the optimizer is effectively disabled for all (sub)expressions containing 922 // redefined built-ins, as there's no way to be sure those values won't change 923 // from one run to the next. 924 // 925 // See the comment for type Program for more details. 926 // 927 // # Example 928 // 929 // var c fmscripts.Compiler 930 // 931 // defs := map[string]any{ 932 // "x": 0, // define `x`, and initialize it to 0 933 // "k": 4.25, // define `k`, and initialize it to 4.25 934 // "b": true, // define `b`, and initialize it to 1.0 935 // "n": -23, // define `n`, and initialize it to -23.0 936 // "pi": 3, // define `pi`, overriding the default constant named `pi` 937 // 938 // "f": numericKernel // type is func ([]float64) float64 939 // "g": otherFunc // type is func (float64) float64 940 // } 941 // 942 // prog, err := c.Compile("log10(k) + f(sqrt(k) * exp(-x), 45, -0.23)", defs) 943 // // ... 944 // 945 // x, _ := prog.Get("x") // Get returns (*float64, bool) 946 // y, _ := prog.Get("y") // a useless pointer, since program doesn't use `y` 947 // // ... 948 // 949 // for i := 0; i < n; i++ { 950 // *x = float64(i)*dx + minx // you update inputs in place using pointers 951 // f := prog.Run() // method Run gives you a float64 back 952 // // ... 953 // } 954 type Compiler struct { 955 maxStack int // the exact stack size the resulting program needs 956 ops []numOp // the program's operations 957 958 vaddr map[string]int // address lookup for values 959 faddr map[string]int // address lookup for functions 960 ftype map[string]any // keeps track of func types during compilation 961 962 values []float64 // variables and constants available to programs 963 funcs []any // funcs available to programs 964 } 965 966 // Compile parses the script given and generates a fast float64-only Program 967 // made only of sequential steps: any custom funcs you provide it can use 968 // their own internal looping and/or conditional logic, of course. 969 func (c *Compiler) Compile(src string, defs map[string]any) (Program, error) { 970 // turn source code into an abstract syntax-tree 971 root, err := parse(src) 972 if err != nil { 973 return Program{}, err 974 } 975 976 // generate operations 977 if err := c.reset(defs); err != nil { 978 return Program{}, err 979 } 980 if err = c.compile(root, 0); err != nil { 981 return Program{}, err 982 } 983 984 // create the resulting program 985 var p Program 986 p.stack = make([]float64, c.maxStack) 987 p.values = make([]float64, len(c.values)) 988 copy(p.values, c.values) 989 p.ops = make([]numOp, len(c.ops)) 990 copy(p.ops, c.ops) 991 p.funcs = make([]any, len(c.ftype)) 992 copy(p.funcs, c.funcs) 993 p.names = make(map[string]int, len(c.vaddr)) 994 995 // give the program's Get method access only to all allocated variables 996 for k, v := range c.vaddr { 997 // avoid exposing numeric constants on the program's name-lookup 998 // table, which would allow users to change literals across runs 999 if _, err := strconv.ParseFloat(k, 64); err == nil { 1000 continue 1001 } 1002 // expose only actual variable names 1003 p.names[k] = v 1004 } 1005 1006 return p, nil 1007 } 1008 1009 // reset prepares a compiler by satisfying internal preconditions func 1010 // compileExpr relies on, when given an abstract syntax-tree 1011 func (c *Compiler) reset(defs map[string]any) error { 1012 // reset the compiler's internal state 1013 c.maxStack = 0 1014 c.ops = c.ops[:0] 1015 c.vaddr = make(map[string]int) 1016 c.faddr = make(map[string]int) 1017 c.ftype = make(map[string]any) 1018 c.values = c.values[:0] 1019 c.funcs = c.funcs[:0] 1020 1021 // allocate vars and funcs 1022 for k, v := range defs { 1023 if err := c.allocEntry(k, v); err != nil { 1024 return err 1025 } 1026 } 1027 return nil 1028 } 1029 1030 // allocEntry simplifies the control-flow of func reset 1031 func (c *Compiler) allocEntry(k string, v any) error { 1032 const ( 1033 maxExactRound = 2 << 52 1034 rangeErrFmt = "%d is outside range of exact float64 integers" 1035 typeErrFmt = "got value of unsupported type %T" 1036 ) 1037 1038 switch v := v.(type) { 1039 case float64: 1040 _, err := c.allocValue(k, v) 1041 return err 1042 1043 case int: 1044 if math.Abs(float64(v)) > maxExactRound { 1045 return fmt.Errorf(rangeErrFmt, v) 1046 } 1047 _, err := c.allocValue(k, float64(v)) 1048 return err 1049 1050 case bool: 1051 _, err := c.allocValue(k, debool(v)) 1052 return err 1053 1054 default: 1055 if isSupportedFunc(v) { 1056 c.ftype[k] = v 1057 _, err := c.allocFunc(k, v) 1058 return err 1059 } 1060 return fmt.Errorf(typeErrFmt, v) 1061 } 1062 } 1063 1064 // allocValue ensures there's a place to match the name given, returning its 1065 // index when successful; only programs using unreasonably many values will 1066 // cause this func to fail 1067 func (c *Compiler) allocValue(s string, f float64) (int, error) { 1068 if len(c.values) >= maxOpIndex { 1069 const fs = "programs can only use up to %d distinct float64 values" 1070 return -1, fmt.Errorf(fs, maxOpIndex+1) 1071 } 1072 1073 i := len(c.values) 1074 c.vaddr[s] = i 1075 c.values = append(c.values, f) 1076 return i, nil 1077 } 1078 1079 // valueIndex returns the index reserved for an allocated value/variable: 1080 // all unallocated values/variables are allocated here on first access 1081 func (c *Compiler) valueIndex(s string, f float64) (int, error) { 1082 // name is found: return the index of the already-allocated var 1083 if i, ok := c.vaddr[s]; ok { 1084 return i, nil 1085 } 1086 1087 // name not found, but it's a known math constant 1088 if f, ok := mathConst[s]; ok { 1089 return c.allocValue(s, f) 1090 } 1091 // name not found, and it's not of a known math constant 1092 return c.allocValue(s, f) 1093 } 1094 1095 // constIndex allocates a constant as a variable named as its own string 1096 // representation, which avoids multiple entries for repeated uses of the 1097 // same constant value 1098 func (c *Compiler) constIndex(f float64) (int, error) { 1099 // constants have no name, so use a canonical string representation 1100 return c.valueIndex(strconv.FormatFloat(f, 'f', 16, 64), f) 1101 } 1102 1103 // funcIndex returns the index reserved for an allocated function: all 1104 // unallocated functions are allocated here on first access 1105 func (c *Compiler) funcIndex(name string) (int, error) { 1106 // check if func was already allocated 1107 if i, ok := c.faddr[name]; ok { 1108 return i, nil 1109 } 1110 1111 // if name is reserved allocate an index for its matching func 1112 if fn, ok := c.ftype[name]; ok { 1113 return c.allocFunc(name, fn) 1114 } 1115 1116 // if name wasn't reserved, see if it's a standard math func name 1117 if fn := c.autoFuncLookup(name); fn != nil { 1118 return c.allocFunc(name, fn) 1119 } 1120 1121 // name isn't even a standard math func's 1122 return -1, fmt.Errorf("function not found") 1123 } 1124 1125 // allocFunc ensures there's a place for the function name given, returning its 1126 // index when successful; only funcs of unsupported types will cause failure 1127 func (c *Compiler) allocFunc(name string, fn any) (int, error) { 1128 if isSupportedFunc(fn) { 1129 i := len(c.funcs) 1130 c.faddr[name] = i 1131 c.ftype[name] = fn 1132 c.funcs = append(c.funcs, fn) 1133 return i, nil 1134 } 1135 return -1, fmt.Errorf("can't use a %T as a number-crunching function", fn) 1136 } 1137 1138 // autoFuncLookup checks built-in deterministic funcs for the name given: its 1139 // result is nil only if there's no match 1140 func (c *Compiler) autoFuncLookup(name string) any { 1141 if fn, ok := determFuncs[name]; ok { 1142 return fn 1143 } 1144 return nil 1145 } 1146 1147 // genOp generates/adds a basic operation to the program, while keeping track 1148 // of the maximum depth the stack can reach 1149 func (c *Compiler) genOp(op numOp, depth int) { 1150 // add 2 defensively to ensure stack space for the inputs of binary ops 1151 n := depth + 2 1152 if c.maxStack < n { 1153 c.maxStack = n 1154 } 1155 c.ops = append(c.ops, op) 1156 } 1157 1158 // compile is a recursive expression evaluator which does the actual compiling 1159 // as it goes along 1160 func (c *Compiler) compile(expr any, depth int) error { 1161 expr = c.optimize(expr) 1162 1163 switch expr := expr.(type) { 1164 case float64: 1165 return c.compileLiteral(expr, depth) 1166 case string: 1167 return c.compileVariable(expr, depth) 1168 case []any: 1169 return c.compileCombo(expr, depth) 1170 case unaryExpr: 1171 return c.compileUnary(expr.op, expr.x, depth) 1172 case binaryExpr: 1173 return c.compileBinary(expr.op, expr.x, expr.y, depth) 1174 case callExpr: 1175 return c.compileCall(expr, depth) 1176 case assignExpr: 1177 return c.compileAssign(expr, depth) 1178 default: 1179 return fmt.Errorf("unsupported expression type %T", expr) 1180 } 1181 } 1182 1183 // compileLiteral generates a load operation for the constant value given 1184 func (c *Compiler) compileLiteral(f float64, depth int) error { 1185 i, err := c.constIndex(f) 1186 if err != nil { 1187 return err 1188 } 1189 c.genOp(numOp{What: load, Index: opindex(i)}, depth) 1190 return nil 1191 } 1192 1193 // compileVariable generates a load operation for the variable name given 1194 func (c *Compiler) compileVariable(name string, depth int) error { 1195 // handle names which aren't defined, but are known math constants 1196 if _, ok := c.vaddr[name]; !ok { 1197 if f, ok := mathConst[name]; ok { 1198 return c.compileLiteral(f, depth) 1199 } 1200 } 1201 1202 // handle actual variables 1203 i, err := c.valueIndex(name, 0) 1204 if err != nil { 1205 return err 1206 } 1207 c.genOp(numOp{What: load, Index: opindex(i)}, depth) 1208 return nil 1209 } 1210 1211 // compileCombo handles a sequence of expressions 1212 func (c *Compiler) compileCombo(exprs []any, depth int) error { 1213 for _, v := range exprs { 1214 err := c.compile(v, depth) 1215 if err != nil { 1216 return err 1217 } 1218 } 1219 return nil 1220 } 1221 1222 // compileUnary handles unary expressions 1223 func (c *Compiler) compileUnary(op string, x any, depth int) error { 1224 err := c.compile(x, depth+1) 1225 if err != nil { 1226 return err 1227 } 1228 1229 if op, ok := unary2op[op]; ok { 1230 c.genOp(numOp{What: op}, depth) 1231 return nil 1232 } 1233 return fmt.Errorf("unary operation %q is unsupported", op) 1234 } 1235 1236 // compileBinary handles binary expressions 1237 func (c *Compiler) compileBinary(op string, x, y any, depth int) error { 1238 switch op { 1239 case "===", "!==": 1240 // handle binary expressions with no matching basic operation, by 1241 // treating them as aliases for function calls 1242 return c.compileCall(callExpr{name: op, args: []any{x, y}}, depth) 1243 } 1244 1245 err := c.compile(x, depth+1) 1246 if err != nil { 1247 return err 1248 } 1249 err = c.compile(y, depth+2) 1250 if err != nil { 1251 return err 1252 } 1253 1254 if op, ok := binary2op[op]; ok { 1255 c.genOp(numOp{What: op}, depth) 1256 return nil 1257 } 1258 return fmt.Errorf("binary operation %q is unsupported", op) 1259 } 1260 1261 // compileCall handles function-call expressions 1262 func (c *Compiler) compileCall(expr callExpr, depth int) error { 1263 // lookup function name 1264 index, err := c.funcIndex(expr.name) 1265 if err != nil { 1266 return fmt.Errorf("%s: %s", expr.name, err.Error()) 1267 } 1268 // get the func value, as its type determines the calling op to emit 1269 v, ok := c.ftype[expr.name] 1270 if !ok { 1271 return fmt.Errorf("%s: function is undefined", expr.name) 1272 } 1273 1274 // figure which type of function operation to use 1275 op, ok := func2op(v) 1276 if !ok { 1277 return fmt.Errorf("%s: unsupported function type %T", expr.name, v) 1278 } 1279 1280 // ensure number of args given to the func makes sense for the func type 1281 err = checkArgCount(func2info[op], expr.name, len(expr.args)) 1282 if err != nil { 1283 return err 1284 } 1285 1286 // generate operations to evaluate all args 1287 for i, v := range expr.args { 1288 err := c.compile(v, depth+i+1) 1289 if err != nil { 1290 return err 1291 } 1292 } 1293 1294 // generate func-call operation 1295 given := len(expr.args) 1296 next := numOp{What: op, NumArgs: opargs(given), Index: opindex(index)} 1297 c.genOp(next, depth) 1298 return nil 1299 } 1300 1301 // checkArgCount does what it says, returning informative errors when arg 1302 // counts are wrong 1303 func checkArgCount(info funcTypeInfo, name string, nargs int) error { 1304 if info.AtLeast < 0 && info.AtMost < 0 { 1305 return nil 1306 } 1307 1308 if info.AtLeast == info.AtMost && nargs != info.AtMost { 1309 const fs = "%s: expected %d args, got %d instead" 1310 return fmt.Errorf(fs, name, info.AtMost, nargs) 1311 } 1312 1313 if info.AtLeast >= 0 && info.AtMost >= 0 { 1314 const fs = "%s: expected between %d and %d args, got %d instead" 1315 if nargs < info.AtLeast || nargs > info.AtMost { 1316 return fmt.Errorf(fs, name, info.AtLeast, info.AtMost, nargs) 1317 } 1318 } 1319 1320 if info.AtLeast >= 0 && nargs < info.AtLeast { 1321 const fs = "%s: expected at least %d args, got %d instead" 1322 return fmt.Errorf(fs, name, info.AtLeast, nargs) 1323 } 1324 1325 if info.AtMost >= 0 && nargs > info.AtMost { 1326 const fs = "%s: expected at most %d args, got %d instead" 1327 return fmt.Errorf(fs, name, info.AtMost, nargs) 1328 } 1329 1330 // all is good 1331 return nil 1332 } 1333 1334 // compileAssign generates a store operation for the expression given 1335 func (c *Compiler) compileAssign(expr assignExpr, depth int) error { 1336 err := c.compile(expr.expr, depth) 1337 if err != nil { 1338 return err 1339 } 1340 1341 i, err := c.allocValue(expr.name, 0) 1342 if err != nil { 1343 return err 1344 } 1345 1346 c.genOp(numOp{What: store, Index: opindex(i)}, depth) 1347 return nil 1348 } 1349 1350 // func sameFunc(x, y any) bool { 1351 // if x == nil || y == nil { 1352 // return false 1353 // } 1354 // return reflect.ValueOf(x).Pointer() == reflect.ValueOf(y).Pointer() 1355 // } 1356 1357 // isSupportedFunc checks if a value's type is a supported func type 1358 func isSupportedFunc(fn any) bool { 1359 _, ok := func2op(fn) 1360 return ok 1361 } 1362 1363 // funcTypeInfo is an entry in the func2info lookup table 1364 type funcTypeInfo struct { 1365 // AtLeast is the minimum number of inputs the func requires: negative 1366 // values are meant to be ignored 1367 AtLeast int 1368 1369 // AtMost is the maximum number of inputs the func requires: negative 1370 // values are meant to be ignored 1371 AtMost int 1372 } 1373 1374 // func2info is a lookup table to check the number of inputs funcs are given 1375 var func2info = map[opcode]funcTypeInfo{ 1376 call0: {AtLeast: +0, AtMost: +0}, 1377 call1: {AtLeast: +1, AtMost: +1}, 1378 call2: {AtLeast: +2, AtMost: +2}, 1379 call3: {AtLeast: +3, AtMost: +3}, 1380 call4: {AtLeast: +4, AtMost: +4}, 1381 call5: {AtLeast: +5, AtMost: +5}, 1382 callv: {AtLeast: -1, AtMost: -1}, 1383 call1v: {AtLeast: +1, AtMost: -1}, 1384 } 1385 1386 // func2op tries to match a func type into a corresponding opcode 1387 func func2op(v any) (op opcode, ok bool) { 1388 switch v.(type) { 1389 case func() float64: 1390 return call0, true 1391 case func(float64) float64: 1392 return call1, true 1393 case func(float64, float64) float64: 1394 return call2, true 1395 case func(float64, float64, float64) float64: 1396 return call3, true 1397 case func(float64, float64, float64, float64) float64: 1398 return call4, true 1399 case func(float64, float64, float64, float64, float64) float64: 1400 return call5, true 1401 case func(...float64) float64: 1402 return callv, true 1403 case func(float64, ...float64) float64: 1404 return call1v, true 1405 default: 1406 return 0, false 1407 } 1408 } 1409 1410 const ( 1411 // the maximum integer a float64 can represent exactly; -maxflint is the 1412 // minimum such integer, since floating-point values allow such symmetries 1413 maxflint = 2 << 52 1414 1415 // base-2 versions of size multipliers 1416 kilobyte = 1024 * 1.0 1417 megabyte = 1024 * kilobyte 1418 gigabyte = 1024 * megabyte 1419 terabyte = 1024 * gigabyte 1420 petabyte = 1024 * terabyte 1421 1422 // unit-conversion multipliers 1423 mi2kmMult = 1 / 0.6213712 1424 nm2kmMult = 1.852 1425 nmi2kmMult = 1.852 1426 yd2mtMult = 1 / 1.093613 1427 ft2mtMult = 1 / 3.28084 1428 in2cmMult = 2.54 1429 lb2kgMult = 0.4535924 1430 ga2ltMult = 1 / 0.2199692 1431 gal2lMult = 1 / 0.2199692 1432 oz2mlMult = 29.5735295625 1433 cup2lMult = 0.2365882365 1434 mpg2kplMult = 0.2199692 / 0.6213712 1435 ton2kgMult = 1 / 907.18474 1436 psi2paMult = 1 / 0.00014503773800722 1437 deg2radMult = math.Pi / 180 1438 rad2degMult = 180 / math.Pi 1439 ) 1440 1441 // default math constants 1442 var mathConst = map[string]float64{ 1443 "e": math.E, 1444 "pi": math.Pi, 1445 "tau": 2 * math.Pi, 1446 "phi": math.Phi, 1447 "nan": math.NaN(), 1448 "inf": math.Inf(+1), 1449 1450 "minint": -float64(maxflint - 1), 1451 "maxint": +float64(maxflint - 1), 1452 "minsafeint": -float64(maxflint - 1), 1453 "maxsafeint": +float64(maxflint - 1), 1454 1455 "false": 0.0, 1456 "true": 1.0, 1457 "f": 0.0, 1458 "t": 1.0, 1459 1460 // conveniently-named multipliers 1461 "femto": 1e-15, 1462 "pico": 1e-12, 1463 "nano": 1e-09, 1464 "micro": 1e-06, 1465 "milli": 1e-03, 1466 "kilo": 1e+03, 1467 "mega": 1e+06, 1468 "giga": 1e+09, 1469 "tera": 1e+12, 1470 "peta": 1e+15, 1471 1472 // unit-conversion multipliers 1473 "mi2km": mi2kmMult, 1474 "nm2km": nm2kmMult, 1475 "nmi2km": nmi2kmMult, 1476 "yd2mt": yd2mtMult, 1477 "ft2mt": ft2mtMult, 1478 "in2cm": in2cmMult, 1479 "lb2kg": lb2kgMult, 1480 "ga2lt": ga2ltMult, 1481 "gal2l": gal2lMult, 1482 "oz2ml": oz2mlMult, 1483 "cup2l": cup2lMult, 1484 "mpg2kpl": mpg2kplMult, 1485 "ton2kg": ton2kgMult, 1486 "psi2pa": psi2paMult, 1487 "deg2rad": deg2radMult, 1488 "rad2deg": rad2degMult, 1489 1490 // base-2 versions of size multipliers 1491 "kb": kilobyte, 1492 "mb": megabyte, 1493 "gb": gigabyte, 1494 "tb": terabyte, 1495 "pb": petabyte, 1496 "binkilo": kilobyte, 1497 "binmega": megabyte, 1498 "bingiga": gigabyte, 1499 "bintera": terabyte, 1500 "binpeta": petabyte, 1501 1502 // physical constants 1503 "c": 299_792_458, // speed of light in m/s 1504 "g": 6.67430e-11, // gravitational constant in N m2/kg2 1505 "h": 6.62607015e-34, // planck constant in J s 1506 "ec": 1.602176634e-19, // elementary charge in C 1507 "e0": 8.8541878128e-12, // vacuum permittivity in C2/(Nm2) 1508 "mu0": 1.25663706212e-6, // vacuum permeability in T m/A 1509 "k": 1.380649e-23, // boltzmann constant in J/K 1510 "mu": 1.66053906660e-27, // atomic mass constant in kg 1511 "me": 9.1093837015e-31, // electron mass in kg 1512 "mp": 1.67262192369e-27, // proton mass in kg 1513 "mn": 1.67492749804e-27, // neutron mass in kg 1514 1515 // float64s can only vaguely approx. avogadro's mole (6.02214076e23) 1516 } 1517 1518 // deterministic math functions lookup-table generated using the command 1519 // 1520 // go doc math | awk '/func/ { gsub(/func |\(.*/, ""); printf("\"%s\": math.%s,\n", tolower($0), $0) }' 1521 // 1522 // then hand-edited to remove funcs, or to use adapter funcs when needed: removed 1523 // funcs either had multiple returns (like SinCos) or dealt with float32 values 1524 var determFuncs = map[string]any{ 1525 "abs": math.Abs, 1526 "acos": math.Acos, 1527 "acosh": math.Acosh, 1528 "asin": math.Asin, 1529 "asinh": math.Asinh, 1530 "atan": math.Atan, 1531 "atan2": math.Atan2, 1532 "atanh": math.Atanh, 1533 "cbrt": math.Cbrt, 1534 "ceil": math.Ceil, 1535 "copysign": math.Copysign, 1536 "cos": math.Cos, 1537 "cosh": math.Cosh, 1538 "dim": math.Dim, 1539 "erf": math.Erf, 1540 "erfc": math.Erfc, 1541 "erfcinv": math.Erfcinv, 1542 "erfinv": math.Erfinv, 1543 "exp": math.Exp, 1544 "exp2": math.Exp2, 1545 "expm1": math.Expm1, 1546 "fma": math.FMA, 1547 "floor": math.Floor, 1548 "gamma": math.Gamma, 1549 "inf": inf, 1550 "isinf": isInf, 1551 "isnan": isNaN, 1552 "j0": math.J0, 1553 "j1": math.J1, 1554 "jn": jn, 1555 "ldexp": ldexp, 1556 "lgamma": lgamma, 1557 "log10": math.Log10, 1558 "log1p": math.Log1p, 1559 "log2": math.Log2, 1560 "logb": math.Logb, 1561 "mod": math.Mod, 1562 "nan": math.NaN, 1563 "nextafter": math.Nextafter, 1564 "pow": math.Pow, 1565 "pow10": pow10, 1566 "remainder": math.Remainder, 1567 "round": math.Round, 1568 "roundtoeven": math.RoundToEven, 1569 "signbit": signbit, 1570 "sin": math.Sin, 1571 "sinh": math.Sinh, 1572 "sqrt": math.Sqrt, 1573 "tan": math.Tan, 1574 "tanh": math.Tanh, 1575 "trunc": math.Trunc, 1576 "y0": math.Y0, 1577 "y1": math.Y1, 1578 "yn": yn, 1579 1580 // a few aliases for the standard math funcs: some of the single-letter 1581 // aliases are named after the ones in `bc`, the basic calculator tool 1582 "a": math.Abs, 1583 "c": math.Cos, 1584 "ceiling": math.Ceil, 1585 "cosine": math.Cos, 1586 "e": math.Exp, 1587 "isinf0": isAnyInf, 1588 "isinfinite": isAnyInf, 1589 "l": math.Log, 1590 "ln": math.Log, 1591 "lg": math.Log2, 1592 "modulus": math.Mod, 1593 "power": math.Pow, 1594 "rem": math.Remainder, 1595 "s": math.Sin, 1596 "sine": math.Sin, 1597 "t": math.Tan, 1598 "tangent": math.Tan, 1599 "truncate": math.Trunc, 1600 "truncated": math.Trunc, 1601 1602 // not from standard math: these custom funcs were added manually 1603 "beta": beta, 1604 "bool": num2bool, 1605 "clamp": clamp, 1606 "cond": cond, // vector-arg if-else chain 1607 "cube": cube, 1608 "cubed": cube, 1609 "degrees": degrees, 1610 "deinf": deInf, 1611 "denan": deNaN, 1612 "factorial": factorial, 1613 "fract": fract, 1614 "horner": polyval, 1615 "hypot": hypot, 1616 "if": ifElse, 1617 "ifelse": ifElse, 1618 "inv": reciprocal, 1619 "isanyinf": isAnyInf, 1620 "isbad": isBad, 1621 "isfin": isFinite, 1622 "isfinite": isFinite, 1623 "isgood": isGood, 1624 "isinteger": isInteger, 1625 "lbeta": lnBeta, 1626 "len": length, 1627 "length": length, 1628 "lnbeta": lnBeta, 1629 "log": math.Log, 1630 "logistic": logistic, 1631 "mag": length, 1632 "max": max, 1633 "min": min, 1634 "neg": negate, 1635 "negate": negate, 1636 "not": notBool, 1637 "polyval": polyval, 1638 "radians": radians, 1639 "range": rangef, // vector-arg max - min 1640 "reciprocal": reciprocal, 1641 "rev": revalue, 1642 "revalue": revalue, 1643 "scale": scale, 1644 "sgn": sign, 1645 "sign": sign, 1646 "sinc": sinc, 1647 "sq": sqr, 1648 "sqmin": solveQuadMin, 1649 "sqmax": solveQuadMax, 1650 "square": sqr, 1651 "squared": sqr, 1652 "unwrap": unwrap, 1653 "wrap": wrap, 1654 1655 // a few aliases for the custom funcs 1656 "deg": degrees, 1657 "isint": isInteger, 1658 "rad": radians, 1659 1660 // a few quicker 2-value versions of vararg funcs: the optimizer depends 1661 // on these to rewrite 2-input uses of their vararg counterparts 1662 "hypot2": math.Hypot, 1663 "max2": math.Max, 1664 "min2": math.Min, 1665 1666 // a few entries to enable custom syntax: the parser depends on these to 1667 // rewrite binary expressions into func calls 1668 "??": deNaN, 1669 "?:": ifElse, 1670 "===": same, 1671 "!==": notSame, 1672 } 1673 1674 // DefineDetFuncs adds more deterministic funcs to the default set. Such funcs 1675 // are considered optimizable, since calling them with the same constant inputs 1676 // is supposed to return the same constant outputs, as the name `deterministic` 1677 // suggests. 1678 // 1679 // Only call this before compiling any scripts, and ensure all funcs given are 1680 // supported and are deterministic. Random-output funcs certainly won't fit 1681 // the bill here. 1682 func DefineDetFuncs(funcs map[string]any) { 1683 for k, v := range funcs { 1684 determFuncs[k] = v 1685 } 1686 } 1687 1688 func sqr(x float64) float64 { 1689 return x * x 1690 } 1691 1692 func cube(x float64) float64 { 1693 return x * x * x 1694 } 1695 1696 func num2bool(x float64) float64 { 1697 if x == 0 { 1698 return 0 1699 } 1700 return 1 1701 } 1702 1703 func logistic(x float64) float64 { 1704 return 1 / (1 + math.Exp(-x)) 1705 } 1706 1707 func sign(x float64) float64 { 1708 if math.IsNaN(x) { 1709 return x 1710 } 1711 if x > 0 { 1712 return +1 1713 } 1714 if x < 0 { 1715 return -1 1716 } 1717 return 0 1718 } 1719 1720 func sinc(x float64) float64 { 1721 if x == 0 { 1722 return 1 1723 } 1724 return math.Sin(x) / x 1725 } 1726 1727 func isInteger(x float64) float64 { 1728 _, frac := math.Modf(x) 1729 if frac == 0 { 1730 return 1 1731 } 1732 return 0 1733 } 1734 1735 func inf(sign float64) float64 { 1736 return math.Inf(int(sign)) 1737 } 1738 1739 func isInf(x float64, sign float64) float64 { 1740 if math.IsInf(x, int(sign)) { 1741 return 1 1742 } 1743 return 0 1744 } 1745 1746 func isAnyInf(x float64) float64 { 1747 if math.IsInf(x, 0) { 1748 return 1 1749 } 1750 return 0 1751 } 1752 1753 func isFinite(x float64) float64 { 1754 if math.IsInf(x, 0) { 1755 return 0 1756 } 1757 return 1 1758 } 1759 1760 func isNaN(x float64) float64 { 1761 if math.IsNaN(x) { 1762 return 1 1763 } 1764 return 0 1765 } 1766 1767 func isGood(x float64) float64 { 1768 if math.IsNaN(x) || math.IsInf(x, 0) { 1769 return 0 1770 } 1771 return 1 1772 } 1773 1774 func isBad(x float64) float64 { 1775 if math.IsNaN(x) || math.IsInf(x, 0) { 1776 return 1 1777 } 1778 return 0 1779 } 1780 1781 func same(x, y float64) float64 { 1782 if math.IsNaN(x) && math.IsNaN(y) { 1783 return 1 1784 } 1785 return debool(x == y) 1786 } 1787 1788 func notSame(x, y float64) float64 { 1789 if math.IsNaN(x) && math.IsNaN(y) { 1790 return 0 1791 } 1792 return debool(x != y) 1793 } 1794 1795 func deNaN(x float64, instead float64) float64 { 1796 if !math.IsNaN(x) { 1797 return x 1798 } 1799 return instead 1800 } 1801 1802 func deInf(x float64, instead float64) float64 { 1803 if !math.IsInf(x, 0) { 1804 return x 1805 } 1806 return instead 1807 } 1808 1809 func revalue(x float64, instead float64) float64 { 1810 if !math.IsNaN(x) && !math.IsInf(x, 0) { 1811 return x 1812 } 1813 return instead 1814 } 1815 1816 func pow10(n float64) float64 { 1817 return math.Pow10(int(n)) 1818 } 1819 1820 func jn(n, x float64) float64 { 1821 return math.Jn(int(n), x) 1822 } 1823 1824 func ldexp(frac, exp float64) float64 { 1825 return math.Ldexp(frac, int(exp)) 1826 } 1827 1828 func lgamma(x float64) float64 { 1829 y, s := math.Lgamma(x) 1830 if s < 0 { 1831 return math.NaN() 1832 } 1833 return y 1834 } 1835 1836 func signbit(x float64) float64 { 1837 if math.Signbit(x) { 1838 return 1 1839 } 1840 return 0 1841 } 1842 1843 func yn(n, x float64) float64 { 1844 return math.Yn(int(n), x) 1845 } 1846 1847 func negate(x float64) float64 { 1848 return -x 1849 } 1850 1851 func reciprocal(x float64) float64 { 1852 return 1 / x 1853 } 1854 1855 func rangef(v ...float64) float64 { 1856 min := math.Inf(+1) 1857 max := math.Inf(-1) 1858 for _, f := range v { 1859 min = math.Min(min, f) 1860 max = math.Max(max, f) 1861 } 1862 return max - min 1863 } 1864 1865 func cond(v ...float64) float64 { 1866 for { 1867 switch len(v) { 1868 case 0: 1869 // either no values are left, or no values were given at all 1870 return math.NaN() 1871 1872 case 1: 1873 // either all previous conditions failed, and this last value is 1874 // automatically chosen, or only 1 value was given to begin with 1875 return v[0] 1876 1877 default: 1878 // check condition: if true (non-0), return the value after it 1879 if v[0] != 0 { 1880 return v[1] 1881 } 1882 // condition was false, so skip the leading pair of values 1883 v = v[2:] 1884 } 1885 } 1886 } 1887 1888 func notBool(x float64) float64 { 1889 if x == 0 { 1890 return 1 1891 } 1892 return 0 1893 } 1894 1895 func ifElse(cond float64, yes, no float64) float64 { 1896 if cond != 0 { 1897 return yes 1898 } 1899 return no 1900 } 1901 1902 func lnGamma(x float64) float64 { 1903 y, s := math.Lgamma(x) 1904 if s < 0 { 1905 return math.NaN() 1906 } 1907 return y 1908 } 1909 1910 func lnBeta(x float64, y float64) float64 { 1911 return lnGamma(x) + lnGamma(y) - lnGamma(x+y) 1912 } 1913 1914 func beta(x float64, y float64) float64 { 1915 return math.Exp(lnBeta(x, y)) 1916 } 1917 1918 func factorial(n float64) float64 { 1919 return math.Round(math.Gamma(n + 1)) 1920 } 1921 1922 func degrees(rad float64) float64 { 1923 return rad2degMult * rad 1924 } 1925 1926 func radians(deg float64) float64 { 1927 return deg2radMult * deg 1928 } 1929 1930 func fract(x float64) float64 { 1931 return x - math.Floor(x) 1932 } 1933 1934 func min(v ...float64) float64 { 1935 min := +math.Inf(+1) 1936 for _, f := range v { 1937 min = math.Min(min, f) 1938 } 1939 return min 1940 } 1941 1942 func max(v ...float64) float64 { 1943 max := +math.Inf(-1) 1944 for _, f := range v { 1945 max = math.Max(max, f) 1946 } 1947 return max 1948 } 1949 1950 func hypot(v ...float64) float64 { 1951 sumsq := 0.0 1952 for _, f := range v { 1953 sumsq += f * f 1954 } 1955 return math.Sqrt(sumsq) 1956 } 1957 1958 func length(v ...float64) float64 { 1959 ss := 0.0 1960 for _, f := range v { 1961 ss += f * f 1962 } 1963 return math.Sqrt(ss) 1964 } 1965 1966 // solveQuadMin finds the lowest solution of a 2nd-degree polynomial, using 1967 // a formula which is more accurate than the textbook one 1968 func solveQuadMin(a, b, c float64) float64 { 1969 disc := math.Sqrt(b*b - 4*a*c) 1970 r1 := 2 * c / (-b - disc) 1971 return r1 1972 } 1973 1974 // solveQuadMax finds the highest solution of a 2nd-degree polynomial, using 1975 // a formula which is more accurate than the textbook one 1976 func solveQuadMax(a, b, c float64) float64 { 1977 disc := math.Sqrt(b*b - 4*a*c) 1978 r2 := 2 * c / (-b + disc) 1979 return r2 1980 } 1981 1982 func wrap(x float64, min, max float64) float64 { 1983 return (x - min) / (max - min) 1984 } 1985 1986 func unwrap(x float64, min, max float64) float64 { 1987 return (max-min)*x + min 1988 } 1989 1990 func clamp(x float64, min, max float64) float64 { 1991 return math.Min(math.Max(x, min), max) 1992 } 1993 1994 func scale(x float64, xmin, xmax, ymin, ymax float64) float64 { 1995 k := (x - xmin) / (xmax - xmin) 1996 return (ymax-ymin)*k + ymin 1997 } 1998 1999 // polyval runs horner's algorithm on a value, with the polynomial coefficients 2000 // given after it, higher-order first 2001 func polyval(x float64, v ...float64) float64 { 2002 if len(v) == 0 { 2003 return 0 2004 } 2005 2006 x0 := x 2007 x = 1.0 2008 y := 0.0 2009 for i := len(v) - 1; i >= 0; i-- { 2010 y += v[i] * x 2011 x *= x0 2012 } 2013 return y 2014 } 2015 2016 // unary2func matches unary operators to funcs the optimizer can use to eval 2017 // constant-input unary expressions into their results 2018 var unary2func = map[string]func(x float64) float64{ 2019 // avoid unary +, since it's a no-op and thus there's no instruction for 2020 // it, which in turn makes unit-tests for these optimization tables fail 2021 // unnecessarily 2022 // "+": func(x float64) float64 { return +x }, 2023 2024 "-": func(x float64) float64 { return -x }, 2025 "!": func(x float64) float64 { return debool(x == 0) }, 2026 "&": func(x float64) float64 { return math.Abs(x) }, 2027 "*": func(x float64) float64 { return x * x }, 2028 "^": func(x float64) float64 { return x * x }, 2029 "/": func(x float64) float64 { return 1 / x }, 2030 } 2031 2032 // binary2func matches binary operators to funcs the optimizer can use to eval 2033 // constant-input binary expressions into their results 2034 var binary2func = map[string]func(x, y float64) float64{ 2035 "+": func(x, y float64) float64 { return x + y }, 2036 "-": func(x, y float64) float64 { return x - y }, 2037 "*": func(x, y float64) float64 { return x * y }, 2038 "/": func(x, y float64) float64 { return x / y }, 2039 "%": func(x, y float64) float64 { return math.Mod(x, y) }, 2040 "&": func(x, y float64) float64 { return debool(x != 0 && y != 0) }, 2041 "&&": func(x, y float64) float64 { return debool(x != 0 && y != 0) }, 2042 "|": func(x, y float64) float64 { return debool(x != 0 || y != 0) }, 2043 "||": func(x, y float64) float64 { return debool(x != 0 || y != 0) }, 2044 "==": func(x, y float64) float64 { return debool(x == y) }, 2045 "!=": func(x, y float64) float64 { return debool(x != y) }, 2046 "<>": func(x, y float64) float64 { return debool(x != y) }, 2047 "<": func(x, y float64) float64 { return debool(x < y) }, 2048 "<=": func(x, y float64) float64 { return debool(x <= y) }, 2049 ">": func(x, y float64) float64 { return debool(x > y) }, 2050 ">=": func(x, y float64) float64 { return debool(x >= y) }, 2051 "**": func(x, y float64) float64 { return math.Pow(x, y) }, 2052 "^": func(x, y float64) float64 { return math.Pow(x, y) }, 2053 } 2054 2055 // func2unary turns built-in func names into built-in unary operators 2056 var func2unary = map[string]string{ 2057 "a": "&", 2058 "abs": "&", 2059 "inv": "/", 2060 "neg": "-", 2061 "negate": "-", 2062 "reciprocal": "/", 2063 "sq": "*", 2064 "square": "*", 2065 "squared": "*", 2066 } 2067 2068 // func2binary turns built-in func names into built-in binary operators 2069 var func2binary = map[string]string{ 2070 "mod": "%", 2071 "modulus": "%", 2072 "pow": "**", 2073 "power": "**", 2074 } 2075 2076 // vararg2func2 matches variable-argument funcs to their 2-argument versions 2077 var vararg2func2 = map[string]string{ 2078 "hypot": "hypot2", 2079 "max": "max2", 2080 "min": "min2", 2081 } 2082 2083 // optimize tries to simplify the expression given as much as possible, by 2084 // simplifying constants whenever possible, and exploiting known built-in 2085 // funcs which are known to behave deterministically 2086 func (c *Compiler) optimize(expr any) any { 2087 switch expr := expr.(type) { 2088 case []any: 2089 return c.optimizeCombo(expr) 2090 2091 case unaryExpr: 2092 return c.optimizeUnaryExpr(expr) 2093 2094 case binaryExpr: 2095 return c.optimizeBinaryExpr(expr) 2096 2097 case callExpr: 2098 return c.optimizeCallExpr(expr) 2099 2100 case assignExpr: 2101 expr.expr = c.optimize(expr.expr) 2102 return expr 2103 2104 default: 2105 f, ok := c.tryConstant(expr) 2106 if ok { 2107 return f 2108 } 2109 return expr 2110 } 2111 } 2112 2113 // optimizeCombo handles a sequence of expressions for the optimizer 2114 func (c *Compiler) optimizeCombo(exprs []any) any { 2115 if len(exprs) == 1 { 2116 return c.optimize(exprs[0]) 2117 } 2118 2119 // count how many expressions are considered useful: these are 2120 // assignments as well as the last expression, since that's what 2121 // determines the script's final result 2122 useful := 0 2123 for i, v := range exprs { 2124 _, ok := v.(assignExpr) 2125 if ok || i == len(exprs)-1 { 2126 useful++ 2127 } 2128 } 2129 2130 // ignore all expressions which are a waste of time, and optimize 2131 // all other expressions 2132 res := make([]any, 0, useful) 2133 for i, v := range exprs { 2134 _, ok := v.(assignExpr) 2135 if ok || i == len(exprs)-1 { 2136 res = append(res, c.optimize(v)) 2137 } 2138 } 2139 return res 2140 } 2141 2142 // optimizeUnaryExpr handles unary expressions for the optimizer 2143 func (c *Compiler) optimizeUnaryExpr(expr unaryExpr) any { 2144 // recursively optimize input 2145 expr.x = c.optimize(expr.x) 2146 2147 // optimize unary ops on a constant into concrete values 2148 if x, ok := expr.x.(float64); ok { 2149 if fn, ok := unary2func[expr.op]; ok { 2150 return fn(x) 2151 } 2152 } 2153 2154 switch expr.op { 2155 case "+": 2156 // unary plus is an identity operation 2157 return expr.x 2158 2159 default: 2160 return expr 2161 } 2162 } 2163 2164 // optimizeBinaryExpr handles binary expressions for the optimizer 2165 func (c *Compiler) optimizeBinaryExpr(expr binaryExpr) any { 2166 // recursively optimize inputs 2167 expr.x = c.optimize(expr.x) 2168 expr.y = c.optimize(expr.y) 2169 2170 // optimize binary ops on 2 constants into concrete values 2171 if x, ok := expr.x.(float64); ok { 2172 if y, ok := expr.y.(float64); ok { 2173 if fn, ok := binary2func[expr.op]; ok { 2174 return fn(x, y) 2175 } 2176 } 2177 } 2178 2179 switch expr.op { 2180 case "+": 2181 if expr.x == 0.0 { 2182 // 0+y -> y 2183 return expr.y 2184 } 2185 if expr.y == 0.0 { 2186 // x+0 -> x 2187 return expr.x 2188 } 2189 2190 case "-": 2191 if expr.x == 0.0 { 2192 // 0-y -> -y 2193 return c.optimizeUnaryExpr(unaryExpr{op: "-", x: expr.y}) 2194 } 2195 if expr.y == 0.0 { 2196 // x-0 -> x 2197 return expr.x 2198 } 2199 2200 case "*": 2201 if expr.x == 0.0 || expr.y == 0.0 { 2202 // 0*y -> 0 2203 // x*0 -> 0 2204 return 0.0 2205 } 2206 if expr.x == 1.0 { 2207 // 1*y -> y 2208 return expr.y 2209 } 2210 if expr.y == 1.0 { 2211 // x*1 -> x 2212 return expr.x 2213 } 2214 if expr.x == -1.0 { 2215 // -1*y -> -y 2216 return c.optimizeUnaryExpr(unaryExpr{op: "-", x: expr.y}) 2217 } 2218 if expr.y == -1.0 { 2219 // x*-1 -> -x 2220 return c.optimizeUnaryExpr(unaryExpr{op: "-", x: expr.x}) 2221 } 2222 2223 case "/": 2224 if expr.x == 1.0 { 2225 // 1/y -> reciprocal of y 2226 return c.optimizeUnaryExpr(unaryExpr{op: "/", x: expr.y}) 2227 } 2228 if expr.y == 1.0 { 2229 // x/1 -> x 2230 return expr.x 2231 } 2232 if expr.y == -1.0 { 2233 // x/-1 -> -x 2234 return c.optimizeUnaryExpr(unaryExpr{op: "-", x: expr.x}) 2235 } 2236 2237 case "**": 2238 switch expr.y { 2239 case -1.0: 2240 // x**-1 -> 1/x, reciprocal of x 2241 return c.optimizeUnaryExpr(unaryExpr{op: "/", x: expr.x}) 2242 case 0.0: 2243 // x**0 -> 1 2244 return 1.0 2245 case 1.0: 2246 // x**1 -> x 2247 return expr.x 2248 case 2.0: 2249 // x**2 -> *x 2250 return c.optimizeUnaryExpr(unaryExpr{op: "*", x: expr.x}) 2251 case 3.0: 2252 // x**3 -> *x*x 2253 sq := unaryExpr{op: "*", x: expr.x} 2254 return c.optimizeBinaryExpr(binaryExpr{op: "*", x: sq, y: expr.x}) 2255 } 2256 2257 case "&", "&&": 2258 if expr.x == 0.0 || expr.y == 0.0 { 2259 // 0 && y -> 0 2260 // x && 0 -> 0 2261 return 0.0 2262 } 2263 } 2264 2265 // no simplifiable patterns were detected 2266 return expr 2267 } 2268 2269 // optimizeCallExpr optimizes special cases of built-in func calls 2270 func (c *Compiler) optimizeCallExpr(call callExpr) any { 2271 // recursively optimize all inputs, and keep track if they're all 2272 // constants in the end 2273 numlit := 0 2274 for i, v := range call.args { 2275 v = c.optimize(v) 2276 call.args[i] = v 2277 if _, ok := v.(float64); ok { 2278 numlit++ 2279 } 2280 } 2281 2282 // if func is overridden, there's no guarantee the new func works the same 2283 if _, ok := determFuncs[call.name]; c.ftype[call.name] != nil || !ok { 2284 return call 2285 } 2286 2287 // from this point on, func is guaranteed to be built-in and deterministic 2288 2289 // handle all-const inputs, by calling func and using its result 2290 if numlit == len(call.args) { 2291 in := make([]float64, 0, len(call.args)) 2292 for _, v := range call.args { 2293 f, _ := v.(float64) 2294 in = append(in, f) 2295 } 2296 2297 if f, ok := tryCall(determFuncs[call.name], in); ok { 2298 return f 2299 } 2300 } 2301 2302 switch len(call.args) { 2303 case 1: 2304 if op, ok := func2unary[call.name]; ok { 2305 expr := unaryExpr{op: op, x: call.args[0]} 2306 return c.optimizeUnaryExpr(expr) 2307 } 2308 return call 2309 2310 case 2: 2311 if op, ok := func2binary[call.name]; ok { 2312 expr := binaryExpr{op: op, x: call.args[0], y: call.args[1]} 2313 return c.optimizeBinaryExpr(expr) 2314 } 2315 if name, ok := vararg2func2[call.name]; ok { 2316 call.name = name 2317 return call 2318 } 2319 return call 2320 2321 default: 2322 return call 2323 } 2324 } 2325 2326 // tryConstant tries to optimize the expression given into a constant 2327 func (c *Compiler) tryConstant(expr any) (value float64, ok bool) { 2328 switch expr := expr.(type) { 2329 case float64: 2330 return expr, true 2331 2332 case string: 2333 if _, ok := c.vaddr[expr]; !ok { 2334 // name isn't explicitly defined 2335 if f, ok := mathConst[expr]; ok { 2336 // and is a known math constant 2337 return f, true 2338 } 2339 } 2340 return 0, false 2341 2342 default: 2343 return 0, false 2344 } 2345 } 2346 2347 // tryCall tries to simplify the function expression given into a constant 2348 func tryCall(fn any, in []float64) (value float64, ok bool) { 2349 switch fn := fn.(type) { 2350 case func(float64) float64: 2351 if len(in) == 1 { 2352 return fn(in[0]), true 2353 } 2354 return 0, false 2355 2356 case func(float64, float64) float64: 2357 if len(in) == 2 { 2358 return fn(in[0], in[1]), true 2359 } 2360 return 0, false 2361 2362 case func(float64, float64, float64) float64: 2363 if len(in) == 3 { 2364 return fn(in[0], in[1], in[2]), true 2365 } 2366 return 0, false 2367 2368 case func(float64, float64, float64, float64) float64: 2369 if len(in) == 4 { 2370 return fn(in[0], in[1], in[2], in[3]), true 2371 } 2372 return 0, false 2373 2374 case func(float64, float64, float64, float64, float64) float64: 2375 if len(in) == 5 { 2376 return fn(in[0], in[1], in[2], in[3], in[4]), true 2377 } 2378 return 0, false 2379 2380 case func(...float64) float64: 2381 return fn(in...), true 2382 2383 case func(float64, ...float64) float64: 2384 if len(in) >= 1 { 2385 return fn(in[0], in[1:]...), true 2386 } 2387 return 0, false 2388 2389 default: 2390 // type isn't a supported func 2391 return 0, false 2392 } 2393 } 2394 2395 // parse turns source code into an expression type interpreters can use. 2396 func parse(src string) (any, error) { 2397 tok := newTokenizer(src) 2398 par, err := newParser(&tok) 2399 if err != nil { 2400 return nil, err 2401 } 2402 2403 v, err := par.parse() 2404 err = par.improveError(err, src) 2405 return v, err 2406 } 2407 2408 // pickLine slices a string, picking one of its lines via a 1-based index: 2409 // func improveError uses it to isolate the line of code an error came from 2410 func pickLine(src string, linenum int) string { 2411 // skip the lines before the target one 2412 for i := 0; i < linenum && len(src) > 0; i++ { 2413 j := strings.IndexByte(src, '\n') 2414 if j < 0 { 2415 break 2416 } 2417 src = src[j+1:] 2418 } 2419 2420 // limit leftover to a single line 2421 i := strings.IndexByte(src, '\n') 2422 if i >= 0 { 2423 return src[:i] 2424 } 2425 return src 2426 } 2427 2428 // unaryExpr is a unary expression 2429 type unaryExpr struct { 2430 op string 2431 x any 2432 } 2433 2434 // binaryExpr is a binary expression 2435 type binaryExpr struct { 2436 op string 2437 x any 2438 y any 2439 } 2440 2441 // callExpr is a function call and its arguments 2442 type callExpr struct { 2443 name string 2444 args []any 2445 } 2446 2447 // assignExpr is a value/variable assignment 2448 type assignExpr struct { 2449 name string 2450 expr any 2451 } 2452 2453 // parser is a parser for JavaScript-like syntax, limited to operations on 2454 // 64-bit floating-point numbers 2455 type parser struct { 2456 tokens []token 2457 line int 2458 pos int 2459 toklen int 2460 } 2461 2462 // newParser is the constructor for type parser 2463 func newParser(t *tokenizer) (parser, error) { 2464 var p parser 2465 2466 // get all tokens from the source code 2467 for { 2468 v, err := t.next() 2469 if v.kind != unknownToken { 2470 p.tokens = append(p.tokens, v) 2471 } 2472 2473 if err == errEOS { 2474 // done scanning/tokenizing 2475 return p, nil 2476 } 2477 2478 if err != nil { 2479 // handle actual errors 2480 return p, err 2481 } 2482 } 2483 } 2484 2485 // improveError adds more info about where exactly in the source code an error 2486 // came from, thus making error messages much more useful 2487 func (p *parser) improveError(err error, src string) error { 2488 if err == nil { 2489 return nil 2490 } 2491 2492 line := pickLine(src, p.line) 2493 if len(line) == 0 || p.pos < 1 { 2494 const fs = "(line %d: pos %d): %w" 2495 return fmt.Errorf(fs, p.line, p.pos, err) 2496 } 2497 2498 ptr := strings.Repeat(" ", p.pos) + "^" 2499 const fs = "(line %d: pos %d): %w\n%s\n%s" 2500 return fmt.Errorf(fs, p.line, p.pos, err, line, ptr) 2501 } 2502 2503 // parse tries to turn tokens into a compilable abstract syntax tree, and is 2504 // the parser's entry point 2505 func (p *parser) parse() (any, error) { 2506 // source codes with no tokens are always errors 2507 if len(p.tokens) == 0 { 2508 const msg = "source code is empty, or has no useful expressions" 2509 return nil, errors.New(msg) 2510 } 2511 2512 // handle optional excel-like leading equal sign 2513 p.acceptSyntax(`=`) 2514 2515 // ignore trailing semicolons 2516 for len(p.tokens) > 0 { 2517 t := p.tokens[len(p.tokens)-1] 2518 if t.kind != syntaxToken || t.value != `;` { 2519 break 2520 } 2521 p.tokens = p.tokens[:len(p.tokens)-1] 2522 } 2523 2524 // handle single expressions as well as multiple semicolon-separated 2525 // expressions: the latter allow value assignments/updates in scripts 2526 var res []any 2527 for keepGoing := true; keepGoing && len(p.tokens) > 0; { 2528 v, err := p.parseExpression() 2529 if err != nil && err != errEOS { 2530 return v, err 2531 } 2532 2533 res = append(res, v) 2534 // handle optional separator/continuation semicolons 2535 _, keepGoing = p.acceptSyntax(`;`) 2536 } 2537 2538 // unexpected unparsed trailing tokens are always errors; any trailing 2539 // semicolons in the original script are already trimmed 2540 if len(p.tokens) > 0 { 2541 const fs = "unexpected %s" 2542 return res, fmt.Errorf(fs, p.tokens[0].value) 2543 } 2544 2545 // make scripts ending in an assignment also load that value, so they're 2546 // useful, as assignments result in no useful value by themselves 2547 assign, ok := res[len(res)-1].(assignExpr) 2548 if ok { 2549 res = append(res, assign.name) 2550 } 2551 2552 // turn 1-item combo expressions into their only expression: some 2553 // unit tests may rely on that for convenience 2554 if len(res) == 1 { 2555 return res[0], nil 2556 } 2557 return res, nil 2558 } 2559 2560 // acceptSyntax advances the parser on the first syntactic string matched: 2561 // notice any number of alternatives/options are allowed, as the syntax 2562 // allows/requires at various points 2563 func (p *parser) acceptSyntax(syntax ...string) (match string, ok bool) { 2564 if len(p.tokens) == 0 { 2565 return "", false 2566 } 2567 2568 t := p.tokens[0] 2569 if t.kind != syntaxToken { 2570 return "", false 2571 } 2572 2573 for _, s := range syntax { 2574 if t.value == s { 2575 p.advance() 2576 return s, true 2577 } 2578 } 2579 return "", false 2580 } 2581 2582 // advance skips the current leading token, if there are still any left 2583 func (p *parser) advance() { 2584 if len(p.tokens) == 0 { 2585 return 2586 } 2587 2588 t := p.tokens[0] 2589 p.tokens = p.tokens[1:] 2590 p.line = t.line 2591 p.pos = t.pos 2592 p.toklen = len(t.value) 2593 } 2594 2595 // acceptNumeric advances the parser on a numeric value, but only if it's 2596 // the leading token: conversely, any other type of token doesn't advance 2597 // the parser; when matches happen the resulting strings need parsing via 2598 // func parseNumber 2599 func (p *parser) acceptNumeric() (numliteral string, ok bool) { 2600 if len(p.tokens) == 0 { 2601 return "", false 2602 } 2603 2604 t := p.tokens[0] 2605 if t.kind == numberToken { 2606 p.advance() 2607 return t.value, true 2608 } 2609 return "", false 2610 } 2611 2612 // demandSyntax imposes a specific syntactic element to follow, or else it's 2613 // an error 2614 func (p *parser) demandSyntax(syntax string) error { 2615 if len(p.tokens) == 0 { 2616 return fmt.Errorf("expected %s instead of the end of source", syntax) 2617 } 2618 2619 first := p.tokens[0] 2620 if first.kind == syntaxToken && first.value == syntax { 2621 p.advance() 2622 return nil 2623 } 2624 return fmt.Errorf("expected %s instead of %s", syntax, first.value) 2625 } 2626 2627 func (p *parser) parseExpression() (any, error) { 2628 x, err := p.parseComparison() 2629 if err != nil { 2630 return x, err 2631 } 2632 2633 // handle assignment statements 2634 if _, ok := p.acceptSyntax(`=`, `:=`); ok { 2635 varname, ok := x.(string) 2636 if !ok { 2637 const fs = "expected a variable name, instead of a %T" 2638 return nil, fmt.Errorf(fs, x) 2639 } 2640 2641 x, err := p.parseExpression() 2642 expr := assignExpr{name: varname, expr: x} 2643 return expr, err 2644 } 2645 2646 // handle and/or logical chains 2647 for { 2648 if op, ok := p.acceptSyntax(`&&`, `||`, `&`, `|`); ok { 2649 y, err := p.parseExpression() 2650 if err != nil { 2651 return y, err 2652 } 2653 x = binaryExpr{op: op, x: x, y: y} 2654 continue 2655 } 2656 break 2657 } 2658 2659 // handle maybe-properties 2660 if _, ok := p.acceptSyntax(`??`); ok { 2661 y, err := p.parseExpression() 2662 expr := callExpr{name: "??", args: []any{x, y}} 2663 return expr, err 2664 } 2665 2666 // handle choice/ternary operator 2667 if _, ok := p.acceptSyntax(`?`); ok { 2668 y, err := p.parseExpression() 2669 if err != nil { 2670 expr := callExpr{name: "?:", args: []any{x, nil, nil}} 2671 return expr, err 2672 } 2673 2674 if _, ok := p.acceptSyntax(`:`); ok { 2675 z, err := p.parseExpression() 2676 expr := callExpr{name: "?:", args: []any{x, y, z}} 2677 return expr, err 2678 } 2679 2680 if len(p.tokens) == 0 { 2681 expr := callExpr{name: "?:", args: []any{x, y, nil}} 2682 return expr, errors.New("expected `:`") 2683 } 2684 2685 s := p.tokens[0].value 2686 expr := callExpr{name: "?:", args: []any{x, y, nil}} 2687 err = fmt.Errorf("expected `:`, but got %q instead", s) 2688 return expr, err 2689 } 2690 2691 // expression was just a comparison, or simpler 2692 return x, nil 2693 } 2694 2695 func (p *parser) parseComparison() (any, error) { 2696 x, err := p.parseTerm() 2697 if err != nil { 2698 return x, err 2699 } 2700 2701 op, ok := p.acceptSyntax(`==`, `!=`, `<`, `>`, `<=`, `>=`, `<>`, `===`, `!==`) 2702 if ok { 2703 y, err := p.parseTerm() 2704 return binaryExpr{op: op, x: x, y: y}, err 2705 } 2706 return x, err 2707 } 2708 2709 // parseBinary handles binary operations, by recursing depth-first on the left 2710 // side of binary expressions; going tail-recursive on these would reverse the 2711 // order of arguments instead, which is obviously wrong 2712 func (p *parser) parseBinary(parse func() (any, error), syntax ...string) (any, error) { 2713 x, err := parse() 2714 if err != nil { 2715 return x, err 2716 } 2717 2718 for { 2719 op, ok := p.acceptSyntax(syntax...) 2720 if !ok { 2721 return x, nil 2722 } 2723 2724 y, err := parse() 2725 x = binaryExpr{op: op, x: x, y: y} 2726 if err != nil { 2727 return x, err 2728 } 2729 } 2730 } 2731 2732 func (p *parser) parseTerm() (any, error) { 2733 return p.parseBinary(func() (any, error) { 2734 return p.parseProduct() 2735 }, `+`, `-`, `^`) 2736 } 2737 2738 func (p *parser) parseProduct() (any, error) { 2739 return p.parseBinary(func() (any, error) { 2740 return p.parsePower() 2741 }, `*`, `/`, `%`) 2742 } 2743 2744 func (p *parser) parsePower() (any, error) { 2745 return p.parseBinary(func() (any, error) { 2746 return p.parseValue() 2747 }, `**`, `^`) 2748 } 2749 2750 func (p *parser) parseValue() (any, error) { 2751 // handle unary operators which can also be considered part of numeric 2752 // literals, and thus should be simplified away 2753 if op, ok := p.acceptSyntax(`+`, `-`); ok { 2754 if s, ok := p.acceptNumeric(); ok { 2755 x, err := strconv.ParseFloat(s, 64) 2756 if err != nil { 2757 return nil, err 2758 } 2759 if simpler, ok := simplifyNumber(op, x); ok { 2760 return simpler, nil 2761 } 2762 return unaryExpr{op: op, x: x}, nil 2763 } 2764 2765 x, err := p.parsePower() 2766 return unaryExpr{op: op, x: x}, err 2767 } 2768 2769 // handle all other unary operators 2770 if op, ok := p.acceptSyntax(`!`, `&`, `*`, `^`); ok { 2771 x, err := p.parsePower() 2772 return unaryExpr{op: op, x: x}, err 2773 } 2774 2775 // handle subexpression in parentheses 2776 if _, ok := p.acceptSyntax(`(`); ok { 2777 x, err := p.parseExpression() 2778 if err != nil { 2779 return x, err 2780 } 2781 2782 if err := p.demandSyntax(`)`); err != nil { 2783 return x, err 2784 } 2785 return p.parseAccessors(x) 2786 } 2787 2788 // handle subexpression in square brackets: it's just an alternative to 2789 // using parentheses for subexpressions 2790 if _, ok := p.acceptSyntax(`[`); ok { 2791 x, err := p.parseExpression() 2792 if err != nil { 2793 return x, err 2794 } 2795 2796 if err := p.demandSyntax(`]`); err != nil { 2797 return x, err 2798 } 2799 return p.parseAccessors(x) 2800 } 2801 2802 // handle all other cases 2803 x, err := p.parseSimpleValue() 2804 if err != nil { 2805 return x, err 2806 } 2807 2808 // handle arbitrarily-long chains of accessors 2809 return p.parseAccessors(x) 2810 } 2811 2812 // parseSimpleValue handles a numeric literal or a variable/func name, also 2813 // known as identifier 2814 func (p *parser) parseSimpleValue() (any, error) { 2815 if len(p.tokens) == 0 { 2816 return nil, errEOS 2817 } 2818 t := p.tokens[0] 2819 2820 switch t.kind { 2821 case identifierToken: 2822 p.advance() 2823 // handle func calls, such as f(...) 2824 if _, ok := p.acceptSyntax(`(`); ok { 2825 args, err := p.parseList(`)`) 2826 expr := callExpr{name: t.value, args: args} 2827 return expr, err 2828 } 2829 // handle func calls, such as f[...] 2830 if _, ok := p.acceptSyntax(`[`); ok { 2831 args, err := p.parseList(`]`) 2832 expr := callExpr{name: t.value, args: args} 2833 return expr, err 2834 } 2835 return t.value, nil 2836 2837 case numberToken: 2838 p.advance() 2839 return strconv.ParseFloat(t.value, 64) 2840 2841 default: 2842 const fs = "unexpected %s (token type %d)" 2843 return nil, fmt.Errorf(fs, t.value, t.kind) 2844 } 2845 } 2846 2847 // parseAccessors handles an arbitrarily-long chain of accessors 2848 func (p *parser) parseAccessors(x any) (any, error) { 2849 for { 2850 s, ok := p.acceptSyntax(`.`, `@`) 2851 if !ok { 2852 // dot-chain is over 2853 return x, nil 2854 } 2855 2856 // handle property/method accessors 2857 v, err := p.parseDot(s, x) 2858 if err != nil { 2859 return v, err 2860 } 2861 x = v 2862 } 2863 } 2864 2865 // parseDot handles what follows a syntactic dot, as opposed to a dot which 2866 // may be part of a numeric literal 2867 func (p *parser) parseDot(after string, x any) (any, error) { 2868 if len(p.tokens) == 0 { 2869 const fs = "unexpected end of source after a %q" 2870 return x, fmt.Errorf(fs, after) 2871 } 2872 2873 t := p.tokens[0] 2874 p.advance() 2875 if t.kind != identifierToken { 2876 const fs = "expected a valid property name, but got %s instead" 2877 return x, fmt.Errorf(fs, t.value) 2878 } 2879 2880 if _, ok := p.acceptSyntax(`(`); ok { 2881 items, err := p.parseList(`)`) 2882 args := append([]any{x}, items...) 2883 return callExpr{name: t.value, args: args}, err 2884 } 2885 2886 if _, ok := p.acceptSyntax(`[`); ok { 2887 items, err := p.parseList(`]`) 2888 args := append([]any{x}, items...) 2889 return callExpr{name: t.value, args: args}, err 2890 } 2891 2892 return callExpr{name: t.value, args: []any{x}}, nil 2893 } 2894 2895 // parseList handles the argument-list following a `(` or a `[` 2896 func (p *parser) parseList(end string) ([]any, error) { 2897 var arr []any 2898 for len(p.tokens) > 0 { 2899 if _, ok := p.acceptSyntax(`,`); ok { 2900 // ensure extra/trailing commas are allowed/ignored 2901 continue 2902 } 2903 2904 if _, ok := p.acceptSyntax(end); ok { 2905 return arr, nil 2906 } 2907 2908 v, err := p.parseExpression() 2909 if err != nil { 2910 return arr, err 2911 } 2912 arr = append(arr, v) 2913 } 2914 2915 // return an appropriate error for the unexpected end of the source 2916 return arr, p.demandSyntax(`)`) 2917 } 2918 2919 // simplifyNumber tries to simplify a few trivial unary operations on 2920 // numeric constants 2921 func simplifyNumber(op string, x any) (v any, ok bool) { 2922 if x, ok := x.(float64); ok { 2923 switch op { 2924 case `+`: 2925 return x, true 2926 case `-`: 2927 return -x, true 2928 default: 2929 return x, false 2930 } 2931 } 2932 return x, false 2933 } 2934 2935 // So far, my apps relying on this package are 2936 // 2937 // - fh, a Function Heatmapper which color-codes f(x, y), seen from above 2938 // - star, a STAtistical Resampler to calculate stats from tabular data 2939 // - waveout, an app to emit sample-by-sample wave-audio by formula 2940 2941 // Quick Notes on Performance 2942 // 2943 // Note: while these microbenchmarks are far from proper benchmarks, Club Pulse 2944 // can still give a general idea of what relative performance to expect. 2945 // 2946 // While funscripts.Interpreter can do anything a Program can do, and much more, 2947 // a fmscripts.Program is way faster for float64-only tasks: typical speed-ups 2948 // are between 20-to-50 times. Various simple benchmarks suggest running speed 2949 // is between 1/2 and 1/4 of native funcs. 2950 // 2951 // Reimplementations of the Club Pulse benchmark, when run 50 million times, 2952 // suggest this package is starting to measure the relative speed of various 2953 // trigonometric func across JIT-ted script runners, running somewhat slower 2954 // than NodeJS/V8, faster than PyPy, and much faster than Python. 2955 // 2956 // Other scripted tests, which aren't as trigonometry/exponential-heavy, hint 2957 // at even higher speed-ups compared to Python and PyPy, as well as being 2958 // almost as fast as NodeJS/V8. 2959 // 2960 // Being so close to Node's V8 (0.6x - 0.8x its speed) is a surprisingly good 2961 // result for the relatively little effort this package took. 2962 2963 const ( 2964 maxArgsLen = 1<<16 - 1 // max length for the []float64 input of callv 2965 maxOpIndex = 1<<32 - 1 // max index for values 2966 ) 2967 2968 // types to keep compiler code the same, while changing sizes of numOp fields 2969 type ( 2970 opcode uint16 2971 opargs uint16 // opcode directly affects max correct value of maxArgsLen 2972 opindex uint32 // opcode directly affects max correct value of maxOpIndex 2973 ) 2974 2975 // numOp is a numeric operation in a program 2976 type numOp struct { 2977 What opcode // what to do 2978 NumArgs opargs // vector-length only used for callv and call1v 2979 Index opindex // index to load a value or call a function 2980 } 2981 2982 // Since the go 1.19 compiler started compiling dense switch statements using 2983 // jump tables, adding more basic ops doesn't slow things down anymore when 2984 // reaching the next power-of-2 thresholds in the number of cases. 2985 2986 const ( 2987 // load float64 value to top of the stack 2988 load opcode = iota 2989 2990 // pop top of stack and store it into value 2991 store opcode = iota 2992 2993 // unary operations 2994 neg opcode = iota 2995 not opcode = iota 2996 abs opcode = iota 2997 square opcode = iota 2998 // 1/x, the reciprocal/inverse value 2999 rec opcode = iota 3000 // x%1, but faster than math.Mod(x, 1) 3001 mod1 opcode = iota 3002 3003 // arithmetic and logic operations: 2 float64 inputs, 1 float64 output 3004 3005 add opcode = iota 3006 sub opcode = iota 3007 mul opcode = iota 3008 div opcode = iota 3009 mod opcode = iota 3010 pow opcode = iota 3011 and opcode = iota 3012 or opcode = iota 3013 3014 // binary comparisons: 2 float64 inputs, 1 booleanized float64 output 3015 3016 equal opcode = iota 3017 notequal opcode = iota 3018 less opcode = iota 3019 lessoreq opcode = iota 3020 more opcode = iota 3021 moreoreq opcode = iota 3022 3023 // function callers with 1..5 float64 inputs and a float64 result 3024 3025 call0 opcode = iota 3026 call1 opcode = iota 3027 call2 opcode = iota 3028 call3 opcode = iota 3029 call4 opcode = iota 3030 call5 opcode = iota 3031 3032 // var-arg input function callers 3033 3034 callv opcode = iota 3035 call1v opcode = iota 3036 ) 3037 3038 // Program runs a sequence of float64 operations, with no explicit control-flow: 3039 // implicit control-flow is available in the float64-only functions you make 3040 // available to such programs. Such custom funcs must all return a float64, and 3041 // either take from 0 to 5 float64 arguments, or a single float64 array. 3042 // 3043 // As the name suggests, don't create such objects directly, but instead use 3044 // Compiler.Compile to create them. Only a Compiler lets you register variables 3045 // and functions, which then become part of your numeric Program. 3046 // 3047 // A Program lets you change values before each run, using pointers from method 3048 // Program.Get: such pointers are guaranteed never to change before or across 3049 // runs. Just ensure you get a variable defined in the Compiler used to make the 3050 // Program, or the pointer will be to a dummy place which has no effect on final 3051 // results. 3052 // 3053 // Expressions using literals are automatically optimized into their results: 3054 // this also applies to func calls with standard math functions and constants. 3055 // 3056 // The only way to limit such optimizations is to redefine common math funcs and 3057 // const values explicitly when compiling: after doing so, all those value-names 3058 // will stand for externally-updatable values which can change from one run to 3059 // the next. Similarly, there's no guarantee the (re)defined functions will be 3060 // deterministic, like the defaults they replaced. 3061 // 3062 // # Example 3063 // 3064 // var c fmscripts.Compiler 3065 // 3066 // defs := map[string]any{ 3067 // "x": 0, // define `x`, and initialize it to 0 3068 // "k": 4.25, // define `k`, and initialize it to 4.25 3069 // "b": true, // define `b`, and initialize it to 1.0 3070 // "n": -23, // define `n`, and initialize it to -23.0 3071 // "pi": 3, // define `pi`, overriding the automatic constant named `pi` 3072 // 3073 // "f": numericKernel // type is func ([]float64) float64 3074 // "g": otherFunc // type is func (float64) float64 3075 // } 3076 // 3077 // prog, err := c.Compile("log10(k) + f(sqrt(k) * exp(-x), 45, -0.23)", defs) 3078 // // ... 3079 // 3080 // x, _ := prog.Get("x") // Get returns (*float64, bool) 3081 // y, _ := prog.Get("y") // a useless pointer, since program doesn't use `y` 3082 // // ... 3083 // 3084 // for i := 0; i < n; i++ { 3085 // *x = float64(i)*dx + minx // you update inputs in place using pointers 3086 // f := prog.Run() // method Run gives you a float64 back 3087 // // ... 3088 // } 3089 type Program struct { 3090 sp int // stack pointer, even though it's an index 3091 3092 stack []float64 // pre-allocated by compiler to max length needed by program 3093 values []float64 // holds all values, whether literals, or variables 3094 3095 ops []numOp // all sequential operations for each run 3096 funcs []any // all funcs used by the program 3097 3098 names map[string]int // variable-name to index lookup 3099 dummy float64 // all pointers for undefined variables point here 3100 3101 // data []float64 3102 } 3103 3104 // Get lets you change parameters/variables before each time you run a program, 3105 // since it doesn't return a value, but a pointer to it, so you can update it 3106 // in place. 3107 // 3108 // If the name given isn't available, the result is a pointer to a dummy place: 3109 // this ensures non-nil pointers, which are always safe to use, even though 3110 // updates to the dummy destination have no effect on program results. 3111 func (p *Program) Get(name string) (ptr *float64, useful bool) { 3112 if i, ok := p.names[name]; ok { 3113 return &p.values[i], true 3114 } 3115 return &p.dummy, false 3116 } 3117 3118 // Clone creates an exact copy of a Program with all values in their current 3119 // state: this is useful when dispatching embarassingly-parallel tasks to 3120 // multiple programs. 3121 func (p Program) Clone() Program { 3122 // can't share the stack nor the values 3123 stack := make([]float64, len(p.stack)) 3124 values := make([]float64, len(p.values)) 3125 copy(stack, p.stack) 3126 copy(values, p.values) 3127 p.stack = stack 3128 p.values = values 3129 3130 // can share everything else as is 3131 return p 3132 } 3133 3134 // Memo: the command to show all bound checks is 3135 // go test -gcflags="-d=ssa/check_bce/debug=1" fmscripts/programs.go 3136 3137 // Discussion about go compiler optimizing switch statements into jump tables 3138 // https://go-review.googlesource.com/c/go/+/357330/ 3139 3140 // Run executes the program once. Before each run, update input values using 3141 // pointers from method Get. 3142 func (p *Program) Run() float64 { 3143 // Check for empty programs: these happen either when a compilation error 3144 // was ignored, or when a program was explicitly declared as a variable. 3145 if len(p.ops) == 0 { 3146 return math.NaN() 3147 } 3148 3149 p.sp = 0 3150 p.runAllOps() 3151 return p.stack[p.sp] 3152 } 3153 3154 type func4 = func(float64, float64, float64, float64) float64 3155 type func5 = func(float64, float64, float64, float64, float64) float64 3156 3157 // runAllOps runs all operations in a loop: when done, the program's result is 3158 // ready as the only item left in the stack 3159 func (p *Program) runAllOps() { 3160 for _, op := range p.ops { 3161 // shortcut for the current stack pointer 3162 sp := p.sp 3163 3164 // Preceding binary ops and func calls by _ = p.stack[i-n] prevents 3165 // an extra bound check for the lhs of assignments, but a quick 3166 // statistical summary of benchmarks doesn't show clear speedups, 3167 // let alone major ones. 3168 // 3169 // Separating different func types into different arrays to avoid 3170 // type-checking at runtime doesn't seem to be worth it either, 3171 // and makes the compiler more complicated. 3172 3173 switch op.What { 3174 case load: 3175 // store above top of stack 3176 p.stack[sp+1] = p.values[op.Index] 3177 p.sp++ 3178 case store: 3179 // store from top of the stack 3180 p.values[op.Index] = p.stack[sp] 3181 p.sp-- 3182 3183 // unary operations 3184 case neg: 3185 p.stack[sp] = -p.stack[sp] 3186 case not: 3187 p.stack[sp] = deboolNot(p.stack[sp]) 3188 case abs: 3189 p.stack[sp] = math.Abs(p.stack[sp]) 3190 case square: 3191 p.stack[sp] *= p.stack[sp] 3192 case rec: 3193 p.stack[sp] = 1 / p.stack[sp] 3194 3195 // binary arithmetic ops 3196 case add: 3197 p.stack[sp-1] += p.stack[sp] 3198 p.sp-- 3199 case sub: 3200 p.stack[sp-1] -= p.stack[sp] 3201 p.sp-- 3202 case mul: 3203 p.stack[sp-1] *= p.stack[sp] 3204 p.sp-- 3205 case div: 3206 p.stack[sp-1] /= p.stack[sp] 3207 p.sp-- 3208 case mod: 3209 p.stack[sp-1] = math.Mod(p.stack[sp-1], p.stack[sp]) 3210 p.sp-- 3211 case pow: 3212 p.stack[sp-1] = math.Pow(p.stack[sp-1], p.stack[sp]) 3213 p.sp-- 3214 3215 // binary boolean ops / binary comparisons 3216 case and: 3217 p.stack[sp-1] = deboolAnd(p.stack[sp-1], p.stack[sp]) 3218 p.sp-- 3219 case or: 3220 p.stack[sp-1] = deboolOr(p.stack[sp-1], p.stack[sp]) 3221 p.sp-- 3222 case equal: 3223 p.stack[sp-1] = debool(p.stack[sp-1] == p.stack[sp]) 3224 p.sp-- 3225 case notequal: 3226 p.stack[sp-1] = debool(p.stack[sp-1] != p.stack[sp]) 3227 p.sp-- 3228 case less: 3229 p.stack[sp-1] = debool(p.stack[sp-1] < p.stack[sp]) 3230 p.sp-- 3231 case lessoreq: 3232 p.stack[sp-1] = debool(p.stack[sp-1] <= p.stack[sp]) 3233 p.sp-- 3234 case more: 3235 p.stack[sp-1] = debool(p.stack[sp-1] > p.stack[sp]) 3236 p.sp-- 3237 case moreoreq: 3238 p.stack[sp-1] = debool(p.stack[sp-1] >= p.stack[sp]) 3239 p.sp-- 3240 3241 // function calls 3242 case call0: 3243 f := p.funcs[op.Index].(func() float64) 3244 // store above top of stack 3245 p.stack[sp+1] = f() 3246 p.sp++ 3247 case call1: 3248 f := p.funcs[op.Index].(func(float64) float64) 3249 p.stack[sp] = f(p.stack[sp]) 3250 case call2: 3251 f := p.funcs[op.Index].(func(float64, float64) float64) 3252 p.stack[sp-1] = f(p.stack[sp-1], p.stack[sp]) 3253 p.sp-- 3254 case call3: 3255 f := p.funcs[op.Index].(func(float64, float64, float64) float64) 3256 p.stack[sp-2] = f(p.stack[sp-2], p.stack[sp-1], p.stack[sp]) 3257 p.sp -= 2 3258 case call4: 3259 f := p.funcs[op.Index].(func4) 3260 st := p.stack 3261 p.stack[sp-3] = f(st[sp-3], st[sp-2], st[sp-1], st[sp]) 3262 p.sp -= 3 3263 case call5: 3264 f := p.funcs[op.Index].(func5) 3265 st := p.stack 3266 p.stack[sp-4] = f(st[sp-4], st[sp-3], st[sp-2], st[sp-1], st[sp]) 3267 p.sp -= 4 3268 case callv: 3269 i := sp - int(op.NumArgs) + 1 3270 f := p.funcs[op.Index].(func(...float64) float64) 3271 p.stack[sp-i+1] = f(p.stack[i : sp+1]...) 3272 p.sp = sp - i + 1 3273 case call1v: 3274 i := sp - int(op.NumArgs) + 1 3275 f := p.funcs[op.Index].(func(float64, ...float64) float64) 3276 p.stack[sp-i+1] = f(p.stack[i], p.stack[i+1:sp+1]...) 3277 p.sp = sp - i + 1 3278 } 3279 } 3280 } 3281 3282 // debool is only used to turn boolean values used in comparison operations 3283 // into float64s, since those are the only type accepted on a program stack 3284 func debool(b bool) float64 { 3285 if b { 3286 return 1 3287 } 3288 return 0 3289 } 3290 3291 // deboolNot runs the basic `not` operation 3292 func deboolNot(x float64) float64 { 3293 return debool(x == 0) 3294 } 3295 3296 // deboolAnd runs the basic `and` operation 3297 func deboolAnd(x, y float64) float64 { 3298 return debool(x != 0 && y != 0) 3299 } 3300 3301 // deboolOr runs the basic `or` operation 3302 func deboolOr(x, y float64) float64 { 3303 return debool(x != 0 || y != 0) 3304 } 3305 3306 // These funcs are kept separate from the larger lookup table so they aren't 3307 // mistaken for deterministic funcs which can be optimized away. 3308 // 3309 // var randomFuncs = map[string]any{ 3310 // "rand": Random, 3311 // "rint": RandomInt, 3312 // "runif": RandomUnif, 3313 // "rexp": RandomExp, 3314 // "rnorm": RandomNorm, 3315 // "rgamma": RandomGamma, 3316 // "rbeta": RandomBeta, 3317 // } 3318 3319 func Random(r *rand.Rand) float64 { 3320 return r.Float64() 3321 } 3322 3323 func RandomInt(r *rand.Rand, min, max float64) float64 { 3324 fmin := math.Trunc(min) 3325 fmax := math.Trunc(max) 3326 if fmin == fmax { 3327 return fmin 3328 } 3329 3330 diff := math.Abs(fmax - fmin) 3331 return float64(r.Intn(int(diff)+1)) + fmin 3332 } 3333 3334 func RandomUnif(r *rand.Rand, min, max float64) float64 { 3335 return (max-min)*r.Float64() + min 3336 } 3337 3338 func RandomExp(r *rand.Rand, scale float64) float64 { 3339 return scale * r.ExpFloat64() 3340 } 3341 3342 func RandomNorm(r *rand.Rand, mu, sigma float64) float64 { 3343 return sigma*r.NormFloat64() + mu 3344 } 3345 3346 // Gamma generates a gamma-distributed real value, using a scale parameter. 3347 // 3348 // The algorithm is from Marsaglia and Tsang, as described in 3349 // 3350 // A simple method for generating gamma variables 3351 // https://dl.acm.org/doi/10.1145/358407.358414 3352 func RandomGamma(r *rand.Rand, scale float64) float64 { 3353 d := scale - 1.0/3.0 3354 c := 1 / math.Sqrt(9/d) 3355 3356 for { 3357 // generate candidate value 3358 var x, v float64 3359 for { 3360 x = r.NormFloat64() 3361 v = 1 + c*x 3362 if v > 0 { 3363 break 3364 } 3365 } 3366 v = v * v * v 3367 3368 // accept or reject candidate value 3369 x2 := x * x 3370 u := r.Float64() 3371 if u < 1-0.0331*x2*x2 { 3372 return d * v 3373 } 3374 if math.Log(u) < 0.5*x2+d*(1-v+math.Log(v)) { 3375 return d * v 3376 } 3377 } 3378 } 3379 3380 // Beta generates a beta-distributed real value. 3381 func RandomBeta(r *rand.Rand, a, b float64) float64 { 3382 return RandomGamma(r, a) / (RandomGamma(r, a) + RandomGamma(r, b)) 3383 } 3384 3385 // tokenType is the specific type of a token; tokens never represent 3386 // whitespace-like text between recognized tokens 3387 type tokenType int 3388 3389 const ( 3390 // default zero-value type for tokens 3391 unknownToken tokenType = iota 3392 3393 // a name 3394 identifierToken tokenType = iota 3395 3396 // a literal numeric value 3397 numberToken tokenType = iota 3398 3399 // any syntactic element made of 1 or more runes 3400 syntaxToken tokenType = iota 3401 ) 3402 3403 // errEOS signals the end of source code, and is the only token-related error 3404 // which the parser should ignore 3405 var errEOS = errors.New(`no more source code`) 3406 3407 // token is either a name, value, or syntactic element coming from a script's 3408 // source code 3409 type token struct { 3410 kind tokenType 3411 value string 3412 line int 3413 pos int 3414 } 3415 3416 // tokenizer splits a string into a stream tokens, via its `next` method 3417 type tokenizer struct { 3418 cur string 3419 linenum int 3420 linepos int 3421 } 3422 3423 // newTokenizer is the constructor for type tokenizer 3424 func newTokenizer(src string) tokenizer { 3425 return tokenizer{ 3426 cur: src, 3427 linenum: 1, 3428 linepos: 1, 3429 } 3430 } 3431 3432 // next advances the tokenizer, giving back a token, unless it's done 3433 func (t *tokenizer) next() (token, error) { 3434 // label to allow looping back after skipping comments, and thus avoid 3435 // an explicit tail-recursion for each commented line 3436 rerun: 3437 3438 // always ignore any whitespace-like source 3439 if err := t.skipWhitespace(); err != nil { 3440 return token{}, err 3441 } 3442 3443 if len(t.cur) == 0 { 3444 return token{}, errEOS 3445 } 3446 3447 // remember starting position, in case of error 3448 line := t.linenum 3449 pos := t.linepos 3450 3451 // use the leading rune to probe what's next 3452 r, _ := utf8.DecodeRuneInString(t.cur) 3453 3454 switch r { 3455 case '0', '1', '2', '3', '4', '5', '6', '7', '8', '9': 3456 s, err := t.scanNumber() 3457 return token{kind: numberToken, value: s, line: line, pos: pos}, err 3458 3459 case '(', ')', '[', ']', '{', '}', ',', '+', '-', '%', '^', '~', '@', ';': 3460 s := t.cur[:1] 3461 t.cur = t.cur[1:] 3462 t.linepos++ 3463 res := token{kind: syntaxToken, value: s, line: line, pos: pos} 3464 return res, t.eos() 3465 3466 case ':': 3467 s, err := t.tryPrefixes(`:=`, `:`) 3468 return token{kind: syntaxToken, value: s, line: line, pos: pos}, err 3469 3470 case '*': 3471 s, err := t.tryPrefixes(`**`, `*`) 3472 return token{kind: syntaxToken, value: s, line: line, pos: pos}, err 3473 3474 case '/': 3475 s, err := t.tryPrefixes(`//`, `/*`, `/`) 3476 // double-slash starts a comment until the end of the line 3477 if s == `//` { 3478 t.skipLine() 3479 goto rerun 3480 } 3481 // handle comments which can span multiple lines 3482 if s == `/*` { 3483 err := t.skipComment() 3484 if err != nil { 3485 return token{}, err 3486 } 3487 goto rerun 3488 } 3489 // handle division 3490 return token{kind: syntaxToken, value: s, line: line, pos: pos}, err 3491 3492 case '#': 3493 // even hash starts a comment until the end of the line, making the 3494 // syntax more Python-like 3495 t.skipLine() 3496 goto rerun 3497 3498 case '&': 3499 s, err := t.tryPrefixes(`&&`, `&`) 3500 return token{kind: syntaxToken, value: s, line: line, pos: pos}, err 3501 3502 case '|': 3503 s, err := t.tryPrefixes(`||`, `|`) 3504 return token{kind: syntaxToken, value: s, line: line, pos: pos}, err 3505 3506 case '?': 3507 s, err := t.tryPrefixes(`??`, `?.`, `?`) 3508 return token{kind: syntaxToken, value: s, line: line, pos: pos}, err 3509 3510 case '.': 3511 r, _ := utf8.DecodeRuneInString(t.cur[1:]) 3512 if '0' <= r && r <= '9' { 3513 s, err := t.scanNumber() 3514 res := token{kind: numberToken, value: s, line: line, pos: pos} 3515 return res, err 3516 } 3517 s, err := t.tryPrefixes(`...`, `..`, `.?`, `.`) 3518 return token{kind: syntaxToken, value: s, line: line, pos: pos}, err 3519 3520 case '=': 3521 // triple-equal makes the syntax even more JavaScript-like 3522 s, err := t.tryPrefixes(`===`, `==`, `=>`, `=`) 3523 return token{kind: syntaxToken, value: s, line: line, pos: pos}, err 3524 3525 case '!': 3526 // not-double-equal makes the syntax even more JavaScript-like 3527 s, err := t.tryPrefixes(`!==`, `!=`, `!`) 3528 return token{kind: syntaxToken, value: s, line: line, pos: pos}, err 3529 3530 case '<': 3531 // the less-more/diamond syntax is a SQL-like way to say not equal 3532 s, err := t.tryPrefixes(`<=`, `<>`, `<`) 3533 return token{kind: syntaxToken, value: s, line: line, pos: pos}, err 3534 3535 case '>': 3536 s, err := t.tryPrefixes(`>=`, `>`) 3537 return token{kind: syntaxToken, value: s, line: line, pos: pos}, err 3538 3539 default: 3540 if isIdentStartRune(r) { 3541 s, err := t.scanIdentifier() 3542 res := token{kind: identifierToken, value: s, line: line, pos: pos} 3543 return res, err 3544 } 3545 const fs = `line %d: pos %d: unexpected symbol %c` 3546 return token{}, fmt.Errorf(fs, t.linenum, t.linepos, r) 3547 } 3548 } 3549 3550 // tryPrefixes tries to greedily match any prefix in the order given: when all 3551 // candidates fail, an empty string is returned; when successful, the tokenizer 3552 // updates its state to account for the matched prefix 3553 func (t *tokenizer) tryPrefixes(prefixes ...string) (string, error) { 3554 for _, pre := range prefixes { 3555 if strings.HasPrefix(t.cur, pre) { 3556 t.linepos += len(pre) 3557 t.cur = t.cur[len(pre):] 3558 return pre, t.eos() 3559 } 3560 } 3561 3562 return ``, t.eos() 3563 } 3564 3565 // skipWhitespace updates the tokenizer to ignore runs of consecutive whitespace 3566 // symbols: these are the likes of space, tab, newline, carriage return, etc. 3567 func (t *tokenizer) skipWhitespace() error { 3568 for len(t.cur) > 0 { 3569 r, size := utf8.DecodeRuneInString(t.cur) 3570 if !unicode.IsSpace(r) { 3571 // no more spaces to skip 3572 return nil 3573 } 3574 3575 t.cur = t.cur[size:] 3576 if r == '\n' { 3577 // reached the next line 3578 t.linenum++ 3579 t.linepos = 1 3580 continue 3581 } 3582 // continuing on the same line 3583 t.linepos++ 3584 } 3585 3586 // source code ended 3587 return errEOS 3588 } 3589 3590 // skipLine updates the tokenizer to the end of the current line, or the end of 3591 // the source code, if it's the last line 3592 func (t *tokenizer) skipLine() error { 3593 for len(t.cur) > 0 { 3594 r, size := utf8.DecodeRuneInString(t.cur) 3595 t.cur = t.cur[size:] 3596 if r == '\n' { 3597 // reached the next line, as expected 3598 t.linenum++ 3599 t.linepos = 1 3600 return nil 3601 } 3602 } 3603 3604 // source code ended 3605 t.linenum++ 3606 t.linepos = 1 3607 return errEOS 3608 } 3609 3610 // skipComment updates the tokenizer to the end of the comment started with a 3611 // `/*` and ending with a `*/`, or to the end of the source code 3612 func (t *tokenizer) skipComment() error { 3613 var prev rune 3614 for len(t.cur) > 0 { 3615 r, size := utf8.DecodeRuneInString(t.cur) 3616 t.cur = t.cur[size:] 3617 3618 if r == '\n' { 3619 t.linenum++ 3620 t.linepos = 1 3621 } else { 3622 t.linepos++ 3623 } 3624 3625 if prev == '*' && r == '/' { 3626 return nil 3627 } 3628 prev = r 3629 } 3630 3631 // source code ended 3632 const msg = "comment not ended with a `*/` sequence" 3633 return errors.New(msg) 3634 } 3635 3636 func (t *tokenizer) scanIdentifier() (string, error) { 3637 end := 0 3638 for len(t.cur) > 0 { 3639 r, size := utf8.DecodeRuneInString(t.cur[end:]) 3640 if (end == 0 && !isIdentStartRune(r)) || !isIdentRestRune(r) { 3641 // identifier ended, and there's more source code after it 3642 name := t.cur[:end] 3643 t.cur = t.cur[end:] 3644 return name, nil 3645 } 3646 end += size 3647 t.linepos++ 3648 } 3649 3650 // source code ended with an identifier name 3651 name := t.cur 3652 t.cur = `` 3653 return name, nil 3654 } 3655 3656 func (t *tokenizer) scanNumber() (string, error) { 3657 dots := 0 3658 end := 0 3659 var prev rune 3660 3661 for len(t.cur) > 0 { 3662 r, size := utf8.DecodeRuneInString(t.cur[end:]) 3663 3664 switch r { 3665 case '0', '1', '2', '3', '4', '5', '6', '7', '8', '9', '_', 'e': 3666 end += size 3667 t.linepos++ 3668 prev = r 3669 3670 case '+', '-': 3671 if end > 0 && prev != 'e' { 3672 // number ended, and there's more source code after it 3673 num := t.cur[:end] 3674 t.cur = t.cur[end:] 3675 return num, nil 3676 } 3677 end += size 3678 t.linepos++ 3679 prev = r 3680 3681 case '.': 3682 nr, _ := utf8.DecodeRuneInString(t.cur[end+size:]) 3683 if dots == 1 || isIdentStartRune(nr) || unicode.IsSpace(nr) || nr == '.' { 3684 // number ended, and there's more source code after it 3685 num := t.cur[:end] 3686 t.cur = t.cur[end:] 3687 return num, nil 3688 } 3689 dots++ 3690 end += size 3691 t.linepos++ 3692 prev = r 3693 3694 default: 3695 // number ended, and there's more source code after it 3696 num := t.cur[:end] 3697 t.cur = t.cur[end:] 3698 return num, nil 3699 } 3700 } 3701 3702 // source code ended with a number 3703 name := t.cur 3704 t.cur = `` 3705 return name, nil 3706 } 3707 3708 // eos checks if the source-code is over: if so, it returns an end-of-file error, 3709 // or a nil error otherwise 3710 func (t *tokenizer) eos() error { 3711 if len(t.cur) == 0 { 3712 return errEOS 3713 } 3714 return nil 3715 } 3716 3717 func isIdentStartRune(r rune) bool { 3718 return ('A' <= r && r <= 'Z') || ('a' <= r && r <= 'z') || r == '_' 3719 } 3720 3721 func isIdentRestRune(r rune) bool { 3722 return isIdentStartRune(r) || ('0' <= r && r <= '9') 3723 } 3724 3725 func Decade(x float64) float64 { return 10 * math.Floor(0.1*x) } 3726 func Century(x float64) float64 { return 100 * math.Floor(0.01*x) } 3727 3728 func Round1(x float64) float64 { return 1e-1 * math.Round(1e+1*x) } 3729 func Round2(x float64) float64 { return 1e-2 * math.Round(1e+2*x) } 3730 func Round3(x float64) float64 { return 1e-3 * math.Round(1e+3*x) } 3731 func Round4(x float64) float64 { return 1e-4 * math.Round(1e+4*x) } 3732 func Round5(x float64) float64 { return 1e-5 * math.Round(1e+5*x) } 3733 func Round6(x float64) float64 { return 1e-6 * math.Round(1e+6*x) } 3734 func Round7(x float64) float64 { return 1e-7 * math.Round(1e+7*x) } 3735 func Round8(x float64) float64 { return 1e-8 * math.Round(1e+8*x) } 3736 func Round9(x float64) float64 { return 1e-9 * math.Round(1e+9*x) } 3737 3738 // Ln calculates the natural logarithm. 3739 func Ln(x float64) float64 { 3740 return math.Log(x) 3741 } 3742 3743 // Log calculates logarithms using the base given. 3744 // func Log(base float64, x float64) float64 { 3745 // return math.Log(x) / math.Log(base) 3746 // } 3747 3748 // Round rounds the number given to the number of decimal places given. 3749 func Round(x float64, decimals int) float64 { 3750 k := math.Pow(10, float64(decimals)) 3751 return math.Round(k*x) / k 3752 } 3753 3754 // RoundBy rounds the number given to the unit-size given 3755 // func RoundBy(x float64, by float64) float64 { 3756 // return math.Round(by*x) / by 3757 // } 3758 3759 // FloorBy quantizes a number downward by the unit-size given 3760 // func FloorBy(x float64, by float64) float64 { 3761 // mod := math.Mod(x, by) 3762 // if x < 0 { 3763 // if mod == 0 { 3764 // return x - mod 3765 // } 3766 // return x - mod - by 3767 // } 3768 // return x - mod 3769 // } 3770 3771 // Wrap linearly interpolates the number given to the range [0...1], 3772 // according to its continuous domain, specified by the limits given 3773 func Wrap(x float64, min, max float64) float64 { 3774 return (x - min) / (max - min) 3775 } 3776 3777 // AnchoredWrap works like Wrap, except when both domain boundaries are positive, 3778 // the minimum becomes 0, and when both domain boundaries are negative, the 3779 // maximum becomes 0. 3780 func AnchoredWrap(x float64, min, max float64) float64 { 3781 if min > 0 && max > 0 { 3782 min = 0 3783 } else if max < 0 && min < 0 { 3784 max = 0 3785 } 3786 return (x - min) / (max - min) 3787 } 3788 3789 // Unwrap undoes Wrap, by turning the [0...1] number given into its equivalent 3790 // in the new domain given. 3791 func Unwrap(x float64, min, max float64) float64 { 3792 return (max-min)*x + min 3793 } 3794 3795 // Clamp constrains the domain of the value given 3796 func Clamp(x float64, min, max float64) float64 { 3797 return math.Min(math.Max(x, min), max) 3798 } 3799 3800 // Scale transforms a number according to its position in [xMin, xMax] into its 3801 // correspondingly-positioned counterpart in [yMin, yMax]: if value isn't in its 3802 // assumed domain, its result will be extrapolated accordingly 3803 func Scale(x float64, xmin, xmax, ymin, ymax float64) float64 { 3804 k := (x - xmin) / (xmax - xmin) 3805 return (ymax-ymin)*k + ymin 3806 } 3807 3808 // AnchoredScale works like Scale, except when both domain boundaries are positive, 3809 // the minimum becomes 0, and when both domain boundaries are negative, the maximum 3810 // becomes 0. This allows for proportionally-correct scaling of quantities, such 3811 // as when showing data visually. 3812 func AnchoredScale(x float64, xmin, xmax, ymin, ymax float64) float64 { 3813 if xmin > 0 && xmax > 0 { 3814 xmin = 0 3815 } else if xmax < 0 && xmin < 0 { 3816 xmax = 0 3817 } 3818 k := (x - xmin) / (xmax - xmin) 3819 return (ymax-ymin)*k + ymin 3820 } 3821 3822 // ForceRange is handy when trying to mold floating-point values into numbers 3823 // valid for JSON, since NaN and Infinity are replaced by the values given; 3824 // the infinity-replacement value is negated for negative infinity. 3825 func ForceRange(x float64, nan, inf float64) float64 { 3826 if math.IsNaN(x) { 3827 return nan 3828 } 3829 if math.IsInf(x, -1) { 3830 return -inf 3831 } 3832 if math.IsInf(x, +1) { 3833 return inf 3834 } 3835 return x 3836 } 3837 3838 // Sign returns the standardized sign of a value, either as -1, 0, or +1: NaN 3839 // values stay as NaN, as is expected when using floating-point values. 3840 func Sign(x float64) float64 { 3841 if x > 0 { 3842 return +1 3843 } 3844 if x < 0 { 3845 return -1 3846 } 3847 if x == 0 { 3848 return 0 3849 } 3850 return math.NaN() 3851 } 3852 3853 // IsInteger checks if a floating-point value is already rounded to an integer 3854 // value. 3855 func IsInteger(x float64) bool { 3856 _, frac := math.Modf(x) 3857 return frac == 0 3858 } 3859 3860 func Square(x float64) float64 { 3861 return x * x 3862 } 3863 3864 func Cube(x float64) float64 { 3865 return x * x * x 3866 } 3867 3868 func RoundTo(x float64, decimals float64) float64 { 3869 k := math.Pow10(int(decimals)) 3870 return math.Round(k*x) / k 3871 } 3872 3873 func RelDiff(x, y float64) float64 { 3874 return (x - y) / y 3875 } 3876 3877 // Bool booleanizes a number: 0 stays 0, and anything else becomes 1. 3878 func Bool(x float64) float64 { 3879 if x == 0 { 3880 return 0 3881 } 3882 return 1 3883 } 3884 3885 // DeNaN replaces NaN with the alternative value given: regular values are 3886 // returned as given. 3887 func DeNaN(x float64, instead float64) float64 { 3888 if !math.IsNaN(x) { 3889 return x 3890 } 3891 return instead 3892 } 3893 3894 // DeInf replaces either infinity with the alternative value given: regular 3895 // values are returned as given. 3896 func DeInf(x float64, instead float64) float64 { 3897 if !math.IsInf(x, 0) { 3898 return x 3899 } 3900 return instead 3901 } 3902 3903 // Revalue replaces NaN and either infinity with the alternative value given: 3904 // regular values are returned as given. 3905 func Revalue(x float64, instead float64) float64 { 3906 if !math.IsNaN(x) && !math.IsInf(x, 0) { 3907 return x 3908 } 3909 return instead 3910 } 3911 3912 // Linear evaluates a linear polynomial, the first arg being the main input, 3913 // followed by the polynomial coefficients in decreasing-power order. 3914 func Linear(x float64, a, b float64) float64 { 3915 return a*x + b 3916 } 3917 3918 // Quadratic evaluates a 2nd-degree polynomial, the first arg being the main 3919 // input, followed by the polynomial coefficients in decreasing-power order. 3920 func Quadratic(x float64, a, b, c float64) float64 { 3921 return (a*x+b)*x + c 3922 } 3923 3924 // Cubic evaluates a cubic polynomial, the first arg being the main input, 3925 // followed by the polynomial coefficients in decreasing-power order. 3926 func Cubic(x float64, a, b, c, d float64) float64 { 3927 return ((a*x+b)*x+c)*x + d 3928 } 3929 3930 func LinearFMA(x float64, a, b float64) float64 { 3931 return math.FMA(x, a, b) 3932 } 3933 3934 func QuadraticFMA(x float64, a, b, c float64) float64 { 3935 lin := math.FMA(x, a, b) 3936 return math.FMA(lin, x, c) 3937 } 3938 3939 func CubicFMA(x float64, a, b, c, d float64) float64 { 3940 lin := math.FMA(x, a, b) 3941 quad := math.FMA(lin, x, c) 3942 return math.FMA(quad, x, d) 3943 } 3944 3945 // Radians converts angular degrees into angular radians: 180 degrees are pi 3946 // pi radians. 3947 func Radians(deg float64) float64 { 3948 const k = math.Pi / 180 3949 return k * deg 3950 } 3951 3952 // Degrees converts angular radians into angular degrees: pi radians are 180 3953 // degrees. 3954 func Degrees(rad float64) float64 { 3955 const k = 180 / math.Pi 3956 return k * rad 3957 } 3958 3959 // Fract calculates the non-integer/fractional part of a number. 3960 func Fract(x float64) float64 { 3961 return x - math.Floor(x) 3962 } 3963 3964 // Mix interpolates 2 numbers using a third number, used as an interpolation 3965 // coefficient. This parameter naturally falls in the range [0, 1], but doesn't 3966 // have to be: when given outside that range, the parameter can extrapolate in 3967 // either direction instead. 3968 func Mix(x, y float64, k float64) float64 { 3969 return (1-k)*(y-x) + x 3970 } 3971 3972 // Step implements a step function with a parametric threshold. 3973 func Step(edge, x float64) float64 { 3974 if x < edge { 3975 return 0 3976 } 3977 return 1 3978 } 3979 3980 // SmoothStep is like the `smoothstep` func found in GLSL, using a cubic 3981 // interpolator in the transition region. 3982 func SmoothStep(edge0, edge1, x float64) float64 { 3983 if x <= edge0 { 3984 return 0 3985 } 3986 if x >= edge1 { 3987 return 1 3988 } 3989 3990 // use the cubic hermite interpolator 3x^2 - 2x^3 in the transition band 3991 return x * x * (3 - 2*x) 3992 } 3993 3994 // Logistic approximates the math func of the same name. 3995 func Logistic(x float64) float64 { 3996 return 1 / (1 + math.Exp(-x)) 3997 } 3998 3999 // Sinc approximates the math func of the same name. 4000 func Sinc(x float64) float64 { 4001 if x != 0 { 4002 return math.Sin(x) / x 4003 } 4004 return 1 4005 } 4006 4007 // Sum adds all the numbers in an array. 4008 func Sum(v ...float64) float64 { 4009 s := 0.0 4010 for _, f := range v { 4011 s += f 4012 } 4013 return s 4014 } 4015 4016 // Product multiplies all the numbers in an array. 4017 func Product(v ...float64) float64 { 4018 p := 1.0 4019 for _, f := range v { 4020 p *= f 4021 } 4022 return p 4023 } 4024 4025 // Length calculates the Euclidean length of the vector given. 4026 func Length(v ...float64) float64 { 4027 ss := 0.0 4028 for _, f := range v { 4029 ss += f * f 4030 } 4031 return math.Sqrt(ss) 4032 } 4033 4034 // Dot calculates the dot product of 2 vectors 4035 func Dot(x []float64, y []float64) float64 { 4036 l := len(x) 4037 if len(y) < l { 4038 l = len(y) 4039 } 4040 4041 dot := 0.0 4042 for i := 0; i < l; i++ { 4043 dot += x[i] * y[i] 4044 } 4045 return dot 4046 } 4047 4048 // Min finds the minimum value from the numbers given. 4049 func Min(v ...float64) float64 { 4050 min := +math.Inf(+1) 4051 for _, f := range v { 4052 min = math.Min(min, f) 4053 } 4054 return min 4055 } 4056 4057 // Max finds the maximum value from the numbers given. 4058 func Max(v ...float64) float64 { 4059 max := +math.Inf(-1) 4060 for _, f := range v { 4061 max = math.Max(max, f) 4062 } 4063 return max 4064 } 4065 4066 // Hypot calculates the Euclidean n-dimensional hypothenuse from the numbers 4067 // given: all numbers can be lengths, or simply positional coordinates. 4068 func Hypot(v ...float64) float64 { 4069 sumsq := 0.0 4070 for _, f := range v { 4071 sumsq += f * f 4072 } 4073 return math.Sqrt(sumsq) 4074 } 4075 4076 // Polyval evaluates a polynomial using Horner's algorithm. The array has all 4077 // the coefficients in textbook order, from the highest power down to the final 4078 // constant. 4079 func Polyval(x float64, v ...float64) float64 { 4080 if len(v) == 0 { 4081 return 0 4082 } 4083 4084 x0 := x 4085 x = 1.0 4086 y := 0.0 4087 for i := len(v) - 1; i >= 0; i-- { 4088 y += v[i] * x 4089 x *= x0 4090 } 4091 return y 4092 } 4093 4094 // LnGamma is a 1-input 1-output version of math.Lgamma from the stdlib. 4095 func LnGamma(x float64) float64 { 4096 y, s := math.Lgamma(x) 4097 if s < 0 { 4098 return math.NaN() 4099 } 4100 return y 4101 } 4102 4103 // LnBeta calculates the natural-logarithm of the beta function. 4104 func LnBeta(x float64, y float64) float64 { 4105 return LnGamma(x) + LnGamma(y) - LnGamma(x+y) 4106 } 4107 4108 // Beta calculates the beta function. 4109 func Beta(x float64, y float64) float64 { 4110 return math.Exp(LnBeta(x, y)) 4111 } 4112 4113 // Factorial calculates the product of all integers in [1, n] 4114 func Factorial(n int) int64 { 4115 return int64(math.Round(math.Gamma(float64(n + 1)))) 4116 } 4117 4118 // IsPrime checks whether an integer is bigger than 1 and can only be fully 4119 // divided by 1 and itself, which is the definition of a prime number. 4120 func IsPrime(n int64) bool { 4121 // prime numbers start at 2 4122 if n < 2 { 4123 return false 4124 } 4125 4126 // 2 is the only even prime 4127 if n%2 == 0 { 4128 return n == 2 4129 } 4130 4131 // no divisor can be more than the square root of the target number: 4132 // this limit makes the loop an O(sqrt(n)) one, instead of O(n); this 4133 // is a major algorithmic speedup both in theory and in practice 4134 max := int64(math.Floor(math.Sqrt(float64(n)))) 4135 4136 // the only possible full-divisors are odd integers 3..sqrt(n), 4137 // since reaching this point guarantees n is odd and n > 2 4138 for d := int64(3); d <= max; d += 2 { 4139 if n%d == 0 { 4140 return false 4141 } 4142 } 4143 return true 4144 } 4145 4146 // LCM finds the least common-multiple of 2 positive integers; when one or 4147 // both inputs aren't positive, this func returns 0. 4148 func LCM(x, y int64) int64 { 4149 if gcd := GCD(x, y); gcd > 0 { 4150 return x * y / gcd 4151 } 4152 return 0 4153 } 4154 4155 // GCD finds the greatest common-divisor of 2 positive integers; when one or 4156 // both inputs aren't positive, this func returns 0. 4157 func GCD(x, y int64) int64 { 4158 if x < 1 || y < 1 { 4159 return 0 4160 } 4161 4162 // the loop below requires a >= b 4163 a, b := x, y 4164 if a < b { 4165 a, b = y, x 4166 } 4167 4168 for b > 0 { 4169 a, b = b, a%b 4170 } 4171 return a 4172 } 4173 4174 // Perm counts the number of all possible permutations from n objects when 4175 // picking k times. When one or both inputs aren't positive, the result is 0. 4176 func Perm(n, k int) int64 { 4177 if n < k || n < 0 || k < 0 { 4178 return 0 4179 } 4180 4181 perm := int64(1) 4182 for i := n - k + 1; i <= n; i++ { 4183 perm *= int64(i) 4184 } 4185 return perm 4186 } 4187 4188 // Choose counts the number of all possible combinations from n objects when 4189 // picking k times. When one or both inputs aren't positive, the result is 0. 4190 func Choose(n, k int) int64 { 4191 if n < k || n < 0 || k < 0 { 4192 return 0 4193 } 4194 4195 // the log trick isn't always more accurate when there's no overflow: 4196 // for those cases calculate using the textbook definition 4197 f := math.Round(float64(Perm(n, k) / Factorial(k))) 4198 if !math.IsInf(f, 0) { 4199 return int64(f) 4200 } 4201 4202 // calculate using the log-factorial of n, k, and n - k 4203 a, _ := math.Lgamma(float64(n + 1)) 4204 b, _ := math.Lgamma(float64(k + 1)) 4205 c, _ := math.Lgamma(float64(n - k + 1)) 4206 return int64(math.Round(math.Exp(a - b - c))) 4207 } 4208 4209 // BinomialMass calculates the probability mass of the binomial random process 4210 // given. When the probability given isn't between 0 and 1, the result is NaN. 4211 func BinomialMass(x, n int, p float64) float64 { 4212 // invalid probability input 4213 if p < 0 || p > 1 { 4214 return math.NaN() 4215 } 4216 // events outside the support are impossible by definition 4217 if x < 0 || x > n { 4218 return 0 4219 } 4220 4221 q := 1 - p 4222 m := n - x 4223 ncx := float64(Choose(n, x)) 4224 return ncx * math.Pow(p, float64(x)) * math.Pow(q, float64(m)) 4225 } 4226 4227 // CumulativeBinomialDensity calculates cumulative probabilities/masses up to 4228 // the value given for the binomial random process given. When the probability 4229 // given isn't between 0 and 1, the result is NaN. 4230 func CumulativeBinomialDensity(x, n int, p float64) float64 { 4231 // invalid probability input 4232 if p < 0 || p > 1 { 4233 return math.NaN() 4234 } 4235 if x < 0 { 4236 return 0 4237 } 4238 if x >= n { 4239 return 1 4240 } 4241 4242 p0 := p 4243 q0 := 1 - p0 4244 q := math.Pow(q0, float64(n)) 4245 4246 pbinom := 0.0 4247 np1 := float64(n + 1) 4248 for k := 0; k < x; k++ { 4249 a, _ := math.Lgamma(np1) 4250 b, _ := math.Lgamma(float64(k + 1)) 4251 c, _ := math.Lgamma(float64(n - k + 1)) 4252 // count all possible combinations for this event 4253 ncomb := math.Round(math.Exp(a - b - c)) 4254 pbinom += ncomb * p * q 4255 p *= p0 4256 q /= q0 4257 } 4258 return pbinom 4259 } 4260 4261 // NormalDensity calculates the density at a point along the normal distribution 4262 // given. 4263 func NormalDensity(x float64, mu, sigma float64) float64 { 4264 z := (x - mu) / sigma 4265 return math.Sqrt(0.5/sigma) * math.Exp(-(z*z)/sigma) 4266 } 4267 4268 // CumulativeNormalDensity calculates the probability of a normal variate of 4269 // being up to the value given. 4270 func CumulativeNormalDensity(x float64, mu, sigma float64) float64 { 4271 z := (x - mu) / sigma 4272 return 0.5 + 0.5*math.Erf(z/math.Sqrt2) 4273 } 4274 4275 // Epanechnikov is a commonly-used kernel function. 4276 func Epanechnikov(x float64) float64 { 4277 if math.Abs(x) > 1 { 4278 // func is 0 ouside -1..+1 4279 return 0 4280 } 4281 return 0.75 * (1 - x*x) 4282 } 4283 4284 // Gauss is the commonly-used Gaussian kernel function. 4285 func Gauss(x float64) float64 { 4286 return math.Exp(-(x * x)) 4287 } 4288 4289 // Tricube is a commonly-used kernel function. 4290 func Tricube(x float64) float64 { 4291 a := math.Abs(x) 4292 if a > 1 { 4293 // func is 0 ouside -1..+1 4294 return 0 4295 } 4296 4297 b := a * a * a 4298 c := 1 - b 4299 return 70.0 / 81.0 * c * c * c 4300 } 4301 4302 // SolveQuad finds the solutions of a 2nd-degree polynomial, using a formula 4303 // which is more accurate than the textbook one. 4304 func SolveQuad(a, b, c float64) (x1 float64, x2 float64) { 4305 div := 2 * c 4306 disc := math.Sqrt(b*b - 4*a*c) 4307 x1 = div / (-b - disc) 4308 x2 = div / (-b + disc) 4309 return x1, x2 4310 }