/* This file contains some of the old nuclear family update routines */
/* shared by modified versions of the LODSCORE, ILINK, LINKMAP, and
MLINK programs*/

Local double msegsex(LINK)
struct LOC_seg *LINK;
{
  int g1, g2, g3, g4, g5, g6, g7, g8, ms, mf, ms1, ms2, mf1, mf2, j, k, f1,
       f2, s1, s2;
  double val, temp2;
  double temp[maxchild];
  int FORLIM;
  gennurec *WITH;
  thetavalues *WITH1;
  int FORLIM1;
  thisarray *WITH2;

  FORLIM = nchild;
  for (k = 0; k < FORLIM; k++)
    temp[k] = 0.0;
  WITH = gennustruct;
  if ((*LINK->p)->male) {
    mf = muthap[LINK->fseg - 1];
    LINK->secondseg = LINK->sseg;
    WITH1 = LINK->secondsex;
    FORLIM = LINK->send;
    for (j = LINK->sstart - 1; j < FORLIM; j++) {
      if (WITH1->segprob[j] == 0.0)
	LINK->secondseg++;
      else {
	temp2 = WITH1->segprob[j];
	s1 = invgenenum1[LINK->secondseg - 1];
	s2 = invgenenum2[LINK->secondseg - 1];
	ms1 = muthap[s1 - 1];
	ms2 = muthap[s2 - 1];
	if (s1 != s2) {
	  g1 = WITH->genenumber[LINK->fseg - 1][s1 - 1];
	  g2 = WITH->genenumber[LINK->fseg - 1][s2 - 1];
	  g3 = WITH->genenumber[LINK->fseg - 1][ms1 - 1];
	  g4 = WITH->genenumber[LINK->fseg - 1][ms2 - 1];
	  g5 = WITH->genenumber[mf - 1][s1 - 1];
	  g6 = WITH->genenumber[mf - 1][s2 - 1];
	  g7 = WITH->genenumber[mf - 1][ms1 - 1];
	  g8 = WITH->genenumber[mf - 1][ms2 - 1];
	  FORLIM1 = nchild;
	  for (k = 0; k < FORLIM1; k++) {
	    WITH2 = thischild[k];
	    if (malechild[k])
	      temp[k] += temp2 * ((1 - LINK->ps) * (WITH2->genarray[s1 - 1] +
				    WITH2->genarray[s2 - 1]) + LINK->ps *
		      (WITH2->genarray[ms1 - 1] + WITH2->genarray[ms2 - 1]));
	    else
	      temp[k] += temp2 * ((1 - LINK->pf) * (1 - LINK->ps) *
		      (WITH2->genarray[g1 - 1] + WITH2->genarray[g2 - 1]) +
		    (1 - LINK->pf) * LINK->ps * (WITH2->genarray[g3 - 1] +
			WITH2->genarray[g4 - 1]) + LINK->pf * (1 - LINK->ps) *
		      (WITH2->genarray[g5 - 1] + WITH2->genarray[g6 - 1]) +
		    LINK->pf * LINK->ps *
		      (WITH2->genarray[g7 - 1] + WITH2->genarray[g8 - 1]));
	  }
	} else {
	  g1 = WITH->genenumber[LINK->fseg - 1][s1 - 1];
	  g3 = WITH->genenumber[LINK->fseg - 1][ms1 - 1];
	  g5 = WITH->genenumber[mf - 1][s1 - 1];
	  g7 = WITH->genenumber[mf - 1][ms1 - 1];
	  FORLIM1 = nchild;
	  for (k = 0; k < FORLIM1; k++) {
	    WITH2 = thischild[k];
	    if (malechild[k])
	      temp[k] += 2.0 * temp2 * ((1 - LINK->ps) * WITH2->genarray[s1 - 1] +
					LINK->ps * WITH2->genarray[ms1 - 1]);
	    else
	      temp[k] += 2.0 * temp2 *
		  ((1 - LINK->pf) * (1 - LINK->ps) * WITH2->genarray[g1 - 1] +
		   LINK->ps * (1 - LINK->pf) * WITH2->genarray[g3 - 1] +
		   LINK->pf * (1 - LINK->ps) * WITH2->genarray[g5 - 1] +
		   LINK->pf * LINK->ps * WITH2->genarray[g7 - 1]);
	  }
	}
	LINK->secondseg++;
      }
    }
  } else {
    LINK->firstseg = LINK->fseg;
    ms = muthap[LINK->sseg - 1];
    WITH1 = LINK->firstsex;
    FORLIM = LINK->fend;
    for (j = LINK->fstart - 1; j < FORLIM; j++) {
      if (WITH1->segprob[j] == 0.0)
	LINK->firstseg++;
      else {
	temp2 = WITH1->segprob[j];
	f1 = invgenenum1[LINK->firstseg - 1];
	f2 = invgenenum2[LINK->firstseg - 1];
	mf1 = muthap[f1 - 1];
	mf2 = muthap[f2 - 1];
	if (f1 != f2) {
	  g1 = WITH->genenumber[LINK->sseg - 1][f1 - 1];
	  g2 = WITH->genenumber[LINK->sseg - 1][f2 - 1];
	  g3 = WITH->genenumber[LINK->sseg - 1][mf1 - 1];
	  g4 = WITH->genenumber[LINK->sseg - 1][mf2 - 1];
	  g5 = WITH->genenumber[ms - 1][f1 - 1];
	  g6 = WITH->genenumber[ms - 1][f2 - 1];
	  g7 = WITH->genenumber[ms - 1][mf1 - 1];
	  g8 = WITH->genenumber[ms - 1][mf2 - 1];
	  FORLIM1 = nchild;
	  for (k = 0; k < FORLIM1; k++) {
	    WITH2 = thischild[k];
	    if (malechild[k])
	      temp[k] += temp2 * ((1 - LINK->pf) * (WITH2->genarray[f1 - 1] +
				    WITH2->genarray[f2 - 1]) + LINK->pf *
		      (WITH2->genarray[mf1 - 1] + WITH2->genarray[mf2 - 1]));
	    else
	      temp[k] += temp2 * ((1 - LINK->pf) * (1 - LINK->ps) *
		      (WITH2->genarray[g1 - 1] + WITH2->genarray[g2 - 1]) +
		    (1 - LINK->ps) * LINK->pf * (WITH2->genarray[g3 - 1] +
			WITH2->genarray[g4 - 1]) + LINK->ps * (1 - LINK->pf) *
		      (WITH2->genarray[g5 - 1] + WITH2->genarray[g6 - 1]) +
		    LINK->pf * LINK->ps *
		      (WITH2->genarray[g7 - 1] + WITH2->genarray[g8 - 1]));
	  }
	} else {
	  g1 = WITH->genenumber[LINK->sseg - 1][f1 - 1];
	  g3 = WITH->genenumber[LINK->sseg - 1][mf1 - 1];
	  g5 = WITH->genenumber[ms - 1][f1 - 1];
	  g7 = WITH->genenumber[ms - 1][mf1 - 1];
	  FORLIM1 = nchild;
	  for (k = 0; k < FORLIM1; k++) {
	    WITH2 = thischild[k];
	    if (malechild[k])
	      temp[k] += 2.0 * temp2 * ((1 - LINK->pf) * WITH2->genarray[f1 - 1] +
					LINK->pf * WITH2->genarray[mf1 - 1]);
	    else
	      temp[k] += 2.0 * temp2 *
		  ((1 - LINK->pf) * (1 - LINK->ps) * WITH2->genarray[g1 - 1] +
		   LINK->pf * (1 - LINK->ps) * WITH2->genarray[g3 - 1] +
		   LINK->ps * (1 - LINK->pf) * WITH2->genarray[g5 - 1] +
		   LINK->ps * LINK->pf * WITH2->genarray[g7 - 1]);
	  }
	}
	LINK->firstseg++;
      }
    }
  }
  val = 1.0;
  FORLIM = nchild;
  for (k = 0; k < FORLIM; k++)
    val *= temp[k];
  return val;
}  /*msegsex*/


Local double msegsexf(LINK)
struct LOC_seg *LINK;
{
  int g1, g2, g3, g4, g5, g6, g7, g8, mf, ms1, ms2, j, k, l, s1, s2, slength;
  double val, temp1;
  double temp[maxchild][maxseg];
  double temp2[maxseg];
  int FORLIM, FORLIM1;
  gennurec *WITH;
  thetavalues *WITH1;
  thisarray *WITH2;

