text:
read("metaMPFR_common.mpl"):

FUNCTION_SERIES := 0:
FUNCTION_SERIES_RATIONAL := %+1:
CONSTANT_SERIES := %+1:

###############################################################################
######################### We can now generate the code ########################
###############################################################################
# This procedure generate code for a straightforward evaluation of the series
# corresponding to a recurrence.
#   rec is a recurrence with inital conditions
#   type can take three values, depending on which series should be evaluated:
#     FUNCTION_SERIES -> produces code for evaluating sum( a(i)*x^i )
#     FUNCTION_SERIES_RATIONAL -> produces code for evaluating sum( a(i)*(p/q)^i )
#     CONSTANT_SERIES -> produces code for evaluating sum( a(i) )
#   name is the name the should be given to the produced procedure.
#   fofx is an optional parameter. It is the function being implemented.
#     -> if provided, this argument will be used to (heuristically) find the limits
#        of the function at +/-oo and find its asymptotical behavior.
#   *IMPORTANT NOTE*: it must be a function of the variable 'x'.
#                     Moreover, neither _f nor x shall be assigned at the time of
#                     calling the function
#
generateStraightforwardAlgo := proc(rec, aofn, type, name, filename, fofx := _f(x))
  local a, b, n, fd, init_cond, nc, d, exponents, formalRec, c, s1, s2, p, q, hardconstant, f0, i0, ri0, i, j, var, var1, var2, var3, guard_bits, required_mulsi, required_divsi, temp, error_counter, error_in_loop, maxofti, additional_error:

  getname(aofn, a, n):

  required_mulsi := {}:
  required_divsi := {}:
  fd := fopen(filename, WRITE):

  # We check that we have a recurrence of the form a(n+d)=r(n)*a(n)
  # and we extract its canonical form:
  #    a(n+d) = c * s1(n)/s1(n+d) * s2(n+d)/s2(n) * p(n)/q(n) * a(n)
  #
  formalRec := checkOneTermRecurrence(rec, a(n)):
  decomposeOneTermRecurrence(formalRec, c, s1, s2, p, q):
  d := formalRec:-order:

  # We keep only non-trivial inital conditions
  init_cond := removeTrivialConditions(rec, a(n));
  nc := nops(init_cond):
  exponents := {1, d}:
  for i from 1 to nc do
    exponents := { op(exponents), init_cond[i][1] }
  od:
  exponents := findFixpointOfDifferences(exponents):
  exponents := exponents minus {0}:

  error_counter := []:

  fprintf(fd, "/* Evaluation by a straightforward algorithm */\n"):
  fprintf(fd, "/* Code automatically generated by metaMPFR. */\n"):

  fprintf(fd, "static int\n"):
  if (type=FUNCTION_SERIES)
    then fprintf(fd, "mpfr_%a (mpfr_ptr res, mpfr_srcptr x, mpfr_rnd_t rnd)\n", name):
  elif (type=FUNCTION_SERIES_RATIONAL)
    then fprintf(fd, "mpfr_%a (mpfr_ptr res, int u, int v, mpfr_rnd_t rnd)\n", name):
  elif (type=CONSTANT_SERIES)
    then fprintf(fd, "mpfr_%a (mpfr_ptr res, mpfr_rnd_t rnd)\n", name):
  fi:
  fprintf(fd, "{\n"):

  ######################################################
  #################### Declarations ####################
  ######################################################

  fprintf(fd, "  MPFR_ZIV_DECL (loop);\n"):
  fprintf(fd, "  MPFR_SAVE_EXPO_DECL (expo);\n"):
  fprintf(fd, "  mpfr_prec_t wprec;             /* working precision */\n"):
  fprintf(fd, "  mpfr_prec_t prec;              /* target precision */\n"):
  fprintf(fd, "  mpfr_prec_t err;               /* used to estimate the evaluation error */\n"):
  fprintf(fd, "  mpfr_prec_t correctBits;       /* estimates the number of correct bits*/\n"):
  fprintf(fd, "  unsigned long int k;\n"):
  fprintf(fd, "  unsigned long int conditionNumber;        /* condition number of the series */\n"):
  fprintf(fd, "  unsigned assumed_exponent;     /* used as a lowerbound of -EXP(f(x)) */\n"):
  fprintf(fd, "  int r;                         /* returned ternary value */\n"):
  fprintf(fd, "  mpfr_t s;                      /* used to store the partial sum */\n"):

  if (whattype(c) = 'fraction') or (whattype(c) = 'integer')
    then hardconstant := 0
    else hardconstant := 1
  fi:
  if (type=FUNCTION_SERIES) then hardconstant := 1 fi:

  if (hardconstant=1)
  then
    if (type=CONSTANT_SERIES)
      then fprintf(fd, "  mpfr_t x%d;                     /* used to store %a */\n", d, c):
    elif (type=FUNCTION_SERIES_RATIONAL)
      then fprintf(fd, "  mpfr_t x%d;                     /* used to store %a */\n", d, c*(u/v)^d):
    elif (type=FUNCTION_SERIES)
      then fprintf(fd, "  mpfr_t x%d;                     /* used to store %a */\n", d, c*x^d):
    fi:
  fi:

  if (type=FUNCTION_SERIES)
  then
    fprintf(fd, "  mpfr_t tmp;\n"):
    if (nops(exponents)-1 >= 1) then fprintf(fd, "  mpfr_t ") fi:
    for i from 1 to nops(exponents)-1 do
      fprintf(fd, "x%d", exponents[i]):
      if (i<nops(exponents)-1) then fprintf(fd, ", ") else fprintf(fd, ";                     /* used to store x^i */\n") fi
    od:
  fi:

  fprintf(fd, "  mpfr_t ");
  for i from 1 to nc do
    fprintf(fd, "tip%d", init_cond[i][1]):
    if (i<nc) then fprintf(fd, ", ") else fprintf(fd, "; /* used to store successive values of t_i */\n") fi:
  od:
  fprintf(fd, "  int "):
  for i from 1 to nc do
    fprintf(fd, "test%d", init_cond[i][1]):
    if (i<nc) then fprintf(fd, ", ") else fprintf(fd, ";\n") fi:
  od:
  fprintf(fd, "  int global_test;               /* used to test when the sum can be stopped */"):

  fprintf(fd, "\n"):
  fprintf(fd, "  /* Logging */\n"):
  if (type=FUNCTION_SERIES)
    then fprintf(fd, "  MPFR_LOG_FUNC ( (\"x[%%#R]=%%R rnd=%%d\", x, x, rnd), (\"res[%%#R]=%%R\", res, res) );\n\n")
  elif (type=FUNCTION_SERIES_RATIONAL)
    then fprintf(fd, "  MPFR_LOG_FUNC ( (\"x=u/v with u=%%d and v=%%d, rnd=%%d\", u, v, rnd), (\"res[%%#R]=%%R\", res, res) );\n\n")
  else fprintf(fd, "  MPFR_LOG_FUNC ( (\"rnd=%%d\", rnd), (\"res[%%#R]=%%R\", res, res) );\n\n")
  fi:

  ######################################################
  #################### Special cases ###################
  ######################################################

  if ( (type=FUNCTION_SERIES) or (type=FUNCTION_SERIES_RATIONAL) )
  then
    fprintf(fd, "  /* Special cases */\n"):

    if (type=FUNCTION_SERIES)
      then fprintf(fd, "  if (MPFR_UNLIKELY (MPFR_IS_NAN (x)))\n")
      else fprintf(fd, "  if (MPFR_UNLIKELY (v==0))\n")
    fi:
    fprintf(fd, "    {\n"):
    fprintf(fd, "      MPFR_SET_NAN (res);\n"):
    fprintf(fd, "      MPFR_RET_NAN;\n"):
    fprintf(fd, "    }\n"):

    if (init_cond[1][1] > 0) then f0 := 0 else f0 := init_cond[1][2] fi:
    if (type=FUNCTION_SERIES)
      then fprintf(fd, "  if (MPFR_UNLIKELY (MPFR_IS_ZERO (x)))\n")
      else fprintf(fd, "  if (MPFR_UNLIKELY (u==0))\n")
    fi:
    if (whattype(f0) = 'integer')
    then
      fprintf(fd, "    {\n"):
      fprintf(fd, "      return mpfr_set_si (res, %a, rnd);\n", f0):
      fprintf(fd, "    }\n")
    else
      printf("You need to provide a function mpfr_%a0(mpfr_t, mpfr_rnd_t) that evaluates %a with correct rounding\n", name, f0):
      fprintf(fd, "    {\n"):
      fprintf(fd, "      return mpfr_%a0 (res, rnd);\n", name):
      fprintf(fd, "    }\n")
    fi:

    if (type=FUNCTION_SERIES)
    then
      fprintf(fd, "  if (MPFR_UNLIKELY (MPFR_IS_INF (x)))\n"):
      fprintf(fd, "    {\n"):
      for i from -1 to 1 by 2 do  # Trick to handle both -oo and +oo
        if (i<0)
          then fprintf(fd, "      if (MPFR_IS_NEG (x))\n"):
          else fprintf(fd, "      else\n")
        fi:
        fprintf(fd, "        {\n"):
        if (i<0) then f0 := "m" else f0 := "p" fi:
        if (fofx <> _f(x))
        then
          f0 := limit(fofx, x=i*infinity):
          if (whattype(f0)='integer')
          then fprintf(fd, "          return mpfr_set_si (res, %a, rnd);\n", f0):
          elif (f0 = infinity) or (f0 = -infinity)
          then
            fprintf(fd, "          MPFR_SET_INF(res);\n"):
            if (f0>0)
            then fprintf(fd, "          MPFR_SET_POS(res);\n")
            else fprintf(fd, "          MPFR_SET_NEG(res);\n")
            fi:
            fprintf(fd, "          MPFR_RET(0);\n"):
          else
            if (i<0) then f0 := "m" else f0 := "p" fi:
          fi:
        fi:
        if ((f0 = "p") or (f0 = "m"))
        then
          printf("You need to provide a function mpfr_%a%sinf(mpfr_t, mpfr_rnd_t) that evaluates lim(f(x), x=", name, f0):
          if (f0 = "m") then printf("-") fi:
          printf("inf) with correct rounding.\n"):
          fprintf(fd, "          return mpfr_%a%sinf (res, rnd);\n", name, f0)
        fi:
        fprintf(fd, "        }\n")
      od:
      fprintf(fd, "    }\n")
    fi:
  fi:

  ######################################################
  ################## Precomputations ###################
  ######################################################

  fprintf(fd, "\n"):
  fprintf(fd, "  /* Save current exponents range */\n"):
  fprintf(fd, "  MPFR_SAVE_EXPO_MARK (expo);\n\n"):

  if ( (type=FUNCTION_SERIES) or (type=FUNCTION_SERIES_RATIONAL) )
  then  fprintf(fd, "  /* FIXME: special case for large values of |x| ? */\n\n")
  fi:

  # Note : prec is the value such that we will try to compute an approximation
  # with relative error smaller than 2^(1-prec).
  # Several things may happen:
  #   1) We do not achieve the intended error: this is because we badly estimated the exponent of the result
  #   2) We achieve the error, but it is not sufficient to decide correct rounding (Ziv's bad case)
  # Against 1), we can try to make our estimate of the exponent better with any heuristic.
  # Against 2), we can consider more guard bits. 11 guard bits seem a good value for the beginning
  #   (statistically, we expect to fail in less than 0.1 % of the cases)
  # wprec is the precision used during the computation, in order to ensure the final relative error 2^(1-prec)
  #
  fprintf(fd, "  /* We begin with 11 guard bits */\n"):
  fprintf(fd, "  prec = MPFR_PREC (res) + 11;\n"):
  fprintf(fd, "  MPFR_ZIV_INIT (loop, prec);\n"):

  # TODO: the value here can be chosen completely heursitically. We could do something much better
  # when fofx is known by using asympt(fofx, x, 1). A clean implementation appears to be complex though.
  # We must catch errors if the development does not exist (e.g. AiryAi(-x));
  # We must find a separation after which the asymptotic behavior is valid (e.g. x>1)
  printf("The code contains a variable assumed_exponent arbitrarily set to 10. You can put any value heuristically chosen. The closer it is to -log_2(|f(x)|), the better it is.\n"):
  fprintf(fd, "  assumed_exponent = 10; /* TIP: You can put any value heuristically chosen. The closer it is to -log_2(|f(x)|), the better it is */\n"):

  # TODO: find a way of putting a rigorous value here.
  # This value *must* be rigorous: the safety of the implementation relies on it.
  # Precisely, we need to have sum(|a(i)*x^i|) <= 2^conditionNumber
  fprintf(fd, "  conditionNumber = xxx; /* FIXME: set a value such that sum(|a(i)*x^i|) <= 2^conditionNumber */\n"):
  printf("The code contains a variable conditionNumber that you must manually set to a suitable value, in order to ensure that sum_{i=0}^{infinity} |a(i)*x^i| <= 2^conditionNumber\n"):
  fprintf(fd, "  wprec = prec + ERRORANALYSISPREC + conditionNumber + assumed_exponent;\n"):

  ######################################################
  ################## Initialisations ###################
  ######################################################

  if (hardconstant=1)
  then
    fprintf(fd, "  mpfr_init (x%d);\n", d):
  fi:
  if (type = FUNCTION_SERIES)
  then
    fprintf(fd, "  mpfr_init (tmp);\n"):
    for i from 1 to nops(exponents)-1 do
      fprintf(fd, "  mpfr_init (x%d);\n", exponents[i])
    od:
  fi:

  for i from 1 to nc do
    fprintf(fd, "  mpfr_init (tip%d);\n", init_cond[i][1])
  od:
  fprintf(fd, "  mpfr_init (s);\n\n"):

  ######################################################
  ########## Ziv' loop: setting the precision ##########
  ######################################################

  fprintf(fd, "  /* ZIV' loop */\n"):
  fprintf(fd, "  for (;;)\n"):
  fprintf(fd, "    {\n"):

  fprintf(fd, "      MPFR_LOG_MSG ((\"Working precision: %%d\\n\", wprec, 0));\n\n"):
  if (hardconstant=1)
  then
    fprintf(fd, "      mpfr_set_prec (x%d, wprec);\n", d):
  fi:
  if (type = FUNCTION_SERIES)
  then
    fprintf(fd, "      mpfr_set_prec (tmp, wprec);\n"):
    fprintf(fd, "      if(mpfr_get_prec (x) > wprec)\n"):
    fprintf(fd, "        mpfr_set_prec (x1, wprec);\n"):
    fprintf(fd, "      else\n"):
    fprintf(fd, "        mpfr_set_prec (x1, mpfr_get_prec (x));\n"):
    for i from 2 to nops(exponents)-1 do
      fprintf(fd, "      mpfr_set_prec (x%d, wprec);\n", exponents[i])
    od:
  fi:

  for i from 1 to nc do
    fprintf(fd, "      mpfr_set_prec (tip%d, wprec);\n", init_cond[i][1])
  od:
  fprintf(fd, "      mpfr_set_prec (s, wprec);\n\n"):

  ######################################################
  ############ Ziv' loop: initial conditions ###########
  ######################################################

  fprintf(fd, "      mpfr_set_ui (s, 0, MPFR_RNDN);\n"):

  if (type = FUNCTION_SERIES)
  then
    fprintf(fd, "      mpfr_set (x1, x, MPFR_RNDN);\n"):
    error_counter := init_error_counter("x1", error_counter):
    for i from 2 to nops(exponents) do
        fprintf(fd, "      mpfr_mul (x%d, x%d, x%d, MPFR_RNDN);\n", exponents[i], exponents[i-1], exponents[i]-exponents[i-1]):
	var1 := sprintf("x%d", exponents[i]):
        var2 := sprintf("x%d", exponents[i-1]):
        var3 := sprintf("x%d", exponents[i]-exponents[i-1]):
	error_counter := error_counter_of_a_multiplication(var1, var2, var3, error_counter):
    od:
  fi:

  for i from 1 to nc do
    i0 := init_cond[i][1]:
    ri0 := init_cond[i][2]:  # We implement t_{i0} <- ri0
    if (whattype(ri0)='integer') or (whattype(ri0)='fraction')
    then
      if (type = FUNCTION_SERIES) and (i0 <> 0)
        then
          var :=  sprintf("      mpfr_mul_si (tip%d, x%d, ", i0, i0):
          var1 := sprintf("tip%d", i0):
          var2 := sprintf("x%d", i0):
          error_counter := error_counter_of_a_multiplication(var1, var2, "", error_counter):
        else
          var := sprintf("      mpfr_set_si (tip%d, ", i0):
          var1 := sprintf("tip%d", i0):
          error_counter := init_error_counter(var1, error_counter):
      fi:
      fprintf(fd, "%s%d, MPFR_RNDN);\n", var, numer(ri0)):
      if (whattype(ri0)='fraction')
      then
        fprintf(fd, "      mpfr_div_si (tip%d, tip%d, %d, MPFR_RNDN);\n", i0, i0, denom(ri0)):
        var1 := sprintf("tip%d", i0):
        error_counter := error_counter_of_a_multiplication(var1, var1, "", error_counter):
      fi:
    else
      printf("You need to provide a function mpfr_%a%d (mpfr_t, mpfr_rnd_t) that evaluates %a with faithful rounding.\n", name, i0, ri0):
      fprintf(fd, "      mpfr_%a%d (tip%d, MPFR_RNDN);\n", name, i0, i0):
      var1 := sprintf("tip%d", i0):
      error_counter := init_error_counter(var1, error_counter):
      if (type = FUNCTION_SERIES) and (i0 <> 0)
      then
        fprintf(fd, "      mpfr_mul (tip%d, tip%d, x%d, MPFR_RNDN);\n", i0, i0, i0):
        var1 := sprintf("tip%d", i0):
        var2 := sprintf("x%d", i0):
        error_counter := error_counter_of_a_multiplication(var1, var1, var2, error_counter):
      fi:
    fi:

    if (type = FUNCTION_SERIES_RATIONAL) and (i0 <> 0)
    then
      var1 := sprintf("tip%d", i0):
      if (i0 = 1)
      then fprintf(fd, "      mpfr_mul_si (tip%d, tip%d, ", i0, i0):
      else
        required_mulsi := { op(required_mulsi), i0 }:
        fprintf(fd, "      mpfr_mul_si%d (tip%d, tip%d, ", i0, i0, i0):
      fi:
      for j from 1 to i0 do
        fprintf(fd, "u, "):
        error_counter := error_counter_of_a_multiplication(var1, var1, "", error_counter):
      od:
      fprintf(fd, "MPFR_RNDN);\n"):
      if (i0 = 1)
      then fprintf(fd, "      mpfr_div_si (tip%d, tip%d, ", i0, i0):
      else
        required_divsi := { op(required_divsi), i0 }:
        fprintf(fd, "      mpfr_div_si%d (tip%d, tip%d, ", i0, i0, i0):
      fi:
      for j from 1 to i0 do
        fprintf(fd, "v, "):
        error_counter := error_counter_of_a_multiplication(var1, var1, "", error_counter):
      od:
      fprintf(fd, "MPFR_RNDN);\n"):
    fi:

    fprintf(fd, "      mpfr_add (s, s, tip%d, MPFR_RNDN);\n\n", i0):
  od:

  if (whattype(c) = 'integer') or (whattype(c) = 'fraction')
  then
    if (type = FUNCTION_SERIES) and (c = -1)
    then
      fprintf(fd, "      MPFR_CHANGE_SIGN (x%d);\n", d):
    elif (type = FUNCTION_SERIES) and (c <> 1)
    then
      fprintf(fd, "      mpfr_mul_si (x%d, x%d, %d, MPFR_RNDN);\n", d, d, numer(c)):
      var1 := sprintf("x%d", d):
      error_counter := error_counter_of_a_multiplication(var1, var1, "", error_counter):
      if (whattype(c) = 'fraction')
      then
        fprintf(fd, "      mpfr_div_si (x%d, x%d, %d, MPFR_RNDN);\n", d, d, denom(c)):
        error_counter := error_counter_of_a_multiplication(var1, var1, "", error_counter):
      fi:
    fi:
  else
    printf("You need to provide a function mpfr_%a_cste (mpfr_t, mpfr_rnd_t) that evaluates %a with faithful rounding.\n", name, c):
    if (type = CONSTANT_SERIES)
    then
      fprintf(fd, "      mpfr_%a_cste (x%d, MPFR_RNDN);\n", name, d):
      var1 := sprintf("x%d", d):
      error_counter := init_error_counter(var1, error_counter):
    elif (type = FUNCTION_SERIES) then
      fprintf(fd, "      mpfr_%a_cste (tmp, MPFR_RNDN);\n", name):
      error_counter := init_error_counter("tmp", error_counter):
      fprintf(fd, "      mpfr_mul (x%d, tmp, x%d, MPFR_RNDN);\n", d, d):
      var1 := sprintf("x%d", d):
      error_counter := error_counter_of_a_multiplication(var1, "tmp", var1, error_counter):
    elif (type = FUNCTION_SERIES_RATIONAL) then
      fprintf(fd, "      mpfr_%a_cste (x%d, MPFR_RNDN);\n", name, d):
      var1 := sprintf("x%d", d):
      error_counter := init_error_counter(var1, error_counter):
      if (d = 1)
      then fprintf(fd, "      mpfr_mul_si (x%d, x%d, ", d, d):
      else
        required_mulsi := { op(required_mulsi), d }:
        fprintf(fd, "      mpfr_mul_si%d (x%d, x%d, ", d, d, d):
      fi:
      for j from 1 to d do
        fprintf(fd, "u, "):
        error_counter := error_counter_of_a_multiplication(var1, var1, "", error_counter):
      od:
      fprintf(fd, "MPFR_RNDN);\n"):
      var1 := sprintf("x%d", d):
      if (d = 1)
      then fprintf(fd, "      mpfr_div_si (x%d, x%d, ", d, d):
      else
        required_divsi := { op(required_divsi), d }:
        fprintf(fd, "      mpfr_div_si%d (x%d, x%d, ", d, d, d):
      fi:
      for j from 1 to d do
        fprintf(fd, "v, "):
        error_counter := error_counter_of_a_multiplication(var1, var1, "", error_counter):
      od:
      fprintf(fd, "MPFR_RNDN);\n"):
    fi:
  fi:

  ######################################################
  ######### Ziv' loop: evaluation of the series ########
  ######################################################

  fprintf(fd, "\n"):
  fprintf(fd, "      /* Evaluation of the series */\n"):
  fprintf(fd, "      k = %d;\n", d):
  fprintf(fd, "      for (;;)\n"):
  fprintf(fd, "        {\n"):
  if (init_cond[1][1] <> 0) then fprintf(fd, "          k += %d;\n", init_cond[1][1]) fi:

  for i from 1 to nc do
    error_in_loop := 0:
    i0 := init_cond[i][1]:
    if (hardconstant = 1)
    then
      fprintf(fd, "          mpfr_mul (tip%d, tip%d, x%d, MPFR_RNDN);\n", i0, i0, d):
      var1 := sprintf("x%d", d):
      error_in_loop := error_in_loop + 1 + find_in_error_counter(var1, error_counter):
    else
      var := sprintf("tip%d", i0):
      temp := generate_multiply_rational(fd, var, var, c, [[var, error_in_loop]], "          "):
      error_in_loop := find_in_error_counter(var, temp):
      if (type = FUNCTION_SERIES_RATIONAL)
      then
        if (d=1)
        then
          fprintf(fd, "          mpfr_mul_si (tip%d, tip%d, u, MPFR_RNDN);\n", i0, i0):
          fprintf(fd, "          mpfr_div_si (tip%d, tip%d, v, MPFR_RNDN);\n", i0, i0):
	  error_in_loop := error_in_loop + 2:
        else
          required_mulsi := { op(required_mulsi), d }:
          fprintf(fd, "          mpfr_mul_si%d (tip%d, tip%d", d, i0, i0):
          for j from 1 to d do fprintf(fd, ", u") od:
          fprintf(fd, ", MPFR_RNDN);\n"):
	  error_in_loop := error_in_loop + d:

          required_divsi := { op(required_divsi), d }:
          fprintf(fd, "          mpfr_div_si%d (tip%d, tip%d", d, i0, i0):
          for j from 1 to d do fprintf(fd, ", v") od:
          fprintf(fd, ", MPFR_RNDN);\n" ):
	  error_in_loop := error_in_loop + d:
        fi
      fi
    fi:
    var := sprintf("tip%d", i0):

    temp := generate_multiply_poly(fd, var, var, subs(n=k-d, p/q), [[var, error_in_loop]], "          "):
    required_mulsi := { op(required_mulsi), op(temp[1]) }:
    required_divsi := { op(required_divsi), op(temp[2]) }:
    error_in_loop := find_in_error_counter(var, temp[3]):
    temp := generate_multiply_poly(fd, "tmp", var, subs(n=k, s2/s1), [[var, error_in_loop]], "          "):
    required_mulsi := { op(required_mulsi), op(temp[1]) }:
    required_divsi := { op(required_divsi), op(temp[2]) }:

    fprintf(fd, "          mpfr_add (s, s, tmp, MPFR_RNDN);\n"):

    if (i<nc) then fprintf(fd, "\n          k += %d;\n", init_cond[i+1][1]-i0)
    else  fprintf(fd, "\n          k += %d;\n", d-i0)
    fi:
  od:

  ######################################################
  ################### Error analysis ###################
  ######################################################

  maxofti := 0: # store the maximum of the error counters of the initial conditions
  for i from 1 to nc do
    var := sprintf("tip%d", init_cond[i][1]):
    if find_in_error_counter(var, error_counter) > maxofti
      then maxofti := find_in_error_counter(var, error_counter):
    fi:
  od:
  additional_error := find_in_error_counter("tmp", temp[3]) - error_in_loop:

  ######################################################
  #### Ziv' loop: stopping criterion for the series ####
  ######################################################

  # The first neglected term is tk, so the remainder is made by
  # tk + t(k+d) + t(k+2d).... and the corresponding series
  # beginning with t(k+1), t(k+2), etc. up to t(k+d-1).
  #
  # We have t(k0+d) = c*s1(k0)/s1(k0+d) * s2(k0+d)/s2(k0)* p(k0)/q(k0) * x^d t(k0)
  # (where x=u/v or x=1 in cases of rational series or constant series)
  # So it suffices that:
  # forall k0>=k-d, |c * s1(k0)/s1(k0+d) * s2(k0+d)/s2(k0) * p(k0)/q(k0) * x^d| <= 1/2     (1)
  #
  # If this is true, |tk| = |c*s1(k-d)/s1(k) * s2(k)/s2(k-d)* p(k-d)/q(k-d) * x^d t(k-d)| <= t(k-d)/2
  # This is also true for larger values of k, so we can bound |tk + t(k+d) + t(k+2d) + ...| by |t(k-d)|.
  # And the same holds for |t(k+1) + ...|, |t(k+2) + ...|, etc. up to |t(k+d-1)+...|.
  #
  # global_test depends on k and we must satisfy:
  #  "if (global_test) then (1) holds".
  #
  # the total remainder is bounded by 2*nc*tk.

  fprintf(fd, "          global_test = xxx; /* FIXME: set the value in order to ensure that, whenever global_test is true, we have: forall k'>=k, |r(k')*x^d| <= 1/2, where r is the fraction such that a(n)=r(n)a(n-d)*/\n"):
  printf("The code contains a variable global_test that you must manually set to a suitable value, in order to ensure that when global_test is true, the following holds:\n"):
  printf("        forall k'>=k, |r(k')*x^d| <= 1/2, where r is the fraction such that a(n)=r(n)a(n-d)\n"):
  guard_bits := 1+1+ceil(log[2](nc)):
  for i from 1 to nc do
    i0 := init_cond[i][1]:
    fprintf(fd, "          test%d = ( (!MPFR_IS_ZERO(s))\n", i0):
    fprintf(fd, "                    && ( MPFR_IS_ZERO(tip%d)\n", i0):
    fprintf(fd, "                         || (MPFR_EXP(tip%d) + (mp_exp_t)prec + %d <= MPFR_EXP(s))\n", i0, guard_bits):
    fprintf(fd, "                         )\n"):
    fprintf(fd, "                    );\n"):
  od:
  fprintf(fd, "          if (");
  for i from 1 to nc do
    fprintf(fd, "test%d && ", init_cond[i][1]):
  od:
  fprintf(fd, "global_test)\n"):
  fprintf(fd, "            break;\n"):
  fprintf(fd, "        }\n\n"):

  ######################################################
  ############## Ziv' loop: testing final ##############
  ######################################################

  fprintf(fd, "      MPFR_LOG_MSG ((\"Truncation rank: %%lu\\n\", k));\n\n"):
  fprintf(fd, "      err = ERRORANALYSISK + conditionNumber - MPFR_GET_EXP (s);\n\n"):
  fprintf(fd, "      /* err is the number of bits lost due to the evaluation error */\n"):
  fprintf(fd, "      /* wprec-(prec+1): number of bits lost due to the approximation error */\n"):
  fprintf(fd, "      MPFR_LOG_MSG ((\"Roundoff error: %%Pu\\n\", err));\n"):
  fprintf(fd, "      MPFR_LOG_MSG ((\"Approxim error: %%Pu\\n\", wprec-prec-1));\n\n"):
  fprintf(fd, "      if (wprec < err+1)\n"):
  fprintf(fd, "        correct_bits=0;\n"):
  fprintf(fd, "      else\n"):
  fprintf(fd, "        {\n"):
  fprintf(fd, "          if (wprec < err+prec+1)\n"):
  fprintf(fd, "            correct_bits =  wprec - err - 1;\n"):
  fprintf(fd, "          else\n"):
  fprintf(fd, "            correct_bits = prec;\n"):
  fprintf(fd, "        }\n\n"):
  fprintf(fd, "      if (MPFR_LIKELY (MPFR_CAN_ROUND (s, correct_bits, MPFR_PREC (y), rnd)))\n"):
  fprintf(fd, "        break;\n\n"):

  fprintf(fd, "      if (correct_bits == 0)\n"):
  fprintf(fd, "        {\n"):
  fprintf(fd, "          assumed_exponent *= 2;\n"):
  fprintf(fd, "          MPFR_LOG_MSG ((\"Not a single bit correct (assumed_exponent=%%lu)\\n\",\n"):
  fprintf(fd, "                         assumed_exponent));\n"):
  fprintf(fd, "          wprec = prec + ERRORANALYSISK + conditionNumber + assumed_exponent;\n"):
  fprintf(fd, "        }\n"):
  fprintf(fd, "      else\n"):
  fprintf(fd, "        {\n"):
  fprintf(fd, "          if (correct_bits < prec)\n"):
  fprintf(fd, "            { /* The precision was badly chosen */\n"):
  fprintf(fd, "              MPFR_LOG_MSG ((\"Bad assumption on the exponent of %s(x)\", 0));\n", name):
  fprintf(fd, "              MPFR_LOG_MSG ((\" (E=%%ld)\\n\", (long) MPFR_GET_EXP (s)));\n"):
  fprintf(fd, "              wprec = prec + err + 1;\n"):
  fprintf(fd, "            }\n"):
  fprintf(fd, "          else\n"):
  fprintf(fd, "            { /* We are really in a bad case of the TMD */\n"):
  fprintf(fd, "              MPFR_ZIV_NEXT (loop, prec);\n\n"):

  fprintf(fd, "              /* We update wprec */\n"):
  fprintf(fd, "              /* We assume that K will not be multiplied by more than 4 */\n"):
  fprintf(fd, "              wprec = prec + ERRORANALYSIS4K + conditionNumber\n"):
  fprintf(fd, "                - MPFR_GET_EXP (s);\n"):
  fprintf(fd, "            }\n"):
  fprintf(fd, "        }\n\n"):

  fprintf(fd, "    } /* End of ZIV loop */\n\n"):
  fprintf(fd, "  MPFR_ZIV_FREE (loop);\n\n"):
  fprintf(fd, "  r = mpfr_set (res, s, rnd);\n\n"):

  ######################################################
  ################ Clearing everything #################
  ######################################################

  fprintf(fd, "  mpfr_clear (s);\n"):
  if (hardconstant=1)
  then
    fprintf(fd, "  mpfr_clear (x%d);\n", d):
  fi:
  if (type = FUNCTION_SERIES)
  then
    fprintf(fd, "  mpfr_clear (tmp);\n"):
    for i from 1 to nops(exponents)-1 do
      fprintf(fd, "  mpfr_clear (x%d);\n", exponents[i]):
    od:
  fi:

  for i from 1 to nc do
    fprintf(fd, "  mpfr_clear (tip%d);\n", init_cond[i][1]):
  od:

  fprintf(fd, "\n"):
  fprintf(fd, "  MPFR_SAVE_EXPO_FREE (expo);\n"):
  fprintf(fd, "  return mpfr_check_range (res, r, rnd);\n"):
  fprintf(fd, "}\n"):

  fclose(fd):

  for i from 1 to nops(required_mulsi) do
    printf("You need to provide a mpfr_mul_si%d function.\n", required_mulsi[i]):
    printf("  -> This can be achieved by a call to generate_muldivsin(\"mul\", %d):\n",  required_mulsi[i]):
  od:
  for i from 1 to nops(required_divsi) do
    printf("You need to provide a mpfr_div_si%d function.\n", required_divsi[i]):
    printf("  -> This can be achieved by a call to generate_muldivsin(\"div\", %d):\n",  required_divsi[i]):
  od:

  ######################################################
  ################## Error analysis ####################
  ######################################################

  printf("\n\n"):
  printf("Before the loop, we have "):
  for i from 1 to nc do
    var := sprintf("tip%d", init_cond[i][1]):
    printf("%s {%d}", var, find_in_error_counter(var, error_counter)):
    if (i <> nc) then printf(", ") else printf("\n") fi:
  od:
  printf("Each step of the loop adds another {%d}\n", error_in_loop):
  if (additional_error <> 0)
  then printf("Moreover, the multiplication by %a adds another {%d} to each term before it is summed.\n", subs(n=k, s2/s1), additional_error)
  fi:
  printf("Finally, we have s = sum_(i=0)^(k-1) ( ti{%d + %dk} )\n", maxofti + additional_error + 1 - error_in_loop, error_in_loop):
  printf("We bound it by {(k+%d)*2^(%d)}\n", ceil( (maxofti + additional_error + 1 - error_in_loop)/error_in_loop), ceil(log[2](error_in_loop))):

  a := ceil( (maxofti + additional_error + 1 - error_in_loop)/error_in_loop):
  b := ceil(log[2](error_in_loop)):

  if (a > 0)
    then var := sprintf("MPFR_INT_CEIL_LOG2 (prec + %d)", a)
    elif (a=0) then  var := sprintf("MPFR_INT_CEIL_LOG2 (prec)")
    else sprintf("MPFR_INT_CEIL_LOG2 (prec - %d)", -a)
  fi:
  if (b > 0) then var := sprintf("%s + %d", var, b+2) fi:
  var := sprintf("sed -n -i 's/ERRORANALYSISPREC/%s/g;p' %s", var, filename):
  system(var):

  if (a > 0)
    then var := sprintf("MPFR_INT_CEIL_LOG2 (k + %d)", a)
    elif (a=0) then  var := sprintf("MPFR_INT_CEIL_LOG2 (k)")
    else sprintf("MPFR_INT_CEIL_LOG2 (k - %d)", -a)
  fi:
  if (b > 0) then var := sprintf("%s + %d", var, b+2) fi:
  var := sprintf("sed -n -i 's/ERRORANALYSISK/%s/g;p' %s", var, filename):
  system(var):

  if (a > 0)
    then var := sprintf("MPFR_INT_CEIL_LOG2 (k + %d)", a)
    elif (a=0) then  var := sprintf("MPFR_INT_CEIL_LOG2 (k)")
    else sprintf("MPFR_INT_CEIL_LOG2 (k - %d)", -a)
  fi:
  var := sprintf("%s + %d", var, b+4):
  var := sprintf("sed -n -i 's/ERRORANALYSIS4K/%s/g;p' %s", var, filename):
  system(var):

