/* The MIT License (MIT) Copyright © 2020-2025 pacman64 Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. */ /* Single-file source-code for waveout. To compile a smaller-sized command-line app, you can use the `go` command as follows: go build -ldflags "-s -w" -trimpath waveout.go */ package main import ( "bufio" "encoding/base64" "encoding/binary" "errors" "fmt" "io" "math" "math/rand" "os" "strconv" "strings" "time" "unicode" "unicode/utf8" ) // aiff header format // // http://paulbourke.net/dataformats/audio/ // // wav header format // // http://soundfile.sapp.org/doc/WaveFormat/ // http://www-mmsp.ece.mcgill.ca/Documents/AudioFormats/WAVE/WAVE.html // https://docs.fileformat.com/audio/wav/ const ( // maxInt helps convert float64 values into int16 ones maxInt = 1<<15 - 1 // wavIntPCM declares integer PCM sound-data in a wav header wavIntPCM = 1 // wavFloatPCM declares floating-point PCM sound-data in a wav header wavFloatPCM = 3 ) // emitInt16LE writes a 16-bit signed integer in little-endian byte order func emitInt16LE(w io.Writer, f float64) { // binary.Write(w, binary.LittleEndian, int16(maxInt*f)) var buf [2]byte binary.LittleEndian.PutUint16(buf[:2], uint16(int16(maxInt*f))) w.Write(buf[:2]) } // emitFloat32LE writes a 32-bit float in little-endian byte order func emitFloat32LE(w io.Writer, f float64) { var buf [4]byte binary.LittleEndian.PutUint32(buf[:4], math.Float32bits(float32(f))) w.Write(buf[:4]) } // emitInt16BE writes a 16-bit signed integer in big-endian byte order func emitInt16BE(w io.Writer, f float64) { // binary.Write(w, binary.BigEndian, int16(maxInt*f)) var buf [2]byte binary.BigEndian.PutUint16(buf[:2], uint16(int16(maxInt*f))) w.Write(buf[:2]) } // emitFloat32BE writes a 32-bit float in big-endian byte order func emitFloat32BE(w io.Writer, f float64) { var buf [4]byte binary.BigEndian.PutUint32(buf[:4], math.Float32bits(float32(f))) w.Write(buf[:4]) } // wavSettings is an item in the type2wavSettings table type wavSettings struct { Type byte BitsPerSample byte } // type2wavSettings encodes values used when emitting wav headers var type2wavSettings = map[sampleFormat]wavSettings{ int16LE: {wavIntPCM, 16}, float32LE: {wavFloatPCM, 32}, } // emitWaveHeader writes the start of a valid .wav file: since it also starts // the wav data section and emits its size, you only need to write all samples // after calling this func func emitWaveHeader(w io.Writer, cfg outputConfig) error { const fmtChunkSize = 16 duration := cfg.MaxTime numchan := uint32(len(cfg.Scripts)) sampleRate := cfg.SampleRate ws, ok := type2wavSettings[cfg.Samples] if !ok { const fs = `internal error: invalid output-format code %d` return fmt.Errorf(fs, cfg.Samples) } kind := uint16(ws.Type) bps := uint32(ws.BitsPerSample) // byte rate br := sampleRate * bps * numchan / 8 // data size in bytes dataSize := uint32(float64(br) * duration) // total file size totalSize := uint32(dataSize + 44) // general descriptor w.Write([]byte(`RIFF`)) binary.Write(w, binary.LittleEndian, uint32(totalSize)) w.Write([]byte(`WAVE`)) // fmt chunk w.Write([]byte(`fmt `)) binary.Write(w, binary.LittleEndian, uint32(fmtChunkSize)) binary.Write(w, binary.LittleEndian, uint16(kind)) binary.Write(w, binary.LittleEndian, uint16(numchan)) binary.Write(w, binary.LittleEndian, uint32(sampleRate)) binary.Write(w, binary.LittleEndian, uint32(br)) binary.Write(w, binary.LittleEndian, uint16(bps*numchan/8)) binary.Write(w, binary.LittleEndian, uint16(bps)) // start data chunk w.Write([]byte(`data`)) binary.Write(w, binary.LittleEndian, uint32(dataSize)) return nil } // config has all the parsed cmd-line options type config struct { // Scripts has the source codes of all scripts for all channels Scripts []string // To is the output format To string // MaxTime is the play duration of the resulting sound MaxTime float64 // SampleRate is the number of samples per second for all channels SampleRate uint } // parseFlags is the constructor for type config func parseFlags(usage string) (config, error) { cfg := config{ To: `wav`, MaxTime: math.NaN(), SampleRate: 48_000, } args := os.Args[1:] if len(args) == 0 { fmt.Fprint(os.Stderr, usage) os.Exit(0) } for _, s := range args { switch s { case `help`, `-h`, `--h`, `-help`, `--help`: fmt.Fprint(os.Stdout, usage) os.Exit(0) } err := cfg.handleArg(s) if err != nil { return cfg, err } } if math.IsNaN(cfg.MaxTime) { cfg.MaxTime = 1 } if cfg.MaxTime < 0 { const fs = `error: given negative duration %f` return cfg, fmt.Errorf(fs, cfg.MaxTime) } return cfg, nil } func (c *config) handleArg(s string) error { switch s { case `44.1k`, `44.1K`: c.SampleRate = 44_100 return nil case `48k`, `48K`: c.SampleRate = 48_000 return nil case `dat`, `DAT`: c.SampleRate = 48_000 return nil case `cd`, `cda`, `CD`, `CDA`: c.SampleRate = 44_100 return nil } // handle output-format names and their aliases if kind, ok := name2type[s]; ok { c.To = kind return nil } // handle time formats, except when they're pure numbers if math.IsNaN(c.MaxTime) { dur, derr := time.ParseDuration(s) if derr == nil { c.MaxTime = float64(dur) / float64(time.Second) return nil } } // handle sample-rate, given either in hertz or kilohertz lc := strings.ToLower(s) if strings.HasSuffix(lc, `khz`) { lc = strings.TrimSuffix(lc, `khz`) khz, err := strconv.ParseFloat(lc, 64) if err != nil || isBadNumber(khz) || khz <= 0 { const fs = `invalid sample-rate frequency %q` return fmt.Errorf(fs, s) } c.SampleRate = uint(1_000 * khz) return nil } else if strings.HasSuffix(lc, `hz`) { lc = strings.TrimSuffix(lc, `hz`) hz, err := strconv.ParseUint(lc, 10, 64) if err != nil { const fs = `invalid sample-rate frequency %q` return fmt.Errorf(fs, s) } c.SampleRate = uint(hz) return nil } c.Scripts = append(c.Scripts, s) return nil } type encoding byte type headerType byte type sampleFormat byte const ( directEncoding encoding = 1 uriEncoding encoding = 2 noHeader headerType = 1 wavHeader headerType = 2 int16BE sampleFormat = 1 int16LE sampleFormat = 2 float32BE sampleFormat = 3 float32LE sampleFormat = 4 ) // name2type normalizes keys used for type2settings var name2type = map[string]string{ `datauri`: `data-uri`, `dataurl`: `data-uri`, `data-uri`: `data-uri`, `data-url`: `data-uri`, `uri`: `data-uri`, `url`: `data-uri`, `raw`: `raw`, `raw16be`: `raw16be`, `raw16le`: `raw16le`, `raw32be`: `raw32be`, `raw32le`: `raw32le`, `audio/x-wav`: `wave-16`, `audio/x-wave`: `wave-16`, `wav`: `wave-16`, `wave`: `wave-16`, `wav16`: `wave-16`, `wave16`: `wave-16`, `wav-16`: `wave-16`, `wave-16`: `wave-16`, `x-wav`: `wave-16`, `x-wave`: `wave-16`, `wav16uri`: `wave-16-uri`, `wave-16-uri`: `wave-16-uri`, `wav32uri`: `wave-32-uri`, `wave-32-uri`: `wave-32-uri`, `wav32`: `wave-32`, `wave32`: `wave-32`, `wav-32`: `wave-32`, `wave-32`: `wave-32`, } // outputSettings are format-specific settings which are controlled by the // output-format option on the cmd-line type outputSettings struct { Encoding encoding Header headerType Samples sampleFormat } // type2settings translates output-format names into the specific settings // these imply var type2settings = map[string]outputSettings{ ``: {directEncoding, wavHeader, int16LE}, `data-uri`: {uriEncoding, wavHeader, int16LE}, `raw`: {directEncoding, noHeader, int16LE}, `raw16be`: {directEncoding, noHeader, int16BE}, `raw16le`: {directEncoding, noHeader, int16LE}, `wave-16`: {directEncoding, wavHeader, int16LE}, `wave-16-uri`: {uriEncoding, wavHeader, int16LE}, `raw32be`: {directEncoding, noHeader, float32BE}, `raw32le`: {directEncoding, noHeader, float32LE}, `wave-32`: {directEncoding, wavHeader, float32LE}, `wave-32-uri`: {uriEncoding, wavHeader, float32LE}, } // outputConfig has all the info the core of this app needs to make sound type outputConfig struct { // Scripts has the source codes of all scripts for all channels Scripts []string // MaxTime is the play duration of the resulting sound MaxTime float64 // SampleRate is the number of samples per second for all channels SampleRate uint32 // all the configuration details needed to emit output outputSettings } // newOutputConfig is the constructor for type outputConfig, translating the // cmd-line info from type config func newOutputConfig(cfg config) (outputConfig, error) { oc := outputConfig{ Scripts: cfg.Scripts, MaxTime: cfg.MaxTime, SampleRate: uint32(cfg.SampleRate), } if len(oc.Scripts) == 0 { return oc, errors.New(`no formulas given`) } outFmt := strings.ToLower(strings.TrimSpace(cfg.To)) if alias, ok := name2type[outFmt]; ok { outFmt = alias } set, ok := type2settings[outFmt] if !ok { const fs = `unsupported output format %q` return oc, fmt.Errorf(fs, cfg.To) } oc.outputSettings = set return oc, nil } // mimeType gives the format's corresponding MIME type, or an empty string // if the type isn't URI-encodable func (oc outputConfig) mimeType() string { if oc.Header == wavHeader { return `audio/x-wav` } return `` } func isBadNumber(f float64) bool { return math.IsNaN(f) || math.IsInf(f, 0) } const usage = ` waveout [options...] [duration...] [formulas...] This app emits wave-sound binary data using the script(s) given. Scripts give you the float64-related functionality you may expect, from numeric operations to several math functions. When given 1 formula, the result is mono; when given 2 formulas (left and right), the result is stereo, and so on. Output is always uncompressed audio: "waveout" can emit that as is, or as a base64-encoded data-URI, which you can use as a "src" attribute value in an HTML audio tag. Output duration is 1 second by default, but you can change that too by using a recognized time format. The first recognized time format is the familiar hh:mm:ss, where the hours are optional, and where seconds can have a decimal part after it. The second recognized time format uses 1-letter shortcuts instead of colons for each time component, each of which is optional: "h" stands for hour, "m" for minutes, and "s" for seconds. Output Formats encoding header samples endian more info wav direct wave int16 little default format wav16 direct wave int16 little alias for "wav" wav32 direct wave float32 little uri data-URI wave int16 little MIME type is audio/x-wav raw direct none int16 little raw16le direct none int16 little alias for "raw" raw32le direct none float32 little raw16be direct none int16 big raw32be direct none float32 big Concrete Examples # low-tones commonly used in club music as beats waveout 2s 'sin(10 * tau * exp(-20 * u)) * exp(-2 * u)' > club-beats.wav # 1 minute and 5 seconds of static-like random noise waveout 1m5s 'rand()' > random-noise.wav # many bell-like clicks in quick succession; can be a cellphone's ringtone waveout 'sin(2048 * tau * t) * exp(-50 * (t%0.1))' > ringtone.wav # similar to the door-opening sound from a dsc powerseries home alarm waveout 'sin(4096 * tau * t) * exp(-10 * (t%0.1))' > home-alarm.wav # watch your ears: quickly increases frequency up to 2khz waveout 'sin(2_000 * t * tau * t)' > frequency-sweep.wav # 1-second 400hz test tone waveout 'sin(400 * tau * t)' > test-tone-400.wav # 2s of a 440hz test tone, also called an A440 sound waveout 2s 'sin(440 * tau * t)' > a440.wav # 1s 400hz test tone with sudden volume drop at the end, to avoid clip waveout 'sin(400 * tau * t) * min(1, exp(-100*(t-0.9)))' > nice-tone.wav # old ringtone used in north america waveout '0.5*sin(350 * tau * t) + 0.5*sin(450 * tau * t)' > na-ringtone.wav # 20 seconds of periodic pings waveout 20s 'sin(800 * tau * u) * exp(-20 * u)' > pings.wav # 2 seconds of a european-style dial-tone waveout 2s '(sin(350 * tau * t) + sin(450 * tau * t)) / 2' > dial-tone.wav # 4 seconds of a north-american-style busy-phone signal waveout 4s '(u < 0.5) * (sin(480*tau * t) + sin(620*tau * t)) / 2' > na-busy.wav # hit the 51st key on a synthetic piano-like instrument waveout 'sin(tau * 440 * 2**((51 - 49)/12) * t) * exp(-10*u)' > piano-key.wav # hit of a synthetic snare-like sound waveout 'random() * exp(-10 * t)' > synth-snare.wav # a stereotypical "laser" sound waveout 'sin(100 * tau * exp(-40 * t))' > laser.wav ` func main() { cfg, err := parseFlags(usage[1:]) if err != nil { fmt.Fprintln(os.Stderr, err.Error()) os.Exit(1) } oc, err := newOutputConfig(cfg) if err != nil { fmt.Fprintln(os.Stderr, err.Error()) os.Exit(1) } addDetermFuncs() if err := run(oc); err != nil { fmt.Fprintln(os.Stderr, err.Error()) os.Exit(1) } } func run(cfg outputConfig) error { // f, err := os.Create(`waveout.prof`) // if err != nil { // return err // } // defer f.Close() // pprof.StartCPUProfile(f) // defer pprof.StopCPUProfile() w := bufio.NewWriterSize(os.Stdout, 64*1024) defer w.Flush() switch cfg.Encoding { case directEncoding: return runDirect(w, cfg) case uriEncoding: mtype := cfg.mimeType() if mtype == `` { return errors.New(`internal error: no MIME type`) } fmt.Fprintf(w, `data:%s;base64,`, mtype) enc := base64.NewEncoder(base64.StdEncoding, w) defer enc.Close() return runDirect(enc, cfg) default: const fs = `internal error: wrong output-encoding code %d` return fmt.Errorf(fs, cfg.Encoding) } } // type2emitter chooses sample-emitter funcs from the format given var type2emitter = map[sampleFormat]func(io.Writer, float64){ int16LE: emitInt16LE, int16BE: emitInt16BE, float32LE: emitFloat32LE, float32BE: emitFloat32BE, } // runDirect emits sound-data bytes: this func can be called with writers // which keep bytes as given, or with re-encoders, such as base64 writers func runDirect(w io.Writer, cfg outputConfig) error { switch cfg.Header { case noHeader: // do nothing, while avoiding error case wavHeader: emitWaveHeader(w, cfg) default: const fs = `internal error: wrong header code %d` return fmt.Errorf(fs, cfg.Header) } emitter, ok := type2emitter[cfg.Samples] if !ok { const fs = `internal error: wrong output-format code %d` return fmt.Errorf(fs, cfg.Samples) } if len(cfg.Scripts) == 1 { return emitMono(w, cfg, emitter) } return emit(w, cfg, emitter) } // makeDefs makes extra funcs and values available to scripts func makeDefs(cfg outputConfig) map[string]any { // copy extra built-in funcs defs := make(map[string]any, len(extras)+6+5) for k, v := range extras { defs[k] = v } // add extra variables defs[`t`] = 0.0 defs[`u`] = 0.0 defs[`d`] = cfg.MaxTime defs[`dur`] = cfg.MaxTime defs[`duration`] = cfg.MaxTime defs[`end`] = cfg.MaxTime // add pseudo-random funcs seed := time.Now().UnixNano() r := rand.New(rand.NewSource(seed)) rand := func() float64 { return random01(r) } randomf := func() float64 { return random(r) } rexpf := func(scale float64) float64 { return rexp(r, scale) } rnormf := func(mu, sigma float64) float64 { return rnorm(r, mu, sigma) } defs[`rand`] = rand defs[`rand01`] = rand defs[`random`] = randomf defs[`rexp`] = rexpf defs[`rnorm`] = rnormf return defs } type emitFunc = func(io.Writer, float64) // emit runs the formulas given to emit all wave samples func emit(w io.Writer, cfg outputConfig, emit emitFunc) error { var c Compiler defs := makeDefs(cfg) programs := make([]Program, 0, len(cfg.Scripts)) tvars := make([]*float64, 0, len(cfg.Scripts)) uvars := make([]*float64, 0, len(cfg.Scripts)) for _, s := range cfg.Scripts { p, err := c.Compile(s, defs) if err != nil { return err } programs = append(programs, p) t, _ := p.Get(`t`) u, _ := p.Get(`u`) tvars = append(tvars, t) uvars = append(uvars, u) } dt := 1.0 / float64(cfg.SampleRate) end := cfg.MaxTime for i := 0.0; true; i++ { now := dt * i if now >= end { return nil } _, u := math.Modf(now) for j, p := range programs { *tvars[j] = now *uvars[j] = u emit(w, p.Run()) } } return nil } // emitMono runs the formula given to emit all single-channel wave samples func emitMono(w io.Writer, cfg outputConfig, emit emitFunc) error { var c Compiler mono, err := c.Compile(cfg.Scripts[0], makeDefs(cfg)) if err != nil { return err } t, _ := mono.Get(`t`) u, needsu := mono.Get(`u`) dt := 1.0 / float64(cfg.SampleRate) end := cfg.MaxTime // update variable `u` only if script uses it: this can speed things // up considerably when that variable isn't used if needsu { for i := 0.0; true; i++ { now := dt * i if now >= end { return nil } *t = now _, *u = math.Modf(now) emit(w, mono.Run()) } return nil } for i := 0.0; true; i++ { now := dt * i if now >= end { return nil } *t = now emit(w, mono.Run()) } return nil } // tau is exactly 1 loop around a circle, which is handy to turn frequencies // into trigonometric angles, since they're measured in radians const tau = 2 * math.Pi // extras has funcs beyond what the script built-ins offer: those built-ins // are for general math calculations, while these are specific for sound // effects, other sound-related calculations, or to make pseudo-random values var extras = map[string]any{ `hihat`: hihat, } // addDetermFuncs does what it says, ensuring these funcs are optimizable when // they're given all-constant expressions as inputs func addDetermFuncs() { DefineDetFuncs(map[string]any{ `ascale`: AnchoredScale, `awrap`: AnchoredWrap, `clamp`: Clamp, `epa`: Epanechnikov, `epanechnikov`: Epanechnikov, `fract`: Fract, `gauss`: Gauss, `horner`: Polyval, `logistic`: Logistic, `mix`: Mix, `polyval`: Polyval, `scale`: Scale, `sign`: Sign, `sinc`: Sinc, `smoothstep`: SmoothStep, `step`: Step, `tricube`: Tricube, `unwrap`: Unwrap, `wrap`: Wrap, `drop`: dropsince, `dropfrom`: dropsince, `dropoff`: dropsince, `dropsince`: dropsince, `kick`: kick, `kicklow`: kicklow, `piano`: piano, `pianokey`: piano, `pickval`: pickval, `pickvalue`: pickval, `sched`: schedule, `schedule`: schedule, `timeval`: timeval, `timevalues`: timeval, }) } // random01 returns a random value in 0 .. 1 func random01(r *rand.Rand) float64 { return r.Float64() } // random returns a random value in -1 .. +1 func random(r *rand.Rand) float64 { return (2 * r.Float64()) - 1 } // rexp returns an exponentially-distributed random value using the scale // (expected value) given func rexp(r *rand.Rand, scale float64) float64 { return scale * r.ExpFloat64() } // rnorm returns a normally-distributed random value using the mean and // standard deviation given func rnorm(r *rand.Rand, mu, sigma float64) float64 { return r.NormFloat64()*sigma + mu } // make sample for a synthetic-drum kick func kick(t float64, f, k float64) float64 { const p = 0.085 return math.Sin(tau*f*math.Pow(p, t)) * math.Exp(-k*t) } // make sample for a heavier-sounding synthetic-drum kick func kicklow(t float64, f, k float64) float64 { const p = 0.08 return math.Sin(tau*f*math.Pow(p, t)) * math.Exp(-k*t) } // make sample for a synthetic hi-hat hit func hihat(t float64, k float64) float64 { return rand.Float64() * math.Exp(-k*t) } // schedule rearranges time, without being a time machine func schedule(t float64, period, delay float64) float64 { v := t + (1 - delay) if v < 0 { return 0 } return math.Mod(v*period, period) } // make sample for a synthetic piano key being hit func piano(t float64, n float64) float64 { p := (math.Floor(n) - 49) / 12 f := 440 * math.Pow(2, p) return math.Sin(tau * f * t) } // multiply rest of a formula with this for a quick volume drop at the end: // this is handy to avoid clips when sounds end playing func dropsince(t float64, start float64) float64 { // return math.Min(1, math.Exp(-100*(t-start))) if t <= start { return 1 } return math.Exp(-100 * (t - start)) } // pickval requires at least 3 args, the first 2 being the current time and // each slot's duration, respectively: these 2 are followed by all the values // to pick for all time slots func pickval(args ...float64) float64 { if len(args) < 3 { return 0 } t := args[0] slotdur := args[1] values := args[2:] u, _ := math.Modf(t / slotdur) n := len(values) i := int(u) % n if 0 <= i && i < n { return values[i] } return 0 } // timeval requires at least 2 args, the first 2 being the current time and // the total looping-period, respectively: these 2 are followed by pairs of // numbers, each consisting of a timestamp and a matching value, in order func timeval(args ...float64) float64 { if len(args) < 2 { return 0 } t := args[0] period := args[1] u, _ := math.Modf(t / period) // find the first value whose periodic timestamp is due for rest := args[2:]; len(rest) >= 2; rest = rest[2:] { if u >= rest[0]/period { return rest[1] } } return 0 } // unary2op turns unary operators into their corresponding basic operations; // some entries are only for the optimizer, and aren't accessible directly // from valid source code var unary2op = map[string]opcode{ "-": neg, "!": not, "&": abs, "*": square, "^": square, "/": rec, "%": mod1, } // binary2op turns binary operators into their corresponding basic operations var binary2op = map[string]opcode{ "+": add, "-": sub, "*": mul, "/": div, "%": mod, "&&": and, "&": and, "||": or, "|": or, "==": equal, "!=": notequal, "<>": notequal, "<": less, "<=": lessoreq, ">": more, ">=": moreoreq, "**": pow, "^": pow, } // Compiler lets you create Program objects, which you can then run. The whole // point of this is to create quicker-to-run numeric scripts. You can even add // variables with initial values, as well as functions for the script to use. // // Common math funcs and constants are automatically detected, and constant // results are optimized away, unless builtins are redefined. In other words, // the optimizer is effectively disabled for all (sub)expressions containing // redefined built-ins, as there's no way to be sure those values won't change // from one run to the next. // // See the comment for type Program for more details. // // # Example // // var c fmscripts.Compiler // // defs := map[string]any{ // "x": 0, // define `x`, and initialize it to 0 // "k": 4.25, // define `k`, and initialize it to 4.25 // "b": true, // define `b`, and initialize it to 1.0 // "n": -23, // define `n`, and initialize it to -23.0 // "pi": 3, // define `pi`, overriding the default constant named `pi` // // "f": numericKernel // type is func ([]float64) float64 // "g": otherFunc // type is func (float64) float64 // } // // prog, err := c.Compile("log10(k) + f(sqrt(k) * exp(-x), 45, -0.23)", defs) // // ... // // x, _ := prog.Get("x") // Get returns (*float64, bool) // y, _ := prog.Get("y") // a useless pointer, since program doesn't use `y` // // ... // // for i := 0; i < n; i++ { // *x = float64(i)*dx + minx // you update inputs in place using pointers // f := prog.Run() // method Run gives you a float64 back // // ... // } type Compiler struct { maxStack int // the exact stack size the resulting program needs ops []numOp // the program's operations vaddr map[string]int // address lookup for values faddr map[string]int // address lookup for functions ftype map[string]any // keeps track of func types during compilation values []float64 // variables and constants available to programs funcs []any // funcs available to programs } // Compile parses the script given and generates a fast float64-only Program // made only of sequential steps: any custom funcs you provide it can use // their own internal looping and/or conditional logic, of course. func (c *Compiler) Compile(src string, defs map[string]any) (Program, error) { // turn source code into an abstract syntax-tree root, err := parse(src) if err != nil { return Program{}, err } // generate operations if err := c.reset(defs); err != nil { return Program{}, err } if err = c.compile(root, 0); err != nil { return Program{}, err } // create the resulting program var p Program p.stack = make([]float64, c.maxStack) p.values = make([]float64, len(c.values)) copy(p.values, c.values) p.ops = make([]numOp, len(c.ops)) copy(p.ops, c.ops) p.funcs = make([]any, len(c.ftype)) copy(p.funcs, c.funcs) p.names = make(map[string]int, len(c.vaddr)) // give the program's Get method access only to all allocated variables for k, v := range c.vaddr { // avoid exposing numeric constants on the program's name-lookup // table, which would allow users to change literals across runs if _, err := strconv.ParseFloat(k, 64); err == nil { continue } // expose only actual variable names p.names[k] = v } return p, nil } // reset prepares a compiler by satisfying internal preconditions func // compileExpr relies on, when given an abstract syntax-tree func (c *Compiler) reset(defs map[string]any) error { // reset the compiler's internal state c.maxStack = 0 c.ops = c.ops[:0] c.vaddr = make(map[string]int) c.faddr = make(map[string]int) c.ftype = make(map[string]any) c.values = c.values[:0] c.funcs = c.funcs[:0] // allocate vars and funcs for k, v := range defs { if err := c.allocEntry(k, v); err != nil { return err } } return nil } // allocEntry simplifies the control-flow of func reset func (c *Compiler) allocEntry(k string, v any) error { const ( maxExactRound = 2 << 52 rangeErrFmt = "%d is outside range of exact float64 integers" typeErrFmt = "got value of unsupported type %T" ) switch v := v.(type) { case float64: _, err := c.allocValue(k, v) return err case int: if math.Abs(float64(v)) > maxExactRound { return fmt.Errorf(rangeErrFmt, v) } _, err := c.allocValue(k, float64(v)) return err case bool: _, err := c.allocValue(k, debool(v)) return err default: if isSupportedFunc(v) { c.ftype[k] = v _, err := c.allocFunc(k, v) return err } return fmt.Errorf(typeErrFmt, v) } } // allocValue ensures there's a place to match the name given, returning its // index when successful; only programs using unreasonably many values will // cause this func to fail func (c *Compiler) allocValue(s string, f float64) (int, error) { if len(c.values) >= maxOpIndex { const fs = "programs can only use up to %d distinct float64 values" return -1, fmt.Errorf(fs, maxOpIndex+1) } i := len(c.values) c.vaddr[s] = i c.values = append(c.values, f) return i, nil } // valueIndex returns the index reserved for an allocated value/variable: // all unallocated values/variables are allocated here on first access func (c *Compiler) valueIndex(s string, f float64) (int, error) { // name is found: return the index of the already-allocated var if i, ok := c.vaddr[s]; ok { return i, nil } // name not found, but it's a known math constant if f, ok := mathConst[s]; ok { return c.allocValue(s, f) } // name not found, and it's not of a known math constant return c.allocValue(s, f) } // constIndex allocates a constant as a variable named as its own string // representation, which avoids multiple entries for repeated uses of the // same constant value func (c *Compiler) constIndex(f float64) (int, error) { // constants have no name, so use a canonical string representation return c.valueIndex(strconv.FormatFloat(f, 'f', 16, 64), f) } // funcIndex returns the index reserved for an allocated function: all // unallocated functions are allocated here on first access func (c *Compiler) funcIndex(name string) (int, error) { // check if func was already allocated if i, ok := c.faddr[name]; ok { return i, nil } // if name is reserved allocate an index for its matching func if fn, ok := c.ftype[name]; ok { return c.allocFunc(name, fn) } // if name wasn't reserved, see if it's a standard math func name if fn := c.autoFuncLookup(name); fn != nil { return c.allocFunc(name, fn) } // name isn't even a standard math func's return -1, fmt.Errorf("function not found") } // allocFunc ensures there's a place for the function name given, returning its // index when successful; only funcs of unsupported types will cause failure func (c *Compiler) allocFunc(name string, fn any) (int, error) { if isSupportedFunc(fn) { i := len(c.funcs) c.faddr[name] = i c.ftype[name] = fn c.funcs = append(c.funcs, fn) return i, nil } return -1, fmt.Errorf("can't use a %T as a number-crunching function", fn) } // autoFuncLookup checks built-in deterministic funcs for the name given: its // result is nil only if there's no match func (c *Compiler) autoFuncLookup(name string) any { if fn, ok := determFuncs[name]; ok { return fn } return nil } // genOp generates/adds a basic operation to the program, while keeping track // of the maximum depth the stack can reach func (c *Compiler) genOp(op numOp, depth int) { // add 2 defensively to ensure stack space for the inputs of binary ops n := depth + 2 if c.maxStack < n { c.maxStack = n } c.ops = append(c.ops, op) } // compile is a recursive expression evaluator which does the actual compiling // as it goes along func (c *Compiler) compile(expr any, depth int) error { expr = c.optimize(expr) switch expr := expr.(type) { case float64: return c.compileLiteral(expr, depth) case string: return c.compileVariable(expr, depth) case []any: return c.compileCombo(expr, depth) case unaryExpr: return c.compileUnary(expr.op, expr.x, depth) case binaryExpr: return c.compileBinary(expr.op, expr.x, expr.y, depth) case callExpr: return c.compileCall(expr, depth) case assignExpr: return c.compileAssign(expr, depth) default: return fmt.Errorf("unsupported expression type %T", expr) } } // compileLiteral generates a load operation for the constant value given func (c *Compiler) compileLiteral(f float64, depth int) error { i, err := c.constIndex(f) if err != nil { return err } c.genOp(numOp{What: load, Index: opindex(i)}, depth) return nil } // compileVariable generates a load operation for the variable name given func (c *Compiler) compileVariable(name string, depth int) error { // handle names which aren't defined, but are known math constants if _, ok := c.vaddr[name]; !ok { if f, ok := mathConst[name]; ok { return c.compileLiteral(f, depth) } } // handle actual variables i, err := c.valueIndex(name, 0) if err != nil { return err } c.genOp(numOp{What: load, Index: opindex(i)}, depth) return nil } // compileCombo handles a sequence of expressions func (c *Compiler) compileCombo(exprs []any, depth int) error { for _, v := range exprs { err := c.compile(v, depth) if err != nil { return err } } return nil } // compileUnary handles unary expressions func (c *Compiler) compileUnary(op string, x any, depth int) error { err := c.compile(x, depth+1) if err != nil { return err } if op, ok := unary2op[op]; ok { c.genOp(numOp{What: op}, depth) return nil } return fmt.Errorf("unary operation %q is unsupported", op) } // compileBinary handles binary expressions func (c *Compiler) compileBinary(op string, x, y any, depth int) error { switch op { case "===", "!==": // handle binary expressions with no matching basic operation, by // treating them as aliases for function calls return c.compileCall(callExpr{name: op, args: []any{x, y}}, depth) } err := c.compile(x, depth+1) if err != nil { return err } err = c.compile(y, depth+2) if err != nil { return err } if op, ok := binary2op[op]; ok { c.genOp(numOp{What: op}, depth) return nil } return fmt.Errorf("binary operation %q is unsupported", op) } // compileCall handles function-call expressions func (c *Compiler) compileCall(expr callExpr, depth int) error { // lookup function name index, err := c.funcIndex(expr.name) if err != nil { return fmt.Errorf("%s: %s", expr.name, err.Error()) } // get the func value, as its type determines the calling op to emit v, ok := c.ftype[expr.name] if !ok { return fmt.Errorf("%s: function is undefined", expr.name) } // figure which type of function operation to use op, ok := func2op(v) if !ok { return fmt.Errorf("%s: unsupported function type %T", expr.name, v) } // ensure number of args given to the func makes sense for the func type err = checkArgCount(func2info[op], expr.name, len(expr.args)) if err != nil { return err } // generate operations to evaluate all args for i, v := range expr.args { err := c.compile(v, depth+i+1) if err != nil { return err } } // generate func-call operation given := len(expr.args) next := numOp{What: op, NumArgs: opargs(given), Index: opindex(index)} c.genOp(next, depth) return nil } // checkArgCount does what it says, returning informative errors when arg // counts are wrong func checkArgCount(info funcTypeInfo, name string, nargs int) error { if info.AtLeast < 0 && info.AtMost < 0 { return nil } if info.AtLeast == info.AtMost && nargs != info.AtMost { const fs = "%s: expected %d args, got %d instead" return fmt.Errorf(fs, name, info.AtMost, nargs) } if info.AtLeast >= 0 && info.AtMost >= 0 { const fs = "%s: expected between %d and %d args, got %d instead" if nargs < info.AtLeast || nargs > info.AtMost { return fmt.Errorf(fs, name, info.AtLeast, info.AtMost, nargs) } } if info.AtLeast >= 0 && nargs < info.AtLeast { const fs = "%s: expected at least %d args, got %d instead" return fmt.Errorf(fs, name, info.AtLeast, nargs) } if info.AtMost >= 0 && nargs > info.AtMost { const fs = "%s: expected at most %d args, got %d instead" return fmt.Errorf(fs, name, info.AtMost, nargs) } // all is good return nil } // compileAssign generates a store operation for the expression given func (c *Compiler) compileAssign(expr assignExpr, depth int) error { err := c.compile(expr.expr, depth) if err != nil { return err } i, err := c.allocValue(expr.name, 0) if err != nil { return err } c.genOp(numOp{What: store, Index: opindex(i)}, depth) return nil } // func sameFunc(x, y any) bool { // if x == nil || y == nil { // return false // } // return reflect.ValueOf(x).Pointer() == reflect.ValueOf(y).Pointer() // } // isSupportedFunc checks if a value's type is a supported func type func isSupportedFunc(fn any) bool { _, ok := func2op(fn) return ok } // funcTypeInfo is an entry in the func2info lookup table type funcTypeInfo struct { // AtLeast is the minimum number of inputs the func requires: negative // values are meant to be ignored AtLeast int // AtMost is the maximum number of inputs the func requires: negative // values are meant to be ignored AtMost int } // func2info is a lookup table to check the number of inputs funcs are given var func2info = map[opcode]funcTypeInfo{ call0: {AtLeast: +0, AtMost: +0}, call1: {AtLeast: +1, AtMost: +1}, call2: {AtLeast: +2, AtMost: +2}, call3: {AtLeast: +3, AtMost: +3}, call4: {AtLeast: +4, AtMost: +4}, call5: {AtLeast: +5, AtMost: +5}, callv: {AtLeast: -1, AtMost: -1}, call1v: {AtLeast: +1, AtMost: -1}, } // func2op tries to match a func type into a corresponding opcode func func2op(v any) (op opcode, ok bool) { switch v.(type) { case func() float64: return call0, true case func(float64) float64: return call1, true case func(float64, float64) float64: return call2, true case func(float64, float64, float64) float64: return call3, true case func(float64, float64, float64, float64) float64: return call4, true case func(float64, float64, float64, float64, float64) float64: return call5, true case func(...float64) float64: return callv, true case func(float64, ...float64) float64: return call1v, true default: return 0, false } } const ( // the maximum integer a float64 can represent exactly; -maxflint is the // minimum such integer, since floating-point values allow such symmetries maxflint = 2 << 52 // base-2 versions of size multipliers kilobyte = 1024 * 1.0 megabyte = 1024 * kilobyte gigabyte = 1024 * megabyte terabyte = 1024 * gigabyte petabyte = 1024 * terabyte // unit-conversion multipliers mi2kmMult = 1 / 0.6213712 nm2kmMult = 1.852 nmi2kmMult = 1.852 yd2mtMult = 1 / 1.093613 ft2mtMult = 1 / 3.28084 in2cmMult = 2.54 lb2kgMult = 0.4535924 ga2ltMult = 1 / 0.2199692 gal2lMult = 1 / 0.2199692 oz2mlMult = 29.5735295625 cup2lMult = 0.2365882365 mpg2kplMult = 0.2199692 / 0.6213712 ton2kgMult = 1 / 907.18474 psi2paMult = 1 / 0.00014503773800722 deg2radMult = math.Pi / 180 rad2degMult = 180 / math.Pi ) // default math constants var mathConst = map[string]float64{ "e": math.E, "pi": math.Pi, "tau": 2 * math.Pi, "phi": math.Phi, "nan": math.NaN(), "inf": math.Inf(+1), "minint": -float64(maxflint - 1), "maxint": +float64(maxflint - 1), "minsafeint": -float64(maxflint - 1), "maxsafeint": +float64(maxflint - 1), "false": 0.0, "true": 1.0, "f": 0.0, "t": 1.0, // conveniently-named multipliers "femto": 1e-15, "pico": 1e-12, "nano": 1e-09, "micro": 1e-06, "milli": 1e-03, "kilo": 1e+03, "mega": 1e+06, "giga": 1e+09, "tera": 1e+12, "peta": 1e+15, // unit-conversion multipliers "mi2km": mi2kmMult, "nm2km": nm2kmMult, "nmi2km": nmi2kmMult, "yd2mt": yd2mtMult, "ft2mt": ft2mtMult, "in2cm": in2cmMult, "lb2kg": lb2kgMult, "ga2lt": ga2ltMult, "gal2l": gal2lMult, "oz2ml": oz2mlMult, "cup2l": cup2lMult, "mpg2kpl": mpg2kplMult, "ton2kg": ton2kgMult, "psi2pa": psi2paMult, "deg2rad": deg2radMult, "rad2deg": rad2degMult, // base-2 versions of size multipliers "kb": kilobyte, "mb": megabyte, "gb": gigabyte, "tb": terabyte, "pb": petabyte, "binkilo": kilobyte, "binmega": megabyte, "bingiga": gigabyte, "bintera": terabyte, "binpeta": petabyte, // physical constants "c": 299_792_458, // speed of light in m/s "g": 6.67430e-11, // gravitational constant in N m2/kg2 "h": 6.62607015e-34, // planck constant in J s "ec": 1.602176634e-19, // elementary charge in C "e0": 8.8541878128e-12, // vacuum permittivity in C2/(Nm2) "mu0": 1.25663706212e-6, // vacuum permeability in T m/A "k": 1.380649e-23, // boltzmann constant in J/K "mu": 1.66053906660e-27, // atomic mass constant in kg "me": 9.1093837015e-31, // electron mass in kg "mp": 1.67262192369e-27, // proton mass in kg "mn": 1.67492749804e-27, // neutron mass in kg // float64s can only vaguely approx. avogadro's mole (6.02214076e23) } // deterministic math functions lookup-table generated using the command // // go doc math | awk '/func/ { gsub(/func |\(.*/, ""); printf("\"%s\": math.%s,\n", tolower($0), $0) }' // // then hand-edited to remove funcs, or to use adapter funcs when needed: removed // funcs either had multiple returns (like SinCos) or dealt with float32 values var determFuncs = map[string]any{ "abs": math.Abs, "acos": math.Acos, "acosh": math.Acosh, "asin": math.Asin, "asinh": math.Asinh, "atan": math.Atan, "atan2": math.Atan2, "atanh": math.Atanh, "cbrt": math.Cbrt, "ceil": math.Ceil, "copysign": math.Copysign, "cos": math.Cos, "cosh": math.Cosh, "dim": math.Dim, "erf": math.Erf, "erfc": math.Erfc, "erfcinv": math.Erfcinv, "erfinv": math.Erfinv, "exp": math.Exp, "exp2": math.Exp2, "expm1": math.Expm1, "fma": math.FMA, "floor": math.Floor, "gamma": math.Gamma, "inf": inf, "isinf": isInf, "isnan": isNaN, "j0": math.J0, "j1": math.J1, "jn": jn, "ldexp": ldexp, "lgamma": lgamma, "log10": math.Log10, "log1p": math.Log1p, "log2": math.Log2, "logb": math.Logb, "mod": math.Mod, "nan": math.NaN, "nextafter": math.Nextafter, "pow": math.Pow, "pow10": pow10, "remainder": math.Remainder, "round": math.Round, "roundtoeven": math.RoundToEven, "signbit": signbit, "sin": math.Sin, "sinh": math.Sinh, "sqrt": math.Sqrt, "tan": math.Tan, "tanh": math.Tanh, "trunc": math.Trunc, "y0": math.Y0, "y1": math.Y1, "yn": yn, // a few aliases for the standard math funcs: some of the single-letter // aliases are named after the ones in `bc`, the basic calculator tool "a": math.Abs, "c": math.Cos, "ceiling": math.Ceil, "cosine": math.Cos, "e": math.Exp, "isinf0": isAnyInf, "isinfinite": isAnyInf, "l": math.Log, "ln": math.Log, "lg": math.Log2, "modulus": math.Mod, "power": math.Pow, "rem": math.Remainder, "s": math.Sin, "sine": math.Sin, "t": math.Tan, "tangent": math.Tan, "truncate": math.Trunc, "truncated": math.Trunc, // not from standard math: these custom funcs were added manually "beta": beta, "bool": num2bool, "clamp": clamp, "cond": cond, // vector-arg if-else chain "cube": cube, "cubed": cube, "degrees": degrees, "deinf": deInf, "denan": deNaN, "factorial": factorial, "fract": fract, "horner": polyval, "hypot": hypot, "if": ifElse, "ifelse": ifElse, "inv": reciprocal, "isanyinf": isAnyInf, "isbad": isBad, "isfin": isFinite, "isfinite": isFinite, "isgood": isGood, "isinteger": isInteger, "lbeta": lnBeta, "len": length, "length": length, "lnbeta": lnBeta, "log": math.Log, "logistic": logistic, "mag": length, "max": max, "min": min, "neg": negate, "negate": negate, "not": notBool, "polyval": polyval, "radians": radians, "range": rangef, // vector-arg max - min "reciprocal": reciprocal, "rev": revalue, "revalue": revalue, "scale": scale, "sgn": sign, "sign": sign, "sinc": sinc, "sq": sqr, "sqmin": solveQuadMin, "sqmax": solveQuadMax, "square": sqr, "squared": sqr, "unwrap": unwrap, "wrap": wrap, // a few aliases for the custom funcs "deg": degrees, "isint": isInteger, "rad": radians, // a few quicker 2-value versions of vararg funcs: the optimizer depends // on these to rewrite 2-input uses of their vararg counterparts "hypot2": math.Hypot, "max2": math.Max, "min2": math.Min, // a few entries to enable custom syntax: the parser depends on these to // rewrite binary expressions into func calls "??": deNaN, "?:": ifElse, "===": same, "!==": notSame, } // DefineDetFuncs adds more deterministic funcs to the default set. Such funcs // are considered optimizable, since calling them with the same constant inputs // is supposed to return the same constant outputs, as the name `deterministic` // suggests. // // Only call this before compiling any scripts, and ensure all funcs given are // supported and are deterministic. Random-output funcs certainly won't fit // the bill here. func DefineDetFuncs(funcs map[string]any) { for k, v := range funcs { determFuncs[k] = v } } func sqr(x float64) float64 { return x * x } func cube(x float64) float64 { return x * x * x } func num2bool(x float64) float64 { if x == 0 { return 0 } return 1 } func logistic(x float64) float64 { return 1 / (1 + math.Exp(-x)) } func sign(x float64) float64 { if math.IsNaN(x) { return x } if x > 0 { return +1 } if x < 0 { return -1 } return 0 } func sinc(x float64) float64 { if x == 0 { return 1 } return math.Sin(x) / x } func isInteger(x float64) float64 { _, frac := math.Modf(x) if frac == 0 { return 1 } return 0 } func inf(sign float64) float64 { return math.Inf(int(sign)) } func isInf(x float64, sign float64) float64 { if math.IsInf(x, int(sign)) { return 1 } return 0 } func isAnyInf(x float64) float64 { if math.IsInf(x, 0) { return 1 } return 0 } func isFinite(x float64) float64 { if math.IsInf(x, 0) { return 0 } return 1 } func isNaN(x float64) float64 { if math.IsNaN(x) { return 1 } return 0 } func isGood(x float64) float64 { if math.IsNaN(x) || math.IsInf(x, 0) { return 0 } return 1 } func isBad(x float64) float64 { if math.IsNaN(x) || math.IsInf(x, 0) { return 1 } return 0 } func same(x, y float64) float64 { if math.IsNaN(x) && math.IsNaN(y) { return 1 } return debool(x == y) } func notSame(x, y float64) float64 { if math.IsNaN(x) && math.IsNaN(y) { return 0 } return debool(x != y) } func deNaN(x float64, instead float64) float64 { if !math.IsNaN(x) { return x } return instead } func deInf(x float64, instead float64) float64 { if !math.IsInf(x, 0) { return x } return instead } func revalue(x float64, instead float64) float64 { if !math.IsNaN(x) && !math.IsInf(x, 0) { return x } return instead } func pow10(n float64) float64 { return math.Pow10(int(n)) } func jn(n, x float64) float64 { return math.Jn(int(n), x) } func ldexp(frac, exp float64) float64 { return math.Ldexp(frac, int(exp)) } func lgamma(x float64) float64 { y, s := math.Lgamma(x) if s < 0 { return math.NaN() } return y } func signbit(x float64) float64 { if math.Signbit(x) { return 1 } return 0 } func yn(n, x float64) float64 { return math.Yn(int(n), x) } func negate(x float64) float64 { return -x } func reciprocal(x float64) float64 { return 1 / x } func rangef(v ...float64) float64 { min := math.Inf(+1) max := math.Inf(-1) for _, f := range v { min = math.Min(min, f) max = math.Max(max, f) } return max - min } func cond(v ...float64) float64 { for { switch len(v) { case 0: // either no values are left, or no values were given at all return math.NaN() case 1: // either all previous conditions failed, and this last value is // automatically chosen, or only 1 value was given to begin with return v[0] default: // check condition: if true (non-0), return the value after it if v[0] != 0 { return v[1] } // condition was false, so skip the leading pair of values v = v[2:] } } } func notBool(x float64) float64 { if x == 0 { return 1 } return 0 } func ifElse(cond float64, yes, no float64) float64 { if cond != 0 { return yes } return no } func lnGamma(x float64) float64 { y, s := math.Lgamma(x) if s < 0 { return math.NaN() } return y } func lnBeta(x float64, y float64) float64 { return lnGamma(x) + lnGamma(y) - lnGamma(x+y) } func beta(x float64, y float64) float64 { return math.Exp(lnBeta(x, y)) } func factorial(n float64) float64 { return math.Round(math.Gamma(n + 1)) } func degrees(rad float64) float64 { return rad2degMult * rad } func radians(deg float64) float64 { return deg2radMult * deg } func fract(x float64) float64 { return x - math.Floor(x) } func min(v ...float64) float64 { min := +math.Inf(+1) for _, f := range v { min = math.Min(min, f) } return min } func max(v ...float64) float64 { max := +math.Inf(-1) for _, f := range v { max = math.Max(max, f) } return max } func hypot(v ...float64) float64 { sumsq := 0.0 for _, f := range v { sumsq += f * f } return math.Sqrt(sumsq) } func length(v ...float64) float64 { ss := 0.0 for _, f := range v { ss += f * f } return math.Sqrt(ss) } // solveQuadMin finds the lowest solution of a 2nd-degree polynomial, using // a formula which is more accurate than the textbook one func solveQuadMin(a, b, c float64) float64 { disc := math.Sqrt(b*b - 4*a*c) r1 := 2 * c / (-b - disc) return r1 } // solveQuadMax finds the highest solution of a 2nd-degree polynomial, using // a formula which is more accurate than the textbook one func solveQuadMax(a, b, c float64) float64 { disc := math.Sqrt(b*b - 4*a*c) r2 := 2 * c / (-b + disc) return r2 } func wrap(x float64, min, max float64) float64 { return (x - min) / (max - min) } func unwrap(x float64, min, max float64) float64 { return (max-min)*x + min } func clamp(x float64, min, max float64) float64 { return math.Min(math.Max(x, min), max) } func scale(x float64, xmin, xmax, ymin, ymax float64) float64 { k := (x - xmin) / (xmax - xmin) return (ymax-ymin)*k + ymin } // polyval runs horner's algorithm on a value, with the polynomial coefficients // given after it, higher-order first func polyval(x float64, v ...float64) float64 { if len(v) == 0 { return 0 } x0 := x x = 1.0 y := 0.0 for i := len(v) - 1; i >= 0; i-- { y += v[i] * x x *= x0 } return y } // unary2func matches unary operators to funcs the optimizer can use to eval // constant-input unary expressions into their results var unary2func = map[string]func(x float64) float64{ // avoid unary +, since it's a no-op and thus there's no instruction for // it, which in turn makes unit-tests for these optimization tables fail // unnecessarily // "+": func(x float64) float64 { return +x }, "-": func(x float64) float64 { return -x }, "!": func(x float64) float64 { return debool(x == 0) }, "&": func(x float64) float64 { return math.Abs(x) }, "*": func(x float64) float64 { return x * x }, "^": func(x float64) float64 { return x * x }, "/": func(x float64) float64 { return 1 / x }, } // binary2func matches binary operators to funcs the optimizer can use to eval // constant-input binary expressions into their results var binary2func = map[string]func(x, y float64) float64{ "+": func(x, y float64) float64 { return x + y }, "-": func(x, y float64) float64 { return x - y }, "*": func(x, y float64) float64 { return x * y }, "/": func(x, y float64) float64 { return x / y }, "%": func(x, y float64) float64 { return math.Mod(x, y) }, "&": func(x, y float64) float64 { return debool(x != 0 && y != 0) }, "&&": func(x, y float64) float64 { return debool(x != 0 && y != 0) }, "|": func(x, y float64) float64 { return debool(x != 0 || y != 0) }, "||": func(x, y float64) float64 { return debool(x != 0 || y != 0) }, "==": func(x, y float64) float64 { return debool(x == y) }, "!=": func(x, y float64) float64 { return debool(x != y) }, "<>": func(x, y float64) float64 { return debool(x != y) }, "<": func(x, y float64) float64 { return debool(x < y) }, "<=": func(x, y float64) float64 { return debool(x <= y) }, ">": func(x, y float64) float64 { return debool(x > y) }, ">=": func(x, y float64) float64 { return debool(x >= y) }, "**": func(x, y float64) float64 { return math.Pow(x, y) }, "^": func(x, y float64) float64 { return math.Pow(x, y) }, } // func2unary turns built-in func names into built-in unary operators var func2unary = map[string]string{ "a": "&", "abs": "&", "inv": "/", "neg": "-", "negate": "-", "reciprocal": "/", "sq": "*", "square": "*", "squared": "*", } // func2binary turns built-in func names into built-in binary operators var func2binary = map[string]string{ "mod": "%", "modulus": "%", "pow": "**", "power": "**", } // vararg2func2 matches variable-argument funcs to their 2-argument versions var vararg2func2 = map[string]string{ "hypot": "hypot2", "max": "max2", "min": "min2", } // optimize tries to simplify the expression given as much as possible, by // simplifying constants whenever possible, and exploiting known built-in // funcs which are known to behave deterministically func (c *Compiler) optimize(expr any) any { switch expr := expr.(type) { case []any: return c.optimizeCombo(expr) case unaryExpr: return c.optimizeUnaryExpr(expr) case binaryExpr: return c.optimizeBinaryExpr(expr) case callExpr: return c.optimizeCallExpr(expr) case assignExpr: expr.expr = c.optimize(expr.expr) return expr default: f, ok := c.tryConstant(expr) if ok { return f } return expr } } // optimizeCombo handles a sequence of expressions for the optimizer func (c *Compiler) optimizeCombo(exprs []any) any { if len(exprs) == 1 { return c.optimize(exprs[0]) } // count how many expressions are considered useful: these are // assignments as well as the last expression, since that's what // determines the script's final result useful := 0 for i, v := range exprs { _, ok := v.(assignExpr) if ok || i == len(exprs)-1 { useful++ } } // ignore all expressions which are a waste of time, and optimize // all other expressions res := make([]any, 0, useful) for i, v := range exprs { _, ok := v.(assignExpr) if ok || i == len(exprs)-1 { res = append(res, c.optimize(v)) } } return res } // optimizeUnaryExpr handles unary expressions for the optimizer func (c *Compiler) optimizeUnaryExpr(expr unaryExpr) any { // recursively optimize input expr.x = c.optimize(expr.x) // optimize unary ops on a constant into concrete values if x, ok := expr.x.(float64); ok { if fn, ok := unary2func[expr.op]; ok { return fn(x) } } switch expr.op { case "+": // unary plus is an identity operation return expr.x default: return expr } } // optimizeBinaryExpr handles binary expressions for the optimizer func (c *Compiler) optimizeBinaryExpr(expr binaryExpr) any { // recursively optimize inputs expr.x = c.optimize(expr.x) expr.y = c.optimize(expr.y) // optimize binary ops on 2 constants into concrete values if x, ok := expr.x.(float64); ok { if y, ok := expr.y.(float64); ok { if fn, ok := binary2func[expr.op]; ok { return fn(x, y) } } } switch expr.op { case "+": if expr.x == 0.0 { // 0+y -> y return expr.y } if expr.y == 0.0 { // x+0 -> x return expr.x } case "-": if expr.x == 0.0 { // 0-y -> -y return c.optimizeUnaryExpr(unaryExpr{op: "-", x: expr.y}) } if expr.y == 0.0 { // x-0 -> x return expr.x } case "*": if expr.x == 0.0 || expr.y == 0.0 { // 0*y -> 0 // x*0 -> 0 return 0.0 } if expr.x == 1.0 { // 1*y -> y return expr.y } if expr.y == 1.0 { // x*1 -> x return expr.x } if expr.x == -1.0 { // -1*y -> -y return c.optimizeUnaryExpr(unaryExpr{op: "-", x: expr.y}) } if expr.y == -1.0 { // x*-1 -> -x return c.optimizeUnaryExpr(unaryExpr{op: "-", x: expr.x}) } case "/": if expr.x == 1.0 { // 1/y -> reciprocal of y return c.optimizeUnaryExpr(unaryExpr{op: "/", x: expr.y}) } if expr.y == 1.0 { // x/1 -> x return expr.x } if expr.y == -1.0 { // x/-1 -> -x return c.optimizeUnaryExpr(unaryExpr{op: "-", x: expr.x}) } case "**": switch expr.y { case -1.0: // x**-1 -> 1/x, reciprocal of x return c.optimizeUnaryExpr(unaryExpr{op: "/", x: expr.x}) case 0.0: // x**0 -> 1 return 1.0 case 1.0: // x**1 -> x return expr.x case 2.0: // x**2 -> *x return c.optimizeUnaryExpr(unaryExpr{op: "*", x: expr.x}) case 3.0: // x**3 -> *x*x sq := unaryExpr{op: "*", x: expr.x} return c.optimizeBinaryExpr(binaryExpr{op: "*", x: sq, y: expr.x}) } case "&", "&&": if expr.x == 0.0 || expr.y == 0.0 { // 0 && y -> 0 // x && 0 -> 0 return 0.0 } } // no simplifiable patterns were detected return expr } // optimizeCallExpr optimizes special cases of built-in func calls func (c *Compiler) optimizeCallExpr(call callExpr) any { // recursively optimize all inputs, and keep track if they're all // constants in the end numlit := 0 for i, v := range call.args { v = c.optimize(v) call.args[i] = v if _, ok := v.(float64); ok { numlit++ } } // if func is overridden, there's no guarantee the new func works the same if _, ok := determFuncs[call.name]; c.ftype[call.name] != nil || !ok { return call } // from this point on, func is guaranteed to be built-in and deterministic // handle all-const inputs, by calling func and using its result if numlit == len(call.args) { in := make([]float64, 0, len(call.args)) for _, v := range call.args { f, _ := v.(float64) in = append(in, f) } if f, ok := tryCall(determFuncs[call.name], in); ok { return f } } switch len(call.args) { case 1: if op, ok := func2unary[call.name]; ok { expr := unaryExpr{op: op, x: call.args[0]} return c.optimizeUnaryExpr(expr) } return call case 2: if op, ok := func2binary[call.name]; ok { expr := binaryExpr{op: op, x: call.args[0], y: call.args[1]} return c.optimizeBinaryExpr(expr) } if name, ok := vararg2func2[call.name]; ok { call.name = name return call } return call default: return call } } // tryConstant tries to optimize the expression given into a constant func (c *Compiler) tryConstant(expr any) (value float64, ok bool) { switch expr := expr.(type) { case float64: return expr, true case string: if _, ok := c.vaddr[expr]; !ok { // name isn't explicitly defined if f, ok := mathConst[expr]; ok { // and is a known math constant return f, true } } return 0, false default: return 0, false } } // tryCall tries to simplify the function expression given into a constant func tryCall(fn any, in []float64) (value float64, ok bool) { switch fn := fn.(type) { case func(float64) float64: if len(in) == 1 { return fn(in[0]), true } return 0, false case func(float64, float64) float64: if len(in) == 2 { return fn(in[0], in[1]), true } return 0, false case func(float64, float64, float64) float64: if len(in) == 3 { return fn(in[0], in[1], in[2]), true } return 0, false case func(float64, float64, float64, float64) float64: if len(in) == 4 { return fn(in[0], in[1], in[2], in[3]), true } return 0, false case func(float64, float64, float64, float64, float64) float64: if len(in) == 5 { return fn(in[0], in[1], in[2], in[3], in[4]), true } return 0, false case func(...float64) float64: return fn(in...), true case func(float64, ...float64) float64: if len(in) >= 1 { return fn(in[0], in[1:]...), true } return 0, false default: // type isn't a supported func return 0, false } } // parse turns source code into an expression type interpreters can use. func parse(src string) (any, error) { tok := newTokenizer(src) par, err := newParser(&tok) if err != nil { return nil, err } v, err := par.parse() err = par.improveError(err, src) return v, err } // pickLine slices a string, picking one of its lines via a 1-based index: // func improveError uses it to isolate the line of code an error came from func pickLine(src string, linenum int) string { // skip the lines before the target one for i := 0; i < linenum && len(src) > 0; i++ { j := strings.IndexByte(src, '\n') if j < 0 { break } src = src[j+1:] } // limit leftover to a single line i := strings.IndexByte(src, '\n') if i >= 0 { return src[:i] } return src } // unaryExpr is a unary expression type unaryExpr struct { op string x any } // binaryExpr is a binary expression type binaryExpr struct { op string x any y any } // callExpr is a function call and its arguments type callExpr struct { name string args []any } // assignExpr is a value/variable assignment type assignExpr struct { name string expr any } // parser is a parser for JavaScript-like syntax, limited to operations on // 64-bit floating-point numbers type parser struct { tokens []token line int pos int toklen int } // newParser is the constructor for type parser func newParser(t *tokenizer) (parser, error) { var p parser // get all tokens from the source code for { v, err := t.next() if v.kind != unknownToken { p.tokens = append(p.tokens, v) } if err == errEOS { // done scanning/tokenizing return p, nil } if err != nil { // handle actual errors return p, err } } } // improveError adds more info about where exactly in the source code an error // came from, thus making error messages much more useful func (p *parser) improveError(err error, src string) error { if err == nil { return nil } line := pickLine(src, p.line) if len(line) == 0 || p.pos < 1 { const fs = "(line %d: pos %d): %w" return fmt.Errorf(fs, p.line, p.pos, err) } ptr := strings.Repeat(" ", p.pos) + "^" const fs = "(line %d: pos %d): %w\n%s\n%s" return fmt.Errorf(fs, p.line, p.pos, err, line, ptr) } // parse tries to turn tokens into a compilable abstract syntax tree, and is // the parser's entry point func (p *parser) parse() (any, error) { // source codes with no tokens are always errors if len(p.tokens) == 0 { const msg = "source code is empty, or has no useful expressions" return nil, errors.New(msg) } // handle optional excel-like leading equal sign p.acceptSyntax(`=`) // ignore trailing semicolons for len(p.tokens) > 0 { t := p.tokens[len(p.tokens)-1] if t.kind != syntaxToken || t.value != `;` { break } p.tokens = p.tokens[:len(p.tokens)-1] } // handle single expressions as well as multiple semicolon-separated // expressions: the latter allow value assignments/updates in scripts var res []any for keepGoing := true; keepGoing && len(p.tokens) > 0; { v, err := p.parseExpression() if err != nil && err != errEOS { return v, err } res = append(res, v) // handle optional separator/continuation semicolons _, keepGoing = p.acceptSyntax(`;`) } // unexpected unparsed trailing tokens are always errors; any trailing // semicolons in the original script are already trimmed if len(p.tokens) > 0 { const fs = "unexpected %s" return res, fmt.Errorf(fs, p.tokens[0].value) } // make scripts ending in an assignment also load that value, so they're // useful, as assignments result in no useful value by themselves assign, ok := res[len(res)-1].(assignExpr) if ok { res = append(res, assign.name) } // turn 1-item combo expressions into their only expression: some // unit tests may rely on that for convenience if len(res) == 1 { return res[0], nil } return res, nil } // acceptSyntax advances the parser on the first syntactic string matched: // notice any number of alternatives/options are allowed, as the syntax // allows/requires at various points func (p *parser) acceptSyntax(syntax ...string) (match string, ok bool) { if len(p.tokens) == 0 { return "", false } t := p.tokens[0] if t.kind != syntaxToken { return "", false } for _, s := range syntax { if t.value == s { p.advance() return s, true } } return "", false } // advance skips the current leading token, if there are still any left func (p *parser) advance() { if len(p.tokens) == 0 { return } t := p.tokens[0] p.tokens = p.tokens[1:] p.line = t.line p.pos = t.pos p.toklen = len(t.value) } // acceptNumeric advances the parser on a numeric value, but only if it's // the leading token: conversely, any other type of token doesn't advance // the parser; when matches happen the resulting strings need parsing via // func parseNumber func (p *parser) acceptNumeric() (numliteral string, ok bool) { if len(p.tokens) == 0 { return "", false } t := p.tokens[0] if t.kind == numberToken { p.advance() return t.value, true } return "", false } // demandSyntax imposes a specific syntactic element to follow, or else it's // an error func (p *parser) demandSyntax(syntax string) error { if len(p.tokens) == 0 { return fmt.Errorf("expected %s instead of the end of source", syntax) } first := p.tokens[0] if first.kind == syntaxToken && first.value == syntax { p.advance() return nil } return fmt.Errorf("expected %s instead of %s", syntax, first.value) } func (p *parser) parseExpression() (any, error) { x, err := p.parseComparison() if err != nil { return x, err } // handle assignment statements if _, ok := p.acceptSyntax(`=`, `:=`); ok { varname, ok := x.(string) if !ok { const fs = "expected a variable name, instead of a %T" return nil, fmt.Errorf(fs, x) } x, err := p.parseExpression() expr := assignExpr{name: varname, expr: x} return expr, err } // handle and/or logical chains for { if op, ok := p.acceptSyntax(`&&`, `||`, `&`, `|`); ok { y, err := p.parseExpression() if err != nil { return y, err } x = binaryExpr{op: op, x: x, y: y} continue } break } // handle maybe-properties if _, ok := p.acceptSyntax(`??`); ok { y, err := p.parseExpression() expr := callExpr{name: "??", args: []any{x, y}} return expr, err } // handle choice/ternary operator if _, ok := p.acceptSyntax(`?`); ok { y, err := p.parseExpression() if err != nil { expr := callExpr{name: "?:", args: []any{x, nil, nil}} return expr, err } if _, ok := p.acceptSyntax(`:`); ok { z, err := p.parseExpression() expr := callExpr{name: "?:", args: []any{x, y, z}} return expr, err } if len(p.tokens) == 0 { expr := callExpr{name: "?:", args: []any{x, y, nil}} return expr, errors.New("expected `:`") } s := p.tokens[0].value expr := callExpr{name: "?:", args: []any{x, y, nil}} err = fmt.Errorf("expected `:`, but got %q instead", s) return expr, err } // expression was just a comparison, or simpler return x, nil } func (p *parser) parseComparison() (any, error) { x, err := p.parseTerm() if err != nil { return x, err } op, ok := p.acceptSyntax(`==`, `!=`, `<`, `>`, `<=`, `>=`, `<>`, `===`, `!==`) if ok { y, err := p.parseTerm() return binaryExpr{op: op, x: x, y: y}, err } return x, err } // parseBinary handles binary operations, by recursing depth-first on the left // side of binary expressions; going tail-recursive on these would reverse the // order of arguments instead, which is obviously wrong func (p *parser) parseBinary(parse func() (any, error), syntax ...string) (any, error) { x, err := parse() if err != nil { return x, err } for { op, ok := p.acceptSyntax(syntax...) if !ok { return x, nil } y, err := parse() x = binaryExpr{op: op, x: x, y: y} if err != nil { return x, err } } } func (p *parser) parseTerm() (any, error) { return p.parseBinary(func() (any, error) { return p.parseProduct() }, `+`, `-`, `^`) } func (p *parser) parseProduct() (any, error) { return p.parseBinary(func() (any, error) { return p.parsePower() }, `*`, `/`, `%`) } func (p *parser) parsePower() (any, error) { return p.parseBinary(func() (any, error) { return p.parseValue() }, `**`, `^`) } func (p *parser) parseValue() (any, error) { // handle unary operators which can also be considered part of numeric // literals, and thus should be simplified away if op, ok := p.acceptSyntax(`+`, `-`); ok { if s, ok := p.acceptNumeric(); ok { x, err := strconv.ParseFloat(s, 64) if err != nil { return nil, err } if simpler, ok := simplifyNumber(op, x); ok { return simpler, nil } return unaryExpr{op: op, x: x}, nil } x, err := p.parsePower() return unaryExpr{op: op, x: x}, err } // handle all other unary operators if op, ok := p.acceptSyntax(`!`, `&`, `*`, `^`); ok { x, err := p.parsePower() return unaryExpr{op: op, x: x}, err } // handle subexpression in parentheses if _, ok := p.acceptSyntax(`(`); ok { x, err := p.parseExpression() if err != nil { return x, err } if err := p.demandSyntax(`)`); err != nil { return x, err } return p.parseAccessors(x) } // handle subexpression in square brackets: it's just an alternative to // using parentheses for subexpressions if _, ok := p.acceptSyntax(`[`); ok { x, err := p.parseExpression() if err != nil { return x, err } if err := p.demandSyntax(`]`); err != nil { return x, err } return p.parseAccessors(x) } // handle all other cases x, err := p.parseSimpleValue() if err != nil { return x, err } // handle arbitrarily-long chains of accessors return p.parseAccessors(x) } // parseSimpleValue handles a numeric literal or a variable/func name, also // known as identifier func (p *parser) parseSimpleValue() (any, error) { if len(p.tokens) == 0 { return nil, errEOS } t := p.tokens[0] switch t.kind { case identifierToken: p.advance() // handle func calls, such as f(...) if _, ok := p.acceptSyntax(`(`); ok { args, err := p.parseList(`)`) expr := callExpr{name: t.value, args: args} return expr, err } // handle func calls, such as f[...] if _, ok := p.acceptSyntax(`[`); ok { args, err := p.parseList(`]`) expr := callExpr{name: t.value, args: args} return expr, err } return t.value, nil case numberToken: p.advance() return strconv.ParseFloat(t.value, 64) default: const fs = "unexpected %s (token type %d)" return nil, fmt.Errorf(fs, t.value, t.kind) } } // parseAccessors handles an arbitrarily-long chain of accessors func (p *parser) parseAccessors(x any) (any, error) { for { s, ok := p.acceptSyntax(`.`, `@`) if !ok { // dot-chain is over return x, nil } // handle property/method accessors v, err := p.parseDot(s, x) if err != nil { return v, err } x = v } } // parseDot handles what follows a syntactic dot, as opposed to a dot which // may be part of a numeric literal func (p *parser) parseDot(after string, x any) (any, error) { if len(p.tokens) == 0 { const fs = "unexpected end of source after a %q" return x, fmt.Errorf(fs, after) } t := p.tokens[0] p.advance() if t.kind != identifierToken { const fs = "expected a valid property name, but got %s instead" return x, fmt.Errorf(fs, t.value) } if _, ok := p.acceptSyntax(`(`); ok { items, err := p.parseList(`)`) args := append([]any{x}, items...) return callExpr{name: t.value, args: args}, err } if _, ok := p.acceptSyntax(`[`); ok { items, err := p.parseList(`]`) args := append([]any{x}, items...) return callExpr{name: t.value, args: args}, err } return callExpr{name: t.value, args: []any{x}}, nil } // parseList handles the argument-list following a `(` or a `[` func (p *parser) parseList(end string) ([]any, error) { var arr []any for len(p.tokens) > 0 { if _, ok := p.acceptSyntax(`,`); ok { // ensure extra/trailing commas are allowed/ignored continue } if _, ok := p.acceptSyntax(end); ok { return arr, nil } v, err := p.parseExpression() if err != nil { return arr, err } arr = append(arr, v) } // return an appropriate error for the unexpected end of the source return arr, p.demandSyntax(`)`) } // simplifyNumber tries to simplify a few trivial unary operations on // numeric constants func simplifyNumber(op string, x any) (v any, ok bool) { if x, ok := x.(float64); ok { switch op { case `+`: return x, true case `-`: return -x, true default: return x, false } } return x, false } // So far, my apps relying on this package are // // - fh, a Function Heatmapper which color-codes f(x, y), seen from above // - star, a STAtistical Resampler to calculate stats from tabular data // - waveout, an app to emit sample-by-sample wave-audio by formula // Quick Notes on Performance // // Note: while these microbenchmarks are far from proper benchmarks, Club Pulse // can still give a general idea of what relative performance to expect. // // While funscripts.Interpreter can do anything a Program can do, and much more, // a fmscripts.Program is way faster for float64-only tasks: typical speed-ups // are between 20-to-50 times. Various simple benchmarks suggest running speed // is between 1/2 and 1/4 of native funcs. // // Reimplementations of the Club Pulse benchmark, when run 50 million times, // suggest this package is starting to measure the relative speed of various // trigonometric func across JIT-ted script runners, running somewhat slower // than NodeJS/V8, faster than PyPy, and much faster than Python. // // Other scripted tests, which aren't as trigonometry/exponential-heavy, hint // at even higher speed-ups compared to Python and PyPy, as well as being // almost as fast as NodeJS/V8. // // Being so close to Node's V8 (0.6x - 0.8x its speed) is a surprisingly good // result for the relatively little effort this package took. const ( maxArgsLen = 1<<16 - 1 // max length for the []float64 input of callv maxOpIndex = 1<<32 - 1 // max index for values ) // types to keep compiler code the same, while changing sizes of numOp fields type ( opcode uint16 opargs uint16 // opcode directly affects max correct value of maxArgsLen opindex uint32 // opcode directly affects max correct value of maxOpIndex ) // numOp is a numeric operation in a program type numOp struct { What opcode // what to do NumArgs opargs // vector-length only used for callv and call1v Index opindex // index to load a value or call a function } // Since the go 1.19 compiler started compiling dense switch statements using // jump tables, adding more basic ops doesn't slow things down anymore when // reaching the next power-of-2 thresholds in the number of cases. const ( // load float64 value to top of the stack load opcode = iota // pop top of stack and store it into value store opcode = iota // unary operations neg opcode = iota not opcode = iota abs opcode = iota square opcode = iota // 1/x, the reciprocal/inverse value rec opcode = iota // x%1, but faster than math.Mod(x, 1) mod1 opcode = iota // arithmetic and logic operations: 2 float64 inputs, 1 float64 output add opcode = iota sub opcode = iota mul opcode = iota div opcode = iota mod opcode = iota pow opcode = iota and opcode = iota or opcode = iota // binary comparisons: 2 float64 inputs, 1 booleanized float64 output equal opcode = iota notequal opcode = iota less opcode = iota lessoreq opcode = iota more opcode = iota moreoreq opcode = iota // function callers with 1..5 float64 inputs and a float64 result call0 opcode = iota call1 opcode = iota call2 opcode = iota call3 opcode = iota call4 opcode = iota call5 opcode = iota // var-arg input function callers callv opcode = iota call1v opcode = iota ) // Program runs a sequence of float64 operations, with no explicit control-flow: // implicit control-flow is available in the float64-only functions you make // available to such programs. Such custom funcs must all return a float64, and // either take from 0 to 5 float64 arguments, or a single float64 array. // // As the name suggests, don't create such objects directly, but instead use // Compiler.Compile to create them. Only a Compiler lets you register variables // and functions, which then become part of your numeric Program. // // A Program lets you change values before each run, using pointers from method // Program.Get: such pointers are guaranteed never to change before or across // runs. Just ensure you get a variable defined in the Compiler used to make the // Program, or the pointer will be to a dummy place which has no effect on final // results. // // Expressions using literals are automatically optimized into their results: // this also applies to func calls with standard math functions and constants. // // The only way to limit such optimizations is to redefine common math funcs and // const values explicitly when compiling: after doing so, all those value-names // will stand for externally-updatable values which can change from one run to // the next. Similarly, there's no guarantee the (re)defined functions will be // deterministic, like the defaults they replaced. // // # Example // // var c fmscripts.Compiler // // defs := map[string]any{ // "x": 0, // define `x`, and initialize it to 0 // "k": 4.25, // define `k`, and initialize it to 4.25 // "b": true, // define `b`, and initialize it to 1.0 // "n": -23, // define `n`, and initialize it to -23.0 // "pi": 3, // define `pi`, overriding the automatic constant named `pi` // // "f": numericKernel // type is func ([]float64) float64 // "g": otherFunc // type is func (float64) float64 // } // // prog, err := c.Compile("log10(k) + f(sqrt(k) * exp(-x), 45, -0.23)", defs) // // ... // // x, _ := prog.Get("x") // Get returns (*float64, bool) // y, _ := prog.Get("y") // a useless pointer, since program doesn't use `y` // // ... // // for i := 0; i < n; i++ { // *x = float64(i)*dx + minx // you update inputs in place using pointers // f := prog.Run() // method Run gives you a float64 back // // ... // } type Program struct { sp int // stack pointer, even though it's an index stack []float64 // pre-allocated by compiler to max length needed by program values []float64 // holds all values, whether literals, or variables ops []numOp // all sequential operations for each run funcs []any // all funcs used by the program names map[string]int // variable-name to index lookup dummy float64 // all pointers for undefined variables point here // data []float64 } // Get lets you change parameters/variables before each time you run a program, // since it doesn't return a value, but a pointer to it, so you can update it // in place. // // If the name given isn't available, the result is a pointer to a dummy place: // this ensures non-nil pointers, which are always safe to use, even though // updates to the dummy destination have no effect on program results. func (p *Program) Get(name string) (ptr *float64, useful bool) { if i, ok := p.names[name]; ok { return &p.values[i], true } return &p.dummy, false } // Clone creates an exact copy of a Program with all values in their current // state: this is useful when dispatching embarassingly-parallel tasks to // multiple programs. func (p Program) Clone() Program { // can't share the stack nor the values stack := make([]float64, len(p.stack)) values := make([]float64, len(p.values)) copy(stack, p.stack) copy(values, p.values) p.stack = stack p.values = values // can share everything else as is return p } // Memo: the command to show all bound checks is // go test -gcflags="-d=ssa/check_bce/debug=1" fmscripts/programs.go // Discussion about go compiler optimizing switch statements into jump tables // https://go-review.googlesource.com/c/go/+/357330/ // Run executes the program once. Before each run, update input values using // pointers from method Get. func (p *Program) Run() float64 { // Check for empty programs: these happen either when a compilation error // was ignored, or when a program was explicitly declared as a variable. if len(p.ops) == 0 { return math.NaN() } p.sp = 0 p.runAllOps() return p.stack[p.sp] } type func4 = func(float64, float64, float64, float64) float64 type func5 = func(float64, float64, float64, float64, float64) float64 // runAllOps runs all operations in a loop: when done, the program's result is // ready as the only item left in the stack func (p *Program) runAllOps() { for _, op := range p.ops { // shortcut for the current stack pointer sp := p.sp // Preceding binary ops and func calls by _ = p.stack[i-n] prevents // an extra bound check for the lhs of assignments, but a quick // statistical summary of benchmarks doesn't show clear speedups, // let alone major ones. // // Separating different func types into different arrays to avoid // type-checking at runtime doesn't seem to be worth it either, // and makes the compiler more complicated. switch op.What { case load: // store above top of stack p.stack[sp+1] = p.values[op.Index] p.sp++ case store: // store from top of the stack p.values[op.Index] = p.stack[sp] p.sp-- // unary operations case neg: p.stack[sp] = -p.stack[sp] case not: p.stack[sp] = deboolNot(p.stack[sp]) case abs: p.stack[sp] = math.Abs(p.stack[sp]) case square: p.stack[sp] *= p.stack[sp] case rec: p.stack[sp] = 1 / p.stack[sp] // binary arithmetic ops case add: p.stack[sp-1] += p.stack[sp] p.sp-- case sub: p.stack[sp-1] -= p.stack[sp] p.sp-- case mul: p.stack[sp-1] *= p.stack[sp] p.sp-- case div: p.stack[sp-1] /= p.stack[sp] p.sp-- case mod: p.stack[sp-1] = math.Mod(p.stack[sp-1], p.stack[sp]) p.sp-- case pow: p.stack[sp-1] = math.Pow(p.stack[sp-1], p.stack[sp]) p.sp-- // binary boolean ops / binary comparisons case and: p.stack[sp-1] = deboolAnd(p.stack[sp-1], p.stack[sp]) p.sp-- case or: p.stack[sp-1] = deboolOr(p.stack[sp-1], p.stack[sp]) p.sp-- case equal: p.stack[sp-1] = debool(p.stack[sp-1] == p.stack[sp]) p.sp-- case notequal: p.stack[sp-1] = debool(p.stack[sp-1] != p.stack[sp]) p.sp-- case less: p.stack[sp-1] = debool(p.stack[sp-1] < p.stack[sp]) p.sp-- case lessoreq: p.stack[sp-1] = debool(p.stack[sp-1] <= p.stack[sp]) p.sp-- case more: p.stack[sp-1] = debool(p.stack[sp-1] > p.stack[sp]) p.sp-- case moreoreq: p.stack[sp-1] = debool(p.stack[sp-1] >= p.stack[sp]) p.sp-- // function calls case call0: f := p.funcs[op.Index].(func() float64) // store above top of stack p.stack[sp+1] = f() p.sp++ case call1: f := p.funcs[op.Index].(func(float64) float64) p.stack[sp] = f(p.stack[sp]) case call2: f := p.funcs[op.Index].(func(float64, float64) float64) p.stack[sp-1] = f(p.stack[sp-1], p.stack[sp]) p.sp-- case call3: f := p.funcs[op.Index].(func(float64, float64, float64) float64) p.stack[sp-2] = f(p.stack[sp-2], p.stack[sp-1], p.stack[sp]) p.sp -= 2 case call4: f := p.funcs[op.Index].(func4) st := p.stack p.stack[sp-3] = f(st[sp-3], st[sp-2], st[sp-1], st[sp]) p.sp -= 3 case call5: f := p.funcs[op.Index].(func5) st := p.stack p.stack[sp-4] = f(st[sp-4], st[sp-3], st[sp-2], st[sp-1], st[sp]) p.sp -= 4 case callv: i := sp - int(op.NumArgs) + 1 f := p.funcs[op.Index].(func(...float64) float64) p.stack[sp-i+1] = f(p.stack[i : sp+1]...) p.sp = sp - i + 1 case call1v: i := sp - int(op.NumArgs) + 1 f := p.funcs[op.Index].(func(float64, ...float64) float64) p.stack[sp-i+1] = f(p.stack[i], p.stack[i+1:sp+1]...) p.sp = sp - i + 1 } } } // debool is only used to turn boolean values used in comparison operations // into float64s, since those are the only type accepted on a program stack func debool(b bool) float64 { if b { return 1 } return 0 } // deboolNot runs the basic `not` operation func deboolNot(x float64) float64 { return debool(x == 0) } // deboolAnd runs the basic `and` operation func deboolAnd(x, y float64) float64 { return debool(x != 0 && y != 0) } // deboolOr runs the basic `or` operation func deboolOr(x, y float64) float64 { return debool(x != 0 || y != 0) } // These funcs are kept separate from the larger lookup table so they aren't // mistaken for deterministic funcs which can be optimized away. // // var randomFuncs = map[string]any{ // "rand": Random, // "rint": RandomInt, // "runif": RandomUnif, // "rexp": RandomExp, // "rnorm": RandomNorm, // "rgamma": RandomGamma, // "rbeta": RandomBeta, // } func Random(r *rand.Rand) float64 { return r.Float64() } func RandomInt(r *rand.Rand, min, max float64) float64 { fmin := math.Trunc(min) fmax := math.Trunc(max) if fmin == fmax { return fmin } diff := math.Abs(fmax - fmin) return float64(r.Intn(int(diff)+1)) + fmin } func RandomUnif(r *rand.Rand, min, max float64) float64 { return (max-min)*r.Float64() + min } func RandomExp(r *rand.Rand, scale float64) float64 { return scale * r.ExpFloat64() } func RandomNorm(r *rand.Rand, mu, sigma float64) float64 { return sigma*r.NormFloat64() + mu } // Gamma generates a gamma-distributed real value, using a scale parameter. // // The algorithm is from Marsaglia and Tsang, as described in // // A simple method for generating gamma variables // https://dl.acm.org/doi/10.1145/358407.358414 func RandomGamma(r *rand.Rand, scale float64) float64 { d := scale - 1.0/3.0 c := 1 / math.Sqrt(9/d) for { // generate candidate value var x, v float64 for { x = r.NormFloat64() v = 1 + c*x if v > 0 { break } } v = v * v * v // accept or reject candidate value x2 := x * x u := r.Float64() if u < 1-0.0331*x2*x2 { return d * v } if math.Log(u) < 0.5*x2+d*(1-v+math.Log(v)) { return d * v } } } // Beta generates a beta-distributed real value. func RandomBeta(r *rand.Rand, a, b float64) float64 { return RandomGamma(r, a) / (RandomGamma(r, a) + RandomGamma(r, b)) } // tokenType is the specific type of a token; tokens never represent // whitespace-like text between recognized tokens type tokenType int const ( // default zero-value type for tokens unknownToken tokenType = iota // a name identifierToken tokenType = iota // a literal numeric value numberToken tokenType = iota // any syntactic element made of 1 or more runes syntaxToken tokenType = iota ) // errEOS signals the end of source code, and is the only token-related error // which the parser should ignore var errEOS = errors.New(`no more source code`) // token is either a name, value, or syntactic element coming from a script's // source code type token struct { kind tokenType value string line int pos int } // tokenizer splits a string into a stream tokens, via its `next` method type tokenizer struct { cur string linenum int linepos int } // newTokenizer is the constructor for type tokenizer func newTokenizer(src string) tokenizer { return tokenizer{ cur: src, linenum: 1, linepos: 1, } } // next advances the tokenizer, giving back a token, unless it's done func (t *tokenizer) next() (token, error) { // label to allow looping back after skipping comments, and thus avoid // an explicit tail-recursion for each commented line rerun: // always ignore any whitespace-like source if err := t.skipWhitespace(); err != nil { return token{}, err } if len(t.cur) == 0 { return token{}, errEOS } // remember starting position, in case of error line := t.linenum pos := t.linepos // use the leading rune to probe what's next r, _ := utf8.DecodeRuneInString(t.cur) switch r { case '0', '1', '2', '3', '4', '5', '6', '7', '8', '9': s, err := t.scanNumber() return token{kind: numberToken, value: s, line: line, pos: pos}, err case '(', ')', '[', ']', '{', '}', ',', '+', '-', '%', '^', '~', '@', ';': s := t.cur[:1] t.cur = t.cur[1:] t.linepos++ res := token{kind: syntaxToken, value: s, line: line, pos: pos} return res, t.eos() case ':': s, err := t.tryPrefixes(`:=`, `:`) return token{kind: syntaxToken, value: s, line: line, pos: pos}, err case '*': s, err := t.tryPrefixes(`**`, `*`) return token{kind: syntaxToken, value: s, line: line, pos: pos}, err case '/': s, err := t.tryPrefixes(`//`, `/*`, `/`) // double-slash starts a comment until the end of the line if s == `//` { t.skipLine() goto rerun } // handle comments which can span multiple lines if s == `/*` { err := t.skipComment() if err != nil { return token{}, err } goto rerun } // handle division return token{kind: syntaxToken, value: s, line: line, pos: pos}, err case '#': // even hash starts a comment until the end of the line, making the // syntax more Python-like t.skipLine() goto rerun case '&': s, err := t.tryPrefixes(`&&`, `&`) return token{kind: syntaxToken, value: s, line: line, pos: pos}, err case '|': s, err := t.tryPrefixes(`||`, `|`) return token{kind: syntaxToken, value: s, line: line, pos: pos}, err case '?': s, err := t.tryPrefixes(`??`, `?.`, `?`) return token{kind: syntaxToken, value: s, line: line, pos: pos}, err case '.': r, _ := utf8.DecodeRuneInString(t.cur[1:]) if '0' <= r && r <= '9' { s, err := t.scanNumber() res := token{kind: numberToken, value: s, line: line, pos: pos} return res, err } s, err := t.tryPrefixes(`...`, `..`, `.?`, `.`) return token{kind: syntaxToken, value: s, line: line, pos: pos}, err case '=': // triple-equal makes the syntax even more JavaScript-like s, err := t.tryPrefixes(`===`, `==`, `=>`, `=`) return token{kind: syntaxToken, value: s, line: line, pos: pos}, err case '!': // not-double-equal makes the syntax even more JavaScript-like s, err := t.tryPrefixes(`!==`, `!=`, `!`) return token{kind: syntaxToken, value: s, line: line, pos: pos}, err case '<': // the less-more/diamond syntax is a SQL-like way to say not equal s, err := t.tryPrefixes(`<=`, `<>`, `<`) return token{kind: syntaxToken, value: s, line: line, pos: pos}, err case '>': s, err := t.tryPrefixes(`>=`, `>`) return token{kind: syntaxToken, value: s, line: line, pos: pos}, err default: if isIdentStartRune(r) { s, err := t.scanIdentifier() res := token{kind: identifierToken, value: s, line: line, pos: pos} return res, err } const fs = `line %d: pos %d: unexpected symbol %c` return token{}, fmt.Errorf(fs, t.linenum, t.linepos, r) } } // tryPrefixes tries to greedily match any prefix in the order given: when all // candidates fail, an empty string is returned; when successful, the tokenizer // updates its state to account for the matched prefix func (t *tokenizer) tryPrefixes(prefixes ...string) (string, error) { for _, pre := range prefixes { if strings.HasPrefix(t.cur, pre) { t.linepos += len(pre) t.cur = t.cur[len(pre):] return pre, t.eos() } } return ``, t.eos() } // skipWhitespace updates the tokenizer to ignore runs of consecutive whitespace // symbols: these are the likes of space, tab, newline, carriage return, etc. func (t *tokenizer) skipWhitespace() error { for len(t.cur) > 0 { r, size := utf8.DecodeRuneInString(t.cur) if !unicode.IsSpace(r) { // no more spaces to skip return nil } t.cur = t.cur[size:] if r == '\n' { // reached the next line t.linenum++ t.linepos = 1 continue } // continuing on the same line t.linepos++ } // source code ended return errEOS } // skipLine updates the tokenizer to the end of the current line, or the end of // the source code, if it's the last line func (t *tokenizer) skipLine() error { for len(t.cur) > 0 { r, size := utf8.DecodeRuneInString(t.cur) t.cur = t.cur[size:] if r == '\n' { // reached the next line, as expected t.linenum++ t.linepos = 1 return nil } } // source code ended t.linenum++ t.linepos = 1 return errEOS } // skipComment updates the tokenizer to the end of the comment started with a // `/*` and ending with a `*/`, or to the end of the source code func (t *tokenizer) skipComment() error { var prev rune for len(t.cur) > 0 { r, size := utf8.DecodeRuneInString(t.cur) t.cur = t.cur[size:] if r == '\n' { t.linenum++ t.linepos = 1 } else { t.linepos++ } if prev == '*' && r == '/' { return nil } prev = r } // source code ended const msg = "comment not ended with a `*/` sequence" return errors.New(msg) } func (t *tokenizer) scanIdentifier() (string, error) { end := 0 for len(t.cur) > 0 { r, size := utf8.DecodeRuneInString(t.cur[end:]) if (end == 0 && !isIdentStartRune(r)) || !isIdentRestRune(r) { // identifier ended, and there's more source code after it name := t.cur[:end] t.cur = t.cur[end:] return name, nil } end += size t.linepos++ } // source code ended with an identifier name name := t.cur t.cur = `` return name, nil } func (t *tokenizer) scanNumber() (string, error) { dots := 0 end := 0 var prev rune for len(t.cur) > 0 { r, size := utf8.DecodeRuneInString(t.cur[end:]) switch r { case '0', '1', '2', '3', '4', '5', '6', '7', '8', '9', '_', 'e': end += size t.linepos++ prev = r case '+', '-': if end > 0 && prev != 'e' { // number ended, and there's more source code after it num := t.cur[:end] t.cur = t.cur[end:] return num, nil } end += size t.linepos++ prev = r case '.': nr, _ := utf8.DecodeRuneInString(t.cur[end+size:]) if dots == 1 || isIdentStartRune(nr) || unicode.IsSpace(nr) || nr == '.' { // number ended, and there's more source code after it num := t.cur[:end] t.cur = t.cur[end:] return num, nil } dots++ end += size t.linepos++ prev = r default: // number ended, and there's more source code after it num := t.cur[:end] t.cur = t.cur[end:] return num, nil } } // source code ended with a number name := t.cur t.cur = `` return name, nil } // eos checks if the source-code is over: if so, it returns an end-of-file error, // or a nil error otherwise func (t *tokenizer) eos() error { if len(t.cur) == 0 { return errEOS } return nil } func isIdentStartRune(r rune) bool { return ('A' <= r && r <= 'Z') || ('a' <= r && r <= 'z') || r == '_' } func isIdentRestRune(r rune) bool { return isIdentStartRune(r) || ('0' <= r && r <= '9') } func Decade(x float64) float64 { return 10 * math.Floor(0.1*x) } func Century(x float64) float64 { return 100 * math.Floor(0.01*x) } func Round1(x float64) float64 { return 1e-1 * math.Round(1e+1*x) } func Round2(x float64) float64 { return 1e-2 * math.Round(1e+2*x) } func Round3(x float64) float64 { return 1e-3 * math.Round(1e+3*x) } func Round4(x float64) float64 { return 1e-4 * math.Round(1e+4*x) } func Round5(x float64) float64 { return 1e-5 * math.Round(1e+5*x) } func Round6(x float64) float64 { return 1e-6 * math.Round(1e+6*x) } func Round7(x float64) float64 { return 1e-7 * math.Round(1e+7*x) } func Round8(x float64) float64 { return 1e-8 * math.Round(1e+8*x) } func Round9(x float64) float64 { return 1e-9 * math.Round(1e+9*x) } // Ln calculates the natural logarithm. func Ln(x float64) float64 { return math.Log(x) } // Log calculates logarithms using the base given. // func Log(base float64, x float64) float64 { // return math.Log(x) / math.Log(base) // } // Round rounds the number given to the number of decimal places given. func Round(x float64, decimals int) float64 { k := math.Pow(10, float64(decimals)) return math.Round(k*x) / k } // RoundBy rounds the number given to the unit-size given // func RoundBy(x float64, by float64) float64 { // return math.Round(by*x) / by // } // FloorBy quantizes a number downward by the unit-size given // func FloorBy(x float64, by float64) float64 { // mod := math.Mod(x, by) // if x < 0 { // if mod == 0 { // return x - mod // } // return x - mod - by // } // return x - mod // } // Wrap linearly interpolates the number given to the range [0...1], // according to its continuous domain, specified by the limits given func Wrap(x float64, min, max float64) float64 { return (x - min) / (max - min) } // AnchoredWrap works like Wrap, except when both domain boundaries are positive, // the minimum becomes 0, and when both domain boundaries are negative, the // maximum becomes 0. func AnchoredWrap(x float64, min, max float64) float64 { if min > 0 && max > 0 { min = 0 } else if max < 0 && min < 0 { max = 0 } return (x - min) / (max - min) } // Unwrap undoes Wrap, by turning the [0...1] number given into its equivalent // in the new domain given. func Unwrap(x float64, min, max float64) float64 { return (max-min)*x + min } // Clamp constrains the domain of the value given func Clamp(x float64, min, max float64) float64 { return math.Min(math.Max(x, min), max) } // Scale transforms a number according to its position in [xMin, xMax] into its // correspondingly-positioned counterpart in [yMin, yMax]: if value isn't in its // assumed domain, its result will be extrapolated accordingly func Scale(x float64, xmin, xmax, ymin, ymax float64) float64 { k := (x - xmin) / (xmax - xmin) return (ymax-ymin)*k + ymin } // AnchoredScale works like Scale, except when both domain boundaries are positive, // the minimum becomes 0, and when both domain boundaries are negative, the maximum // becomes 0. This allows for proportionally-correct scaling of quantities, such // as when showing data visually. func AnchoredScale(x float64, xmin, xmax, ymin, ymax float64) float64 { if xmin > 0 && xmax > 0 { xmin = 0 } else if xmax < 0 && xmin < 0 { xmax = 0 } k := (x - xmin) / (xmax - xmin) return (ymax-ymin)*k + ymin } // ForceRange is handy when trying to mold floating-point values into numbers // valid for JSON, since NaN and Infinity are replaced by the values given; // the infinity-replacement value is negated for negative infinity. func ForceRange(x float64, nan, inf float64) float64 { if math.IsNaN(x) { return nan } if math.IsInf(x, -1) { return -inf } if math.IsInf(x, +1) { return inf } return x } // Sign returns the standardized sign of a value, either as -1, 0, or +1: NaN // values stay as NaN, as is expected when using floating-point values. func Sign(x float64) float64 { if x > 0 { return +1 } if x < 0 { return -1 } if x == 0 { return 0 } return math.NaN() } // IsInteger checks if a floating-point value is already rounded to an integer // value. func IsInteger(x float64) bool { _, frac := math.Modf(x) return frac == 0 } func Square(x float64) float64 { return x * x } func Cube(x float64) float64 { return x * x * x } func RoundTo(x float64, decimals float64) float64 { k := math.Pow10(int(decimals)) return math.Round(k*x) / k } func RelDiff(x, y float64) float64 { return (x - y) / y } // Bool booleanizes a number: 0 stays 0, and anything else becomes 1. func Bool(x float64) float64 { if x == 0 { return 0 } return 1 } // DeNaN replaces NaN with the alternative value given: regular values are // returned as given. func DeNaN(x float64, instead float64) float64 { if !math.IsNaN(x) { return x } return instead } // DeInf replaces either infinity with the alternative value given: regular // values are returned as given. func DeInf(x float64, instead float64) float64 { if !math.IsInf(x, 0) { return x } return instead } // Revalue replaces NaN and either infinity with the alternative value given: // regular values are returned as given. func Revalue(x float64, instead float64) float64 { if !math.IsNaN(x) && !math.IsInf(x, 0) { return x } return instead } // Linear evaluates a linear polynomial, the first arg being the main input, // followed by the polynomial coefficients in decreasing-power order. func Linear(x float64, a, b float64) float64 { return a*x + b } // Quadratic evaluates a 2nd-degree polynomial, the first arg being the main // input, followed by the polynomial coefficients in decreasing-power order. func Quadratic(x float64, a, b, c float64) float64 { return (a*x+b)*x + c } // Cubic evaluates a cubic polynomial, the first arg being the main input, // followed by the polynomial coefficients in decreasing-power order. func Cubic(x float64, a, b, c, d float64) float64 { return ((a*x+b)*x+c)*x + d } func LinearFMA(x float64, a, b float64) float64 { return math.FMA(x, a, b) } func QuadraticFMA(x float64, a, b, c float64) float64 { lin := math.FMA(x, a, b) return math.FMA(lin, x, c) } func CubicFMA(x float64, a, b, c, d float64) float64 { lin := math.FMA(x, a, b) quad := math.FMA(lin, x, c) return math.FMA(quad, x, d) } // Radians converts angular degrees into angular radians: 180 degrees are pi // pi radians. func Radians(deg float64) float64 { const k = math.Pi / 180 return k * deg } // Degrees converts angular radians into angular degrees: pi radians are 180 // degrees. func Degrees(rad float64) float64 { const k = 180 / math.Pi return k * rad } // Fract calculates the non-integer/fractional part of a number. func Fract(x float64) float64 { return x - math.Floor(x) } // Mix interpolates 2 numbers using a third number, used as an interpolation // coefficient. This parameter naturally falls in the range [0, 1], but doesn't // have to be: when given outside that range, the parameter can extrapolate in // either direction instead. func Mix(x, y float64, k float64) float64 { return (1-k)*(y-x) + x } // Step implements a step function with a parametric threshold. func Step(edge, x float64) float64 { if x < edge { return 0 } return 1 } // SmoothStep is like the `smoothstep` func found in GLSL, using a cubic // interpolator in the transition region. func SmoothStep(edge0, edge1, x float64) float64 { if x <= edge0 { return 0 } if x >= edge1 { return 1 } // use the cubic hermite interpolator 3x^2 - 2x^3 in the transition band return x * x * (3 - 2*x) } // Logistic approximates the math func of the same name. func Logistic(x float64) float64 { return 1 / (1 + math.Exp(-x)) } // Sinc approximates the math func of the same name. func Sinc(x float64) float64 { if x != 0 { return math.Sin(x) / x } return 1 } // Sum adds all the numbers in an array. func Sum(v ...float64) float64 { s := 0.0 for _, f := range v { s += f } return s } // Product multiplies all the numbers in an array. func Product(v ...float64) float64 { p := 1.0 for _, f := range v { p *= f } return p } // Length calculates the Euclidean length of the vector given. func Length(v ...float64) float64 { ss := 0.0 for _, f := range v { ss += f * f } return math.Sqrt(ss) } // Dot calculates the dot product of 2 vectors func Dot(x []float64, y []float64) float64 { l := len(x) if len(y) < l { l = len(y) } dot := 0.0 for i := 0; i < l; i++ { dot += x[i] * y[i] } return dot } // Min finds the minimum value from the numbers given. func Min(v ...float64) float64 { min := +math.Inf(+1) for _, f := range v { min = math.Min(min, f) } return min } // Max finds the maximum value from the numbers given. func Max(v ...float64) float64 { max := +math.Inf(-1) for _, f := range v { max = math.Max(max, f) } return max } // Hypot calculates the Euclidean n-dimensional hypothenuse from the numbers // given: all numbers can be lengths, or simply positional coordinates. func Hypot(v ...float64) float64 { sumsq := 0.0 for _, f := range v { sumsq += f * f } return math.Sqrt(sumsq) } // Polyval evaluates a polynomial using Horner's algorithm. The array has all // the coefficients in textbook order, from the highest power down to the final // constant. func Polyval(x float64, v ...float64) float64 { if len(v) == 0 { return 0 } x0 := x x = 1.0 y := 0.0 for i := len(v) - 1; i >= 0; i-- { y += v[i] * x x *= x0 } return y } // LnGamma is a 1-input 1-output version of math.Lgamma from the stdlib. func LnGamma(x float64) float64 { y, s := math.Lgamma(x) if s < 0 { return math.NaN() } return y } // LnBeta calculates the natural-logarithm of the beta function. func LnBeta(x float64, y float64) float64 { return LnGamma(x) + LnGamma(y) - LnGamma(x+y) } // Beta calculates the beta function. func Beta(x float64, y float64) float64 { return math.Exp(LnBeta(x, y)) } // Factorial calculates the product of all integers in [1, n] func Factorial(n int) int64 { return int64(math.Round(math.Gamma(float64(n + 1)))) } // IsPrime checks whether an integer is bigger than 1 and can only be fully // divided by 1 and itself, which is the definition of a prime number. func IsPrime(n int64) bool { // prime numbers start at 2 if n < 2 { return false } // 2 is the only even prime if n%2 == 0 { return n == 2 } // no divisor can be more than the square root of the target number: // this limit makes the loop an O(sqrt(n)) one, instead of O(n); this // is a major algorithmic speedup both in theory and in practice max := int64(math.Floor(math.Sqrt(float64(n)))) // the only possible full-divisors are odd integers 3..sqrt(n), // since reaching this point guarantees n is odd and n > 2 for d := int64(3); d <= max; d += 2 { if n%d == 0 { return false } } return true } // LCM finds the least common-multiple of 2 positive integers; when one or // both inputs aren't positive, this func returns 0. func LCM(x, y int64) int64 { if gcd := GCD(x, y); gcd > 0 { return x * y / gcd } return 0 } // GCD finds the greatest common-divisor of 2 positive integers; when one or // both inputs aren't positive, this func returns 0. func GCD(x, y int64) int64 { if x < 1 || y < 1 { return 0 } // the loop below requires a >= b a, b := x, y if a < b { a, b = y, x } for b > 0 { a, b = b, a%b } return a } // Perm counts the number of all possible permutations from n objects when // picking k times. When one or both inputs aren't positive, the result is 0. func Perm(n, k int) int64 { if n < k || n < 0 || k < 0 { return 0 } perm := int64(1) for i := n - k + 1; i <= n; i++ { perm *= int64(i) } return perm } // Choose counts the number of all possible combinations from n objects when // picking k times. When one or both inputs aren't positive, the result is 0. func Choose(n, k int) int64 { if n < k || n < 0 || k < 0 { return 0 } // the log trick isn't always more accurate when there's no overflow: // for those cases calculate using the textbook definition f := math.Round(float64(Perm(n, k) / Factorial(k))) if !math.IsInf(f, 0) { return int64(f) } // calculate using the log-factorial of n, k, and n - k a, _ := math.Lgamma(float64(n + 1)) b, _ := math.Lgamma(float64(k + 1)) c, _ := math.Lgamma(float64(n - k + 1)) return int64(math.Round(math.Exp(a - b - c))) } // BinomialMass calculates the probability mass of the binomial random process // given. When the probability given isn't between 0 and 1, the result is NaN. func BinomialMass(x, n int, p float64) float64 { // invalid probability input if p < 0 || p > 1 { return math.NaN() } // events outside the support are impossible by definition if x < 0 || x > n { return 0 } q := 1 - p m := n - x ncx := float64(Choose(n, x)) return ncx * math.Pow(p, float64(x)) * math.Pow(q, float64(m)) } // CumulativeBinomialDensity calculates cumulative probabilities/masses up to // the value given for the binomial random process given. When the probability // given isn't between 0 and 1, the result is NaN. func CumulativeBinomialDensity(x, n int, p float64) float64 { // invalid probability input if p < 0 || p > 1 { return math.NaN() } if x < 0 { return 0 } if x >= n { return 1 } p0 := p q0 := 1 - p0 q := math.Pow(q0, float64(n)) pbinom := 0.0 np1 := float64(n + 1) for k := 0; k < x; k++ { a, _ := math.Lgamma(np1) b, _ := math.Lgamma(float64(k + 1)) c, _ := math.Lgamma(float64(n - k + 1)) // count all possible combinations for this event ncomb := math.Round(math.Exp(a - b - c)) pbinom += ncomb * p * q p *= p0 q /= q0 } return pbinom } // NormalDensity calculates the density at a point along the normal distribution // given. func NormalDensity(x float64, mu, sigma float64) float64 { z := (x - mu) / sigma return math.Sqrt(0.5/sigma) * math.Exp(-(z*z)/sigma) } // CumulativeNormalDensity calculates the probability of a normal variate of // being up to the value given. func CumulativeNormalDensity(x float64, mu, sigma float64) float64 { z := (x - mu) / sigma return 0.5 + 0.5*math.Erf(z/math.Sqrt2) } // Epanechnikov is a commonly-used kernel function. func Epanechnikov(x float64) float64 { if math.Abs(x) > 1 { // func is 0 ouside -1..+1 return 0 } return 0.75 * (1 - x*x) } // Gauss is the commonly-used Gaussian kernel function. func Gauss(x float64) float64 { return math.Exp(-(x * x)) } // Tricube is a commonly-used kernel function. func Tricube(x float64) float64 { a := math.Abs(x) if a > 1 { // func is 0 ouside -1..+1 return 0 } b := a * a * a c := 1 - b return 70.0 / 81.0 * c * c * c } // SolveQuad finds the solutions of a 2nd-degree polynomial, using a formula // which is more accurate than the textbook one. func SolveQuad(a, b, c float64) (x1 float64, x2 float64) { div := 2 * c disc := math.Sqrt(b*b - 4*a*c) x1 = div / (-b - disc) x2 = div / (-b + disc) return x1, x2 }