  mf = muthap[LINK->fseg - 1];
  slength = LINK->send - LINK->sstart + 1;
  FORLIM = nchild;
  for (k = 0; k < FORLIM; k++) {
    for (l = 0; l < slength; l++)
      temp[k][l] = 0.0;
  }
  LINK->secondseg = LINK->sseg;
  WITH = gennustruct;
  WITH1 = LINK->secondsex;
  FORLIM = LINK->send;
  for (j = LINK->sstart - 1; j < FORLIM; j++) {
    for (l = 0; l < slength; l++)
      temp2[l] = WITH1->segprob[j + l * slength];
    s1 = invgenenum1[LINK->secondseg - 1];
    s2 = invgenenum2[LINK->secondseg - 1];
    ms1 = muthap[s1 - 1];
    ms2 = muthap[s2 - 1];
    if (s1 != s2) {
      g1 = WITH->genenumber[LINK->fseg - 1][s1 - 1];
      g2 = WITH->genenumber[LINK->fseg - 1][s2 - 1];
      g3 = WITH->genenumber[LINK->fseg - 1][ms1 - 1];
      g4 = WITH->genenumber[LINK->fseg - 1][ms2 - 1];
      g5 = WITH->genenumber[mf - 1][s1 - 1];
      g6 = WITH->genenumber[mf - 1][s2 - 1];
      g7 = WITH->genenumber[mf - 1][ms1 - 1];
      g8 = WITH->genenumber[mf - 1][ms2 - 1];
      FORLIM1 = nchild;
      for (k = 0; k < FORLIM1; k++) {
	WITH2 = thischild[k];
	if (malechild[k])
	  val = (1 - LINK->ps) * (WITH2->genarray[s1 - 1] + WITH2->genarray[s2 - 1]) +
		LINK->ps * (WITH2->genarray[ms1 - 1] + WITH2->genarray[ms2 - 1]);
	else
	  val = (1 - LINK->pf) * (1 - LINK->ps) * (WITH2->genarray[g1 - 1] +
		  WITH2->genarray[g2 - 1]) + (1 - LINK->pf) * LINK->ps *
		(WITH2->genarray[g3 - 1] + WITH2->genarray[g4 - 1]) +
	      LINK->pf * (1 - LINK->ps) * (WITH2->genarray[g5 - 1] + WITH2->
		    genarray[g6 - 1]) + LINK->pf * LINK->ps *
		(WITH2->genarray[g7 - 1] + WITH2->genarray[g8 - 1]);
	for (l = 0; l < slength; l++)
	  temp[k][l] += temp2[l] * val;
      }
    } else {
      g1 = WITH->genenumber[LINK->fseg - 1][s1 - 1];
      g3 = WITH->genenumber[LINK->fseg - 1][ms1 - 1];
      g5 = WITH->genenumber[mf - 1][s1 - 1];
      g7 = WITH->genenumber[mf - 1][ms1 - 1];
      FORLIM1 = nchild;
      for (k = 0; k < FORLIM1; k++) {
	WITH2 = thischild[k];
	if (malechild[k])
	  val = 2.0 * ((1 - LINK->ps) * WITH2->genarray[s1 - 1] +
		       LINK->ps * WITH2->genarray[ms1 - 1]);
	else
	  val = 2.0 * ((1 - LINK->pf) * (1 - LINK->ps) * WITH2->genarray[g1 - 1] +
		       LINK->ps * (1 - LINK->pf) * WITH2->genarray[g3 - 1] +
		       LINK->pf * (1 - LINK->ps) * WITH2->genarray[g5 - 1] +
		       LINK->ps * LINK->pf * WITH2->genarray[g7 - 1]);
	for (l = 0; l < slength; l++)
	  temp[k][l] += temp2[l] * val;
      }
    }
    LINK->secondseg++;
  }
  temp1 = 0.0;
  for (l = 0; l < slength; l++) {
    val = 1.0;
    FORLIM1 = nchild;
    for (k = 0; k < FORLIM1; k++)
      val *= temp[k][l];
    temp1 += val;
  }
  return temp1;
}  /*msegsexf*/


Local double segsex(LINK)
struct LOC_seg *LINK;
{
  int g1, g2, j, k, f1, f2, s1, s2;
  double val, temp2;
  double temp[maxchild];
  int FORLIM;
  gennurec *WITH;
  thetavalues *WITH1;
  int FORLIM1;
  thisarray *WITH2;

  FORLIM = nchild;
  for (k = 0; k < FORLIM; k++)
    temp[k] = 0.0;
  WITH = gennustruct;
  if ((*LINK->p)->male) {
    LINK->secondseg = LINK->sseg;
    WITH1 = LINK->secondsex;
    FORLIM = LINK->send;
    for (j = LINK->sstart - 1; j < FORLIM; j++) {
      if (WITH1->segprob[j] == 0.0)
	LINK->secondseg++;
      else {
	temp2 = WITH1->segprob[j];
	s1 = invgenenum1[LINK->secondseg - 1];
	s2 = invgenenum2[LINK->secondseg - 1];
	if (s1 != s2) {
	  g1 = WITH->genenumber[LINK->fseg - 1][s1 - 1];
	  g2 = WITH->genenumber[LINK->fseg - 1][s2 - 1];
	  FORLIM1 = nchild;
	  for (k = 0; k < FORLIM1; k++) {
	    WITH2 = thischild[k];
	    if (malechild[k])
	      temp[k] += temp2 *
			 (WITH2->genarray[s1 - 1] + WITH2->genarray[s2 - 1]);
	    else
	      temp[k] += temp2 *
			 (WITH2->genarray[g1 - 1] + WITH2->genarray[g2 - 1]);
	  }
	} else {
	  g1 = WITH->genenumber[LINK->fseg - 1][s1 - 1];
	  FORLIM1 = nchild;
	  for (k = 0; k < FORLIM1; k++) {
	    WITH2 = thischild[k];
	    if (malechild[k])
	      temp[k] += 2.0 * temp2 * WITH2->genarray[s1 - 1];
	    else
	      temp[k] += 2.0 * temp2 * WITH2->genarray[g1 - 1];
	  }
	}
	LINK->secondseg++;
      }
    }
  } else {
    LINK->firstseg = LINK->fseg;
    WITH1 = LINK->firstsex;
    FORLIM = LINK->fend;
    for (j = LINK->fstart - 1; j < FORLIM; j++) {
      if (WITH1->segprob[j] == 0.0)
	LINK->firstseg++;
      else {
	temp2 = WITH1->segprob[j];
	f1 = invgenenum1[LINK->firstseg - 1];
	f2 = invgenenum2[LINK->firstseg - 1];
	if (f1 != f2) {
	  g1 = WITH->genenumber[LINK->sseg - 1][f1 - 1];
	  g2 = WITH->genenumber[LINK->sseg - 1][f2 - 1];
	  FORLIM1 = nchild;
	  for (k = 0; k < FORLIM1; k++) {
	    WITH2 = thischild[k];
	    if (malechild[k])
	      temp[k] += temp2 *
			 (WITH2->genarray[f1 - 1] + WITH2->genarray[f2 - 1]);
	    else
	      temp[k] += temp2 *
			 (WITH2->genarray[g1 - 1] + WITH2->genarray[g2 - 1]);
	  }
	} else {
	  g1 = WITH->genenumber[LINK->sseg - 1][f1 - 1];
	  FORLIM1 = nchild;
	  for (k = 0; k < FORLIM1; k++) {
	    WITH2 = thischild[k];
	    if (malechild[k])
	      temp[k] += 2.0 * temp2 * WITH2->genarray[f1 - 1];
	    else
	      temp[k] += 2.0 * temp2 * WITH2->genarray[g1 - 1];
	  }
	}
	LINK->firstseg++;
      }
    }
  }
  val = 1.0;
  FORLIM = nchild;
  for (k = 0; k < FORLIM; k++)
    val *= temp[k];
  return val;
}  /*segsex*/

Local double segsexf(LINK)
struct LOC_seg *LINK;
{
  int g1, g2, j, k, l, s1, s2, slength;
  double val, temp1;
  double temp[maxchild][maxseg];
  double temp2[maxseg];
  int FORLIM, FORLIM1;
  gennurec *WITH;
  thetavalues *WITH1;
  thisarray *WITH2;