end proc:


text:
with(DifferentialAlgebra):
ring_diff := DifferentialRing(blocks = [[s, x2, i, r, x1], [y1, y2] ], derivations = [t]):
sys := [
diff(s(t), t) - (-b1*b0*s(t)*x1(t)*i(t) - b0*s(t)*i(t) - s(t)*mu + mu + g*r(t)),
diff(x2(t), t) - (M*x1(t)),
diff(i(t), t) - (-nu*i(t) + b1*b0*s(t)*x1(t)*i(t) + b0*s(t)*i(t) - mu*i(t)),
diff(r(t), t) - (nu*i(t) - mu*r(t) - g*r(t)),
diff(x1(t), t) - (-M*x2(t)),
y1(t) - (i(t)),
y2(t) - (r(t))
];
res := CodeTools[CPUTime](RosenfeldGroebner(sys, ring_diff, singsol=none));

ext: mpl


text:
with(DifferentialAlgebra):
ring_diff := DifferentialRing(blocks = [[x1, x2], [y1] ], derivations = [t]):
sys := [
diff(x1(t), t) - (a*x1(t) + b*x1(t) - x2(t)*c*x1(t)),
diff(x2(t), t) - (-a*b*x2(t) + d*x2(t)*x1(t)),
y1(t) - (x1(t))
];
res := CodeTools[CPUTime](RosenfeldGroebner(sys, ring_diff, singsol=none));


