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