  slength = LINK->send - LINK->sstart + 1;
  FORLIM = nchild;
  for (k = 0; k < FORLIM; k++) {
    for (l = 0; l < slength; l++)
      temp[k][l] = 0.0;
  }
  LINK->secondseg = LINK->sseg;
  WITH = gennustruct;
  WITH1 = LINK->secondsex;
  FORLIM = LINK->send;
  for (j = LINK->sstart - 1; j < FORLIM; j++) {
    for (l = 0; l < slength; l++)
      temp2[l] = WITH1->segprob[j + l * slength];
    s1 = invgenenum1[LINK->secondseg - 1];
    s2 = invgenenum2[LINK->secondseg - 1];
    if (s1 != s2) {
      g1 = WITH->genenumber[LINK->fseg - 1][s1 - 1];
      g2 = WITH->genenumber[LINK->fseg - 1][s2 - 1];
      FORLIM1 = nchild;
      for (k = 0; k < FORLIM1; k++) {
	WITH2 = thischild[k];
	if (malechild[k])
	  val = WITH2->genarray[s1 - 1] + WITH2->genarray[s2 - 1];
	else
	  val = WITH2->genarray[g1 - 1] + WITH2->genarray[g2 - 1];
	if (val != 0.0) {
	  for (l = 0; l < slength; l++)
	    temp[k][l] += temp2[l] * val;
	}
      }
    } else {
      g1 = WITH->genenumber[LINK->fseg - 1][s1 - 1];
      FORLIM1 = nchild;
      for (k = 0; k < FORLIM1; k++) {
	WITH2 = thischild[k];
	if (malechild[k])
	  val = 2.0 * WITH2->genarray[s1 - 1];
	else
	  val = 2.0 * WITH2->genarray[g1 - 1];
	if (val != 0.0) {
	  for (l = 0; l < slength; l++)
	    temp[k][l] += temp2[l] * val;
	}
      }
    }
    LINK->secondseg++;
  }
  temp1 = 0.0;
  for (l = 0; l < slength; l++) {
    val = 1.0;
    FORLIM1 = nchild;
    for (k = 0; k < FORLIM1; k++)
      val *= temp[k][l];
    temp1 += val;
  }
  return temp1;
}

Local double segfun(LINK)
struct LOC_seg *LINK;
{
  int g1, g2, g3, g4, i, j, k, f1, f2, s1, s2;
  double val, temp1, temp2;
  double temp[maxchild];
  int FORLIM;
  gennurec *WITH;
  thetavalues *WITH1, *WITH2;
  int FORLIM1, FORLIM2;
  thisarray *WITH3;

  LINK->firstseg = LINK->fseg;
  FORLIM = nchild;
  for (k = 0; k < FORLIM; k++)
    temp[k] = 0.0;
  WITH = gennustruct;
  WITH1 = LINK->firstsex;
  FORLIM = LINK->fend;
  for (i = LINK->fstart - 1; i < FORLIM; i++) {
    if (WITH1->segprob[i] == 0.0)
      LINK->firstseg++;
    else {
      temp1 = WITH1->segprob[i];
      f1 = invgenenum1[LINK->firstseg - 1];
      f2 = invgenenum2[LINK->firstseg - 1];
      LINK->secondseg = LINK->sseg;
      WITH2 = LINK->secondsex;
      if (f1 != f2) {
	FORLIM1 = LINK->send;
	for (j = LINK->sstart - 1; j < FORLIM1; j++) {
	  if (WITH2->segprob[j] == 0.0)
	    LINK->secondseg++;
	  else {
	    temp2 = temp1 * WITH2->segprob[j];
	    s1 = invgenenum1[LINK->secondseg - 1];
	    s2 = invgenenum2[LINK->secondseg - 1];
	    if (s1 != s2) {
	      g1 = WITH->genenumber[f1 - 1][s1 - 1];
	      g2 = WITH->genenumber[f1 - 1][s2 - 1];
	      g3 = WITH->genenumber[f2 - 1][s1 - 1];
	      g4 = WITH->genenumber[f2 - 1][s2 - 1];
	      FORLIM2 = nchild;
	      for (k = 0; k < FORLIM2; k++) {
		WITH3 = thischild[k];
		temp[k] +=
		  temp2 * (WITH3->genarray[g1 - 1] + WITH3->genarray[g2 - 1] +
			   WITH3->genarray[g3 - 1] + WITH3->genarray[g4 - 1]);
	      }
	    } else {
	      g1 = WITH->genenumber[f1 - 1][s1 - 1];
	      g3 = WITH->genenumber[f2 - 1][s1 - 1];
	      FORLIM2 = nchild;
	      for (k = 0; k < FORLIM2; k++) {
		WITH3 = thischild[k];
		temp[k] += 2 * temp2 *
			   (WITH3->genarray[g1 - 1] + WITH3->genarray[g3 - 1]);
	      }
	    }
	    LINK->secondseg++;
	  }
	}
      } else {
	FORLIM1 = LINK->send;
	for (j = LINK->sstart - 1; j < FORLIM1; j++) {
	  if (WITH2->segprob[j] == 0.0)
	    LINK->secondseg++;
	  else {
	    temp2 = temp1 * WITH2->segprob[j];
	    s1 = invgenenum1[LINK->secondseg - 1];
	    s2 = invgenenum2[LINK->secondseg - 1];
	    if (s1 != s2) {
	      g1 = WITH->genenumber[f1 - 1][s1 - 1];
	      g2 = WITH->genenumber[f1 - 1][s2 - 1];
	      FORLIM2 = nchild;
	      for (k = 0; k < FORLIM2; k++) {
		WITH3 = thischild[k];
		temp[k] += 2 * temp2 *
			   (WITH3->genarray[g1 - 1] + WITH3->genarray[g2 - 1]);
	      }
	    } else {
	      g1 = WITH->genenumber[f1 - 1][s1 - 1];
	      FORLIM2 = nchild;
	      for (k = 0; k < FORLIM2; k++) {
		WITH3 = thischild[k];
		temp[k] += 4 * temp2 * WITH3->genarray[g1 - 1];
	      }
	    }
	    LINK->secondseg++;
	  }
	}
      }
      LINK->firstseg++;
    }
  }
  val = 1.0;
  FORLIM = nchild;
  for (k = 0; k < FORLIM; k++)
    val *= temp[k];
  return val;
}

Local double msegfast(LINK)
struct LOC_seg *LINK;
{
  int g1, g2, g3, g4, g5, g6, g7, g8, g9, g10, g11, g12, g13, g14, g15, g16,
       i, j, k, l, f1, f2, s1, s2, ms1, ms2, mf1, mf2, slength;
  double val, temp1;
  double temp[maxchild][maxseg];
  double temp2[maxseg];
  int FORLIM, FORLIM1;
  gennurec *WITH;
  thetavalues *WITH1, *WITH2;
  int FORLIM2;
  thisarray *WITH3;