text:
//impulse: /give @p spawn_egg 1 0 {display:{Name:"Redirecting Cube Spawner"},EntityTag:{id:"minecraft:chicken",CustomName:"MMH_RedirectingCubeSpawner",Silent:1,NoGravity:1}}
//MMH_loadChunks()
///clone ~1 ~1 ~1 ~3 ~4 ~3 20 1 20

#MMH
repeat process MMH_discouragementRedirectionCubeSpawner {
if: /testfor @e[type=Chicken,name=MMH_RedirectingCubeSpawner]
then {
  MMH_loadChunks()
  /execute @e[type=Chicken,name=MMH_RedirectingCubeSpawner] ~ ~ ~ clone 20 1 20 22 4 22 ~-1 ~ ~-1 masked
  /execute @e[type=Chicken,name=MMH_RedirectingCubeSpawner] ~ ~ ~ fill ~ ~ ~ ~ ~ ~ air 0 replace pumpkin
  /tp @e[type=Chicken,name=MMH_RedirectingCubeSpawner] ~ -100 ~
}
}

//repeat: /testfor @e[type=armor_stand,name=ACV_Cube,tag=1]
///blockdata ${-1} {SuccessCount:1}
//conditional: /summon armor_stand ~1 ~-2 ~-1 {CustomName:"ACV_Cube",DisabledSlots:2035728,ArmorItems:[{},{},{},{Count:1,id:"dispenser"}],NoGravity:1b,Invisible:1b,Invulnerable:1b,Tags:["ACV_Redirecting","1"]}


text:
single_case_Partition := module()
  uses Domain, Domain_Type;

  export SimplName  := "Single case partition";
  export SimplOrder := 11;

    export ModuleApply := proc(vs :: DomBound, sh :: DomShape, $)
        subsindets(sh, DomSplit, x->do_simp(op(x)));
    end proc;

    local do_simp := proc(p:: Partition,$)::DomShape;
        local r := Partition:-Simpl:-single_nonzero_piece_cps(
            proc(c,v) if v::DomConstrain then DConstrain(conv_bool(c),op(v)) else p end if
            end proc,p,_testzero=(x->x=DSum()));
        if r :: Partition then DSplit(r) else r end if;
    end proc;

    local conv_bool := proc(r, $)
        if r :: {specfunc(`And`), `and`} then
            op(map(conv_bool,r))
        else
            r
        end if;
    end proc;
end module;

Partition_simpl := module()
  uses Domain, Domain_Type;

  export SimplName  := "Partition simpl";
  export SimplOrder := (10+1/2);

  local TRY := proc(bnds, kb, kb_rn, as, pr)
    local r, ns;
    r := subs(kb_rn, op(1,pr));

    ns := op(1,bnds);
    ns := select(b->op(3,b) in {`Int`,`Ints`}, ns);
    ns := map(curry(op,1), ns);

    r := Partition:-Simpl(r, kb, _name_cands=ns, _testequal=`=`) assuming op(as);
    if not r :: Partition then r else DSplit(r) end if;
  end proc;

  export ModuleApply := proc(vs :: DomBound, sh :: DomShape, $)
    local as := Domain:-Bound:-toConstraints(vs, 'bound_types');
    subsindets(sh, DomSplit, curry(TRY,vs, op(Domain:-Bound:-toKB(vs)),as));
  end proc;
end module;

max_stars_repo_name: vmchale/hakaru

max_issues_repo_name: vmchale/hakaru

max_forks_repo_name: vmchale/hakaru


text:
program sample29;
var unused1 : integer;
    UnusedArrayForTest : array[200] of char;
procedure gcmlcm(m, n, gc, lc : integer);
var a, b, r : integer;
    begin
        a := m;
        b := n;
        while b <> 0 do begin
            r := a - (a div b) * b;
            a := b;
            b := r
        end;
        gc := a;
        lc := (m div gc) * n
    end;
procedure abs(a, b : integer);
        if a < 0 then
            b := -a
        else
            b := a
    end;
procedure gcm(a, b, gc : integer);
var lc, aa, bb : integer;
        if (a = 0) or (b = 0) then
            gc := 1
        else
                call abs(a, aa);
                call abs(b, bb);
                call gcmlcm(aa, bb, gc, lc)
            end
    end;
