#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <ctype.h>
#include <string.h>
#include <limits.h>

#include "ViennaRNA/fold_vars.h"
#include "ViennaRNA/energy_par.h"
#include "ViennaRNA/utils.h"
#include "ViennaRNA/constraints.h"
#include "ViennaRNA/gquad.h"
#include "ViennaRNA/structured_domains.h"
#include "ViennaRNA/unstructured_domains.h"
#include "ViennaRNA/exterior_loops.h"

#ifdef __GNUC__
# define INLINE inline
#else
# define INLINE
#endif

struct default_data {
  int                       *idx;
  char                      *mx;
  int                       cp;
  int                       *hc_up;
  void                      *hc_dat;
  vrna_callback_hc_evaluate *hc_f;
};


/*
 #################################
 # PRIVATE FUNCTION DECLARATIONS #
 #################################
 */

PRIVATE FLT_OR_DBL
exp_E_ext_fast(vrna_fold_compound_t *vc,
               int                  i,
               int                  j,
               vrna_mx_pf_aux_el_t  *aux_mx);


PRIVATE FLT_OR_DBL
exp_E_ext_fast_comparative(vrna_fold_compound_t *vc,
                           int                  i,
                           int                  j,
                           vrna_mx_pf_aux_el_t  *aux_mx);


PRIVATE char
hc_default(int  i,
           int  j,
           int  k,
           int  l,
           char d,
           void *data);


PRIVATE char
hc_default_user(int   i,
                int   j,
                int   k,
                int   l,
                char  d,
                void  *data);


/*
 #################################
 # BEGIN OF FUNCTION DEFINITIONS #
 #################################
 */
PUBLIC int
E_ext_loop(int                  i,
           int                  j,
           vrna_fold_compound_t *vc)
{
  char                      *ptype, *hard_constraints;
  short                     *S;
  int                       ij, en, e, type, cp, *idx;
  vrna_param_t              *P;
  vrna_md_t                 *md;
  vrna_sc_t                 *sc;
  vrna_callback_hc_evaluate *evaluate;
  struct default_data       hc_dat_local;

  cp                = vc->cutpoint;
  S                 = vc->sequence_encoding;
  idx               = vc->jindx;
  ptype             = vc->ptype;
  P                 = vc->params;
  md                = &(P->model_details);
  hard_constraints  = vc->hc->matrix;
  sc                = vc->sc;

  hc_dat_local.idx  = idx;
  hc_dat_local.mx   = hard_constraints;
  hc_dat_local.cp   = cp;

  if (vc->hc->f) {
    evaluate            = &hc_default_user;
    hc_dat_local.hc_f   = vc->hc->f;
    hc_dat_local.hc_dat = vc->hc->data;
  } else {
    evaluate = &hc_default;
  }


  e     = INF;
  ij    = idx[j] + i;
  type  = ptype[ij];

  if ((cp < 0) || (((i) >= cp) || ((j) < cp))) {
    /* regular exterior loop */
    if (evaluate(i, j, i, j, VRNA_DECOMP_EXT_STEM, &hc_dat_local)) {
      if (type == 0)
        type = 7;

      switch (md->dangles) {
        case 2:
          e = E_ExtLoop(type, S[i - 1], S[j + 1], P);
          break;

        case 0:
        /* fall through */

        default:
          e = E_ExtLoop(type, -1, -1, P);
          break;
      }
      if (sc)
        if (sc->f)
          e += sc->f(i, j, i, j, VRNA_DECOMP_EXT_STEM, sc->data);
    }

    if (md->dangles % 2) {
      ij = idx[j - 1] + i;
      if (evaluate(i, j, i, j - 1, VRNA_DECOMP_EXT_STEM, &hc_dat_local)) {
        type = vc->ptype[ij];

        if (type == 0)
          type = 7;

        en = E_ExtLoop(type, -1, S[j], P);
        if (sc)
          if (sc->f)
            en += sc->f(i, j, i, j - 1, VRNA_DECOMP_EXT_STEM, sc->data);
        e = MIN2(e, en);
      }

      ij = idx[j] + i + 1;
      if (evaluate(i, j, i + 1, j, VRNA_DECOMP_EXT_STEM, &hc_dat_local)) {
        type = vc->ptype[ij];

        if (type == 0)
          type = 7;

        en = E_ExtLoop(type, S[i], -1, P);
        if (sc)
          if (sc->f)
            en += sc->f(i, j, i + 1, j, VRNA_DECOMP_EXT_STEM, sc->data);
        e = MIN2(e, en);
      }
    }
  }

  return e;
}


