hyperrogue/util.cpp

// Hyperbolic Rogue -- basic utility functions
// Copyright (C) 2011-2018 Zeno Rogue, see 'hyper.cpp' for details

/** \file util.cpp
 *  \brief basic utility functions: maths, parsing expressions
 */

#include "hyper.h"
namespace hr {

#if CAP_TIMEOFDAY
#if !CAP_SDL
int lastusec;
int uticks;

EX int SDL_GetTicks() {
  struct timeval tim;
  gettimeofday(&tim, NULL);
  int newusec = tim.tv_usec;
  uticks += newusec - lastusec;
  if(newusec <= lastusec)
    uticks += 1000000;
  lastusec = newusec;
  return uticks / 1000;
  }

#endif
#endif

EX long double sqr(long double x) { return x*x; }

EX ld round_nearest(ld x) { if(x > 0) return int(x+.5); else return -int(.5-x); }
EX ld round_nearest(ld x, ld multiple_of) { return multiple_of * round_nearest(x / multiple_of); }

EX int gcd(int i, int j) {
  return i ? gcd(j%i, i) : j;
  }

EX int gmod(int i, int j) {
  i %= j; if(i<0) i += j;
  return i;
  }

EX int zgmod(int a, int b) { return b ? gmod(a, b) : a; }

EX int szgmod(int a, int b) { 
  if(!b) return a;
  a = gmod(a, b);
  if(2*a >= b) return a - b;
  return a;
  }

EX int gdiv(int i, int j) {
  return (i - gmod(i, j)) / j;
  }

EX ld frac(ld x) {
  x -= int(x);
  if(x < 0) x++;
  return x;
  }

EX ld lerp(ld a0, ld a1, ld x) {
  return a0 + (a1-a0) * x;
  }

EX cld lerp(cld a0, cld a1, ld x) {
  return a0 + (a1-a0) * x;
  }

EX ld ilerp(ld a0, ld a1, ld x) {
  return (x-a0) / (a1-a0);
  }

// debug utilities

#if CAP_PROFILING

#define FRAMES 64
#define CATS 16

long long proftable[16][FRAMES];
int pframeid;

void profile_frame() { 
  pframeid++; pframeid %=  FRAMES;
  for(int t=0; t<16; t++) proftable[t][pframeid] = 0;
  }

EX void profile_start(int t) { proftable[t][pframeid] -= getms(); }
EX void profile_stop(int t) { proftable[t][pframeid] += getms(); }

void profile_info() {
  for(int t=0; t<16; t++) {
    sort(proftable[t], proftable[t]+FRAMES);
    if(proftable[t][FRAMES-1] == 0) continue;
    long long sum = 0;
    for(int f=0; f<FRAMES; f++) sum += proftable[t][f];
    printf("Category %d: avg = %Ld, %Ld..%Ld..%Ld..%Ld..%Ld\n",
      t, sum / FRAMES, proftable[t][0], proftable[t][16], proftable[t][32],
      proftable[t][48], proftable[t][63]);
    }
  }
#endif

#if !CAP_PROFILING
#if HDR
#define profile_frame()
#define profile_start(t)
#define profile_stop(t)
#define profile_info()
#endif
#endif

EX purehookset hooks_tests;

EX string simplify(const string& s) {
  string res;
  for(char c: s) if(isalnum(c)) res += c;
  return res;
  }

EX bool appears(const string& haystack, const string& needle) {
  return simplify(haystack).find(simplify(needle)) != string::npos;
  }

#if HDR
struct hr_parse_exception : hr_exception {
  string s;
  hr_parse_exception(const string& z) : s(z) {}
  ~hr_parse_exception() noexcept(true) {}
  };

struct exp_parser {
  string s;
  int at;
  int line_number, last_line;
  exp_parser() { at = 0; line_number = 1; last_line = 0; }
  
  string where() { 
    if(s.find('\n')) return "(line " + its(line_number) + ", pos " + its(at-last_line) + ")";
    else return "(pos " + its(at) + ")";
    }
  