procedure lcm(a, b, lc : integer);
var gc, aa, bb : integer;
        if (a = 0) or (b = 0) then
            lc := 1
        else
                call abs(a, aa);
                call abs(b, bb);
                call gcmlcm(aa, bb, gc, lc)
            end
    end;
var unusedchar : char;
procedure reduce(a1, a2 : integer);
var gc : integer;
        if a1 = 0 then
                a2 := 1;
                return
            end;
        if a2 = 0 then
                a1 := 1;
                return
            end;
        if a2 < 0 then
                a1 := -a1;
                a2 := -a2
            end;
        call gcm(a1, a2, gc);
        a1 := a1 div gc;
        a2 := a2 div gc
    end;
procedure sum(x1, x2, y1, y2 : integer);
var lc, y11 : integer;
        call lcm(x2, y2, lc);
        x1 := x1 * (lc div x2);
        y11 := y1 * (lc div y2);
        x1 := x1 + y11;
        x2 := lc;
        call reduce(x1, x2)
    end;
procedure sub(x1, x2, y1, y2 : integer);
var lc, y11 : integer;
        call sum(x1, x2, -y1, y2)
    end;
procedure mult(x1, x2, y1, y2 : integer);
var gc, y22, y11 : integer;
        call gcm(x1, y2, gc);
        x1 := x1 div gc;
        y22 := y2 div gc;
        call gcm(x2, y1, gc);
        x2 := x2 div gc;
        y11 := y1 div gc;
        x1 := x1 * y11;
        x2 := x2 * y22;
        call reduce(x1, x2)
    end;
procedure divide(x1, x2, y1, y2 : integer);
        call mult(x1, x2, y2, y1)
    end;
var unusedarray : array[100] of char;
procedure printfinal(a, b : integer);
        if a = 0 then
            writeln('Final Result =', a)
        else if b = 1 then
            writeln('Final Result =', a)
        else
            writeln('Final Result =', a, '/', b)
    end;
procedure printtemp(a, b : integer);
        if a = 0 then
            writeln('Temporary Result =', a)
        else if b = 1 then
            writeln('Temporary Result =', a)
        else
            writeln('Temporary Result =', a, '/', b)
    end;
var x1, x2, y1, y2 : integer;
var com : char;
    endflag : boolean;
    writeln('   *** Calculator -- h for help ***');
    x1 := 0;
    x2 := 1;
    endflag := false;
    while notendflag do begin
        writeln(' Please input command :');
        readln(com,y1);
        y2 := 1;
        if (com = 'c') or (com = 'C') then
                x1 := y1;
                x2 := y2
            end
        else if com = '+' then
            call sum(x1, x2, y1, y2)
        else if com = '-' then
            call sub(x1, x2, y1, y2)
        else if com = '*' then
            call mult(x1, x2, y1, y2)
        else if com = '/' then
            call divide(x1, x2, y1, y2)
        else if (com = 'o') or (com = 'O') then
            endflag := true
        else
                writeln;
                writeln('Calculator Usage:');
                writeln('  c number : clear & set it');
                writeln('  + number : add it');
                writeln('  - number : subtract it');
                writeln('  * number : multiply it');
                writeln('  / number : divide by it');
                writeln('  o        : off(terminate execution)');
                writeln
            end;
        if endflag then
            call printfinal(x1, x2)
        else
            call printtemp(x1, x2)
    end
end.

max_stars_repo_name: Hiroya-W/lang-processing

max_issues_repo_name: Hiroya-W/lang-processing

max_forks_repo_name: Hiroya-W/lang-processing


text:
standard_colour_list :=
 map(u -> ColorTools[Color](u /~ 255.),[
  [255,  0,  0], [229,122,  0], [ 76,115,  0], [ 64,255, 89],
  [ 38,145,153], [  0, 41,153], [213,128,255], [230,115,161],
  [ 89,  0,  0], [ 89, 58, 22], [204,197,102], [ 51,102, 58],
  [  0,204,255], [115,130,230], [ 96, 57,115], [242,  0, 65],
  [191, 67, 48], [178,137, 89], [255,238,  0], [ 45,179, 98],
  [  0, 68,128], [ 36,  0, 89], [204,  0,163], [127, 64, 72],
  [255,162,128], [102, 82,  0], [195,255,128], [128,255,246],
  [ 77,117,153], [136,  0,204], [102,  0, 54]
 ]):

standard_colour := proc(i)
 local n;
 n := nops(standard_colour_list);
 standard_colour_list[modp(i-1,n)+1];
end:


text:
(* type: work_gga_x *)
(* prefix:
  gga_x_b88_params *params;

  assert(p->params != NULL);
  params = (gga_x_b88_params * )(p->params);
*)

f0 := x -> 1.0 + params_a_beta/X_FACTOR_C*x^2/(1.0 + params_a_gamma*params_a_beta*x*arcsinh(x)):

f := x -> f0(x):

max_stars_repo_name: rdietric/lsms

max_issues_repo_name: rdietric/lsms

max_forks_repo_name: rdietric/lsms


text:
input := FileTools:-Text:-ReadFile("AoC-2021-20-input.txt"):
with(StringTools):

split:=Split(Trim(input), "\n"):
code := split[1]; length(code);

bin2dec := proc(s)
     if length(s) <> 9 then error "length of %1 is not 9, it is %2", s, length(s); end if;
     local l := Explode(s);
     add(ifelse(l[i]=".", 0, ifelse(l[i]="#", 1, X))*2^(9-i), i=1..9);
end proc:

bin2dec("...#...#.") = 34; code[34+1];

image := split[3..-1]: {Explode(Join(image,""))[]};

iwidth := length(image[1]);
ilength := nops(image);
ipad := 2;
p := code[bin2dec(cat(code[1]$9))+1];

image := [ cat(p$(iwidth+2*ipad))$ipad,
    seq(cat(p$ipad,image[i],p$ipad), i=1..iwidth),
    cat(p$(iwidth+2*ipad))$ipad
    ]:

{map(length, image)[]};

viewimage := proc(image)
  for local line in image do
     printf(cat(line,"\n"));
  end do;
end proc:

viewimage(image);
# walk the image
nimage := image:

for gen to 50 do
    nwidth := length(nimage[1]);
    nlength := nops(nimage);
    nbuf := StringBuffer():
    for i from 2 to nlength-1 do
        for j from 2 to nwidth-1 do
            nbhd := cat( seq(nimage[i+k][j-1..j+1], k=-1..1));
            pxidx := bin2dec(nbhd);
            nbuf:-append(code[pxidx+1]);
        end do;
    end do;
    tmp :=  nbuf:-value();
    nimage := [ seq(tmp[(i-1)*(nwidth-2)+1..i*(nwidth-2)], i=1..nlength-2) ];
    nwidth := length(nimage[1]);
    nlength := nops(nimage);
    ipad := 2;
    p := ifelse(gen mod 2 = 1, code[1], code[bin2dec(cat(code[1]$9))+1]);
    nimage := [ cat(p$(nwidth+2*ipad))$ipad, seq(cat(p$ipad,nimage[i],p$ipad), i=1..nwidth), cat(p$(nwidth+2*ipad))$ipad ]:
end do:

viewimage(nimage);

answer2 := add(CountCharacterOccurrences(nimage[i], "#"), i=1..nops(nimage));


text: aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa


text:
`find_bases/homology` := proc(H)
 local dM,dZM,nC,nZ,nB,nH,d,d_max,i,j,r,LM,ZM,BZM,BM,HZM,HM,
  S,M0,M1,M2,M3,pivots,non_pivots;

 dM := H[differential_matrix]:
 nC := eval(H['chain_number']);
 d_max := H['dim'];
 for d from 0 to d_max do
  if d = 0 then
   LM[0] := Matrix(nC[0],0):
   ZM[0] := IdentityMatrix(nC[0]):
   nZ[0] := nC[0];
  else
   ZM[d] := [op(NullSpace(dM[d]))];
   nZ[d] := nops(ZM[d]);
   ZM[d] := Matrix(nC[d],nZ[d],ZM[d]);
  fi;

  if d = d_max then
   BZM[d] := Matrix(nZ[d],0);
   BM[d] := Matrix(nC[d],0);
   nB[d] := 0;
   HZM[d] := IdentityMatrix(nZ[d]);
   HM[d] := ZM[d];
   nH[d] := nZ[d];
  else
   dZM[d+1] := LinearSolve(ZM[d],dM[d+1]);
   S := LUDecomposition(Transpose(dZM[d+1]),output=['P','L','U1','R']):
   M0,M1 := Transpose(S[4]),Transpose(S[1].S[2].S[3]):
   pivots := {}:
   j := 1:
   i := 1:
   while i <= nC[d+1] do
    while j <= nZ[d] and M0[j,i] = 0 do j := j+1; od;
    if j <= nZ[d] then pivots := {op(pivots),j}; fi;
    i := i+1;
   od:
   r := nops(pivots);
   nB[d] := r;
   nH[d] := nZ[d]-r;

   non_pivots := {seq(j,j=1..nZ[d])} minus pivots;

   M2 := Matrix(nZ[d],nZ[d]-r);
   for i from 1 to nZ[d]-r do M2[non_pivots[i],i] := 1; od:
   M3 := SubMatrix(M0,1..nZ[d],1..r);
   HM[d] := ZM[d] . M2;
   BM[d] := ZM[d] . M3;
   ZM[d] := <HM[d]|BM[d]>;
   HZM[d] := <IdentityMatrix(nH[d]),Matrix(nZ[d]-nH[d],nH[d])>;
   BZM[d] := <Matrix(nZ[d]-nB[d],nB[d]),IdentityMatrix(nB[d])>;
   LM[d+1] := (1/M1).<IdentityMatrix(r),Matrix(nC[d+1]-r,r)>;
  fi;
 od:

 H['cycle_matrix'] := eval(ZM);
 H['boundary_matrix'] := eval(BM);
 H['boundary_cycle_matrix'] := eval(BZM);
 H['boundary_lift_matrix'] := eval(LM);
 H['homology_matrix'] := eval(HM);
 H['homology_cycle_matrix'] := eval(HZM);

 H['cycle_number'] := eval(nZ);
 H['boundary_number'] := eval(nB);
 H['homology_number'] := eval(nH);

 return NULL;
end:

######################################################################

`check_bases/homology` := proc(H)
 local dM,dZM,nC,nZ,nB,nH,d,d_max,i,j,r,LM,ZM,BZM,BM,HZM,HM,
  checks;

 nC  := eval(H['chain_number']);
 nZ  := eval(H['cycle_number']);
 nB  := eval(H['boundary_number']);
 nH  := eval(H['homology_number']);
 ZM  := eval(H['cycle_matrix']);
 BM  := eval(H['boundary_matrix']);
 BZM := eval(H['boundary_cycle_matrix']);
 LM  := eval(H['boundary_lift_matrix']);
 HM  := eval(H['homology_matrix']);
 HZM := eval(H['homology_cycle_matrix']);
 d_max := H['dim'];
 dM := eval(H['differential_matrix']):

 for d from 0 to d_max do
  checks := [
   [Dimension( ZM[d])] = [nC[d],nZ[d]],
   [Dimension( BM[d])] = [nC[d],nB[d]],
   [Dimension(BZM[d])] = [nZ[d],nB[d]],
   [Dimension( LM[d])] = [nC[d],`if`(d>0,nB[d-1],0)],
   [Dimension( HM[d])] = [nC[d],nH[d]],
   [Dimension(HZM[d])] = [nZ[d],nH[d]],
   nC[d] = nH[d] + nB[d] + `if`(d>0,nB[d-1],0),
   nH[d] = H['betti_number'][d]
  ];

  if not(`and`(op(map(evalb,checks)))) then
   print([d,checks]);
   return false;
  fi;
 od:

 for d from 1 to d_max do
  checks := [
   Equal(dM[d] . ZM[d],Matrix(nC[d-1],nZ[d])),
   Equal(dM[d] . BM[d],Matrix(nC[d-1],nB[d])),
   Equal(dM[d] . LM[d],BM[d-1]),
   Rank(dM[d]) = nB[d-1],
   Rank(ZM[d]) = nZ[d],
   Rank(BM[d]) = nB[d],
   Rank(HM[d]) = nH[d],
   Rank(BZM[d]) = nB[d],
   Rank(HZM[d]) = nH[d]
  ];
  if not(`and`(op(map(evalb,checks)))) then
   print([d,checks]);
   return false;
  fi;
 od:
 return true;
end:

######################################################################

`vec/chains` := (H) -> (d::nonnegint) -> proc(u)
 local v,i,uc,uv;

 if u = 0 then
  return Vector(H['chain_number'][d]);
 elif type(u,`+`) then
  # Slightly convoluted because of funny behaviour of + on vectors
  return `+`(op(map(`vec/chains`(H)(d),[op(u)])));
 elif type(u,`*`) then
  uc,uv := selectremove(type,u,rational);
  if uc <> 1 then
   return uc * `vec/chains`(H)(d)(uv);
  fi;
 elif type(u,specfunc(chain)) then
  v := Vector(H['chain_number'][d]);
  i := H['chain_index'][d][u];
  v[i] := 1;
  return v;
 fi;
 return 'procname'(args);
end:

`vec/cycles/vec` := (H) -> (d::nonnegint) -> proc(v)
 LinearSolve(H['cycle_matrix'][d],v);
end:

`vec/cycles` := (H) -> (d::nonnegint) -> proc(u)
 `vec/cycles/vec`(H)(d)(`vec/chains`(H)(d)(u));
end:

`vec/homology/vec` := (H) -> (d::nonnegint) -> proc(v)
 SubVector(LinearSolve(H['cycle_matrix'][d],v),1..H['homology_number'][d]);
end:

`vec/homology` := (H) -> (d::nonnegint) -> proc(u)
 `vec/homology/vec`(H)(d)(`vec/chains`(H)(d)(u));
end:

`basis/homology` := (H) -> proc(d::nonnegint)
 convert(Transpose(H['homology_matrix'][d]) . Vector(H['chains'][d]),list);
end:

######################################################################