PUBLIC void
E_ext_loop_5(vrna_fold_compound_t *vc)
{
  char                      *ptype, *hc;
  short                     *S;
  int                       en, i, j, ij, type, length, *indx, *hc_up, *f5, *c, dangle_model,
                            *ggg, with_gquad, turn, k, u, with_ud;
  vrna_sc_t                 *sc;
  vrna_param_t              *P;
  vrna_ud_t                 *domains_up;
  vrna_callback_hc_evaluate *evaluate;
  struct default_data       hc_dat_local;

  length        = (int)vc->length;
  ptype         = vc->ptype;
  S             = vc->sequence_encoding;
  indx          = vc->jindx;
  hc            = vc->hc->matrix;
  hc_up         = vc->hc->up_ext;
  sc            = vc->sc;
  f5            = vc->matrices->f5;
  c             = vc->matrices->c;
  P             = vc->params;
  dangle_model  = P->model_details.dangles;
  ggg           = vc->matrices->ggg;
  with_gquad    = P->model_details.gquad;
  turn          = P->model_details.min_loop_size;
  domains_up    = vc->domains_up;
  with_ud       = (domains_up && domains_up->energy_cb) ? 1 : 0;

  hc_dat_local.idx    = indx;
  hc_dat_local.mx     = hc;
  hc_dat_local.hc_up  = hc_up;
  hc_dat_local.cp     = vc->cutpoint;

  if (vc->hc->f) {
    evaluate            = &hc_default_user;
    hc_dat_local.hc_f   = vc->hc->f;
    hc_dat_local.hc_dat = vc->hc->data;
  } else {
    evaluate = &hc_default;
  }

  f5[0] = 0;
  for (i = 1; i <= turn + 1; i++) {
    if (f5[i - 1] != INF) {
      if (evaluate(1, i, 1, i - 1, VRNA_DECOMP_EXT_EXT, &hc_dat_local)) {
        f5[i] = f5[i - 1];
        if (sc) {
          if (sc->energy_up)
            f5[i] += sc->energy_up[i][1];
          if (sc->f)
            f5[i] += sc->f(1, i, 1, i - 1, VRNA_DECOMP_EXT_EXT, sc->data);
        }
      } else {
        f5[i] = INF;
      }
    } else {
      f5[i] = INF;
    }
  }

  if (with_ud) {
    /* do we include ligand binding? */
    /*  construct all possible combinations of
     *  f[i-1] + L[i,j] with j <= turn + 1
     */
    for (i = 1; i <= turn + 1; i++) {
      if (f5[i - 1] != INF) {
        for (k = 0; k < domains_up->uniq_motif_count; k++) {
          u = domains_up->uniq_motif_size[k];
          j = i + u - 1;
          if (j <= turn + 1) {
            if (evaluate(1, j, i - 1, i, VRNA_DECOMP_EXT_EXT, &hc_dat_local)) {
              en = f5[i - 1]
                   + domains_up->energy_cb(vc,
                                           i, j,
                                           VRNA_UNSTRUCTURED_DOMAIN_EXT_LOOP,
                                           domains_up->data);
              if (sc) {
                if (sc->energy_up)
                  en += sc->energy_up[i][u];
                if (sc->f)
                  en += sc->f(1, j, 1, j - u, VRNA_DECOMP_EXT_EXT, sc->data);
              }
              f5[j] = MIN2(f5[j], en);
            }
          }
        }
      }
    }
  }

  /* duplicated code may be faster than conditions inside loop ;) */
  switch (dangle_model) {
    /* dont use dangling end and mismatch contributions at all */
    case 0:
      for (j = turn + 2; j <= length; j++) {
        /* initialize with INF */
        f5[j] = INF;

        /* check for 3' extension with one unpaired nucleotide */
        if (f5[j - 1] != INF) {
          if (evaluate(1, j, 1, j - 1, VRNA_DECOMP_EXT_EXT, &hc_dat_local)) {
            f5[j] = f5[j - 1];
            if (sc) {
              if (sc->energy_up)
                f5[j] += sc->energy_up[j][1];
              if (sc->f)
                f5[j] += sc->f(1, j, 1, j - 1, VRNA_DECOMP_EXT_EXT, sc->data);
            }
          }
        }

        if (with_ud) {
          for (k = 0; k < domains_up->uniq_motif_count; k++) {
            u = domains_up->uniq_motif_size[k];
            if ((j - u >= 0) && (f5[j - u] != INF)) {
              if (evaluate(1, j, 1, j - u, VRNA_DECOMP_EXT_EXT, &hc_dat_local)) {
                en = f5[j - u]
                     + domains_up->energy_cb(vc,
                                             j - u + 1, j,
                                             VRNA_UNSTRUCTURED_DOMAIN_EXT_LOOP | VRNA_UNSTRUCTURED_DOMAIN_MOTIF,
                                             domains_up->data);
                if (sc) {
                  if (sc->energy_up)
                    en += sc->energy_up[j - u + 1][u];
                  if (sc->f)
                    en += sc->f(1, j, 1, j - u, VRNA_DECOMP_EXT_EXT, sc->data);
                }

                f5[j] = MIN2(f5[j], en);
              }
            }
          }
        }

        /* check for possible stems branching off the exterior loop */
        if (sc && sc->f) {
          for (i = j - turn - 1; i > 1; i--) {
            if (f5[i - 1] != INF) {
              ij = indx[j] + i;

              if (with_gquad)
                f5[j] = MIN2(f5[j], f5[i - 1] + ggg[ij]);

              if (c[ij] != INF) {
                if (evaluate(1, j, i - 1, i, VRNA_DECOMP_EXT_EXT_STEM, &hc_dat_local)) {
                  type = ptype[ij];

                  if (type == 0)
                    type = 7;

                  en    = f5[i - 1] + c[ij] + E_ExtLoop(type, -1, -1, P);
                  en    += sc->f(1, j, i - 1, i, VRNA_DECOMP_EXT_EXT_STEM, sc->data);
                  f5[j] = MIN2(f5[j], en);
                }
              }
            }
          }
        } else {
          for (i = j - turn - 1; i > 1; i--) {
            if (f5[i - 1] != INF) {
              ij = indx[j] + i;

              if (with_gquad)
                f5[j] = MIN2(f5[j], f5[i - 1] + ggg[ij]);

              if (c[ij] != INF) {
                if (evaluate(1, j, i - 1, i, VRNA_DECOMP_EXT_EXT_STEM, &hc_dat_local)) {
                  type = ptype[ij];

                  if (type == 0)
                    type = 7;

                  en    = f5[i - 1] + c[ij] + E_ExtLoop(type, -1, -1, P);
                  f5[j] = MIN2(f5[j], en);
                }
              }
            }
          }
        }
        ij = indx[j] + 1;

        if (with_gquad)
          f5[j] = MIN2(f5[j], ggg[ij]);

        if (c[ij] != INF) {
          if (evaluate(1, j, 1, j, VRNA_DECOMP_EXT_STEM, &hc_dat_local)) {
            type = ptype[ij];

            if (type == 0)
              type = 7;

            en = c[ij] + E_ExtLoop(type, -1, -1, P);
            if (sc)
              if (sc->f)
                en += sc->f(1, j, 1, j, VRNA_DECOMP_EXT_STEM, sc->data);
            f5[j] = MIN2(f5[j], en);
          }
        }
      }
      break;

    /* always use dangles on both sides */
    case 2:
      for (j = turn + 2; j < length; j++) {
        f5[j] = INF;

        if (f5[j - 1] != INF) {
          if (evaluate(1, j, 1, j - 1, VRNA_DECOMP_EXT_EXT, &hc_dat_local)) {
            f5[j] = f5[j - 1];
            if (sc) {
              if (sc->energy_up)
                f5[j] += sc->energy_up[j][1];
              if (sc->f)
                f5[j] += sc->f(1, j, 1, j - 1, VRNA_DECOMP_EXT_EXT, sc->data);
            }
          }
        }

        if (with_ud) {
          for (k = 0; k < domains_up->uniq_motif_count; k++) {
            u = domains_up->uniq_motif_size[k];
            if ((j - u >= 0) && (f5[j - u] != INF)) {
              if (evaluate(1, j, 1, j - u, VRNA_DECOMP_EXT_EXT, &hc_dat_local)) {
                en = f5[j - u]
                     + domains_up->energy_cb(vc,
                                             j - u + 1, j,
                                             VRNA_UNSTRUCTURED_DOMAIN_EXT_LOOP | VRNA_UNSTRUCTURED_DOMAIN_MOTIF,
                                             domains_up->data);
                if (sc) {
                  if (sc->energy_up)
                    en += sc->energy_up[j - u + 1][u];
                  if (sc->f)
                    en += sc->f(1, j, 1, j - u, VRNA_DECOMP_EXT_EXT, sc->data);
                }

                f5[j] = MIN2(f5[j], en);
              }
            }
          }
        }

        if (sc && sc->f) {
          for (i = j - turn - 1; i > 1; i--) {
            if (f5[i - 1] != INF) {
              ij = indx[j] + i;

              if (with_gquad)
                f5[j] = MIN2(f5[j], f5[i - 1] + ggg[ij]);

              if (c[ij] != INF) {
                if (evaluate(1, j, i - 1, i, VRNA_DECOMP_EXT_EXT_STEM, &hc_dat_local)) {
                  type = ptype[ij];

                  if (type == 0)
                    type = 7;

                  en    = f5[i - 1] + c[ij] + E_ExtLoop(type, S[i - 1], S[j + 1], P);
                  en    += sc->f(1, j, i - 1, i, VRNA_DECOMP_EXT_EXT_STEM, sc->data);
                  f5[j] = MIN2(f5[j], en);
                }
              }
            }
          }
        } else {
          for (i = j - turn - 1; i > 1; i--) {
            if (f5[i - 1] != INF) {
              ij = indx[j] + i;

              if (with_gquad)
                f5[j] = MIN2(f5[j], f5[i - 1] + ggg[ij]);

              if (c[ij] != INF) {
                if (evaluate(1, j, i - 1, i, VRNA_DECOMP_EXT_EXT_STEM, &hc_dat_local)) {
                  type = ptype[ij];

                  if (type == 0)
                    type = 7;

                  en    = f5[i - 1] + c[ij] + E_ExtLoop(type, S[i - 1], S[j + 1], P);
                  f5[j] = MIN2(f5[j], en);
                }
              }
            }
          }
        }
        ij = indx[j] + 1;

        if (with_gquad)
          f5[j] = MIN2(f5[j], ggg[ij]);

        if (c[ij] != INF) {
          if (evaluate(1, j, 1, j, VRNA_DECOMP_EXT_STEM, &hc_dat_local)) {
            type = ptype[ij];

            if (type == 0)
              type = 7;

            en = c[ij] + E_ExtLoop(type, -1, S[j + 1], P);
            if (sc)
              if (sc->f)
                en += sc->f(1, j, 1, j, VRNA_DECOMP_EXT_STEM, sc->data);
            f5[j] = MIN2(f5[j], en);
          }
        }
      }

      f5[length] = INF;
      if (f5[length - 1] != INF) {
        if (evaluate(1, length, 1, length - 1, VRNA_DECOMP_EXT_EXT, &hc_dat_local)) {
          f5[length] = f5[length - 1];
          if (sc) {
            if (sc->energy_up)
              f5[length] += sc->energy_up[length][1];
            if (sc->f)
              f5[length] += sc->f(1, length, 1, length - 1, VRNA_DECOMP_EXT_EXT, sc->data);
          }
        }
      }

      if (with_ud) {
        for (k = 0; k < domains_up->uniq_motif_count; k++) {
          u = domains_up->uniq_motif_size[k];
          if ((length - u >= 0) && (f5[length - u] != INF)) {
            if (evaluate(1, length, 1, length - u, VRNA_DECOMP_EXT_EXT, &hc_dat_local)) {
              en = f5[length - u]
                   + domains_up->energy_cb(vc,
                                           length - u + 1, length,
                                           VRNA_UNSTRUCTURED_DOMAIN_EXT_LOOP | VRNA_UNSTRUCTURED_DOMAIN_MOTIF,
                                           domains_up->data);
              if (sc) {
                if (sc->energy_up)
                  en += sc->energy_up[length - u + 1][u];
                if (sc->f)
                  en += sc->f(1, length, 1, length - u, VRNA_DECOMP_EXT_EXT, sc->data);
              }

              f5[length] = MIN2(f5[length], en);
            }
          }
        }
      }

      if (sc && sc->f) {
        for (i = length - turn - 1; i > 1; i--) {
          if (f5[i - 1] != INF) {
            ij = indx[length] + i;

            if (with_gquad)
              f5[length] = MIN2(f5[length], f5[i - 1] + ggg[ij]);

            if (c[ij] != INF) {
              if (evaluate(1, length, i - 1, i, VRNA_DECOMP_EXT_EXT_STEM, &hc_dat_local)) {
                type = ptype[ij];

                if (type == 0)
                  type = 7;

                en          = f5[i - 1] + c[ij] + E_ExtLoop(type, S[i - 1], -1, P);
                en          += sc->f(1, length, i - 1, i, VRNA_DECOMP_EXT_EXT_STEM, sc->data);
                f5[length]  = MIN2(f5[length], en);
              }
            }
          }
        }
      } else {
        for (i = length - turn - 1; i > 1; i--) {
          if (f5[i - 1] != INF) {
            ij = indx[length] + i;

            if (with_gquad)
              f5[length] = MIN2(f5[length], f5[i - 1] + ggg[ij]);

            if (c[ij] != INF) {
              if (evaluate(1, length, i - 1, i, VRNA_DECOMP_EXT_EXT_STEM, &hc_dat_local)) {
                type = ptype[ij];

                if (type == 0)
                  type = 7;

                en          = f5[i - 1] + c[ij] + E_ExtLoop(type, S[i - 1], -1, P);
                f5[length]  = MIN2(f5[length], en);
              }
            }
          }
        }
      }
      ij = indx[length] + 1;

      if (with_gquad)
        f5[length] = MIN2(f5[length], ggg[ij]);

      if (c[ij] != INF) {
        if (evaluate(1, length, 1, length, VRNA_DECOMP_EXT_STEM, &hc_dat_local)) {
          type = ptype[ij];

          if (type == 0)
            type = 7;

          en = c[ij] + E_ExtLoop(type, -1, -1, P);
          if (sc)
            if (sc->f)
              en += sc->f(1, length, 1, length, VRNA_DECOMP_EXT_STEM, sc->data);
          f5[length] = MIN2(f5[length], en);
        }
      }
      break;

    /* normal dangles, aka dangle_model = 1 || 3 */
    default:
      for (j = turn + 2; j <= length; j++) {
        f5[j] = INF;
        if (f5[j - 1] != INF) {
          if (evaluate(1, j, 1, j - 1, VRNA_DECOMP_EXT_EXT, &hc_dat_local)) {
            f5[j] = f5[j - 1];
            if (sc) {
              if (sc->energy_up)
                f5[j] += sc->energy_up[j][1];
              if (sc->f)
                f5[j] += sc->f(1, j, 1, j - 1, VRNA_DECOMP_EXT_EXT, sc->data);
            }
          }
        }

        if (with_ud) {
          for (k = 0; k < domains_up->uniq_motif_count; k++) {
            u = domains_up->uniq_motif_size[k];
            if ((j - u >= 0) && (f5[j - u] != INF)) {
              if (evaluate(1, j, 1, j - u, VRNA_DECOMP_EXT_EXT, &hc_dat_local)) {
                en = f5[j - u]
                     + domains_up->energy_cb(vc,
                                             j - u + 1, j,
                                             VRNA_UNSTRUCTURED_DOMAIN_EXT_LOOP | VRNA_UNSTRUCTURED_DOMAIN_MOTIF,
                                             domains_up->data);
                if (sc) {
                  if (sc->energy_up)
                    en += sc->energy_up[j - u + 1][u];
                  if (sc->f)
                    en += sc->f(1, j, 1, j - u, VRNA_DECOMP_EXT_EXT, sc->data);
                }

                f5[j] = MIN2(f5[j], en);
              }
            }
          }
        }

        for (i = j - turn - 1; i > 1; i--) {
          ij = indx[j] + i;
          if (f5[i - 1] != INF) {
            if (with_gquad)
              f5[j] = MIN2(f5[j], f5[i - 1] + ggg[ij]);

            if (c[ij] != INF) {
              if (evaluate(1, j, i - 1, i, VRNA_DECOMP_EXT_EXT_STEM, &hc_dat_local)) {
                type = ptype[ij];

                if (type == 0)
                  type = 7;

                en = f5[i - 1] + c[ij] + E_ExtLoop(type, -1, -1, P);
                if (sc)
                  if (sc->f)
                    en += sc->f(1, j, i - 1, i, VRNA_DECOMP_EXT_EXT_STEM, sc->data);
                f5[j] = MIN2(f5[j], en);
              }
            }
          }

          if ((f5[i - 2] != INF) && c[ij] != INF) {
            if (evaluate(1, j, i - 2, i, VRNA_DECOMP_EXT_EXT_STEM, &hc_dat_local)) {
              type = ptype[ij];

              if (type == 0)
                type = 7;

              en = f5[i - 2] + c[ij] + E_ExtLoop(type, S[i - 1], -1, P);

              if (sc) {
                if (sc->energy_up)
                  en += sc->energy_up[i - 1][1];
                if (sc->f)
                  en += sc->f(1, j, i - 2, i, VRNA_DECOMP_EXT_EXT_STEM, sc->data);
              }
              f5[j] = MIN2(f5[j], en);
            }
          }

          ij = indx[j - 1] + i;
          if (c[ij] != INF) {
            if (f5[i - 1] != INF) {
              if (evaluate(1, j, i - 1, i, VRNA_DECOMP_EXT_EXT_STEM1, &hc_dat_local)) {
                type = ptype[ij];

                if (type == 0)
                  type = 7;

                en = f5[i - 1] + c[ij] + E_ExtLoop(type, -1, S[j], P);

                if (sc) {
                  if (sc->energy_up)
                    en += sc->energy_up[j][1];
                  if (sc->f)
                    en += sc->f(1, j, i - 1, i, VRNA_DECOMP_EXT_EXT_STEM1, sc->data);
                }
                f5[j] = MIN2(f5[j], en);
              }
            }

            if (f5[i - 2] != INF) {
              if (evaluate(1, j, i - 2, i, VRNA_DECOMP_EXT_EXT_STEM1, &hc_dat_local)) {
                type = ptype[ij];

                if (type == 0)
                  type = 7;

                en = f5[i - 2] + c[ij] + E_ExtLoop(type, S[i - 1], S[j], P);

                if (sc) {
                  if (sc->energy_up)
                    en += sc->energy_up[i - 1][1] + sc->energy_up[j][1];
                  if (sc->f)
                    en += sc->f(1, j, i - 2, i, VRNA_DECOMP_EXT_EXT_STEM1, sc->data);
                }
                f5[j] = MIN2(f5[j], en);
              }
            }
          }
        }

        ij = indx[j] + 1;

        if (with_gquad)
          f5[j] = MIN2(f5[j], ggg[ij]);

        if (c[ij] != INF) {
          if (evaluate(1, j, 1, j, VRNA_DECOMP_EXT_STEM, &hc_dat_local)) {
            type = ptype[ij];

            if (type == 0)
              type = 7;

            en = c[ij] + E_ExtLoop(type, -1, -1, P);
            if (sc)
              if (sc->f)
                en += sc->f(1, j, 1, j, VRNA_DECOMP_EXT_STEM, sc->data);
            f5[j] = MIN2(f5[j], en);
          }
        }
        ij = indx[j - 1] + 1;
        if (c[ij] != INF) {
          if (evaluate(1, j, 1, j - 1, VRNA_DECOMP_EXT_STEM, &hc_dat_local)) {
            type = ptype[ij];

            if (type == 0)
              type = 7;

            en = c[ij] + E_ExtLoop(type, -1, S[j], P);

            if (sc) {
              if (sc->energy_up)
                en += sc->energy_up[j][1];
              if (sc->f)
                en += sc->f(1, j, 1, j - 1, VRNA_DECOMP_EXT_STEM, sc->data);
            }
            f5[j] = MIN2(f5[j], en);
          }
        }
      }         /* end for j... */
      break;
  }
}