  map<string, cld> extra_params;

  bool ok() { return at == isize(s); }
  char next(int step=0) { if(at >= isize(s)-step) return 0; else return s[at+step]; }
  
  bool eat(const char *c) {
    int orig_at = at;
    while(*c && *c == next()) at++, c++;
    if(*c == 0) return true;
    else at = orig_at;
    return false;
    }

  void skip_white();

  string next_token();

  char snext(int step=0) { skip_white(); return next(step); }

  cld parse(int prio = 0);

  ld rparse(int prio = 0) { return validate_real(parse(prio)); }
  int iparse(int prio = 0) { return int(floor(rparse(prio) + .5)); }

  cld parsepar() {
    cld res = parse();
    force_eat(")");
    return res;
    }

  ld validate_real(cld x) {
    if(kz(imag(x))) throw hr_parse_exception("expected real number but " + lalign(-1, x) + " found at " + where());
    return real(x);
    }
  
  void force_eat(const char *c) {
    skip_white();
    if(!eat(c)) throw hr_parse_exception("expected: " + string(c) + " at " + where());
    }

  };
#endif

void exp_parser::skip_white() {
  while(next() == ' ' || next() == '\n' || next() == '\r' || next() == '\t') {
    if(next() == '\r') last_line++;
    if(next() == '\n') {
      line_number++, last_line = at;
      }
    at++;
    }
  }

string exp_parser::next_token() {
  skip_white();
  string token;
  while(true) {
    char c = next();
    if((c >= '0' && c <= '9') || (c == '.' && next(1) != '.') || (c >= 'a' && c <= 'z') || (c >= 'A' && c <= 'Z') || c == '_')
      token += c, at++;
    else break;
    }
  return token;
  }

cld exp_parser::parse(int prio) {
  cld res;
  skip_white();
  if(eat("sin(")) res = sin(parsepar());
  else if(eat("cos(")) res = cos(parsepar());
  else if(eat("sinh(")) res = sinh(parsepar());
  else if(eat("cosh(")) res = cosh(parsepar());
  else if(eat("asin(")) res = asin(parsepar());
  else if(eat("acos(")) res = acos(parsepar());
  else if(eat("asinh(")) res = asinh(parsepar());
  else if(eat("acosh(")) res = acosh(parsepar());
  else if(eat("exp(")) res = exp(parsepar());
  else if(eat("sqrt(")) res = sqrt(parsepar());
  else if(eat("log(")) res = log(parsepar());
  else if(eat("tan(")) res = tan(parsepar());
  else if(eat("tanh(")) res = tanh(parsepar());
  else if(eat("atan(")) res = atan(parsepar());
  else if(eat("atanh(")) res = atanh(parsepar());
  else if(eat("abs(")) res = abs(parsepar());
  else if(eat("re(")) res = real(parsepar());
  else if(eat("im(")) res = imag(parsepar());
  else if(eat("conj(")) res = std::conj(parsepar());
  else if(eat("floor(")) res = floor(validate_real(parsepar()));
  else if(eat("frac(")) { res = parsepar(); res = res - floor(validate_real(res)); }
  else if(eat("to01(")) { res = parsepar(); return atan(res) / ld(M_PI) + ld(0.5); }
  else if(eat("edge(")) {
    ld a = rparse(0);
    force_eat(",");
    ld b = rparse(0);
    force_eat(")");
    res = edge_of_triangle_with_angles(2*M_PI/a, M_PI/b, M_PI/b);
    }
  else if(eat("edge_angles(")) {
    cld a = rparse(0);
    force_eat(",");
    cld b = rparse(0);
    force_eat(",");
    cld c = rparse(0);
    force_eat(")");

    if(extra_params.count("angleunit")) {    
      
      a *= extra_params["angleunit"];
      b *= extra_params["angleunit"];
      c *= extra_params["angleunit"];
      }

    return edge_of_triangle_with_angles(real(a), real(b), real(c));
    }
  else if(eat("regradius(")) {
    ld a = rparse(0);
    force_eat(",");
    ld b = rparse(0);
    force_eat(")");
    res = edge_of_triangle_with_angles(M_PI/2, M_PI/a, M_PI/b);
    }
  else if(eat("arcmedge(")) {
    vector<int> vals;
    vals.push_back(iparse(0));
    while(true) {
      skip_white();
      if(eat(",")) vals.push_back(iparse(0));
      else break;
      }
    force_eat(")");
    arcm::archimedean_tiling test;
    test.faces = vals;
    test.compute_sum();
    test.compute_geometry();
    res = test.edgelength;
    if(extra_params.count("distunit"))
      res /= extra_params["distunit"];
    }
  else if(eat("regangle(")) {
    cld edgelen = parse(0);
    if(extra_params.count("distunit")) {
      edgelen = edgelen * extra_params["distunit"];
      }
    
    force_eat(",");
    ld edges = rparse(0);
    force_eat(")");
    ld alpha = M_PI / edges;
    ld c = asin_auto(sin_auto(validate_real(edgelen)/2) / sin(alpha));
    hyperpoint h = xpush(c) * spin(M_PI - 2*alpha) * xpush0(c);
    ld result = 2 * atan2(h);
    if(result < 0) result = -result;
    while(result > 2 * M_PI) result -= 2 * M_PI;
    if(result > M_PI) result = 2 * M_PI - result;
    
    if(arb::legacy) {
      res = M_PI - result;
      if(extra_params.count("angleofs"))
        res -= extra_params["angleofs"];
      }
    else
      res = result;

    if(extra_params.count("angleunit"))
      res /= extra_params["angleunit"];    
    }
  else if(eat("test(")) {
    res = parsepar();
    println(hlog, "res = ", res, ": ", fts(real(res), 10), ",", fts(imag(res), 10));
    }
  else if(eat("ifp(")) {
    cld cond = parse(0);
    force_eat(",");
    cld yes = parse(0);
    force_eat(",");
    cld no = parsepar();
    res = real(cond) > 0 ? yes : no;
    }  
  else if(eat("wallif(")) {
    cld val0 = parse(0);
    force_eat(",");
    cld val1 = parsepar();
    if(real(extra_params["p"]) >= 3.5) res = val0;
    else res = val1;
    }
  else if(eat("rgb(")) {     
    cld val0 = parse(0);
    force_eat(",");
    cld val1 = parse(0);
    force_eat(",");
    cld val2 = parsepar();
    switch(int(real(extra_params["p"]) + .5)) {
      case 1: res = val0; break;
      case 2: res = val1; break;
      case 3: res = val2; break;
      default: res = 0;
      }
    }
  else if(eat("let(")) {
    string name = next_token();
    force_eat("=");
    cld val = parse(0);
    force_eat(",");
    dynamicval<cld> d(extra_params[name], val);
    res = parsepar();
    }
  #if CAP_TEXTURE
  else if(eat("txp(")) {
    cld val = parsepar();
    res = texture::get_txp(real(val), imag(val), int(real(extra_params["p"]) + .5)-1);
    }
  #endif
  else if(next() == '(') at++, res = parsepar(); 
  else {
    string number = next_token();
    if(extra_params.count(number)) res = extra_params[number];
    else if(params.count(number)) res = params.at(number);
    else if(number == "e") res = exp(1);
    else if(number == "i") res = cld(0, 1);
    else if(number == "p" || number == "pi") res = M_PI;
    else if(number == "" && next() == '-') { at++; res = -parse(prio); }
    else if(number == "") throw hr_parse_exception("number missing, " + where());
    else if(number == "s") res = ticks / 1000.;
    else if(number == "ms") res = ticks;
    else if(number[0] == '0' && number[1] == 'x') res = strtoll(number.c_str()+2, NULL, 16);
    else if(number == "mousex") res = mousex;
    else if(number == "deg") res = degree;
    else if(number == "ultra_mirror_dist") res = cgi.ultra_mirror_dist;
    else if(number == "step") res = hdist0(tC0(currentmap->adj(cwt.at, 0)));
    else if(number == "mousey") res = mousey;
    else if(number == "random") res = randd();
    else if(number == "mousez") res = cld(mousex - current_display->xcenter, mousey - current_display->ycenter) / cld(current_display->radius, 0);
    else if(number == "shot") res = inHighQual ? 1 : 0;
    else if(number[0] >= 'a' && number[0] <= 'z') throw hr_parse_exception("unknown value: " + number);
    else { std::stringstream ss; res = 0; ss << number; ss >> res; }
    }
  while(true) {
    skip_white();
    #if CAP_ANIMATIONS
    if(next() == '.' && next(1) == '.' && prio == 0) {
      static const cld NO_DERIVATIVE(3.1, 2.5);
      vector<array<cld, 4>> rest = { make_array(res, NO_DERIVATIVE, res, NO_DERIVATIVE) };
      bool second = true;
      while(next() == '.' && next(1) == '.') {
        /* spline interpolation */
        if(next(2) == '/') {
          at += 3;
          rest.back()[second ? 3 : 1] = parse(10);
          continue;
          }
        /* sharp end */
        else if(next(2) == '|') {
          at += 3;
          rest.back()[2] = parse(10);
          rest.back()[3] = NO_DERIVATIVE;
          second = true;
          continue;
          }
        at += 2; 
        auto val = parse(10);
        rest.emplace_back(make_array(val, NO_DERIVATIVE, val, NO_DERIVATIVE));
        second = false;
        }
      ld v = ticks * (isize(rest)-1.) / anims::period;
      int vf = v;
      v -= vf;
      vf %= (isize(rest)-1);
      auto& lft = rest[vf];
      auto& rgt = rest[vf+1];
      if(lft[3] == NO_DERIVATIVE && rgt[1] == NO_DERIVATIVE)
        res = lerp(lft[2], rgt[0], v);
      else if(rgt[1] == NO_DERIVATIVE)
        res = lerp(lft[2] + lft[3] * v, rgt[0], v*v);
      else if(lft[3] == NO_DERIVATIVE)
        res = lerp(lft[2], rgt[0] + rgt[1] * (v-1), (2-v)*v);
      else {
        res = lerp(lft[2] + lft[3] * v, rgt[0] + rgt[1] * (v-1), v*v*(3-2*v));
        }
      return res;
      }
    else 
    #endif
    if(next() == '+' && prio <= 10) at++, res = res + parse(20);
    else if(next() == '-' && prio <= 10) at++, res = res - parse(20);
    else if(next() == '*' && prio <= 20) at++, res = res * parse(30);
    else if(next() == '/' && prio <= 20) at++, res = res / parse(30);
    else if(next() == '^') at++, res = pow(res, parse(40));
    else break;
    }
  return res;
  }

EX ld parseld(const string& s) {
  exp_parser ep;
  ep.s = s;
  return real(ep.parse());
  }

EX string parser_help() {
  return XLAT("Functions available: %1", 
    "(a)sin(h), (a)cos(h), (a)tan(h), exp, log, abs, re, im, conj, let(t=...,...t...), floor, frac, e, i, pi, s, ms, mousex, mousey, mousez, shot [1 if taking screenshot/animation], sqrt, to01, random, edge(7,3), regradius(7,3), ifp(a,v,w) [if positive]");
  }

#if HDR
struct bignum {
  static const int BASE = 1000000000;
  static const long long BASE2 = BASE * (long long)BASE;
  vector<int> digits;
  bignum() {}
  bignum(int i) : digits() { digits.push_back(i); }
  void be(int i) { digits.resize(1); digits[0] = i; }
  bignum& operator +=(const bignum& b);
  void addmul(const bignum& b, int factor);
  string get_str(int max_length) const;
  bignum(ld d);
  