is_homology_sphere := proc(H,d_)
 local d,d1,b,i;
 d1 := H['dim'];
 if nargs > 1 then
  d := d_;
 else
  d := d1;
 fi;
 b := eval(H['betti_number']);

 if d = -1 then # The -1 sphere is empty
  return evalb(d1 = -1);
 fi;

 for i from 0 to d1 do
  if i <> 0 and i <> d and b[i] <> 0 then
   return false;
  fi;
 od;

 if d = 0 then
  return evalb(b[0] = 2);
 else
  return evalb(b[0] = 1 and b[d] = 1);
 fi;
end:

######################################################################

is_homology_bouquet := proc(H,d,r)
 local d1,b,i;
 d1 := H['dim'];
 b := eval(H['betti_number']);

 for i from 0 to d1 do
  if i <> 0 and i <> d and b[i] <> 0 then
   return false;
  fi;
 od;

 if d = 0 then
  return evalb(b[0] = r+1);
 else
  return evalb(b[0] = 1 and b[d] = r);
 fi;
end:

######################################################################

is_homology_acyclic := proc(H)
 local d1,b,i;
 d1 := H['dim'];
 if d1 = -1 then return false; fi;

 b := eval(H['betti_number']);
 if b[0] <> 1 then return false; fi;

 for i from 1 to d1 do
  if b[i] <> 0 then
   return false;
  fi;
 od;

 return true;
end:


text:
program Write;	{sample24}begin
                                writeln  (  'It''s OK?'  )end.

max_issues_repo_name: shouth/LanguageProcessing


text:
######################################################################

`is_element/preord` := (A::set) -> proc(R)
 global reason;

 if not `is_element/autorel`(A)(R) then
  reason := [convert(procname,string),"R is not a relation on A",R,A,reason];
  return false;
 fi;

 if not(`is_reflexive/autorel`(A)(R)) then
  reason := [convert(procname,string),"R is not reflexive",R];
  return false;
 fi;

 if not(`is_transitive/autorel`(A)(R)) then
  reason := [convert(procname,string),"R is not transitive",R];
  return false;
 fi;

 return true;
end;

######################################################################

`is_equal/preord`     := eval(`is_equal/autorel`);
`is_separated/preord` := eval(`is_separated/autorel`);

# We order the set of preorders by reverse inclusion, so that
# separated preorders are maximal.

`is_leq/preord`       := (A::set) -> proc(R,S)
 evalb(S minus R = {});
end;

`op/preord` := eval(`op/autorel`);

######################################################################

`random_element/preord` := (A::set) -> proc(R)
 local pi,n,r,i,a,b,N,RN,A2,RA;

 pi := `random_element/partitions`(A)();
 n := nops(pi);
 r := table();
 for i from 1 to n do
  for a in pi[i] do
   r[a] := i;
  od:
 od:
 N := {seq(i,i=1..n)};
 RN := `random_element/partord`(N)();
 A2 := {seq(seq([a,b],b in A),a in A)};
 RA := select(u -> member([r[u[1]],r[u[2]]],RN),A2);

 return RA;
end:

######################################################################

`is_total/preord` := (A::set) -> proc(R)
 local U,a,b;

 U := R union `op/preord`(A)(R);
 U := {seq(seq([a,b],b in A),a in A)} minus U;

 return evalb(U = {});
end;

######################################################################

`strictify/preord` := (A::set) -> proc(R)
 return R minus `op/preord`(A)(R);
end:

######################################################################

`res/preord` := (A::set,B::set) -> (R) ->
 select(xy -> member(xy[1],B) and member(xy[2],B),R);

######################################################################

`block_partition/preord` := (A::set) -> proc(R)
 local E,pi,X,x,B;

 E := R intersect map(ab -> [ab[2],ab[1]],R);
 pi := NULL;
 X := A;

 while X <> {} do
  x := X[1];
  B := map(ab -> ab[2],select(ab -> ab[1] = x,E));
  pi := pi,B;
  X := X minus B;
 od:

 pi := {pi};

 return pi;
end:

######################################################################
# Note that the function below is not actually a rank function, but
# it is an ingredient in certain rank functions to be defined
# elsewhere.

`rank/preord` := (A::set) -> proc(R)
 nops(`block_partition/preord`(A)(R));
end:

######################################################################

`bump/preord` := (A::set) -> proc(R)
 local pi,B,a,b,B0,R0;

 pi := `block_partition/preord`(A)(R);

 for B in pi do
  if nops(B) > 1 then
   a := B[1];
   B0 := B minus {a};
   R0 := R minus {seq([b,a],b in B0)};
   return R0;
  fi;
 od;

 return FAIL;
end:

######################################################################

`rank_table/preord` := (A::set) -> proc(R)
 local pi,pi0,a,b,r,k,B;

 pi := `block_partition/preord`(A)(R);
 pi0 := map(B -> B[1],pi);
 r := table():
 for B in pi do
  k := nops(select(a -> member([a,B[1]],R),pi0)) - 1;
  for b in B do r[b] := k; od;
 od;

 return eval(r);
end;

######################################################################

`separated_quotient/preord` := (A::set) -> proc(R)
 local A0,R0,a0,b0;

 A0 := `block_partition/preord`(A)(R);
 R0 := NULL;
 for a0 in A0 do
  for b0 in A0 do
   if member([a0[1],b0[1]],R) then
    R0 := R0,[a0,b0];
   fi;
  od;
 od;
 R0 := {R0};

 return [A0,R0];
end:

######################################################################

`preord_is_leq` := (A::set) -> (R) -> (a,b) -> member([a,b],R);

######################################################################

`preord_is_equiv` := (A::set) -> (R) -> (a,b) ->
  member([a,b],R) and member([b,a],R);

######################################################################

`preord_is_less` := (A::set) -> (R) -> (a,b) ->
  member([a,b],R) and not(member([b,a],R));


text:
with(networks):

# encode an array as a single number.
encode := (x,F) -> ( F[output](x[1])*100
	           + F[output](x[2])*10
	           + F[output](x[3])*1 ):

# Define the switch function.
switch := proc(a,b,c,F)
       local x3,x1;
       x1 := F[input](1);
       x3 := F[`+`]( x1, F[`+`]( x1,x1 ) ) ;
       F[`-`]( F[`*`]( x3, F[`*`](a,b) ),c );
end:

# Define my three cluster functions.
f:= (x,F) -> [x[1],x[2],switch(x[1],x[2],x[3],F)]:
g:= (x,F) -> [x[1],switch(x[1],x[3],x[2],F),x[3]]:
h:= (x,F) -> [switch(x[2],x[3],x[1],F),x[2],x[3]]:

# Give me the graph neighbors of x.
N := x -> neighbors(x,G):

# Given a prime F, return a graph G(M,F)
makegraph := proc(F,firstelement)
   local g;
   g := graph({encode(firstelement,F)},{}):
   bleg(g, firstelement, 1, F):
   eval(g);
end:

# useful function
bleg := proc(G,x,depth,F)
  if depth < 100 then
    step(G,x,f(x,F),depth,F);
    step(G,x,g(x,F),depth,F);
    step(G,x,h(x,F),depth,F);
  fi;
end:

# useful function
step := proc(G,x,y,depth,F)
  if (encode(y,F) in vertices(G)) then
    addedge({encode(x,F),encode(y,F)},G):
  else
  #  if markoff(x,F) then
      addvertex(encode(y,F),G):
      addedge({encode(x,F),encode(y,F)},G):
      bleg(G,y,depth+1,F);
  #  fi
  fi
end:

# Usefull functions.
# is this a markoff solution?  For debuging
markoff := (x,F) -> markoff2(x[1],x[2],x[3],F):
markoff2 := proc(x,y,z,F)
    local lhs,rhs,x1,x3;
    x1 := F[input](1);
    x3 := F[`+`](x1,x1,x1) ;
    lhs := F[`+`]( F[`^`](x,2), F[`^`](y,2), F[`^`](z,2));
    rhs := F[`*`]( x3, x,y,z);
    if lhs = rhs then true else false fi;
    #if (x^2 + y^2 + z^2 - 3*x*y*z) mod p = 0 then true else false fi;
end:

decode2 := proc(x)
    local n,one,ten,hun;
    n := 10;
    one := x mod n;
    ten := ((x-one)/n) mod n;
    hun := (x-one-(n*ten))/n^2;
    [hun,ten,one];
end:
decode := proc(x,F)
    local y;
    y := decode2(x);
    [F[input](y[1]), F[input](y[2]), F[input](y[3])];
end:

# Enumerate ALL of the solutions
enumerate := proc(p,k,F)
	local a,b,c,A,B,C,sols;
	sols := {}:
	for a from 0 to p^k-1 do
	    for b from 0 to p^k-1 do
		for c from 0 to p^k-1 do
		    A := F[input](a):
		    B := F[input](b):
		    C := F[input](c):
		    if markoff([A,B,C],F) then
			sols := sols union {[A,B,C]};
		    fi
		od
	    od
	od:
	sols;
end:

# here is where things begin to happen.
#let our field F be
p:=2:
k:=2:
F := GF(p,k):
#seq(F[input](x),x=0..p^k-1);
#G := makegraph(F,[F[input](1),F[input](1),F[input](1)]):
#nops(vertices(G));
#vertices(G);
#draw3d(G);
#draw(G);
sols := enumerate(p,k,F):
nops(sols);
sols;

#found solutions.
u := {}:
# set of connected components.
Gees := {}:
for x in sols do
    if not(x in u)  then
        G := makegraph(F, x);
        for y in vertices(G) do
            u := u union {decode(y,F)};
        od:
        Gees := Gees union {eval(G)};
    fi
od:
nops(Gees);
draw3d(Gees[1]);
seq(nops(vertices(Gees[i])), i = 1..nops(Gees));
#seq((vertices(Gees[i])), i = 1..nops(Gees));


text:
print 5;
print 6;
print 5 foo;
print 7;
print 8;

max_stars_repo_name: jgke/miniplc

max_issues_repo_name: jgke/miniplc

max_forks_repo_name: jgke/miniplc


text:
kernelopts(assertlevel=1):
assert_count := 0;
assert_verbosely := false;
assert_stop_on_fail := true;

#@ _ASSERT
_ASSERT := proc(p,s::string)
 global assert_count,assert_verbosely,assert_stop_on_fail,assert_fail_data;
 assert_count := assert_count+1;
 if assert_verbosely then printf("Checking: %s\n",s); fi;
 if not p then
  assert_fail_data := args[3..-1];
  if assert_stop_on_fail then
   error cat("Assertion failed: ",s);
  else
   printf(cat("Assertion failed: ",s,"\n"));
  fi;
 fi;
end:

######################################################################

# checklist is a list of checking functions.  check_all() runs all
# the checking functions in the list.  add_check() adds a function
# to the list.

checklist := []: #@ checklist

check_all := proc()
 global assert_count;
 local C;

 assert_count := 0;

 for C in checklist do C(); od:
end:

add_check := proc(C)
 global checklist;

 checklist := [op(checklist),C];
end:

######################################################################

num_steps := 100;

check_set_generic := proc(S::string,params := [])
 local
  small_context,params0,
  is_element,random_element,list_elements,count_elements,is_equal,is_leq,
  `has/is_element`,`has/random_element`,`has/list_elements`,
  `has/count_elements`,`has/is_equal`,`has/is_leq`,
  p,x,y,t,i,j,L,L_full,L_rand,L_rand_ok,L_rand_bad,found;

 printf("\nChecking %s\n",S);

 is_element := eval(convert(cat("is_small_element/",S),name));

 if type(is_element,procedure) then
  small_context := true;
  is_element      := eval(convert(cat("is_small_element/"    ,S),name));
  random_element  := eval(convert(cat("random_small_element/",S),name));
  list_elements   := eval(convert(cat("list_small_elements/" ,S),name));
  count_elements  := eval(convert(cat("count_small_elements/",S),name));
  is_equal        := eval(convert(cat("is_equal/"            ,S),name));
  is_leq          := eval(convert(cat("is_leq/"              ,S),name));
  printf("Elements are filtered by size\n");
 else
  small_context := false;
  is_element      := eval(convert(cat("is_element/"    ,S),name));
  random_element  := eval(convert(cat("random_element/",S),name));
  list_elements   := eval(convert(cat("list_elements/" ,S),name));
  count_elements  := eval(convert(cat("count_elements/",S),name));
  is_equal        := eval(convert(cat("is_equal/"      ,S),name));
  is_leq          := eval(convert(cat("is_leq/"        ,S),name));
 fi;

 `has/is_element`     := type(is_element,procedure);
 `has/random_element` := type(random_element,procedure);
 `has/list_elements`  := type(list_elements,procedure);
 `has/count_elements` := type(count_elements,procedure);
 `has/is_equal`       := type(is_equal,procedure);
 `has/is_leq`         := type(is_leq,procedure);

 if `has/is_element` then
  for p in params do is_element := eval(is_element(op(p))); od;
 elif is_element <> NULL then
  printf("Warning: no is_element function has been defined.\n");
 fi;

 if `has/random_element` then
  for p in params do random_element := eval(random_element(op(p))); od;
 elif random_element <> NULL then
  printf("Warning: no random_element function has been defined.\n");
 fi;

 if `has/list_elements` then
  for p in params do list_elements := eval(list_elements(op(p))); od;
 elif list_elements <> NULL then
  printf("Warning: no list_elements function has been defined.\n");
 fi;

 if `has/count_elements` then
  for p in params do count_elements := eval(count_elements(op(p))); od;
 elif count_elements <> NULL then
  printf("Warning: no count_elements function has been defined.\n");
 fi;

 params0 := params;
 if small_context then
  params0 := [op(1..-2,params)];
 fi;

 if `has/is_equal` then
  for p in params0 do is_equal := eval(is_equal(op(p))); od;
 else
  is_equal := (x,y) -> evalb(x = y);
 fi;

 if `has/is_leq` then
  for p in params0 do is_leq := eval(is_leq(op(p))); od;
 elif is_leq <> NULL then
  printf("Warning: no is_leq function has been defined.\n");
 fi;

 L := NULL;

 if `has/random_element` and `has/is_element` then
  printf("Checking is_element on random elements: ");
  L_rand := [seq(random_element(),i=1..10)];
  L_rand_ok,L_rand_bad := selectremove(x -> (x <> FAIL),L_rand);
  t := `and`(true,op(map(is_element,L_rand_ok)));
  if t then
   L := L_rand_ok;
   if L_rand_bad = [] then
    printf("OK\n");
   elif L_rand_ok <> [] then
    printf("OK (but some examples could not be generated)\n");
   else
    printf("Could not generate any examples\n");
   fi;
  else
   printf("FAIL\n");
  fi;
 fi;

 if `has/list_elements` and `has/is_element` then
  printf("Checking is_element on listed elements: ");
  L_full := list_elements;
  if type([L_full],[list]) then
   L := L_full;
   t := `and`(op(map(is_element,L_full)));
   if t then
    printf("OK (%d elements)\n",nops(L_full));
   else
    printf("FAIL\n");
   fi;
  else
   printf("Element list is not a list\n");
   L := [];
   L_full := [];
  fi;
 fi;

 if `has/random_element` and `has/list_elements` and `has/is_element`
     and L_rand_ok <> [] then
  printf("Comparing random elements and listed elements: ");
  t := true;
  for x in L_rand_ok do
   found := false;
   for y in L_full do
    if is_equal(x,y) then
     found := true;
     break;
    fi;
   od;
   if not found then
    t := false;
    break;
   fi;
  od;
  printf(`if`(t,"OK\n","FAIL\n"));
 fi;

 if L <> NULL and `has/is_equal` then
  printf("Checking is_equal: ");
  t := true;
  for i from 1 to nops(L) do
   t := t and is_equal(L[i],L[i]);
   if not t then break; fi;
   if `has/list_elements` then
    for j from i+1 to nops(L) do
     t := t and not(is_equal(L[i],L[j]));
     if not t then break; fi;
    od;
   fi;
   if not t then break; fi;
  od;
  printf(`if`(t,"OK\n","FAIL\n"));
 fi;

 if `has/list_elements` and `has/count_elements` then
  printf("Checking count_elements (%A) vs listed elements: ",count_elements);
  t := evalb(count_elements = nops(L_full));
  printf(`if`(t,"OK\n","FAIL\n"));
 fi;