PUBLIC int
E_Stem(int          type,
       int          si1,
       int          sj1,
       int          extLoop,
       vrna_param_t *P)
{
  int energy  = 0;
  int d5      = (si1 >= 0) ? P->dangle5[type][si1] : 0;
  int d3      = (sj1 >= 0) ? P->dangle3[type][sj1] : 0;

  if (type > 2)
    energy += P->TerminalAU;

  if (si1 >= 0 && sj1 >= 0)
    energy += (extLoop) ? P->mismatchExt[type][si1][sj1] : P->mismatchM[type][si1][sj1];
  else
    energy += d5 + d3;

  if (!extLoop)
    energy += P->MLintern[type];

  return energy;
}


PUBLIC int
E_ExtLoop(int           type,
          int           si1,
          int           sj1,
          vrna_param_t  *P)
{
  int energy = 0;

  if (si1 >= 0 && sj1 >= 0)
    energy += P->mismatchExt[type][si1][sj1];
  else if (si1 >= 0)
    energy += P->dangle5[type][si1];
  else if (sj1 >= 0)
    energy += P->dangle3[type][sj1];

  if (type > 2)
    energy += P->TerminalAU;

  return energy;
}


PUBLIC FLT_OR_DBL
exp_E_Stem(int              type,
           int              si1,
           int              sj1,
           int              extLoop,
           vrna_exp_param_t *P)
{
  double  energy  = 1.0;
  double  d5      = (si1 >= 0) ? P->expdangle5[type][si1] : 1.;
  double  d3      = (sj1 >= 0) ? P->expdangle3[type][sj1] : 1.;

  if (si1 >= 0 && sj1 >= 0)
    energy = (extLoop) ? P->expmismatchExt[type][si1][sj1] : P->expmismatchM[type][si1][sj1];
  else
    energy = d5 * d3;

  if (type > 2)
    energy *= P->expTermAU;

  if (!extLoop)
    energy *= P->expMLintern[type];

  return (FLT_OR_DBL)energy;
}