  bool operator < (const bignum&) const;
  bool operator > (const bignum& b) const { return b < self; }

  ld leading() const {
    switch(isize(digits)) {
      case 0:
        return 0;
      case 1:
        return digits.back();
      default:
        return digits.back() + ld(digits[isize(digits)-2]) / BASE;
      }
    }

  ld approx() const {
    return leading() * pow(BASE, isize(digits) - 1);
    }
  
  ld log_approx() const {
    return log(leading()) * log(BASE) * (isize(digits) - 1);
    }
  
  ld approx_div(const bignum& b) const {
    return leading() / b.leading() * pow(BASE, isize(digits) - isize(b.digits));
    }
  
  int approx_int() const {
    if(isize(digits) > 1) return BASE;
    if(digits.empty()) return 0;
    return digits[0];
    }
  
  bool nonzero() { return approx_ld() != 0; }
  
  bignum randomized_div(int x) const;

  ld approx_ld() const {
    ld res = 0;
    for(int i=0; i<isize(digits); i++) res += digits[i] * pow(BASE, i);
    return res;
    }
    
  long long approx_ll() const {
    if(isize(digits) > 2) return BASE2;
    if(digits.empty()) return 0;
    if(isize(digits) == 1) return digits[0];
    return digits[0] + digits[1] * (long long) BASE;
    }
  