  LINK->firstseg = LINK->fseg;
  slength = LINK->send - LINK->sstart + 1;
  FORLIM = nchild;
  for (k = 0; k < FORLIM; k++) {
    for (l = 0; l < slength; l++)
      temp[k][l] = 0.0;
  }
  WITH = gennustruct;
  WITH1 = LINK->firstsex;
  FORLIM = LINK->fend;
  for (i = LINK->fstart - 1; i < FORLIM; i++) {
    if (WITH1->segprob[i] == 0.0)
      LINK->firstseg++;
    else {
      temp1 = WITH1->segprob[i];
      f1 = invgenenum1[LINK->firstseg - 1];
      f2 = invgenenum2[LINK->firstseg - 1];
      mf1 = muthap[f1 - 1];
      mf2 = muthap[f2 - 1];
      LINK->secondseg = LINK->sseg;
      WITH2 = LINK->secondsex;
      if (f1 != f2) {
	FORLIM1 = LINK->send;
	for (j = LINK->sstart - 1; j < FORLIM1; j++) {
	  for (l = 0; l < slength; l++)
	    temp2[l] = temp1 * WITH2->segprob[j + l * slength];
	  s1 = invgenenum1[LINK->secondseg - 1];
	  s2 = invgenenum2[LINK->secondseg - 1];
	  ms1 = muthap[s1 - 1];
	  ms2 = muthap[s2 - 1];
	  if (s1 != s2) {
	    g1 = WITH->genenumber[f1 - 1][s1 - 1];
	    g2 = WITH->genenumber[f1 - 1][s2 - 1];
	    g3 = WITH->genenumber[f2 - 1][s1 - 1];
	    g4 = WITH->genenumber[f2 - 1][s2 - 1];
	    g5 = WITH->genenumber[f1 - 1][ms1 - 1];
	    g6 = WITH->genenumber[f1 - 1][ms2 - 1];
	    g7 = WITH->genenumber[f2 - 1][ms1 - 1];
	    g8 = WITH->genenumber[f2 - 1][ms2 - 1];
	    g9 = WITH->genenumber[mf1 - 1][s1 - 1];
	    g10 = WITH->genenumber[mf1 - 1][s2 - 1];
	    g11 = WITH->genenumber[mf2 - 1][s1 - 1];
	    g12 = WITH->genenumber[mf2 - 1][s2 - 1];
	    g13 = WITH->genenumber[mf1 - 1][ms1 - 1];
	    g14 = WITH->genenumber[mf1 - 1][ms2 - 1];
	    g15 = WITH->genenumber[mf2 - 1][ms1 - 1];
	    g16 = WITH->genenumber[mf2 - 1][ms2 - 1];
	    FORLIM2 = nchild;
	    for (k = 0; k < FORLIM2; k++) {
	      WITH3 = thischild[k];
	      val = (1 - LINK->ps) * (1 - LINK->pf) *
		    (WITH3->genarray[g1 - 1] + WITH3->genarray[g2 - 1] +
		      WITH3->genarray[g3 - 1] + WITH3->genarray[g4 - 1]) +
		  LINK->ps * (1 - LINK->pf) * (WITH3->genarray[g5 - 1] +
		      WITH3->genarray[g6 - 1] + WITH3->genarray[g7 - 1] +
		      WITH3->genarray[g8 - 1]) + LINK->pf * (1 - LINK->ps) *
		    (WITH3->genarray[g9 - 1] + WITH3->genarray[g10 - 1] +
		      WITH3->genarray[g11 - 1] + WITH3->genarray[g12 - 1]) +
		  LINK->pf * LINK->ps *
		    (WITH3->genarray[g13 - 1] + WITH3->genarray[g14 - 1] +
		      WITH3->genarray[g15 - 1] + WITH3->genarray[g16 - 1]);
	      if (val != 0.0) {
		for (l = 0; l < slength; l++)
		  temp[k][l] += temp2[l] * val;
	      }
	    }
	  } else {
	    g1 = WITH->genenumber[f1 - 1][s1 - 1];
	    g3 = WITH->genenumber[f2 - 1][s1 - 1];
	    g5 = WITH->genenumber[f1 - 1][ms1 - 1];
	    g7 = WITH->genenumber[f2 - 1][ms1 - 1];
	    g9 = WITH->genenumber[mf1 - 1][s1 - 1];
	    g11 = WITH->genenumber[mf2 - 1][s1 - 1];
	    g13 = WITH->genenumber[mf1 - 1][ms1 - 1];
	    g15 = WITH->genenumber[mf2 - 1][ms1 - 1];
	    FORLIM2 = nchild;
	    for (k = 0; k < FORLIM2; k++) {
	      WITH3 = thischild[k];
	      val = 2 * ((1 - LINK->ps) * (1 - LINK->pf) *
		      (WITH3->genarray[g1 - 1] + WITH3->genarray[g3 - 1]) +
		    LINK->ps * (1 - LINK->pf) * (WITH3->genarray[g5 - 1] +
			WITH3->genarray[g7 - 1]) + LINK->pf * (1 - LINK->ps) *
		      (WITH3->genarray[g9 - 1] + WITH3->genarray[g11 - 1]) +
		    LINK->pf * LINK->ps *
		      (WITH3->genarray[g13 - 1] + WITH3->genarray[g15 - 1]));
	      for (l = 0; l < slength; l++)
		temp[k][l] += temp2[l] * val;
	    }
	  }
	  LINK->secondseg++;
	}
      } else {
	FORLIM1 = LINK->send;
	for (j = LINK->sstart - 1; j < FORLIM1; j++) {
	  for (l = 0; l < slength; l++)
	    temp2[l] = temp1 * WITH2->segprob[j + l * slength];
	  s1 = invgenenum1[LINK->secondseg - 1];
	  s2 = invgenenum2[LINK->secondseg - 1];
	  ms1 = muthap[s1 - 1];
	  ms2 = muthap[s2 - 1];
	  if (s1 != s2) {
	    g1 = WITH->genenumber[f1 - 1][s1 - 1];
	    g2 = WITH->genenumber[f1 - 1][s2 - 1];
	    g5 = WITH->genenumber[f1 - 1][ms1 - 1];
	    g6 = WITH->genenumber[f1 - 1][ms2 - 1];
	    g9 = WITH->genenumber[mf1 - 1][s1 - 1];
	    g10 = WITH->genenumber[mf1 - 1][s2 - 1];
	    g13 = WITH->genenumber[mf1 - 1][ms1 - 1];
	    g14 = WITH->genenumber[mf1 - 1][ms2 - 1];
	    FORLIM2 = nchild;
	    for (k = 0; k < FORLIM2; k++) {
	      WITH3 = thischild[k];
	      val = 2 * ((1 - LINK->ps) * (1 - LINK->pf) *
		      (WITH3->genarray[g1 - 1] + WITH3->genarray[g2 - 1]) +
		    LINK->ps * (1 - LINK->pf) * (WITH3->genarray[g5 - 1] +
			WITH3->genarray[g6 - 1]) + LINK->pf * (1 - LINK->ps) *
		      (WITH3->genarray[g9 - 1] + WITH3->genarray[g10 - 1]) +
		    LINK->pf * LINK->ps *
		      (WITH3->genarray[g13 - 1] + WITH3->genarray[g14 - 1]));
	      if (val != 0.0) {
		for (l = 0; l < slength; l++)
		  temp[k][l] += temp2[l] * val;
	      }
	    }
	  } else {
	    g1 = WITH->genenumber[f1 - 1][s1 - 1];
	    g5 = WITH->genenumber[f1 - 1][ms1 - 1];
	    g9 = WITH->genenumber[mf1 - 1][s1 - 1];
	    g13 = WITH->genenumber[mf1 - 1][ms1 - 1];
	    FORLIM2 = nchild;
	    for (k = 0; k < FORLIM2; k++) {
	      WITH3 = thischild[k];
	      val = 4 * ((1 - LINK->ps) * (1 - LINK->pf) * WITH3->genarray[g1 - 1] +
			 LINK->ps * (1 - LINK->pf) * WITH3->genarray[g5 - 1] +
			 LINK->pf * (1 - LINK->ps) * WITH3->genarray[g9 - 1] +
			 LINK->pf * LINK->ps * WITH3->genarray[g13 - 1]);
	      if (val != 0.0) {
		for (l = 0; l < slength; l++)
		  temp[k][l] += temp2[l] * val;
	      }
	    }
	  }
	  LINK->secondseg++;
	}
      }
      LINK->firstseg++;
    }
  }
  temp1 = 0.0;
  for (l = 0; l < slength; l++) {
    val = 1.0;
    FORLIM1 = nchild;
    for (k = 0; k < FORLIM1; k++)
      val *= temp[k][l];
    temp1 += val;
  }
  return temp1;
}  /*msegfast*/


Local double msegfun(LINK)
struct LOC_seg *LINK;
{
  int g1, g2, g3, g4, g5, g6, g7, g8, g9, g10, g11, g12, g13, g14, g15, g16,
       i, j, k, f1, f2, s1, s2, ms1, ms2, mf1, mf2;
  double val, temp1, temp2;
  double temp[maxchild];
  int FORLIM;
  gennurec *WITH;
  thetavalues *WITH1, *WITH2;
  int FORLIM1, FORLIM2;
  thisarray *WITH3;