end:

######################################################################

check_map_generic := proc(fs::string,X::string,Y::string,params_f,params_X,params_Y)
 global reason;
 local P_f,P_X,P_Y,list_elements_X,random_element_X,is_element_X,is_element_Y,p,x,y,L,f,ok,i;

 printf("Checking %s : %s -> %s\n",fs,X,Y);

 P_f := `if`(nargs > 3,params_f,[]);
 P_X := `if`(nargs > 4,params_X,P_f);
 P_Y := `if`(nargs > 5,params_Y,P_X);

 list_elements_X := eval(convert(cat("list_elements/",X),name));
 if list_elements_X <> NULL then
  for p in P_X do
   list_elements_X := eval(list_elements_X(op(p)));
  od;
  L := list_elements_X;
 else
  random_element_X := eval(convert(cat("random_element/",X),name));
  if random_element_X = NULL then
   error sprintf("Neither list_elements/%s nor random_element/%s is defined",X,X);
  fi;
  for p in P_X do
   random_element_X := eval(random_element_X(op(p)));
  od;
  L := [seq(random_element_X(),i=1..num_steps)];
 fi;

 f := eval(convert(fs,name));
 for p in P_f do f := eval(f(op(p))); od;

 is_element_X := eval(convert(cat("is_element/",X),name));
 for p in P_X do
  is_element_X := eval(is_element_X(op(p)));
 od;

 is_element_Y := eval(convert(cat("is_element/",Y),name));
 for p in P_Y do
  is_element_Y := eval(is_element_Y(op(p)));
 od;

 ok := true;
 for x in L do
  if not(is_element_X(x)) then
   ok := false;
   reason := ["x is not in X",x];
   break;
  fi;

  y := f(x);

  if not(is_element_Y(y)) then
   ok := false;
   reason := ["f(x) = y is not in Y",x,y];
   break;
  fi;
 od:

 printf(`if`(ok,"OK\n","FAIL\n"));
end:

######################################################################

check_binop_generic :=
 proc(fs::string,X::string,Y::string,Z::string,params_f,params_X,params_Y,params_Z)
 global reason;
 local P_f,P_X,P_Y,P_Z,
       list_elements_X,random_element_X,is_element_X,
       list_elements_Y,random_element_Y,is_element_Y,
       is_element_Z,p,x,y,z,L_X,L_Y,f,ok,ns,curry,i;

 printf("Checking %s : %s x %s -> %s\n",fs,X,Y,Z);

 P_f := `if`(nargs > 3,params_f,[]);
 P_X := `if`(nargs > 4,params_X,P_f);
 P_Y := `if`(nargs > 5,params_Y,P_X);
 P_Z := `if`(nargs > 5,params_Z,P_Y);
 curry := false;

 if member("curry",[_rest]) then
  curry := true;
 fi;

 ns := ceil(sqrt(evalf(num_steps)));

 list_elements_X := eval(convert(cat("list_elements/",X),name));
 if list_elements_X <> NULL then
  for p in P_X do
   list_elements_X := eval(list_elements_X(op(p)));
  od;
  L_X := list_elements_X;
 else
  random_element_X := eval(convert(cat("random_element/",X),name));
  if random_element_X = NULL then
   error sprintf("Neither list_elements/%s nor random_element/%s is defined",X,X);
  fi;
  for p in P_X do
   random_element_X := eval(random_element_X(op(p)));
  od;
  L_X := [seq(random_element_X(),i=1..ns)];
 fi;

 list_elements_Y := eval(convert(cat("list_elements/",Y),name));
 if list_elements_Y <> NULL then
  for p in P_Y do
   list_elements_Y := eval(list_elements_Y(op(p)));
  od;
  L_Y := list_elements_Y;
 else
  random_element_Y := eval(convert(cat("random_element/",Y),name));
  if random_element_Y = NULL then
   error sprintf("Neither list_elements/%s nor random_element/%s is defined",Y,Y);
  fi;
  for p in P_Y do
   random_element_Y := eval(random_element_Y(op(p)));
  od;
  L_Y := [seq(random_element_Y(),i=1..ns)];
 fi;

 f := eval(convert(fs,name));
 for p in P_f do f := eval(f(op(p))); od;

 is_element_X := eval(convert(cat("is_element/",X),name));
 for p in P_X do
  is_element_X := eval(is_element_X(op(p)));
 od;

 is_element_Y := eval(convert(cat("is_element/",Y),name));
 for p in P_Y do
  is_element_Y := eval(is_element_Y(op(p)));
 od;

 is_element_Z := eval(convert(cat("is_element/",Z),name));
 for p in P_Z do
  is_element_Z := eval(is_element_Z(op(p)));
 od;

 ok := true;
 for x in L_X do
  if not(is_element_X(x)) then
   ok := false;
   reason := ["x is not in X",x];
   break;
  fi;

  for y in L_Y do
   if not(is_element_Y(y)) then
    ok := false;
    reason := ["y is not in Y",y];
    break;
   fi;

   if curry then
    z := f(x,y);
   else
    z := f(x)(y);
   fi;

   if not(is_element_Z(z)) then
    ok := false;
    reason := ["f(x,y) = z is not in Z",x,y,z];
    break;
   fi;
  od;
 od:

 printf(`if`(ok,"OK\n","FAIL\n"));
end:

######################################################################

check_operad := proc(O::string,params,A,B,C)
 global reason;
 local is_element,is_equal,random_element,eta,gamma,
       x,i,p,q,qp,Fp,Fq,Fqp,U,V,W,ok,c,b;

 printf("Checking operad structure for %s: ",O);

 is_element     := eval(convert(cat("is_element/",O),name));
 is_equal       := eval(convert(cat("is_equal/",O),name));
 random_element := eval(convert(cat("random_element/",O),name));
 eta            := eval(convert(cat("eta/",O),name));
 gamma          := eval(convert(cat("gamma/",O),name));

 for x in params do
  is_element     := eval(is_element(op(x)));
  is_equal       := eval(is_equal(op(x)));
  random_element := eval(random_element(op(x)));
  eta            := eval(eta(op(x)));
  gamma          := eval(gamma(op(x)));
 od:

 for i from 1 to num_steps do
  p := `random_element/epi`(A,B)();
  q := `random_element/epi`(B,C)();
  qp := `compose/maps`(A,B,C)(p,q);
  Fp := fibres(A,B)(p);
  Fq := fibres(B,C)(q);
  Fqp := fibres(A,C)(qp);
  U := random_element(C)();
  V := table({seq(c = random_element(Fq[c])(),c in C)});
  W := table({seq(b = random_element(Fp[b])(),b in B)});
  ok := `check_gamma_axiom`(A,B,C,p,q,U,V,W,
   is_element,is_equal,gamma);
  if not ok then printf("FAIL\n");
   reason := ["gamma axiom failed",A,B,C,eval(p),eval(q),eval(U),eval(V),eval(W),reason];
   return false;
  fi;
 od:

 printf("OK\n");
end:

check_operad_morphism := proc(f::string,f_params,
                              C::string,C_params,
                              D::string,D_params,
                              A,B)
 global reason;
 local is_el_CC,is_equal_CC,gamma_CC,
       is_el_DD,is_equal_DD,gamma_DD,
       rand_CC,phi,x,i,p,Fp,U,V,ok,b;

 printf("Checking operad morphism %s : %s -> %s\n",f,C,D);

 is_el_CC    := eval(convert(cat("is_element/",C),name));
 is_equal_CC := eval(convert(cat("is_equal/",C),name));
 gamma_CC    := eval(convert(cat("gamma/",C),name));
 rand_CC     := eval(convert(cat("random_element/",C),name));

 is_el_DD    := eval(convert(cat("is_element/",D),name));
 is_equal_DD := eval(convert(cat("is_equal/",D),name));
 gamma_DD    := eval(convert(cat("gamma/",D),name));

 phi         := eval(convert(f,name));

 for x in f_params do
  phi := eval(phi(op(x)));
 od:

 for x in C_params do
  is_el_CC    := eval(is_el_CC(op(x)));
  is_equal_CC := eval(is_equal_CC(op(x)));
  gamma_CC    := eval(gamma_CC(op(x)));
  rand_CC     := eval(rand_CC(op(x)));
 od;

 for x in D_params do
  is_el_DD    := eval(is_el_DD(op(x)));
  is_equal_DD := eval(is_equal_DD(op(x)));
  gamma_DD    := eval(gamma_DD(op(x)));
 od:

 for i from 1 to num_steps do
  p := `random_element/epi`(A,B)();
  Fp := fibres(A,B)(p);
  U := rand_CC(B)();
  V := table({seq(b = rand_CC(Fp[b])(),b in B)});
  ok := `check_operad_morphism_axiom`(A,B,p,U,V,
      is_el_CC,is_equal_CC,gamma_CC,
      is_el_DD,is_equal_DD,gamma_DD,
      phi);
  if not ok then printf("FAIL\n");
   reason := ["operad morphism axiom failed",A,B,eval(p),eval(U),eval(V),reason];
   return false;
  fi;
 od:

 printf("OK\n");
end:

######################################################################

check_semioperad := proc(O::string,params,L)
 global reason;
 local is_element,is_equal,random_element,gamma,
       x,J,L_,i,k,jj,ii,H,U,V,W,ok,j0,k0;

 printf("Checking operad structure for %s: ",O);

 is_element     := eval(convert(cat("is_element/",O),name));
 is_equal       := eval(convert(cat("is_equal/",O),name));
 random_element := eval(convert(cat("random_element/",O),name));
 gamma          := eval(convert(cat("gamma/",O),name));

 for x in params do
  is_element     := eval(is_element(op(x)));
  is_equal       := eval(is_equal(op(x)));
  random_element := eval(random_element(op(x)));
  gamma          := eval(gamma(op(x)));
 od:

 J := map(nops,L);
 L_ := map(op,L);
 k := nops(J);
 jj := `+`(op(J));
 ii := `+`(op(L_));
 H := [seq(`+`(op(L[k0])),k0=1..k)];

 for i from 1 to num_steps do
  U := random_element(k)();
  V := [seq(random_element(J[k0])(),k0 = 1..k)];
  W := [seq([seq(random_element(L[k0][j0])(),j0=1..J[k0])],k0=1..k)];
  ok := `check_semi_gamma_axiom`(L,U,V,W,
   is_element,is_equal,gamma);
  if not ok then printf("FAIL\n");
   reason := ["gamma axiom failed",L,eval(U),eval(V),eval(W),reason];
   return false;
  fi;
 od:

 printf("OK\n");
end:

######################################################################

check_semi_circ := proc(O::string,params,l,m,n)
 global reason;
 local is_element,is_equal,random_element,circ,
       x,c,i,j,k,U,V,W,ok;

 printf("Checking operad structure for %s: ",O);

 is_element     := eval(convert(cat("is_element/",O),name));
 is_equal       := eval(convert(cat("is_equal/",O),name));
 random_element := eval(convert(cat("random_element/",O),name));
 circ           := eval(convert(cat("circ/",O),name));

 for x in params do
  is_element     := eval(is_element(op(x)));
  is_equal       := eval(is_equal(op(x)));
  random_element := eval(random_element(op(x)));
  circ           := eval(circ(op(x)));
 od:

 for c from 1 to num_steps do
  U := random_element(l)();
  V := random_element(m)();
  W := random_element(n)();
  i := rand(1..l)();
  j := rand(1..m)();
  ok := `check_semi_circ_axiom_nested`(i,j,l,m,n,U,V,W,
   is_element,is_equal,circ);
  if not ok then printf("FAIL\n");
   reason := ["nested circ axiom failed",i,j,l,m,n,eval(U),eval(V),eval(W),reason];
   return false;
  fi;

  i := rand(1..(l-1))();
  k := rand((i+1)..l)();
  ok := `check_semi_circ_axiom_disjoint`(i,k,l,m,n,U,V,W,
   is_element,is_equal,circ);
  if not ok then printf("FAIL\n");
   reason := ["disjoint circ axiom failed",i,j,l,m,n,eval(U),eval(V),eval(W),reason];
   return false;
  fi;
 od:

 printf("OK\n");
