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