  LINK->firstseg = LINK->fseg;
  FORLIM = nchild;
  for (k = 0; k < FORLIM; k++)
    temp[k] = 0.0;
  WITH = gennustruct;
  WITH1 = LINK->firstsex;
  FORLIM = LINK->fend;
  for (i = LINK->fstart - 1; i < FORLIM; i++) {
    if (WITH1->segprob[i] == 0.0)
      LINK->firstseg++;
    else {
      temp1 = WITH1->segprob[i];
      f1 = invgenenum1[LINK->firstseg - 1];
      f2 = invgenenum2[LINK->firstseg - 1];
      mf1 = muthap[f1 - 1];
      mf2 = muthap[f2 - 1];
      LINK->secondseg = LINK->sseg;
      WITH2 = LINK->secondsex;
      if (f1 != f2) {
	FORLIM1 = LINK->send;
	for (j = LINK->sstart - 1; j < FORLIM1; j++) {
	  if (WITH2->segprob[j] == 0.0)
	    LINK->secondseg++;
	  else {
	    temp2 = temp1 * WITH2->segprob[j];
	    s1 = invgenenum1[LINK->secondseg - 1];
	    s2 = invgenenum2[LINK->secondseg - 1];
	    ms1 = muthap[s1 - 1];
	    ms2 = muthap[s2 - 1];
	    if (s1 != s2) {
	      g1 = WITH->genenumber[f1 - 1][s1 - 1];
	      g2 = WITH->genenumber[f1 - 1][s2 - 1];
	      g3 = WITH->genenumber[f2 - 1][s1 - 1];
	      g4 = WITH->genenumber[f2 - 1][s2 - 1];
	      g5 = WITH->genenumber[f1 - 1][ms1 - 1];
	      g6 = WITH->genenumber[f1 - 1][ms2 - 1];
	      g7 = WITH->genenumber[f2 - 1][ms1 - 1];
	      g8 = WITH->genenumber[f2 - 1][ms2 - 1];
	      g9 = WITH->genenumber[mf1 - 1][s1 - 1];
	      g10 = WITH->genenumber[mf1 - 1][s2 - 1];
	      g11 = WITH->genenumber[mf2 - 1][s1 - 1];
	      g12 = WITH->genenumber[mf2 - 1][s2 - 1];
	      g13 = WITH->genenumber[mf1 - 1][ms1 - 1];
	      g14 = WITH->genenumber[mf1 - 1][ms2 - 1];
	      g15 = WITH->genenumber[mf2 - 1][ms1 - 1];
	      g16 = WITH->genenumber[mf2 - 1][ms2 - 1];
	      FORLIM2 = nchild;
	      for (k = 0; k < FORLIM2; k++) {
		WITH3 = thischild[k];
		temp[k] += temp2 * ((1 - LINK->ps) * (1 - LINK->pf) *
			(WITH3->genarray[g1 - 1] + WITH3->genarray[g2 - 1] +
			  WITH3->genarray[g3 - 1] + WITH3->genarray[g4 - 1]) +
		      LINK->ps * (1 - LINK->pf) * (WITH3->genarray[g5 - 1] +
			  WITH3->genarray[g6 - 1] + WITH3->genarray[g7 - 1] +
			  WITH3->genarray[g8 - 1]) + LINK->pf *
			(1 - LINK->ps) * (WITH3->genarray[g9 - 1] + WITH3->
			    genarray[g10 - 1] + WITH3->genarray[g11 - 1] +
			  WITH3->genarray[g12 - 1]) + LINK->pf * LINK->ps *
			(WITH3->genarray[g13 - 1] + WITH3->genarray[g14 - 1] +
			  WITH3->genarray[g15 - 1] +
			  WITH3->genarray[g16 - 1]));
	      }
	    } else {
	      g1 = WITH->genenumber[f1 - 1][s1 - 1];
	      g3 = WITH->genenumber[f2 - 1][s1 - 1];
	      g5 = WITH->genenumber[f1 - 1][ms1 - 1];
	      g7 = WITH->genenumber[f2 - 1][ms1 - 1];
	      g9 = WITH->genenumber[mf1 - 1][s1 - 1];
	      g11 = WITH->genenumber[mf2 - 1][s1 - 1];
	      g13 = WITH->genenumber[mf1 - 1][ms1 - 1];
	      g15 = WITH->genenumber[mf2 - 1][ms1 - 1];
	      FORLIM2 = nchild;
	      for (k = 0; k < FORLIM2; k++) {
		WITH3 = thischild[k];
		temp[k] += 2 * temp2 * ((1 - LINK->ps) * (1 - LINK->pf) *
			(WITH3->genarray[g1 - 1] + WITH3->genarray[g3 - 1]) +
		      LINK->ps * (1 - LINK->pf) *
			(WITH3->genarray[g5 - 1] + WITH3->genarray[g7 - 1]) +
		      LINK->pf * (1 - LINK->ps) *
			(WITH3->genarray[g9 - 1] + WITH3->genarray[g11 - 1]) +
		      LINK->pf * LINK->ps * (WITH3->genarray[g13 - 1] +
			  WITH3->genarray[g15 - 1]));
	      }
	    }
	    LINK->secondseg++;
	  }
	}
      } else {
	FORLIM1 = LINK->send;
	for (j = LINK->sstart - 1; j < FORLIM1; j++) {
	  if (WITH2->segprob[j] == 0.0)
	    LINK->secondseg++;
	  else {
	    temp2 = temp1 * WITH2->segprob[j];
	    s1 = invgenenum1[LINK->secondseg - 1];
	    s2 = invgenenum2[LINK->secondseg - 1];
	    ms1 = muthap[s1 - 1];
	    ms2 = muthap[s2 - 1];
	    if (s1 != s2) {
	      g1 = WITH->genenumber[f1 - 1][s1 - 1];
	      g2 = WITH->genenumber[f1 - 1][s2 - 1];
	      g5 = WITH->genenumber[f1 - 1][ms1 - 1];
	      g6 = WITH->genenumber[f1 - 1][ms2 - 1];
	      g9 = WITH->genenumber[mf1 - 1][s1 - 1];
	      g10 = WITH->genenumber[mf1 - 1][s2 - 1];
	      g13 = WITH->genenumber[mf1 - 1][ms1 - 1];
	      g14 = WITH->genenumber[mf1 - 1][ms2 - 1];
	      FORLIM2 = nchild;
	      for (k = 0; k < FORLIM2; k++) {
		WITH3 = thischild[k];
		temp[k] += 2 * temp2 * ((1 - LINK->ps) * (1 - LINK->pf) *
			(WITH3->genarray[g1 - 1] + WITH3->genarray[g2 - 1]) +
		      LINK->ps * (1 - LINK->pf) *
			(WITH3->genarray[g5 - 1] + WITH3->genarray[g6 - 1]) +
		      LINK->pf * (1 - LINK->ps) *
			(WITH3->genarray[g9 - 1] + WITH3->genarray[g10 - 1]) +
		      LINK->pf * LINK->ps * (WITH3->genarray[g13 - 1] +
			  WITH3->genarray[g14 - 1]));
	      }
	    } else {
	      g1 = WITH->genenumber[f1 - 1][s1 - 1];
	      g5 = WITH->genenumber[f1 - 1][ms1 - 1];
	      g9 = WITH->genenumber[mf1 - 1][s1 - 1];
	      g13 = WITH->genenumber[mf1 - 1][ms1 - 1];
	      FORLIM2 = nchild;
	      for (k = 0; k < FORLIM2; k++) {
		WITH3 = thischild[k];
		temp[k] += 4 * temp2 * ((1 - LINK->ps) * (1 - LINK->pf) *
			WITH3->genarray[g1 - 1] +
		      LINK->ps * (1 - LINK->pf) * WITH3->genarray[g5 - 1] +
		      LINK->pf * (1 - LINK->ps) * WITH3->genarray[g9 - 1] +
		      LINK->pf * LINK->ps * WITH3->genarray[g13 - 1]);
	      }
	    }
	    LINK->secondseg++;
	  }
	}
      }
      LINK->firstseg++;
    }
  }
  val = 1.0;
  FORLIM = nchild;
  for (k = 0; k < FORLIM; k++)
    val *= temp[k];
  return val;
}  /*msegfun*/

Local double segfast(LINK)
struct LOC_seg *LINK;
{
  int g1, g2, g3, g4, i, j, k, l, f1, f2, s1, s2, slength;
  double val, temp1;
  double temp[maxchild][maxseg];
  double temp2[maxseg];
  int FORLIM, FORLIM1;
  gennurec *WITH;
  thetavalues *WITH1, *WITH2;
  int FORLIM2;
  thisarray *WITH3;

  LINK->firstseg = LINK->fseg;
  slength = LINK->send - LINK->sstart + 1;
  FORLIM = nchild;
  for (k = 0; k < FORLIM; k++) {
    for (l = 0; l < slength; l++)
      temp[k][l] = 0.0;
  }
  WITH = gennustruct;
  WITH1 = LINK->firstsex;
  FORLIM = LINK->fend;
  for (i = LINK->fstart - 1; i < FORLIM; i++) {
    if (WITH1->segprob[i] == 0.0)
      LINK->firstseg++;
    else {
      temp1 = WITH1->segprob[i];
      f1 = invgenenum1[LINK->firstseg - 1];
      f2 = invgenenum2[LINK->firstseg - 1];
      LINK->secondseg = LINK->sseg;
      WITH2 = LINK->secondsex;
      if (f1 != f2) {
	FORLIM1 = LINK->send;
	for (j = LINK->sstart - 1; j < FORLIM1; j++) {
	  for (l = 0; l < slength; l++)
	    temp2[l] = temp1 * WITH2->segprob[j + l * slength];
	  s1 = invgenenum1[LINK->secondseg - 1];
	  s2 = invgenenum2[LINK->secondseg - 1];
	  if (s1 != s2) {
	    g1 = WITH->genenumber[f1 - 1][s1 - 1];
	    g2 = WITH->genenumber[f1 - 1][s2 - 1];
	    g3 = WITH->genenumber[f2 - 1][s1 - 1];
	    g4 = WITH->genenumber[f2 - 1][s2 - 1];
	    FORLIM2 = nchild;
	    for (k = 0; k < FORLIM2; k++) {
	      WITH3 = thischild[k];
	      val = WITH3->genarray[g1 - 1] + WITH3->genarray[g2 - 1] +
		    WITH3->genarray[g3 - 1] + WITH3->genarray[g4 - 1];
	      if (val != 0.0) {
		for (l = 0; l < slength; l++)
		  temp[k][l] += temp2[l] * val;
	      }
	    }
	  } else {
	    g1 = WITH->genenumber[f1 - 1][s1 - 1];
	    g3 = WITH->genenumber[f2 - 1][s1 - 1];
	    FORLIM2 = nchild;
	    for (k = 0; k < FORLIM2; k++) {
	      WITH3 = thischild[k];
	      val = 2.0 * (WITH3->genarray[g1 - 1] + WITH3->genarray[g3 - 1]);
	      if (val != 0.0) {
		for (l = 0; l < slength; l++)
		  temp[k][l] += temp2[l] * val;
	      }
	    }
	  }
	  LINK->secondseg++;
	}
      } else {
	FORLIM1 = LINK->send;
	for (j = LINK->sstart - 1; j < FORLIM1; j++) {
	  for (l = 0; l < slength; l++)
	    temp2[l] = temp1 * WITH2->segprob[j + l * slength];
	  s1 = invgenenum1[LINK->secondseg - 1];
	  s2 = invgenenum2[LINK->secondseg - 1];
	  if (s1 != s2) {
	    g1 = WITH->genenumber[f1 - 1][s1 - 1];
	    g2 = WITH->genenumber[f1 - 1][s2 - 1];
	    FORLIM2 = nchild;
	    for (k = 0; k < FORLIM2; k++) {
	      WITH3 = thischild[k];
	      val = 2.0 * (WITH3->genarray[g1 - 1] + WITH3->genarray[g2 - 1]);
	      if (val != 0.0) {
		for (l = 0; l < slength; l++)
		  temp[k][l] += temp2[l] * val;
	      }
	    }
	  } else {
	    g1 = WITH->genenumber[f1 - 1][s1 - 1];
	    FORLIM2 = nchild;
	    for (k = 0; k < FORLIM2; k++) {
	      WITH3 = thischild[k];
	      val = 4.0 * WITH3->genarray[g1 - 1];
	      if (val != 0.0) {
		for (l = 0; l < slength; l++)
		  temp[k][l] += temp2[l] * val;
	      }
	    }
	  }
	  LINK->secondseg++;
	}
      }
      LINK->firstseg++;
    }
  }
  temp1 = 0.0;
  for (l = 0; l < slength; l++) {
    val = 1.0;
    FORLIM1 = nchild;
    for (k = 0; k < FORLIM1; k++)
      val *= temp[k][l];
    temp1 += val;
  }
  return temp1;
}  /*segfast*/