end:


text:
######################################################################

`is_element/cycord` := (A::set) -> proc(C)
 global reason;

 if not type(C,list) then
  reason := ["is_element/cycord","C is not a list",C];
  return false;
 fi;

 if nops(C) <> nops(A) then
  reason := ["is_element/cycord","C does not the same length as A",C,A];
  return false;
 fi;

 if {op(C)} <> A then
  reason := ["is_element/cycord","C is not an enumeration of A",C,A];
  return false;
 fi;

 return true;
end;

######################################################################

`is_equal/cycord` := (A::set) -> proc(C0,C1)
 local n,a,i,C2;

 n := nops(A);
 if n <= 2 then
  return true;
 fi;

 a := C0[1];
 i := 1;
 while i < n and C1[i] <> a do i := i+1; od;

 C2 := [op(i..n,C1),op(1..i-1,C1)];
 return evalb(C0 = C2);
end;

`is_leq/cycord` := NULL;

######################################################################

`op/cycord` := (A::set) -> proc(C)
 local n,i;
 n := nops(A);
 return [seq(C[n-i],i=0..n-1)];
end;

`res/cycord` := (A::set,B::set) -> proc(C)
 return select(a -> member(a,B),C);
end:

######################################################################

`random_element/cycord` := (A::set) -> proc()
 return combinat[randperm](A);
end:

`list_elements/cycord` := proc(A::set)
 local a,B;

 a := A[1];
 B := A minus {a};
 return map(C -> [a,op(C)],combinat[permute](B));
end:

`count_elements/cycord` := proc(A::set)
 return max(0,nops(A) - 1) !;
end:

######################################################################

`cut/cycord` := (A::set) -> (C) -> proc(a)
 local n,i;
 n := nops(C);
 if n <= 1 then return C; fi;
 i := 1;
 while i <= n and C[i] <> a do i := i+1; od;
 if i > n then error("a is not in C"); fi;
 return [op(i..n,C),op(1..i-1,C)];
end:

######################################################################

`cycord_is_strictly_cyclic` := (A::set) -> (C) -> proc(x::list)
 local n,S,a,i,y;

 n := nops(x);
 S := {op(x)};
 if n < 2 then return true; fi;
 if nops(S) < n then return false; fi;
 if S minus A <> {} then return false; fi;
 y := `cut/cycord`(A)(C)(x[1]);
 y := select(i -> member(i,S),y);
 return evalb(x = y);
end:

######################################################################

`cycord_is_cyclic` := (A::set) -> (C) -> proc(x)
 local n,y,i;
 n := nops(x);
 y := NULL;
 for i from 1 to n do
  if (i = n and x[n] <> x[1]) or (i < n and x[i] <> x[i+1]) then
   y := y,x[i];
  fi;
 od;
 y := [y];
 return `cycord_is_strictly_cyclic`(A)(C)(y);
end:

######################################################################

`open_interval/cycord` := (A::set) -> (C) -> proc(a,b)
 local n,i,J,C0;
 if a = b then return []; fi;
 C0 := `cut/cycord`(A)(C)(a);
 J := NULL;
 i := 2;
 n := nops(C0);
 while i <= n and C0[i] <> b do i := i+1; od;
 return [op(2..i-1,C0)];
end:

`half_open_interval/cycord` := (A::set) -> (C) -> proc(a,b)
 return [a,op(`open_interval/cycord`(A)(C)(a,b))];
end:

`closed_interval/cycord` := (A::set) -> (C) -> proc(a,b)
 return [a,op(`open_interval/cycord`(A)(C)(a,b)),b];
end:

`are_crossed/cycord` := (A::set) -> (C) -> proc(U,V)
 local C0,n,i,j,k;
 if U intersect V <> {} then return true; fi;
 if U = {} or V = {} then return false; fi;
 C0 := `cut/cycord`(A)(C)(U[1]);
 n := nops(C0);
 i := 2;
 while i <= n and not(member(C0[i],V)) do i := i+1; od;
 if i > n then return false; fi;
 j := i+1;
 while j <= n and not(member(C0[j],U)) do j := j+1; od;
 if j > n then return false; fi;
 k := j+1;
 while k <= n and not(member(C0[k],V)) do k := k+1; od;
 if k > n then return false; fi;
 return true;
end:

`is_convex/cycord` := (A::set) -> (C) -> proc(W)
 return not(`are_crossed/cycord`(A)(C)(W,A minus W));
end;


text:
read "../IdentifiabilityODE.mpl";

sys := [
diff(EGF_EGFR(t), t) = reaction_1_k1*EGF_EGFR(t) - reaction_9_k1*EGF_EGFR(t) - reaction_1_k2*EGF_EGFR(t),
diff(EGFR(t), t) = -reaction_1_k1*EGF_EGFR(t) + reaction_1_k2*EGF_EGFR(t) - EGFR(t)*EGFR_turnover + EGFR_turnover*pro_EGFR(t),
diff(pEGFR(t), t) = -reaction_4_k1*pEGFR(t) + reaction_9_k1*EGF_EGFR(t) - pEGFR(t)*Akt(t)*reaction_2_k1 + reaction_3_k1*pEGFR_Akt(t) + pEGFR_Akt(t)*reaction_2_k2,
diff(pAkt_S6(t), t) = pAkt(t)*reaction_5_k1*S6(t) - reaction_6_k1*pAkt_S6(t) - pAkt_S6(t)*reaction_5_k2,
diff(pEGFR_Akt(t), t) = pEGFR(t)*Akt(t)*reaction_2_k1 - reaction_3_k1*pEGFR_Akt(t) - pEGFR_Akt(t)*reaction_2_k2,
diff(pAkt(t), t) = -pAkt(t)*reaction_7_k1 - pAkt(t)*reaction_5_k1*S6(t) + reaction_6_k1*pAkt_S6(t) + reaction_3_k1*pEGFR_Akt(t) + pAkt_S6(t)*reaction_5_k2,
diff(S6(t), t) = pS6(t)*reaction_8_k1 - pAkt(t)*reaction_5_k1*S6(t) + pAkt_S6(t)*reaction_5_k2,
diff(pS6(t), t) = -pS6(t)*reaction_8_k1 + reaction_6_k1*pAkt_S6(t),
diff(Akt(t), t) = pAkt(t)*reaction_7_k1 - pEGFR(t)*Akt(t)*reaction_2_k1 + pEGFR_Akt(t)*reaction_2_k2,
y1(t) = pEGFR(t)*a1 + a1*pEGFR_Akt(t),
y2(t) = a2*pAkt(t) + a2*pAkt_S6(t),
y3(t) = pS6(t)*a3
];
CodeTools[CPUTime](IdentifiabilityODE(sys, GetParameters(sys)));


text:
(*
 Copyright (C) 2017 M.A.L. Marques

 This Source Code Form is subject to the terms of the Mozilla Public
 License, v. 2.0. If a copy of the MPL was not distributed with this
 file, You can obtain one at http://mozilla.org/MPL/2.0/.
*)

(* type: mgga_exc *)
(* prefix:
  mgga_x_tau_hcth_params *params;

  assert(p->params != NULL);
  params = (mgga_x_tau_hcth_params * ) (p->params);
*)

hcth_coeff_a := [0, 1, 0, -2, 0, 1]:

(* Equation (29) *)
hcth_gamX := 0.004:
hcth_ux   := x -> hcth_gamX*x^2/(1 + hcth_gamX*x^2):

hcth_gxl  := x -> add(params_a_cx_local [i]*hcth_ux(x)^(i-1), i=1..4):
hcth_gxnl := x -> add(params_a_cx_nlocal[i]*hcth_ux(x)^(i-1), i=1..4):

hcth_f    := (x, u, t) -> hcth_gxl(x) + hcth_gxnl(x)*mgga_series_w(hcth_coeff_a, 6, t):

f := (rs, z, xt, xs0, xs1, u0, u1, t0, t1) -> mgga_exchange(hcth_f, rs, z, xs0, xs1, u0, u1, t0, t1):


text:
TITLE

    DakotaModel_deterministic;
    ! Objective value 4165.0

INDEX

    PRODUCTS := (Desk, Table, Chair);
    RESOURCES := (Lumber, Finishing, Carpentry);

DATA

    Price[PRODUCTS] := (60, 40, 10);
    Cost[RESOURCES] := (2, 4, 5.2);
    Consumption[RESOURCES, PRODUCTS] := (8, 6, 1
					4, 2, 1.5
					2, 1.5, .5);
    Demand[PRODUCTS] := (150, 125, 300);

BINARY VARIABLES

   BinVec[PRODUCTS];

VARIABLES

   Purchase[RESOURCES];
   Production[PRODUCTS];
   BinaryVar;

MIN

    Profit = -sum(PRODUCTS: Price * Production) + sum(RESOURCES: Cost * Purchase) - BinaryVar - sum(PRODUCTS: BinVec);

SUBJECT TO

    ResourceConstraint[RESOURCES]: sum(PRODUCTS: Consumption * Production) <= Purchase;
    DemandConstraint[PRODUCTS]: Production <= Demand;

BOUNDS

    DeskBound: Production[PRODUCTS=Desk] <= 160 ;
    TableBound: 130 >= Production[PRODUCTS=Table] ;

FREE

    Purchase[RESOURCES=Lumber];

BINARY

   BinaryVar;

INTEGER

   Production[PRODUCTS];

END

max_stars_repo_name: pashtetGP/optconvert

max_issues_repo_name: pashtetGP/optconvert

max_forks_repo_name: pashtetGP/optconvert


text:
read "../IdentifiabilityODE.mpl";

sys := [
diff(KS01(t), t) = c0001*KS00(t) + FS11(t)*gamma1101 - F(t)*alpha01*S01(t) + FS01(t)*beta01 - K(t)*S01(t)*a01 + KS01(t)*b01,
diff(FS10(t), t) = -c0001*KS00(t) - b00*KS00(t) + S00(t)*K(t)*a00 - c0011*KS00(t) - c0010*KS00(t),
diff(S11(t), t) = FS11(t)*beta11 + FS11(t)*gamma1100 + FS11(t)*gamma1101 + FS11(t)*gamma1110 - F(t)*alpha10*S10(t) - F(t)*alpha11*S11(t) - F(t)*alpha01*S01(t) + gamma1000*FS10(t) + FS01(t)*beta01 + FS01(t)*gamma0100 + beta10*FS10(t),
diff(S00(t), t) = F(t)*alpha10*S10(t) - gamma1000*FS10(t) - beta10*FS10(t),
diff(S10(t), t) = c0001*KS00(t) - a10*K(t)*S10(t) + KS10(t)*b10 + KS10(t)*c1011 + b00*KS00(t) + c0111*KS01(t) - S00(t)*K(t)*a00 - K(t)*S01(t)*a01 + c0011*KS00(t) + c0010*KS00(t) + KS01(t)*b01,
diff(F(t), t) = F(t)*alpha01*S01(t) - FS01(t)*beta01 - FS01(t)*gamma0100,
diff(KS00(t), t) = FS11(t)*gamma1100 + gamma1000*FS10(t) + FS01(t)*gamma0100 + b00*KS00(t) - S00(t)*K(t)*a00,
diff(S01(t), t) = -FS11(t)*beta11 - FS11(t)*gamma1100 - FS11(t)*gamma1101 - FS11(t)*gamma1110 + F(t)*alpha11*S11(t),
diff(KS10(t), t) = -a10*K(t)*S10(t) + FS11(t)*gamma1110 - F(t)*alpha10*S10(t) + KS10(t)*b10 + beta10*FS10(t) + c0010*KS00(t),
diff(FS01(t), t) = FS11(t)*beta11 - F(t)*alpha11*S11(t) + KS10(t)*c1011 + c0111*KS01(t) + c0011*KS00(t),
diff(FS11(t), t) = -c0111*KS01(t) + K(t)*S01(t)*a01 - KS01(t)*b01,
diff(K(t), t) = a10*K(t)*S10(t) - KS10(t)*b10 - KS10(t)*c1011,
y1(t) = F(t),
y5(t) = S11(t),
y3(t) = S01(t),
y0(t) = K(t),
y4(t) = S10(t),
y2(t) = S00(t)
];
CodeTools[CPUTime](IdentifiabilityODE(sys, GetParameters(sys)));


text:
# Step 1 of 3: from Hakaru to Maple LO (linear operator)

toLO := module()
  export
  ModuleApply := proc(m, $)
    local h;
    h := gensym('h');
    _Env_HakaruSolve := false;
    LO(h, integrate(m, h, []))
  end proc;

  export
  integrate := proc(m, h, loops :: list(name = range) := [], $)
    local x, n, i, res, l, br, m0;

    if m :: known_continuous then
      integrate_known(Int, Ints, 'xx', m, h, loops)
    elif m :: known_discrete then
      integrate_known(Sum, Sums, 'kk', m, h, loops)
    elif m :: 'Ret(anything)' then
      applyintegrand(h, mk_ary(op(1,m), loops))
    elif m :: 'Bind(anything, name, anything)' then
      res := eval(op(3,m), op(2,m) = mk_idx(op(2,m), loops));
      res := eval(Integrand(op(2,m), 'integrate'(res, x, loops)), x=h);
      integrate(op(1,m), res, loops);
    elif m :: 'specfunc(Msum)' and (nops(loops) = 0 or nops(m) = 1) then
      `+`(op(map(integrate, [op(m)], h, loops)))
    elif m :: 'Weight(anything, anything)' then
      foldl(product, op(1,m), op(loops)) * integrate(op(2,m), h, loops)
    elif m :: t_pw_or_part then
      Partition:-AppPartOrPw
                 ( proc(p)
                       if not Partition:-ConditionsDepend(p, lhs~(loops)) then
                           Partition:-Pmap(z->integrate(z,h,loops),p)
                       end if;
                   end proc
                 , m);
    elif m :: t_case and not depends(op(1,m), map(lhs, loops)) then
      subsop(2=map(proc(b :: Branch(anything, anything))
                     eval(subsop(2='integrate'(op(2,b), x, loops),b), x=h)
                   end proc,
                   op(2,m)),
             m);
    elif m :: 'LO(name, anything)' then
      eval(op(2,m), op(1,m) = h)
    elif m :: 'Plate(nonnegint, name, anything)' then
      # Unroll Plate when the bound is known. We unroll Plate here (instead
      # of unrolling Ints in reduce, for example) because we have no other
      # way to integrate certain Plates, namely those whose bodies' control
      # flows depend on the index.
      x := mk_sym('pp', h);
      x := [seq(cat(x,i), i=0..op(1,m)-1)];
      if op(1,m) = 0 then
        res := undefined;
      else
        if op(1,m) = 1 then
          res := op(1,x);
        else
           res := idx(x,op(2,m));