PUBLIC FLT_OR_DBL
exp_E_ExtLoop(int               type,
              int               si1,
              int               sj1,
              vrna_exp_param_t  *P)
{
  double energy = 1.0;

  if (si1 >= 0 && sj1 >= 0)
    energy = P->expmismatchExt[type][si1][sj1];
  else if (si1 >= 0)
    energy = P->expdangle5[type][si1];
  else if (sj1 >= 0)
    energy = P->expdangle3[type][sj1];

  if (type > 2)
    energy *= P->expTermAU;

  return (FLT_OR_DBL)energy;
}


PUBLIC int
vrna_BT_ext_loop_f5(vrna_fold_compound_t  *vc,
                    int                   *k,
                    int                   *i,
                    int                   *j,
                    vrna_bp_stack_t       *bp_stack,
                    int                   *stack_count)
{
  unsigned char             type;
  char                      *ptype;
  short                     mm5, mm3, *S1;
  unsigned int              *sn;
  int                       length, fij, fi, jj, u, en, e, *my_f5, *my_c, *my_ggg, *idx,
                            dangle_model, turn, with_gquad, cnt, ii, with_ud;
  vrna_param_t              *P;
  vrna_md_t                 *md;
  vrna_hc_t                 *hc;
  vrna_sc_t                 *sc;
  vrna_ud_t                 *domains_up;
  vrna_callback_hc_evaluate *evaluate;
  struct default_data       hc_dat_local;

  length        = vc->length;
  P             = vc->params;
  md            = &(P->model_details);
  sn            = vc->strand_number;
  hc            = vc->hc;
  sc            = vc->sc;
  my_f5         = vc->matrices->f5;
  my_c          = vc->matrices->c;
  my_ggg        = vc->matrices->ggg;
  domains_up    = vc->domains_up;
  idx           = vc->jindx;
  ptype         = vc->ptype;
  S1            = vc->sequence_encoding;
  dangle_model  = md->dangles;
  turn          = md->min_loop_size;
  with_gquad    = md->gquad;
  with_ud       = (domains_up && domains_up->energy_cb) ? 1 : 0;

  hc_dat_local.idx    = idx;
  hc_dat_local.mx     = hc->matrix;
  hc_dat_local.hc_up  = hc->up_ext;
  hc_dat_local.cp     = vc->cutpoint;

  if (vc->hc->f) {
    evaluate            = &hc_default_user;
    hc_dat_local.hc_f   = hc->f;
    hc_dat_local.hc_dat = hc->data;
  } else {
    evaluate = &hc_default;
  }

  jj = *k;

  /* nibble off unpaired 3' stretches harboring bound ligands (interspersed with unpaired nucleotides) */
  if (with_ud) {
    do {
      fij = my_f5[jj];
      fi  = INF;

      /* try nibble off one unpaired nucleotide first */
      if (evaluate(1, jj, 1, jj - 1, VRNA_DECOMP_EXT_EXT, &hc_dat_local)) {
        fi = my_f5[jj - 1];

        if (sc) {
          if (sc->energy_up)
            fi += sc->energy_up[jj][1];
          if (sc->f)
            fi += sc->f(1, jj, 1, jj - 1, VRNA_DECOMP_EXT_EXT, sc->data);
        }

        if (jj == 1) {
          /* no more pairs */
          *i  = *j = -1;
          *k  = 0;
          return 1;
        }

        if (fij == fi) {
          jj--;
          continue;
        }
      }

      /* next, try nibble off a ligand */
      for (cnt = 0; cnt < domains_up->uniq_motif_count; cnt++) {
        u   = domains_up->uniq_motif_size[cnt];
        ii  = jj - u + 1;
        if ((ii > 0) && evaluate(1, jj, 1, jj - u, VRNA_DECOMP_EXT_EXT, &hc_dat_local)) {
          en = domains_up->energy_cb(vc,
                                     ii, jj,
                                     VRNA_UNSTRUCTURED_DOMAIN_EXT_LOOP | VRNA_UNSTRUCTURED_DOMAIN_MOTIF,
                                     domains_up->data);
          if (sc) {
            if (sc->energy_up)
              en += sc->energy_up[ii][u];
            if (sc->f)
              en += sc->f(1, jj, 1, jj - u, VRNA_DECOMP_EXT_EXT, sc->data);
          }
          fi  = my_f5[ii - 1];
          fi  += en;

          if (fij == fi) {
            /* skip remaining motifs after first hit */
            jj = ii - 1;
            break;
          }
        }
      }

      if (jj == 0) {
        /* no more pairs */
        *i  = *j = -1;
        *k  = 0;
        return 1;
      }
    } while (fij == fi);
  } else {
    /* nibble off unpaired 3' bases */
    do {
      fij = my_f5[jj];
      fi  = INF;

      if (evaluate(1, jj, 1, jj - 1, VRNA_DECOMP_EXT_EXT, &hc_dat_local)) {
        fi = my_f5[jj - 1];

        if (sc) {
          if (sc->energy_up)
            fi += sc->energy_up[jj][1];
          if (sc->f)
            fi += sc->f(1, jj, 1, jj - 1, VRNA_DECOMP_EXT_EXT, sc->data);
        }
      }

      if (--jj == 0)
        break;
    } while (fij == fi);
    jj++;
  }

  if (jj < turn + 2) {
    /* no more pairs */
    *i  = *j = -1;
    *k  = 0;
    return 1;
  }

  /* must have found a decomposition */
  switch (dangle_model) {
    case 0:   /* j is paired. Find pairing partner */
      for (u = jj - turn - 1; u >= 1; u--) {
        if (with_gquad) {
          if (fij == my_f5[u - 1] + my_ggg[idx[jj] + u]) {
            *i  = *j = -1;
            *k  = u - 1;
            vrna_BT_gquad_mfe(vc, u, jj, bp_stack, stack_count);
            return 1;
          }
        }

        if (evaluate(1, jj, u - 1, u, VRNA_DECOMP_EXT_EXT_STEM, &hc_dat_local)) {
          type = (unsigned char)ptype[idx[jj] + u];

          if (type == 0)
            type = 7;

          en = my_c[idx[jj] + u];
          if (sc)
            if (sc->f)
              en += sc->f(1, jj, u - 1, u, VRNA_DECOMP_EXT_EXT_STEM, sc->data);
          if (sn[jj] != sn[u])
            en += P->DuplexInit;
          if (fij == E_ExtLoop(type, -1, -1, P) + en + my_f5[u - 1]) {
            *i                            = u;
            *j                            = jj;
            *k                            = u - 1;
            bp_stack[++(*stack_count)].i  = u;
            bp_stack[(*stack_count)].j    = jj;
            return 1;
          }
        }
      }
      break;

    case 2:
      mm3 = ((jj < length) && (sn[jj + 1] == sn[jj])) ? S1[jj + 1] : -1;
      for (u = jj - turn - 1; u >= 1; u--) {
        if (with_gquad) {
          if (fij == my_f5[u - 1] + my_ggg[idx[jj] + u]) {
            *i  = *j = -1;
            *k  = u - 1;
            vrna_BT_gquad_mfe(vc, u, jj, bp_stack, stack_count);
            return 1;
          }
        }

        if (evaluate(1, jj, u - 1, u, VRNA_DECOMP_EXT_EXT_STEM, &hc_dat_local)) {
          mm5   = ((u > 1) && (sn[u] == sn[u - 1])) ? S1[u - 1] : -1;
          type  = (unsigned char)ptype[idx[jj] + u];

          if (type == 0)
            type = 7;

          en = my_c[idx[jj] + u];
          if (sc)
            if (sc->f)
              en += sc->f(1, jj, u - 1, u, VRNA_DECOMP_EXT_EXT_STEM, sc->data);
          if (sn[jj] != sn[u])
            en += P->DuplexInit;
          if (fij == E_ExtLoop(type, mm5, mm3, P) + en + my_f5[u - 1]) {
            *i                            = u;
            *j                            = jj;
            *k                            = u - 1;
            bp_stack[++(*stack_count)].i  = u;
            bp_stack[(*stack_count)].j    = jj;
            return 1;
          }
        }
      }
      break;

    default:
      if (with_gquad) {
        if (fij == my_ggg[idx[jj] + 1]) {
          *i  = *j = -1;
          *k  = 0;
          vrna_BT_gquad_mfe(vc, 1, jj, bp_stack, stack_count);
          return 1;
        }
      }

      if (evaluate(1, jj, 1, jj, VRNA_DECOMP_EXT_STEM, &hc_dat_local)) {
        type = (unsigned char)ptype[idx[jj] + 1];

        if (type == 0)
          type = 7;

        en = my_c[idx[jj] + 1];
        if (sc)
          if (sc->f)
            en += sc->f(1, jj, 1, jj, VRNA_DECOMP_EXT_STEM, sc->data);
        if (sn[jj] != sn[1])
          en += P->DuplexInit;
        if (fij == en + E_ExtLoop(type, -1, -1, P)) {
          *i                            = 1;
          *j                            = jj;
          *k                            = 0;
          bp_stack[++(*stack_count)].i  = 1;
          bp_stack[(*stack_count)].j    = jj;
          return 1;
        }
      }

      if (evaluate(1, jj, 1, jj - 1, VRNA_DECOMP_EXT_STEM, &hc_dat_local)) {
        if (sn[jj] == sn[jj - 1]) {
          mm3   = S1[jj];
          type  = (unsigned char)ptype[idx[jj - 1] + 1];

          if (type == 0)
            type = 7;

          en = my_c[idx[jj - 1] + 1];
          if (sc) {
            if (sc->energy_up)
              en += sc->energy_up[jj][1];
            if (sc->f)
              en += sc->f(1, jj, 1, jj - 1, VRNA_DECOMP_EXT_STEM, sc->data);
          }
          if (sn[jj - 1] != sn[1])
            en += P->DuplexInit;

          if (fij == en + E_ExtLoop(type, -1, mm3, P)) {
            *i                            = 1;
            *j                            = jj - 1;
            *k                            = 0;
            bp_stack[++(*stack_count)].i  = 1;
            bp_stack[(*stack_count)].j    = jj - 1;
            return 1;
          }
        }
      }

      for (u = jj - turn - 1; u > 1; u--) {
        if (with_gquad) {
          if (fij == my_f5[u - 1] + my_ggg[idx[jj] + u]) {
            *i  = *j = -1;
            *k  = u - 1;
            vrna_BT_gquad_mfe(vc, u, jj, bp_stack, stack_count);
            return 1;
          }
        }

        type = (unsigned char)ptype[idx[jj] + u];
        if (type == 0)
          type = 7;

        en = my_c[idx[jj] + u];
        if (sn[jj] != sn[u])
          en += P->DuplexInit;

        if (evaluate(1, jj, u - 1, u, VRNA_DECOMP_EXT_EXT_STEM, &hc_dat_local)) {
          e = my_f5[u - 1] + en + E_ExtLoop(type, -1, -1, P);
          if (sc)
            if (sc->f)
              e += sc->f(1, jj, u - 1, u, VRNA_DECOMP_EXT_EXT_STEM, sc->data);
          if (fij == e) {
            *i                            = u;
            *j                            = jj;
            *k                            = u - 1;
            bp_stack[++(*stack_count)].i  = u;
            bp_stack[(*stack_count)].j    = jj;
            return 1;
          }
        }

        if (evaluate(1, jj, u - 2, u, VRNA_DECOMP_EXT_EXT_STEM, &hc_dat_local)) {
          if (sn[u] == sn[u - 1]) {
            mm5 = S1[u - 1];
            e   = my_f5[u - 2] + en + E_ExtLoop(type, mm5, -1, P);
            if (sc) {
              if (sc->energy_up)
                e += sc->energy_up[u - 1][1];
              if (sc->f)
                e += sc->f(1, jj, u - 2, u, VRNA_DECOMP_EXT_EXT_STEM, sc->data);
            }
            if (fij == e) {
              *i                            = u;
              *j                            = jj;
              *k                            = u - 2;
              bp_stack[++(*stack_count)].i  = u;
              bp_stack[(*stack_count)].j    = jj;
              return 1;
            }
          }
        }

        type = (unsigned char)ptype[idx[jj - 1] + u];
        if (type == 0)
          type = 7;

        en = my_c[idx[jj - 1] + u];
        if (sn[jj - 1] != sn[u])
          en += P->DuplexInit;

        mm5 = (sn[u] == sn[u - 1]) ? S1[u - 1] : -1;
        mm3 = (sn[jj] == sn[jj - 1]) ? S1[jj] : -1;

        if (evaluate(1, jj, u - 1, u, VRNA_DECOMP_EXT_EXT_STEM1, &hc_dat_local)) {
          e = my_f5[u - 1] + en + E_ExtLoop(type, -1, mm3, P);

          if (sc) {
            if (sc->energy_up)
              e += sc->energy_up[jj][1];
            if (sc->f)
              e += sc->f(1, jj, u - 1, u, VRNA_DECOMP_EXT_EXT_STEM1, sc->data);
          }

          if (fij == e) {
            *i                            = u;
            *j                            = jj - 1;
            *k                            = u - 1;
            bp_stack[++(*stack_count)].i  = u;
            bp_stack[(*stack_count)].j    = jj - 1;
            return 1;
          }
        }

        if (evaluate(1, jj, u - 2, u, VRNA_DECOMP_EXT_EXT_STEM1, &hc_dat_local)) {
          e = my_f5[u - 2] + en + E_ExtLoop(type, mm5, mm3, P);
          if (sc) {
            if (sc->energy_up)
              e += sc->energy_up[jj][1]
                   + sc->energy_up[u - 1][1];
            if (sc->f)
              e += sc->f(1, jj, u - 2, u, VRNA_DECOMP_EXT_EXT_STEM1, sc->data);
          }
          if (fij == e) {
            *i                            = u;
            *j                            = jj - 1;
            *k                            = u - 2;
            bp_stack[++(*stack_count)].i  = u;
            bp_stack[(*stack_count)].j    = jj - 1;
            return 1;
          }
        }
      }
      break;
  }

  return 0;
}