Void segsexctop(LINK)
struct LOC_seg *LINK;
{
  double segval;
  int first, second, sstop;
  boolean thiscensor, thisrare;
  censorrec *WITH;
  thisarray *WITH1;
  int FORLIM;
  thisarray *WITH2;

  initseg(LINK);
  WITH = censorstruct;
  if ((*LINK->p)->male) {
    WITH1 = (*LINK->p)->gen;
    FORLIM = LINK->nfirst;
    for (first = 0; first < FORLIM; first++) {
      if (WITH1->genarray[first] != 0.0) {
	segval = 0.0;
	LINK->fseg = first + 1;
	thisrare = rare[first];
	second = 1;
	if (thisrare) {
	  WITH2 = (*LINK->q)->gen;
	  do {
	    LINK->sstart = probstart[second - 1];
	    LINK->send = probend[second - 1];
	    LINK->sseg = segstart[second - 1];
	    sstop = second + LINK->send - LINK->sstart + 1;
	    if (thisc < maxcensor) {
	      thisc++;
	      thiscensor = WITH->censor[thisc - minint];
	    } else
	      thiscensor = (WITH2->genarray[second - 1] == 0.0) ||
			    (rare[second - 1] = true);
	    if (!thiscensor) {
	      if (mutsys != 0)
		segval += WITH2->genarray[second - 1] * msegsexf(LINK);
	      else
		segval += WITH2->genarray[second - 1] * segsexf(LINK);
	    }
	    second = sstop;
	  } while (second <= LINK->nsecond);
	} else {
	  WITH2 = (*LINK->q)->gen;
	  do {
	    LINK->sstart = probstart[second - 1];
	    LINK->send = probend[second - 1];
	    LINK->sseg = segstart[second - 1];
	    sstop = second + LINK->send - LINK->sstart + 1;
	    if (thisc < maxcensor) {
	      thisc++;
	      thiscensor = WITH->censor[thisc - minint];
	    } else
	      thiscensor = (WITH2->genarray[second - 1] == 0.0);
	    if (!thiscensor) {
	      if (mutsys != 0)
		segval += WITH2->genarray[second - 1] * msegsexf(LINK);
	      else
		segval += WITH2->genarray[second - 1] * segsexf(LINK);
	    }
	    second = sstop;
	  } while (second <= LINK->nsecond);
	}
	WITH1->genarray[first] *= segval * segscale;
      }
    }
  } else {
    WITH1 = (*LINK->p)->gen;
    FORLIM = LINK->nfirst;
    for (first = 0; first < FORLIM; first++) {
      if (WITH1->genarray[first] != 0.0) {
	segval = 0.0;
	LINK->fstart = probstart[first];
	LINK->fend = probend[first];
	LINK->fseg = segstart[first];
	thisrare = rare[first];
	second = 1;
	if (thisrare) {
	  WITH2 = (*LINK->q)->gen;
	  do {
	    LINK->sseg = second;
	    if (thisc < maxcensor) {
	      thisc++;
	      thiscensor = WITH->censor[thisc - minint];
	    } else
	      thiscensor = (WITH2->genarray[second - 1] == 0.0) ||
			   (rare[second - 1] = true);
	    if (!thiscensor) {
	      if (mutsys != 0)
		segval += WITH2->genarray[second - 1] * msegsex(LINK);
	      else
		segval += WITH2->genarray[second - 1] * segsex(LINK);
	    }
	    second++;
	  } while (second <= LINK->nsecond);
	} else {
	  WITH2 = (*LINK->q)->gen;
	  do {
	    LINK->sseg = second;
	    if (thisc < maxcensor) {
	      thisc++;
	      thiscensor = WITH->censor[thisc - minint];
	    } else
	      thiscensor = (WITH2->genarray[second - 1] == 0.0);
	    if (!thiscensor) {
	      if (mutsys != 0)
		segval += WITH2->genarray[second - 1] * msegsex(LINK);
	      else
		segval += WITH2->genarray[second - 1] * segsex(LINK);
	    }
	    second++;
	  } while (second <= LINK->nsecond);
	}
	WITH1->genarray[first] *= segval * segscale;
      }
    }
  }
  cleanup(LINK->q, LINK->LINK);
  exitseg(LINK);
}  /*segsexctop*/


Void segsextop(LINK)
struct LOC_seg *LINK;
{
  double segval, val;
  int first, second, sstop;
  censorrec *WITH;
  thisarray *WITH1;
  int FORLIM;
  thisarray *WITH2;

  initseg(LINK);
  WITH = censorstruct;
  if ((*LINK->p)->male) {
    WITH1 = (*LINK->p)->gen;
    FORLIM = LINK->nfirst;
    for (first = 1; first <= FORLIM; first++) {
      if (WITH1->genarray[first - 1] != 0.0) {
	segval = 0.0;
	LINK->fseg = first;
	second = 1;
	WITH2 = (*LINK->q)->gen;
	do {
	  LINK->sstart = probstart[second - 1];
	  LINK->send = probend[second - 1];
	  LINK->sseg = segstart[second - 1];
	  sstop = second + LINK->send - LINK->sstart + 1;
	  if (thisc <= maxcensor)
	    thisc++;
	  if (WITH2->genarray[second - 1] == 0.0) {
	    second = sstop;
	    if (thisc <= maxcensor)
	      WITH->censor[thisc - minint] = true;
	  } else {
	    if (mutsys != 0)
	      val = msegsexf(LINK);
	    else
	      val = segsexf(LINK);
	    if (thisc <= maxcensor) {
	      WITH->censor[thisc - minint] = (val==0.0);
	    }
	    segval += WITH2->genarray[second - 1] * val;
	    second = sstop;
	  }
	} while (second <= LINK->nsecond);
	WITH1->genarray[first - 1] *= segval * segscale;
      }
    }
  } else {
    WITH1 = (*LINK->p)->gen;
    FORLIM = LINK->nfirst;
    for (first = 0; first < FORLIM; first++) {
      if (WITH1->genarray[first] != 0.0) {
	segval = 0.0;
	LINK->fstart = probstart[first];
	LINK->fend = probend[first];
	LINK->fseg = segstart[first];
	second = 1;
	WITH2 = (*LINK->q)->gen;
	do {
	  LINK->sseg = second;
	  if (thisc <= maxcensor)
	    thisc++;
	  if (WITH2->genarray[second - 1] == 0.0) {
	    second++;
	    if (thisc <= maxcensor)
	      WITH->censor[thisc - minint] = true;
	  } else {
	    if (mutsys != 0)
	      val = msegsex(LINK);
	    else
	      val = segsex(LINK);
	    if (thisc <= maxcensor) {
	      WITH->censor[thisc - minint] = (val == 0.0);
	    }
	    segval += WITH2->genarray[second - 1] * val;
	    second++;
	  }
	} while (second <= LINK->nsecond);
	WITH1->genarray[first] *= segval * segscale;
      }
    }
  }
  cleanup(LINK->q, LINK->LINK);
  exitseg(LINK);
}  /*segsextop*/