#          res := piecewise(seq(op([op(2,m)=i-1, op(i,x)]), i=2..op(1,m)),
#                         op(1,x));
        end if;
        res := mk_idx(res, loops);
      end if;
      res := applyintegrand(h, mk_ary('ary'(op(1,m), op(2,m), res), loops));
      for i from op(1,m) to 1 by -1 do
        res := integrate(eval(op(3,m), op(2,m)=i-1),
                         Integrand(op(i,x), res), loops);
      end do;
      res
    elif m :: 'Plate(anything, name, anything)' then
      integrate(op(3,m), h, [op(2,m)=0..op(1,m)-1, op(loops)])
    elif m :: 'Context(anything, anything)' then
      applyop(integrate, 2, m, h, loops)
    elif h :: appliable then
      x := gensym('xa');
      'integrate'(m, Integrand(x, h(x)), loops)
    else
      'procname(_passed)'
    end if
  end proc;

  local
  integrate_known := proc(make, makes, var, m, h, loops :: list(name=range), $)
    local x, dens, bds;
    x := mk_sym(var, h);
    dens := density[op(0,m)](op(m));
    bds := bounds[op(0,m)](op(m));
    if loops = [] then
      make(dens(x) * applyintegrand(h, x), x = bds);
    else
      makes(foldl(product, dens(mk_idx(x,loops)), op(loops))
              * applyintegrand(h, x),
            x, bds, loops)
    end if;
  end proc;
end module; #toLO

max_stars_repo_name: zaxtax/hakaru

max_issues_repo_name: zaxtax/hakaru

max_forks_repo_name: zaxtax/hakaru


text:
(*
 Copyright (C) 2019 Susi Lehtola

 This Source Code Form is subject to the terms of the Mozilla Public
 License, v. 2.0. If a copy of the MPL was not distributed with this
 file, You can obtain one at http://mozilla.org/MPL/2.0/.
*)

(* type: gga_exc *)
(* prefix:
  gga_x_ncap_params *params;

  assert(p->params != NULL);
  params = (gga_x_ncap_params * )(p->params);
*)

ncap_f0 := s -> 1 + params_a_mu*tanh(s)*arcsinh(s)*( 1 + params_a_alpha*((1-params_a_zeta)*s*log(1+s) + params_a_zeta*s))/(1 + params_a_beta*tanh(s)*arcsinh(s)):
ncap_f  := x -> ncap_f0(X2S*x):

f := (rs, zeta, xt, xs0, xs1) -> gga_exchange(ncap_f, rs, zeta, xs0, xs1):


text:
######################################################################
# prime_simplex_boundary(A) is the set of maps x : A -> [0,1] with
# min = 0 and max = 1.

`is_element/prime_simplex_boundary` := (A::set) -> proc(x)
 local a,u,v,w;
 global reason;

 if not type(x,table) then
  reason := [convert(procname,string),"x is not a table",x];
  return false;
 fi;

 if map(op,{indices(x)}) <> A then
  reason := [convert(procname,string),"x is not indexed by A",x,A];
  return false;
 fi;

 v := 0;
 w := 1;

 for a in A do
  u := x[a];

  if not (`is_element/RR`(u) and u >= 0 and u <= 1) then
   reason := [convert(procname,string),"x[a] is not in the unit interval",a,u];
   return false;
  fi;

  v := max(u,v);
  w := min(u,w);
 od;

 if simplify(v - 1) <> 0 then
  reason := [convert(procname,string),"max(x[a] : a in A) <> 1",v];
  return false;
 fi;

 if simplify(w) <> 0 then
  reason := [convert(procname,string),"min(x[a] : a in A) <> 0",w];
  return false;
 fi;

 return true;
end;

######################################################################

`is_equal/prime_simplex_boundary` := (A::set) -> proc(x,y)
 local a;
 global reason;

 for a in A do
  if simplify(x[a] - y[a]) <> 0 then
   reason := [convert(procname,string),"x[a] <> y[a]",a,x[a],y[a]];
   return false;
  fi;
 od;

 return true;
end;

######################################################################

`is_leq/prime_simplex_boundary` := NULL;

######################################################################

`random_element/prime_simplex_boundary` := (A::set) -> proc(d::posint := 5)
 local x,a,r1,r2,u,v;

 if nops(A) < 2 then return FAIL; fi;

 r1 := rand(0..d);
 r2 := rand(1..d);
 x := table();
 u := 0;
 v := 0;

 while u = v do
  u := 0;
  v := infinity;
  for a in A do
   x[a] := r1()/r2();
   u := max(u,x[a]);
   v := min(v,x[a]);
  od;
 od;

 for a in A do x[a] := (x[a]-v)/(u-v); od;

 return eval(x);
end;

######################################################################

`list_elements/prime_simplex_boundary` := NULL;
`count_elements/prime_simplex_boundary` := NULL;

######################################################################

`phi/proper_nonempty_subsets/prime_simplex_boundary` := (A::set) -> proc(U)
 local x,a;

 x := table();
 for a in A do x[a] := 0; od;
 for a in U do x[a] := 1; od;

 return eval(x);
end;


text:
# ClusterComplexRoots.mpl
#
# clusteredList := ClusterComplexRoots( rootList, tolerance=s, MaxOrder=M );
#
# Take a list rootList of complex numbers and cluster them
# according to apparent multiplicity.
#
# This routine is heuristic (the problem is hard) and is intended
# to work only up to the given maximum order (default MaxOrder is 3)
#
# The routine uses distances and a fuzz factor F[m] and only considers
# neighbouring points up to about |r[i]-r[j]| <= F[m]*s^(1/m) to be
# possible candidates for a cluster; it also uses approximate
# angle symmetry, looking for angles between adjacent roots to be
# approximately 2*Pi/m .
#
# The returned clusteredList will be a list of pairs [r,m] where
# r is the mean of the original cluster of m roots.
# -----------------------------------------------------------------------------
# (c) Robert M. Corless and Leili Rafiee Sevyeri 2020 (MIT Release licence)
#
# 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.
# -------------------------------------------------------------------------------

ClusterComplexRoots := proc( rootList::list,
                            {tolerance::numeric := Float(1,2-Digits)},
                            {MaxOrder::nonnegint := 3}, $ )
    local anglesequal,
          c, closeenough, closeNeighbours, clustered, combos, countem,
          DistanceMatrix, exactangle, finished_combos, Fuzz,
          i, j, k, m, max_radius, mean,
          n, neighbours, notInAcluster, nroots,
          phasedifference, phases, phaseorder,
          r, radius, trialcluster, trialindices, unseen_m;

    # Safety factors for distance computations.
    # They essentially replace the factor m!*B(r)/(D@@m)(f)(r)
    # which connects the root semimetric with the coefficient
    # metric; in this implementation, no account is taken of
    # the problem in this way (but it can be done in a post-
    # processing validation step).
    Fuzz := Array(2..MaxOrder);
    for m from 2 to MaxOrder do
      Fuzz[m] := evalf(2^(1/m));
    od;

    # Process small rootlists
    nroots := numelems( rootList );
    if nroots=0 then
      error "cannot cluster an empty list of roots";
    elif nroots=1 then
      # Return a cluster of 1
      return [[rootList[1],1]];
    end if;

    # nroots>=2 from now on
    # We will need distances between roots
    # -- to identify "close" roots
    # -- to ignore "far" roots
    DistanceMatrix := Matrix( nroots, nroots, 'shape'='symmetric' );
    for i to nroots do
      for j to nroots do
        DistanceMatrix[i,j] := evalf( abs( rootList[i] - rootList[j]) ) ;
      od;
    od;

    # We order neighbours from nearest to farthest
    neighbours := Array(1..nroots);
    for i to nroots do
      neighbours[i] := sort( [seq(DistanceMatrix[i,j],j=1..nroots)],
                              'output'='permutation');
      # neighbours[i] should be a list of integers such that
      # the nearest root to rootList[i] should be
      # rootList[neighbours[i][1]] (itself), the next nearest should be
      # rootList[neighbours[i][2]], and so on
    od;

    # Now begin the clustering.  We start with no clustered roots at all.
    clustered := NULL;
    # The notInAcluster roots are, at the start, all of them.
    # This states that we *assume* these are not in any m-cluster
    # for m > MaxOrder.
    notInAcluster := seq(k,k=1..nroots);
    for m from MaxOrder by -1 to 2 do
      # The sequence "notInAcluster" can shrink inside this loop
      # The list unseen_m tracks which roots we haven't yet
      # considered for any cluster at level m
      unseen_m := [notInAcluster]; # never any duplicates
      notInAcluster := NULL; # means known not to be part of any m-cluster
      # Assert nroots in clustered, notInAcluster,unseen_m +0(i)
      while numelems(unseen_m) >= m do
        # Checking for m-cycles over all elements of unseen_m
        i := unseen_m[1];
        # print("considering root", i, "or", rootList[i], "for cycles of length", m);
        unseen_m := unseen_m[2..numelems(unseen_m)]; # Assert nroots-1 + 1(i)
        # the root under consideration is rootList[i]
        # we try to find an m-cluster with this in it
        # No need to consider roots farther away than
        # twice the distance to the centre of any m-cluster
        # to this tolerance
        k := nroots;
        while k > 0 and
              DistanceMatrix[i,neighbours[i][k]] > 2*Fuzz[m]*tolerance^(1/m) do
          k := k-1;
        od;
        closeNeighbours := {seq(neighbours[i][j], j=2..k)}; # might be empty
        closeNeighbours := closeNeighbours intersect convert(unseen_m,set);
        closeNeighbours := convert(closeNeighbours,list);
        # examine all sets of m-1 close neighbours
        n := numelems(closeNeighbours); # <= k-1
        if n>=m-1 then
          # Now look at all possible choices of m-clusters
          # from the set of close neighbours
          combos := Iterator:-Combination( n, m-1 ); # at least 1
          finished_combos := false;
          for c in combos do
            # Iterator:-Combination indexes from 0 but lists index from 1
            trialindices := [i,seq(closeNeighbours[c[k]+1],k=1..m-1)];
            trialcluster := seq( rootList[k], k in trialindices );
            mean := evalf( add( r, r in [trialcluster] )/m );
            # Check the distances
            closeenough := true;
            max_radius := 0;
            for r in trialcluster do
              radius := abs(r-mean);
              max_radius := max( max_radius, radius );
              if radius > Fuzz[m]*tolerance^(1/m) then
                closeenough := false;
                break;
              end if;
            od;
            if closeenough then
              # Now check the angles
              # May be nonsense if max_radius=0
              anglesequal := true;
              if max_radius >= Float(1,1-Digits) then
                phases := [seq( argument( trialcluster[k]-mean), k=1..m) ];
                phaseorder := sort( phases, 'output'='permutation' );
                trialcluster := [seq( trialcluster[phaseorder[k]], k=1..m)];
                phases := [seq( phases[phaseorder[k]], k=1..m)];
                # Now check if all the angles are approximately equal
                exactangle := evalf(2*Pi/m);
                for k from 2 to m do
                  phasedifference := phases[k]-phases[k-1];
                  # Heuristic: what constitutes approximately equal angle?
                  if 0.75*exactangle > phasedifference or
                      phasedifference > 1.25*exactangle then
                    anglesequal := false;
                    break; # no need to consider this trial cluster further
                  end if;
                end do; # exit this loop on break for unequal angles
              end if;
              # closeenough is true, remember
              if anglesequal then
                # We have found a cluster of m roots
                clustered := [mean,m,trialcluster],clustered;
                # print([clustered], unseen_m);
                unseen_m := remove( k->member(k,convert(trialindices,set)),
                                    unseen_m );
                # print("removed",trialindices,"got",unseen_m);
                # Assert nroots + 1(i)
                break; # stop considering rootList[i]
              end if;
            # else Not close enough to be a cluster.
            end if;
            finished_combos := Rank(combos)=Number(combos);
            # Go look at the other choices of m-1 elements
          end do; # end for c in combos do; exit this loop via break for found cluster
          # stop looking at rootList[i] but keep looking for m-clusters
          if finished_combos then
            # Loop finished without finding a cluster
            notInAcluster := i,notInAcluster;
            # else We did find a cluster and have seen i
          end if;
          # print(m,"unseen_m",unseen_m);
        else
          # Not enough close neighbours to form an m-cycle
          notInAcluster := i, notInAcluster;
        end if;
      end do; # end while unseen_m has more than m elts do
      # Now look at clusters of smaller order
      for k in unseen_m do
        notInAcluster := k,notInAcluster;
      od;
      unseen_m := [];
    end do;
    # Everything left over is a 1-cluster
    notInAcluster := [notInAcluster];
    return [clustered, seq( [rootList[notInAcluster[k]],1],
                            k=1..numelems(notInAcluster)) ];
  end proc:

max_stars_repo_name: rcorless/ClusterComplexRootsMaple

max_issues_repo_name: rcorless/ClusterComplexRootsMaple

max_forks_repo_name: rcorless/ClusterComplexRootsMaple