PUBLIC vrna_mx_pf_aux_el_t *
vrna_exp_E_ext_fast_init(vrna_fold_compound_t *vc)
{
  vrna_mx_pf_aux_el_t *aux_mx = NULL;

  if (vc) {
    char                      *hc;
    unsigned int              u, s;
    int                       i, j, d, n, turn, ij, *idx, *iidx, *hc_up;
    FLT_OR_DBL                *q, *scale;
    vrna_callback_hc_evaluate *evaluate;
    struct default_data       hc_dat_local;

    n     = (int)vc->length;
    idx   = vc->jindx;
    iidx  = vc->iindx;
    turn  = vc->exp_params->model_details.min_loop_size;
    q     = vc->exp_matrices->q;
    scale = vc->exp_matrices->scale;
    hc    = vc->hc->matrix;
    hc_up = vc->hc->up_ext;

    hc_dat_local.idx    = idx;
    hc_dat_local.mx     = hc;
    hc_dat_local.hc_up  = hc_up;
    hc_dat_local.cp     = vc->cutpoint;

    if (vc->hc->f) {
      evaluate            = &hc_default_user;
      hc_dat_local.hc_f   = vc->hc->f;
      hc_dat_local.hc_dat = vc->hc->data;
    } else {
      evaluate = &hc_default;
    }


    /* allocate memory for helper arrays */
    aux_mx            = (vrna_mx_pf_aux_el_t *)vrna_alloc(sizeof(vrna_mx_pf_aux_el_t));
    aux_mx->qq        = (FLT_OR_DBL *)vrna_alloc(sizeof(FLT_OR_DBL) * (n + 2));
    aux_mx->qq1       = (FLT_OR_DBL *)vrna_alloc(sizeof(FLT_OR_DBL) * (n + 2));
    aux_mx->qqu_size  = 0;
    aux_mx->qqu       = NULL;

    if (vc->type == VRNA_FC_TYPE_SINGLE) {
      vrna_sc_t *sc         = vc->sc;
      vrna_ud_t *domains_up = vc->domains_up;
      int       with_ud     = (domains_up && domains_up->exp_energy_cb);

      /* pre-processing ligand binding production rule(s) and auxiliary memory */
      if (with_ud) {
        int ud_max_size = 0;
        for (u = 0; u < domains_up->uniq_motif_count; u++)
          if (ud_max_size < domains_up->uniq_motif_size[u])
            ud_max_size = domains_up->uniq_motif_size[u];

        aux_mx->qqu_size  = ud_max_size;
        aux_mx->qqu       = (FLT_OR_DBL **)vrna_alloc(sizeof(FLT_OR_DBL *) * (ud_max_size + 1));

        for (u = 0; u <= ud_max_size; u++)
          aux_mx->qqu[u] = (FLT_OR_DBL *)vrna_alloc(sizeof(FLT_OR_DBL) * (n + 2));
      }

      for (d = 0; d <= turn; d++)
        for (i = 1; i <= n - d; i++) {
          j   = i + d;
          ij  = iidx[i] - j;

          if (j > n)
            continue;

          if (evaluate(i, j, i, j, VRNA_DECOMP_EXT_UP, &hc_dat_local)) {
            q[ij] = 1.0 * scale[d + 1];

            if (sc) {
              if (sc->exp_energy_up)
                q[ij] *= sc->exp_energy_up[i][d + 1];
              if (sc->exp_f)
                q[ij] *= sc->exp_f(i, j, i, j, VRNA_DECOMP_EXT_UP, sc->data);
            }

            if (with_ud) {
              q[ij] += q[ij] * domains_up->exp_energy_cb(vc,
                                                         i, j,
                                                         VRNA_UNSTRUCTURED_DOMAIN_EXT_LOOP,
                                                         domains_up->data);
            }
          } else {
            q[ij] = 0.;
          }
        }
    } else if (vc->type == VRNA_FC_TYPE_COMPARATIVE) {
      vrna_sc_t       **scs = vc->scs;
      unsigned short  **a2s = vc->a2s;
      for (d = 0; d <= turn; d++)
        for (i = 1; i <= n - d; i++) {
          j   = i + d;
          ij  = iidx[i] - j;
          if (evaluate(i, j, i, j, VRNA_DECOMP_EXT_UP, &hc_dat_local)) {
            q[ij] = 1.0 * scale[d + 1];

            if (scs) {
              for (s = 0; s < vc->n_seq; s++)
                if (scs[s]) {
                  u = d + 1 /* a2s[s][j] - a2s[s][i] + 1 */;
                  if (scs[s]->exp_energy_up)
                    q[ij] *= scs[s]->exp_energy_up[a2s[s][i]][u];
                }
            }
          } else {
            q[ij] = 0.;
          }
        }
    }
  }

  return aux_mx;
}