Void segctop(LINK)
struct LOC_seg *LINK;
{
  double segval;
  int first, second, sstop;
  boolean thiscensor, thisrare;
  censorrec *WITH;
  thisarray *WITH1;
  int FORLIM;
  thisarray *WITH2;

  initseg(LINK);
  WITH = censorstruct;
  WITH1 = (*LINK->p)->gen;
  FORLIM = fgeno;
  for (first = 0; first < FORLIM; first++) {
    if (WITH1->genarray[first] != 0.0) {
      segval = 0.0;
      LINK->fstart = probstart[first];
      LINK->fend = probend[first];
      LINK->fseg = segstart[first];
      thisrare = rare[first];
      second = 1;
      if (thisrare) {
	WITH2 = (*LINK->q)->gen;
	do {
	  LINK->sstart = probstart[second - 1];
	  LINK->send = probend[second - 1];
	  LINK->sseg = segstart[second - 1];
	  sstop = second + LINK->send - LINK->sstart + 1;
	  if (thisc < maxcensor) {
	    thisc++;
	    thiscensor = WITH->censor[thisc - minint];
	  } else
	    thiscensor = (WITH2->genarray[second - 1] == 0.0) ||
			 rare[second - 1];
	  if (thiscensor)
	    second = sstop;
	  else {
	    if (mutsys != 0)
	      segval += WITH2->genarray[second - 1] * msegfast(LINK);
	    else
	      segval += WITH2->genarray[second - 1] * segfast(LINK);
	    second = sstop;
	  }
	} while (second <= fgeno);
      } else {
	WITH2 = (*LINK->q)->gen;
	do {
	  LINK->sstart = probstart[second - 1];
	  LINK->send = probend[second - 1];
	  LINK->sseg = segstart[second - 1];
	  sstop = second + LINK->send - LINK->sstart + 1;
	  if (thisc < maxcensor) {
	    thisc++;
	    thiscensor = WITH->censor[thisc - minint];
	  } else
	    thiscensor = (WITH2->genarray[second - 1] == 0.0);
	  if (thiscensor)
	    second = sstop;
	  else {
	    if (mutsys != 0)
	      segval += WITH2->genarray[second - 1] * msegfast(LINK);
	    else
	      segval += WITH2->genarray[second - 1] * segfast(LINK);
	    second = sstop;
	  }
	} while (second <= fgeno);
      }
      WITH1->genarray[first] *= segval * segscale;
    }
  }
  cleanup(LINK->q, LINK->LINK);
  exitseg(LINK);
  if (approximate && firstapprox && ((*LINK->p)->pa==NULL))
    getapprox(LINK);
}  /*segctop*/

Void segtop(LINK)
struct LOC_seg *LINK;
{
  double segval, val;
  int first, second, sstop;
  boolean thisrare;
  censorrec *WITH;
  thisarray *WITH1;
  int FORLIM;
  thisarray *WITH2;

  initseg(LINK);
  WITH = censorstruct;
  WITH1 = (*LINK->p)->gen;
  FORLIM = fgeno;
  for (first = 0; first < FORLIM; first++) {
    if (WITH1->genarray[first] != 0.0) {
      thisrare = rare[first];
      segval = 0.0;
      LINK->fstart = probstart[first];
      LINK->fend = probend[first];
      LINK->fseg = segstart[first];
      second = 1;
      if (thisrare) {
	WITH2 = (*LINK->q)->gen;
	do {
	  LINK->sstart = probstart[second - 1];
	  LINK->send = probend[second - 1];
	  LINK->sseg = segstart[second - 1];
	  sstop = second + LINK->send - LINK->sstart + 1;
	  if (thisc <= maxcensor)
	    thisc++;
	  if ((WITH2->genarray[second - 1] == 0.0) || (rare[second - 1])) {
	    second = sstop;
	    if (thisc <= maxcensor)
	      WITH->censor[thisc - minint] = true;
	  } else {
	    if (mutsys != 0)
	      val = msegfast(LINK);
	    else
	      val = segfast(LINK);
	    if (thisc <= maxcensor) {
	      WITH->censor[thisc - minint] = (val == 0.0);
	    }
	    segval += WITH2->genarray[second - 1] * val;
	    second = sstop;
	  }
	} while (second <= fgeno);
      } else {
	WITH2 = (*LINK->q)->gen;
	do {
	  LINK->sstart = probstart[second - 1];
	  LINK->send = probend[second - 1];
	  LINK->sseg = segstart[second - 1];
	  sstop = second + LINK->send - LINK->sstart + 1;
	  if (thisc <= maxcensor)
	    thisc++;
	  if (WITH2->genarray[second - 1] == 0.0) {
	    second = sstop;
	    if (thisc <= maxcensor)
	      WITH->censor[thisc - minint] = true;
	  } else {
	    if (mutsys != 0)
	      val = msegfast(LINK);
	    else
	      val = segfast(LINK);
	    if (thisc <= maxcensor) {
	      WITH->censor[thisc - minint] = (val == 0.0);
	    }
	    segval += WITH2->genarray[second - 1] * val;
	    second = sstop;
	  }
	} while (second <= fgeno);
      }
      WITH1->genarray[first] *= segval * segscale;
    }
  }
  cleanup(LINK->q, LINK->LINK);
  exitseg(LINK);
  if (approximate && firstapprox && ((*LINK->p)->pa==NULL))
    getapprox(LINK);
}  /*segtop*/


Void segcapprox(LINK)
struct LOC_seg *LINK;
{
  double segval;
  int first, second, sstop;
  boolean thiscensor;
  censorrec *WITH;
  approxrec *WITH1;
  thisarray *WITH2;
  int FORLIM;
  thisarray *WITH3;

  initseg(LINK);
  WITH = censorstruct;
  WITH1 = approxstruct;
  WITH2 = (*LINK->p)->gen;
  FORLIM = fgeno;
  for (first = 0; first < FORLIM; first++) {
    if (WITH2->genarray[first] != 0.0) {
      segval = 0.0;
      LINK->fstart = probstart[first];
      LINK->fend = probend[first];
      LINK->fseg = segstart[first];
      second = 1;
      WITH3 = (*LINK->q)->gen;
      do {
	LINK->sstart = probstart[second - 1];
	LINK->send = probend[second - 1];
	LINK->sseg = segstart[second - 1];
	sstop = second + LINK->send - LINK->sstart + 1;
	if (thisc < maxcensor) {
	  thisc++;
	  thiscensor = WITH->censor[thisc - minint];
	} else
	  thiscensor = (WITH3->genarray[second - 1] == 0.0);
	if (thiscensor || 
	    (!(WITH1->approxarray[LINK->LINK->thisped - 1][first])))
	  second = sstop;
	else {
	  if (mutsys != 0)
	    segval += WITH3->genarray[second - 1] * msegfast(LINK);
	  else
	    segval += WITH3->genarray[second - 1] * segfast(LINK);
	  second = sstop;
	}
      } while (second <= fgeno);
      WITH2->genarray[first] *= segval * segscale;
    }
  }
  cleanup(LINK->q, LINK->LINK);
  exitseg(LINK);
}  /*segcapprox*/

#include "oldsegup.c"

Void msegsexdown(LINK)
struct LOC_seg *LINK;
{
  short here;
  genotype gene;
  double val, temp2;
  int ms2, ms1, mf, j, first, second;
  short FORLIM;
  gennurec *WITH;
  censorrec *WITH1;
  thisarray *WITH2;
  int FORLIM1;
  thisarray *WITH3;
  int FORLIM2;
  thetavalues *WITH4;
  int FORLIM3;