  #if CAP_GMP
  mpq_class as_mpq() const {
    string s = get_str(999999);
    string t;
    for(char c: s) if(c != ' ') t += c;
    return mpq_class(t);
    }
  #endif
  
  friend inline bignum operator +(bignum a, const bignum& b) { a.addmul(b, 1); return a; }
  friend inline bignum operator -(bignum a, const bignum& b) { a.addmul(b, -1); return a; }
  };
#endif

bignum& bignum::operator +=(const bignum& b) {
  int K = isize(b.digits);
  if(K > isize(digits)) digits.resize(K);
  int carry = 0;
  for(int i=0; i<K || carry; i++) {
    if(i >= isize(digits)) digits.push_back(0);
    digits[i] += carry;
    if(i < K) digits[i] += b.digits[i];
    if(digits[i] >= BASE) {
      digits[i] -= BASE;
      carry = 1;
      }
    else carry = 0;
    }
  return *this;
  }

bool bignum::operator < (const bignum& b) const {
  if(isize(digits) != isize(b.digits))
    return isize(digits) < isize(b.digits);
  for(int i = isize(digits)-1; i>=0; i--)
    if(digits[i] != b.digits[i])
      return digits[i] < b.digits[i];
  return false;
  }

bignum bignum::randomized_div(int x) const {
  bignum res = self;
  long long carry = 0;
  int K = isize(res.digits);
  for(int i=K-1; i>=0; i--) {
    carry *= BASE;
    carry += digits[i];
    // strange compiler buug:
    // if I do / and %, function 'divmod' is called, and it complains on launch that divmod is unimplemented
    res.digits[i] = int(carry / x);
    carry -= res.digits[i] * (long long)(x);
    }
  while(isize(res.digits) && res.digits.back() == 0) res.digits.pop_back();
  if(rand() % x < carry) res += 1;
  return res;
  }

void bignum::addmul(const bignum& b, int factor) {
  int K = isize(b.digits);
  if(K > isize(digits)) digits.resize(K);
  int carry = 0;
  for(int i=0; i<K || (carry > 0 || carry < -1) || (carry == -1 && i < isize(digits)); i++) {
    if(i >= isize(digits)) digits.push_back(0);
    long long l = digits[i];
    l += carry;
    if(i < K) l += b.digits[i] * factor;
    carry = 0;
    if(l >= BASE) carry = int(l / BASE);
    if(l < 0) carry = -int((BASE-1-l) / BASE);
    l -= carry * BASE;
    digits[i] = int(l);
    }
  if(carry < 0) digits.back() -= BASE;
  while(isize(digits) && digits.back() == 0) digits.pop_back();
  }

EX bignum hrand(bignum b) {
  bignum res;
  int d = isize(b.digits);
  while(true) {
    res.digits.resize(d);
    for(int i=0; i<d-1; i++) res.digits[i] = hrand(bignum::BASE);
    res.digits.back() = hrand(b.digits.back() + 1);
    if(res < b) return res;
    }  
  }

EX void operator ++(bignum &b, int) {
  int i = 0;
  while(true) {
    if(isize(b.digits) == i) { b.digits.push_back(1); break; }
    else if(b.digits[i] == bignum::BASE-1) {
      b.digits[i] = 0;
      i++;
      }
    else {
      b.digits[i]++;
      break;
      }      
    }
  }

EX void operator --(bignum &b, int) {
  int i = 0;
  while(true) {
    if(isize(b.digits) == i) { b.digits.push_back(bignum::BASE-1); break; }
    else if(b.digits[i] == 0) {
      b.digits[i] = bignum::BASE-1;
      i++;
      }
    else {
      b.digits[i]--;
      break;
      }      
    }
  }

string bignum::get_str(int max_length) const {
  if(digits.empty()) return "0";
  string ret = its(digits.back());
  for(int i=isize(digits)-2; i>=0; i--) {
    if(isize(ret) > max_length && i) {
      ret += XLAT(" (%1 more digits)", its(9 * (i+1)));
      return ret;
      }

    ret += " ";
    string val = its(digits[i]);
    while(isize(val) < 9) val = "0" + val;
    ret += val;
    }
  return ret;
  }

EX string short_form(bignum b) {
  if(b < 0) return "-" + short_form(0-b);
  else if(b < 100000) return its(b.approx_int());
  else {
    long long val;
    int q;
    if(isize(b.digits) >= 2) {
      q = max(isize(b.digits) - 2, 0);
      val = b.digits[q] + (long long)(bignum::BASE) * b.digits[q+1];
      }
    else {
      q = 0;
      val = b.digits[0];
      }

    int digits = q * 9;
    while(val >= 1000) { val /= 10; digits++; }
    string str = its(int(val)) + "E" + its(digits + 2);
    str.insert(1, ".");
    return str;
    }
  }

bignum::bignum(ld d) {
  if(d == 0) return;
  int n = 1;
  while(d > BASE) d /= BASE, n++;
  digits.resize(n);
  n--;
  while(n >= 0) { digits[n] = int(d); d -= digits[n]; d *= BASE; n--; }
  }

#if CAP_ZLIB
/* compression/decompression */

EX string compress_string(string s) {
  z_stream strm;
  strm.zalloc = Z_NULL;
  strm.zfree = Z_NULL;
  strm.opaque = Z_NULL;
  println(hlog, "pre init");
  auto ret = deflateInit(&strm, 9);
  if(ret != Z_OK) throw "z-error";
  println(hlog, "init ok");
  strm.avail_in = isize(s);
  strm.next_in = (Bytef*) &s[0];
  vector<char> buf(1000000, 0);
  strm.avail_out = 1000000;
  strm.next_out = (Bytef*) &buf[0];
  if(deflate(&strm, Z_FINISH) != Z_STREAM_END) throw "z-error-2";
  println(hlog, "deflate ok");
  string out(&buf[0], (char*)(strm.next_out) - &buf[0]);
  println(hlog, isize(s), " -> ", isize(out));
  return out;
  }

EX string decompress_string(string s) {
  z_stream strm;
  strm.zalloc = Z_NULL;
  strm.zfree = Z_NULL;
  strm.opaque = Z_NULL;
  auto ret = inflateInit(&strm);
  if(ret != Z_OK) throw "z-error";
  strm.avail_in = isize(s);
  strm.next_in = (Bytef*) &s[0];
  vector<char> buf(1000000, 0);
  strm.avail_out = 1000000;
  strm.next_out = (Bytef*) &buf[0];
  if(inflate(&strm, Z_FINISH) != Z_STREAM_END) throw "z-error-2";
  string out(&buf[0], (char*)(strm.next_out) - &buf[0]);
  println(hlog, isize(s), " -> ", isize(out));
  return out;
  }
#endif

}