PUBLIC void
vrna_exp_E_ext_fast_rotate(vrna_fold_compound_t *vc,
                           vrna_mx_pf_aux_el_t  *aux_mx)
{
  if (vc && aux_mx) {
    int         u;
    FLT_OR_DBL  *tmp;

    tmp         = aux_mx->qq1;
    aux_mx->qq1 = aux_mx->qq;
    aux_mx->qq  = tmp;

    /* rotate auxiliary arrays for unstructured domains */
    if (aux_mx->qqu) {
      tmp = aux_mx->qqu[aux_mx->qqu_size];
      for (u = aux_mx->qqu_size; u > 0; u--)
        aux_mx->qqu[u] = aux_mx->qqu[u - 1];
      aux_mx->qqu[0] = tmp;
    }
  }
}


PUBLIC void
vrna_exp_E_ext_fast_free(vrna_fold_compound_t *vc,
                         vrna_mx_pf_aux_el_t  *aux_mx)
{
  if (vc && aux_mx) {
    int u;

    free(aux_mx->qq);
    free(aux_mx->qq1);

    if (aux_mx->qqu) {
      for (u = 0; u <= aux_mx->qqu_size; u++)
        free(aux_mx->qqu[u]);

      free(aux_mx->qqu);
    }
    free(aux_mx);
  }
}