  initseg(LINK);
  FORLIM = fgeno;
  for (here = 0; here < FORLIM; here++)
    gene[here] = 0.0;
  WITH = gennustruct;
  WITH1 = censorstruct;
  WITH2 = (*LINK->p)->gen;
  FORLIM1 = LINK->nfirst;
  for (first = 0; first < FORLIM1; first++) {
    if (WITH2->genarray[first] != 0.0) {
      LINK->fseg = first + 1;
      second = 1;
      WITH3 = (*LINK->q)->gen;
      FORLIM2 = LINK->nsecond;
      for (second = 0; second < FORLIM2; second++) {
	if (WITH3->genarray[second] != 0.0) {
	  val = WITH3->genarray[second] * (*LINK->p)->gen->genarray[first] *
		segscale;
	  LINK->sstart = probstart[second];
	  LINK->send = probend[second];
	  LINK->sseg = segstart[second];
	  if (nchild != 0)
	    val *= msegsex(LINK);
	  if (val != 0.0) {
	    mf = muthap[invgenenum1[LINK->fseg - 1] - 1];
	    LINK->secondseg = LINK->sseg;
	    WITH4 = femaletheta;
	    FORLIM3 = LINK->send;
	    for (j = LINK->sstart - 1; j < FORLIM3; j++) {
	      ms1 = muthap[invgenenum1[LINK->secondseg - 1] - 1];
	      ms2 = muthap[invgenenum2[LINK->secondseg - 1] - 1];
	      temp2 = WITH4->segprob[j];
	      if (temp2 != 0.0) {
		if ((*LINK->r)->male) {
		  here = invgenenum1[LINK->secondseg - 1];
		  gene[here - 1] += (1 - LINK->ps) * temp2 * val;
		  here = invgenenum2[LINK->secondseg - 1];
		  gene[here - 1] += (1 - LINK->ps) * temp2 * val;
		  here = ms1;
		  gene[here - 1] += LINK->ps * temp2 * val;
		  here = ms2;
		  gene[here - 1] += LINK->ps * temp2 * val;
		} else {
		  here = WITH->genenumber[invgenenum1[LINK->secondseg - 1] - 1]
		    [first];
		  gene[here - 1] += (1 - LINK->pf) * (1 - LINK->ps) * temp2 * val;
		  here = WITH->genenumber[invgenenum2[LINK->secondseg - 1] - 1]
		    [first];
		  gene[here - 1] += (1 - LINK->pf) * (1 - LINK->ps) * temp2 * val;
		  here = WITH->genenumber[invgenenum1[LINK->secondseg - 1] - 1]
		    [mf - 1];
		  gene[here - 1] += LINK->pf * (1 - LINK->ps) * temp2 * val;
		  here = WITH->genenumber[invgenenum2[LINK->secondseg - 1] - 1]
		    [mf - 1];
		  gene[here - 1] += LINK->pf * (1 - LINK->ps) * temp2 * val;
		  here = WITH->genenumber[ms1 - 1][first];
		  gene[here - 1] += (1 - LINK->pf) * LINK->ps * temp2 * val;
		  here = WITH->genenumber[ms2 - 1][first];
		  gene[here - 1] += (1 - LINK->pf) * LINK->ps * temp2 * val;
		  here = WITH->genenumber[ms1 - 1][mf - 1];
		  gene[here - 1] += LINK->pf * LINK->ps * temp2 * val;
		  here = WITH->genenumber[ms2 - 1][mf - 1];
		  gene[here - 1] += LINK->pf * LINK->ps * temp2 * val;
		}
	      }
	      LINK->secondseg++;
	    }
	  }
	}  /*second*/
      }
    }  /*first*/
  }
  memcpy((*LINK->p)->gen->genarray, (*LINK->r)->gen->genarray,
	 sizeof(genotype));
  memcpy((*LINK->r)->gen->genarray, gene, sizeof(genotype));
  WITH2 = (*LINK->r)->gen;
  FORLIM1 = fgeno;
  for (first = 0; first < FORLIM1; first++)
    WITH2->genarray[first] *= (*LINK->p)->gen->genarray[first];
  cleanup(LINK->p, LINK->LINK);
  cleanup(LINK->q, LINK->LINK);
  exitseg(LINK);
}  /*msegsexdown*/


Void msegdown(LINK)
struct LOC_seg *LINK;
{
  short here;
  genotype gene;
  double val, temp, temp1, temp2;
  int i, j, first, second, f1, f2, s1, s2, mf1, mf2, ms1, ms2;
  short FORLIM;
  gennurec *WITH;
  censorrec *WITH1;
  thisarray *WITH2;
  int FORLIM1;
  thisarray *WITH3;
  int FORLIM2;
  thetavalues *WITH4;
  int FORLIM3;
  thetavalues *WITH5;
  int FORLIM4;

  initseg(LINK);
  FORLIM = fgeno;
  for (here = 0; here < FORLIM; here++)
    gene[here] = 0.0;
  WITH = gennustruct;
  WITH1 = censorstruct;
  WITH2 = (*LINK->p)->gen;
  FORLIM1 = fgeno;
  for (first = 0; first < FORLIM1; first++) {
    if (WITH2->genarray[first] != 0.0) {
      LINK->fstart = probstart[first];
      LINK->fend = probend[first];
      LINK->fseg = segstart[first];
      WITH3 = (*LINK->q)->gen;
      FORLIM2 = fgeno;
      for (second = 0; second < FORLIM2; second++) {
	if (WITH3->genarray[second] != 0.0) {
	  LINK->sstart = probstart[second];
	  LINK->send = probend[second];
	  LINK->sseg = segstart[second];
	  val = WITH3->genarray[second] * (*LINK->p)->gen->genarray[first] *
		segscale;
	  if (nchild != 0)
	    val = msegfun(LINK) * val;
	  if (val != 0.0) {
	    LINK->firstseg = LINK->fseg;
	    WITH4 = maletheta;
	    FORLIM3 = LINK->fend;
	    for (i = LINK->fstart - 1; i < FORLIM3; i++) {
	      temp1 = WITH4->segprob[i];
	      LINK->secondseg = LINK->sseg;
	      if (temp1 != 0.0) {
		WITH5 = femaletheta;
		FORLIM4 = LINK->send;
		for (j = LINK->sstart - 1; j < FORLIM4; j++) {
		  if (WITH5->segprob[j] == 0.0)
		    LINK->secondseg++;
		  else {
		    temp2 = WITH5->segprob[j];
		    f1 = invgenenum1[LINK->firstseg - 1];
		    f2 = invgenenum2[LINK->firstseg - 1];
		    s1 = invgenenum1[LINK->secondseg - 1];
		    s2 = invgenenum2[LINK->secondseg - 1];
		    temp = (1 - LINK->pf) * (1 - LINK->ps) * temp1 * temp2 * val;
		    here = WITH->genenumber[s1 - 1][f1 - 1];
		    gene[here - 1] += temp;
		    here = WITH->genenumber[s1 - 1][f2 - 1];
		    gene[here - 1] += temp;
		    here = WITH->genenumber[s2 - 1][f1 - 1];
		    gene[here - 1] += temp;
		    here = WITH->genenumber[s2 - 1][f2 - 1];
		    gene[here - 1] += temp;
		    ms1 = muthap[s1 - 1];
		    ms2 = muthap[s2 - 1];
		    temp = (1 - LINK->pf) * LINK->ps * temp1 * temp2 * val;
		    here = WITH->genenumber[ms1 - 1][f1 - 1];
		    gene[here - 1] += temp;
		    here = WITH->genenumber[ms1 - 1][f2 - 1];
		    gene[here - 1] += temp;
		    here = WITH->genenumber[ms2 - 1][f1 - 1];
		    gene[here - 1] += temp;
		    here = WITH->genenumber[ms2 - 1][f2 - 1];
		    gene[here - 1] += temp;
		    mf1 = muthap[f1 - 1];
		    mf2 = muthap[f2 - 1];
		    temp = LINK->pf * (1 - LINK->ps) * temp1 * temp2 * val;
		    here = WITH->genenumber[mf1 - 1][s1 - 1];
		    gene[here - 1] += temp;
		    here = WITH->genenumber[mf1 - 1][s2 - 1];
		    gene[here - 1] += temp;
		    here = WITH->genenumber[mf2 - 1][s1 - 1];
		    gene[here - 1] += temp;
		    here = WITH->genenumber[mf2 - 1][s2 - 1];
		    gene[here - 1] += temp;
		    temp = LINK->pf * LINK->ps * temp1 * temp2 * val;
		    here = WITH->genenumber[mf1 - 1][ms1 - 1];
		    gene[here - 1] += temp;
		    here = WITH->genenumber[mf1 - 1][ms2 - 1];
		    gene[here - 1] += temp;
		    here = WITH->genenumber[mf2 - 1][ms1 - 1];
		    gene[here - 1] += temp;
		    here = WITH->genenumber[mf2 - 1][ms2 - 1];
		    gene[here - 1] += temp;
		    LINK->secondseg++;
		  }
		}
	      }
	      LINK->firstseg++;
	    }
	  }
	}  /*second*/
      }
    }  /*first*/
  }
  memcpy((*LINK->p)->gen->genarray, (*LINK->r)->gen->genarray,
	 sizeof(genotype));
  memcpy((*LINK->r)->gen->genarray, gene, sizeof(genotype));
  WITH2 = (*LINK->r)->gen;
  FORLIM1 = fgeno;
  for (first = 0; first < FORLIM1; first++)
    WITH2->genarray[first] *= (*LINK->p)->gen->genarray[first];
  cleanup(LINK->p, LINK->LINK);
  cleanup(LINK->q, LINK->LINK);
  exitseg(LINK);


}  /*msegdown*/