PUBLIC FLT_OR_DBL
vrna_exp_E_ext_fast(vrna_fold_compound_t  *vc,
                    int                   i,
                    int                   j,
                    vrna_mx_pf_aux_el_t   *aux_mx)
{
  if (vc) {
    switch (vc->type) {
      case VRNA_FC_TYPE_SINGLE:
        return exp_E_ext_fast(vc, i, j, aux_mx);
        break;

      case VRNA_FC_TYPE_COMPARATIVE:
        return exp_E_ext_fast_comparative(vc, i, j, aux_mx);
        break;

      default:
        vrna_message_warning("vrna_exp_E_ext_fast@exterior_loops.c: Unknown fold_compound type");
        return 0.;
        break;
    }
  } else {
    return 0.;
  }
}


PRIVATE FLT_OR_DBL
exp_E_ext_fast(vrna_fold_compound_t *vc,
               int                  i,
               int                  j,
               vrna_mx_pf_aux_el_t  *aux_mx)
{
  short                     *S1;
  unsigned char             type;
  int                       n, *iidx, k, ij, kl, with_ud, u, circular, with_gquad;
  FLT_OR_DBL                qbt1, *q, *qb, *qq, *qq1, **qqu, q_temp, *scale, q_temp2, *G;
  vrna_md_t                 *md;
  vrna_exp_param_t          *pf_params;
  vrna_ud_t                 *domains_up;
  vrna_hc_t                 *hc;
  vrna_sc_t                 *sc;
  vrna_callback_hc_evaluate *evaluate;
  struct default_data       hc_dat_local;

  n                   = (int)vc->length;
  iidx                = vc->iindx;
  ij                  = iidx[i] - j;
  qq                  = aux_mx->qq;
  qq1                 = aux_mx->qq1;
  qqu                 = aux_mx->qqu;
  q                   = vc->exp_matrices->q;
  qb                  = vc->exp_matrices->qb;
  G                   = vc->exp_matrices->G;
  scale               = vc->exp_matrices->scale;
  pf_params           = vc->exp_params;
  md                  = &(pf_params->model_details);
  hc                  = vc->hc;
  sc                  = vc->sc;
  domains_up          = vc->domains_up;
  circular            = md->circ;
  with_gquad          = md->gquad;
  with_ud             = (domains_up && domains_up->exp_energy_cb);
  hc_dat_local.idx    = vc->jindx;
  hc_dat_local.mx     = hc->matrix;
  hc_dat_local.hc_up  = hc->up_ext;
  hc_dat_local.cp     = vc->cutpoint;

  if (hc->f) {
    evaluate            = &hc_default_user;
    hc_dat_local.hc_f   = hc->f;
    hc_dat_local.hc_dat = hc->data;
  } else {
    evaluate = &hc_default;
  }

  qbt1 = 0.;

  /* all exterior loop parts [i, j] with exactly one stem (i, u) i < u < j */
  if (evaluate(i, j, i, j - 1, VRNA_DECOMP_EXT_EXT, &hc_dat_local)) {
    q_temp = qq1[i] * scale[1];

    if (sc) {
      if (sc->exp_energy_up)
        q_temp *= sc->exp_energy_up[j][1];

      if (sc->exp_f)
        q_temp *= sc->exp_f(i, j, i, j - 1, VRNA_DECOMP_EXT_EXT, sc->data);
    }

    if (with_ud) {
      int cnt;
      for (cnt = 0; cnt < domains_up->uniq_motif_count; cnt++) {
        u = domains_up->uniq_motif_size[cnt];
        if (j - u >= i) {
          if (evaluate(i, j, i, j - u, VRNA_DECOMP_EXT_EXT, &hc_dat_local)) {
            q_temp2 = qqu[u][i]
                      * domains_up->exp_energy_cb(vc,
                                                  j - u + 1, j,
                                                  VRNA_UNSTRUCTURED_DOMAIN_EXT_LOOP | VRNA_UNSTRUCTURED_DOMAIN_MOTIF,
                                                  domains_up->data)
                      * scale[u];

            if (sc) {
              if (sc->exp_energy_up)
                q_temp2 *= sc->exp_energy_up[j - u + 1][u];
              if (sc->exp_f)
                q_temp2 *= sc->exp_f(i, j, i, j - u, VRNA_DECOMP_EXT_EXT, sc->data);
            }

            q_temp += q_temp2;
          }
        }
      }
    }

    qbt1 += q_temp;
  }

  /* exterior loop part with stem (i, j) */
  if (evaluate(i, j, i, j, VRNA_DECOMP_EXT_STEM, &hc_dat_local)) {
    S1    = vc->sequence_encoding;
    type  = md->pair[S1[i]][S1[j]];
    if (type == 0)
      type = 7;

    q_temp = qb[ij]
             * exp_E_ExtLoop(type, ((i > 1) || circular) ? S1[i - 1] : -1, ((j < n) || circular) ? S1[j + 1] : -1, pf_params);

    if (sc)
      if (sc->exp_f)
        q_temp *= sc->exp_f(i, j, i, j, VRNA_DECOMP_EXT_STEM, sc->data);
    qbt1 += q_temp;
  }

  if (with_gquad)
    qbt1 += G[ij];

  qq[i] = qbt1;

  if (with_ud)
    qqu[0][i] = qbt1;

  /* the entire stretch [i,j] is unpaired */
  if (evaluate(i, j, i, j, VRNA_DECOMP_EXT_UP, &hc_dat_local)) {
    u       = j - i + 1;
    q_temp  = 1.0 * scale[u];

    if (sc) {
      if (sc->exp_energy_up)
        q_temp *= sc->exp_energy_up[i][u];

      if (sc->exp_f)
        q_temp *= sc->exp_f(i, j, i, j, VRNA_DECOMP_EXT_UP, sc->data);
    }

    qbt1 += q_temp;

    if (with_ud) {
      qbt1 += q_temp * domains_up->exp_energy_cb(vc,
                                                 i, j,
                                                 VRNA_UNSTRUCTURED_DOMAIN_EXT_LOOP,
                                                 domains_up->data);
    }
  }

  kl = iidx[i] - j + 1;
  if (sc && sc->exp_f) {
    for (k = j; k > i; k--, kl++) {
      q_temp  = q[kl] * qq[k];
      q_temp  *= sc->exp_f(i, j, k - 1, k, VRNA_DECOMP_EXT_EXT_EXT, sc->data);
      qbt1    += q_temp;
    }
  } else {
    for (k = j; k > i; k--, kl++)
      qbt1 += q[kl] * qq[k];
  }

  return qbt1;
}


PRIVATE FLT_OR_DBL
exp_E_ext_fast_comparative(vrna_fold_compound_t *vc,
                           int                  i,
                           int                  j,
                           vrna_mx_pf_aux_el_t  *aux_mx)
{
  int                       n, s, n_seq, *iidx, k, ij, kl, u, circular, type;
  unsigned short            **a2s;
  short                     **S, **S5, **S3;
  FLT_OR_DBL                qbt1, *q, *qb, *qq, *qq1, q_temp, *scale;
  vrna_md_t                 *md;
  vrna_exp_param_t          *pf_params;
  vrna_hc_t                 *hc;
  vrna_sc_t                 **scs;
  vrna_callback_hc_evaluate *evaluate;
  struct default_data       hc_dat_local;

  n                   = (int)vc->length;
  n_seq               = vc->n_seq;
  iidx                = vc->iindx;
  ij                  = iidx[i] - j;
  S                   = vc->S;
  S5                  = vc->S5;     /* S5[s][i] holds next base 5' of i in sequence s */
  S3                  = vc->S3;     /* Sl[s][i] holds next base 3' of i in sequence s */
  a2s                 = vc->a2s;
  qq                  = aux_mx->qq;
  qq1                 = aux_mx->qq1;
  q                   = vc->exp_matrices->q;
  qb                  = vc->exp_matrices->qb;
  scale               = vc->exp_matrices->scale;
  pf_params           = vc->exp_params;
  md                  = &(pf_params->model_details);
  hc                  = vc->hc;
  scs                 = vc->scs;
  circular            = md->circ;
  hc_dat_local.idx    = vc->jindx;
  hc_dat_local.mx     = hc->matrix;
  hc_dat_local.hc_up  = hc->up_ext;
  hc_dat_local.cp     = vc->cutpoint;

  if (hc->f) {
    evaluate            = &hc_default_user;
    hc_dat_local.hc_f   = hc->f;
    hc_dat_local.hc_dat = hc->data;
  } else {
    evaluate = &hc_default;
  }

  qbt1 = 0.;

  /* all exterior loop parts [i, j] with exactly one stem (i, u) i < u < j */
  if (evaluate(i, j, i, j - 1, VRNA_DECOMP_EXT_EXT, &hc_dat_local)) {
    q_temp = qq1[i] * scale[1];

    if (scs) {
      for (s = 0; s < n_seq; s++) {
        if (scs[s])
          if (scs[s]->exp_energy_up)
            q_temp *= scs[s]->exp_energy_up[a2s[s][j]][1];
      }
    }

    qbt1 += q_temp;
  }

  /* exterior loop part with stem (i, j) */
  if (evaluate(i, j, i, j, VRNA_DECOMP_EXT_STEM, &hc_dat_local)) {
    q_temp = qb[ij];

    for (s = 0; s < n_seq; s++) {
      type = md->pair[S[s][i]][S[s][j]];
      if (type == 0)
        type = 7;

      q_temp *= exp_E_ExtLoop(type, ((i > 1) || circular) ? S5[s][i] : -1, ((j < n) || circular) ? S3[s][j] : -1, pf_params);
    }

    qbt1 += q_temp;
  }

  qq[i] = qbt1;

  /* the entire stretch [i,j] is unpaired */
  if (evaluate(i, j, i, j, VRNA_DECOMP_EXT_UP, &hc_dat_local)) {
    u       = j - i + 1;
    q_temp  = 1.0 * scale[u];

    if (scs) {
      for (s = 0; s < n_seq; s++) {
        if (scs[s])
          if (scs[s]->exp_energy_up)
            q_temp *= scs[s]->exp_energy_up[a2s[s][i]][a2s[s][j] - a2s[s][i] + 1];
      }
    }

    qbt1 += q_temp;
  }

  kl = iidx[i] - j + 1;
  for (k = j; k > i; k--, kl++)
    qbt1 += q[kl] * qq[k];

  return qbt1;
}


PRIVATE char
hc_default(int  i,
           int  j,
           int  k,
           int  l,
           char d,
           void *data)
{
  int                 kl, di, dj;
  char                eval;
  struct default_data *dat = (struct default_data *)data;

  eval  = (char)0;
  di    = k - i;
  dj    = j - l;
  switch (d) {
    case VRNA_DECOMP_EXT_EXT_STEM:
      kl = dat->idx[j] + l;
      if (dat->mx[kl] & VRNA_CONSTRAINT_CONTEXT_EXT_LOOP) {
        eval = (char)1;
        if (i != l) {
          /* otherwise, stem spans from i to j */
          di = l - k + 1;
          if ((di != 0) && (dat->hc_up[k + 1] < di))
            eval = (char)0;
        }
      }
      break;

    case VRNA_DECOMP_EXT_EXT_STEM1:
      kl = dat->idx[j - 1] + l;
      if (dat->mx[kl] & VRNA_CONSTRAINT_CONTEXT_EXT_LOOP) {
        eval = (char)1;
        if (i != l) {
          /* otherwise, stem spans from i to j - 1 */
          di = l - k + 1;
          if (dat->hc_up[j] == 0)
            eval = (char)0;
          if ((di != 0) && (dat->hc_up[k + 1] < di))
            eval = (char)0;
        }
      }
      break;

    case VRNA_DECOMP_EXT_STEM:
      kl = dat->idx[l] + k;
      if (dat->mx[kl] & VRNA_CONSTRAINT_CONTEXT_EXT_LOOP) {
        eval = (char)1;
        if ((di != 0) && (dat->hc_up[i] < di))
          eval = (char)0;
        if ((dj != 0) && (dat->hc_up[l + 1] < dj))
          eval = (char)0;
      }
      break;

    case VRNA_DECOMP_EXT_EXT:
      eval = (char)1;
      if ((di != 0) && (dat->hc_up[i] < di))
        eval = (char)0;
      if ((dj != 0) && (dat->hc_up[l + 1] < dj))
        eval = (char)0;
      break;

    case VRNA_DECOMP_EXT_UP:
      di    = j - i + 1;
      eval  = (dat->hc_up[i] >= di) ? (char)1 : (char)0;
      break;

    default:
      nrerror("wtf");
  }
  return eval;
}


PRIVATE char
hc_default_user(int   i,
                int   j,
                int   k,
                int   l,
                char  d,
                void  *data)
{
  char                eval;
  struct default_data *dat = (struct default_data *)data;

  eval  = hc_default(i, j, k, l, d, data);
  eval  = (dat->hc_f(i, j, k, l, d, dat->hc_dat)) ? eval : (char)0;

  return eval;
}