/* Array mapping functions for APL-Scheme.
   Copyright (C) 1994 Radey Shouman.
   
   This program is free software; you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation; either version 1, or (at your option)
   any later version.
   
   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.
   
   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the Free Software
   Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
   */

#include "scm.h"

#ifdef ARRAYS

typedef struct
{
  char *name;
  SCM sproc;
  int (*vproc) ();
}
ra_iproc;

#define BVE_REF(a, i) ((VELTS(a)[(i)/LONG_BIT] & (1L<<((i)%LONG_BIT))) ? 1 : 0)
#define BVE_SET(a, i) (VELTS(a)[(i)/LONG_BIT] |= (1L<<((i)%LONG_BIT)))
#define BVE_CLR(a, i) (VELTS(a)[(i)/LONG_BIT] &= ~(1L<<((i)%LONG_BIT)))
/* Fast, recycling scm_vector ref */
#define RVREF(ra, i, e) (e = scm_cvref(ra, i, e))
/* #define RVREF(ra, i, e) (scm_cvref(ra, i, SCM_UNDEFINED)) to turn off */

/* IVDEP means "ignore scm_vector dependencies", meaning we guarantee that
   elements of scm_vector operands are not aliased */
#ifdef _UNICOS
#define IVDEP(test, line) if (test) {_Pragma("ivdep"); line} else {line}
#else
#define IVDEP(test, line) line
#endif

/* inds must be a uvect or ivect, no check. */
static sizet 
cind (ra, inds)
     SCM ra, inds;
{
  sizet i;
  int k;
  long *ve = VELTS (inds);
  if (!ARRAYP (ra))
    return *ve;
  i = ARRAY_BASE (ra);
  for (k = 0; k < ARRAY_NDIM (ra); k++)
    i += (ve[k] - ARRAY_DIMS (ra)[k].lbnd) * ARRAY_DIMS (ra)[k].inc;
  return i;
}

/* Checker for scm_array mapping functions:
   return values: 4 --> shapes, increments, and bases are the same;
   3 --> shapes and increments are the same;
   2 --> shapes are the same;
   1 --> ras are at least as big as ra0;
   0 --> no match.
   */
int 
scm_ra_matchp (ra0, ras)
     SCM ra0, ras;
{
  SCM ra1;
  scm_array_dim dims;
  scm_array_dim *s0 = &dims;
  scm_array_dim *s1;
  sizet bas0 = 0;
  int i, ndim = 1;
  int exact = 2			/* 4 */ ;	/* Don't care about values >2 (yet?) */
  if IMP
    (ra0) return 0;
  switch TYP7
    (ra0)
      {
      default:
	return 0;
      case tc7_vector:
      case tc7_string:
      case tc7_bvect:
      case tc7_uvect:
      case tc7_ivect:
      case tc7_fvect:
      case tc7_dvect:
      case tc7_cvect:
	s0->lbnd = 0;
	s0->inc = 1;
	s0->ubnd = (long) LENGTH (ra0) - 1;
	break;
      case tc7_smob:
	if (!ARRAYP (ra0))
	  return 0;
	ndim = ARRAY_NDIM (ra0);
	s0 = ARRAY_DIMS (ra0);
	bas0 = ARRAY_BASE (ra0);
	break;
      }
  while NIMP
    (ras)
      {
	ra1 = CAR (ras);
	if IMP
	  (ra1) return 0;
	switch TYP7
	  (ra1)
	    {
	    default:
	      return 0;
	    case tc7_vector:
	    case tc7_string:
	    case tc7_bvect:
	    case tc7_uvect:
	    case tc7_ivect:
	    case tc7_fvect:
	    case tc7_dvect:
	    case tc7_cvect:
	      if (1 != ndim)
		return 0;
	      switch (exact)
		{
		case 4:
		  if (0 != bas0)
		    exact = 3;
		case 3:
		  if (1 != s0->inc)
		    exact = 2;
		case 2:
		  if ((0 == s0->lbnd) && (s0->ubnd == LENGTH (ra1) - 1))
		    break;
		  exact = 1;
		case 1:
		  if (s0->lbnd < 0 || s0->ubnd >= LENGTH (ra1))
		    return 0;
		}
	      break;
	    case tc7_smob:
	      if (!ARRAYP (ra1) || ndim != ARRAY_NDIM (ra1))
		return 0;
	      s1 = ARRAY_DIMS (ra1);
	      if (bas0 != ARRAY_BASE (ra1))
		exact = 3;
	      for (i = 0; i < ndim; i++)
		switch (exact)
		  {
		  case 4:
		  case 3:
		    if (s0[i].inc != s1[i].inc)
		      exact = 2;
		  case 2:
		    if (s0[i].lbnd == s1[i].lbnd && s0[i].ubnd == s1[i].ubnd)
		      break;
		    exact = 1;
		  default:
		    if (s0[i].lbnd < s1[i].lbnd || s0[i].ubnd > s1[i].ubnd)
		      return (s0[i].lbnd <= s0[i].ubnd ? 0 : 1);
		  }
	      break;
	    }
	ras = CDR (ras);
      }
  return exact;
}

static char s_ra_mismatch[] = "array shape mismatch";
int 
scm_ramapc (cproc, data, ra0, lra, what)
     int (*cproc) ();
     SCM data, ra0, lra;
     char *what;
{
  SCM inds, z;
  SCM vra0, ra1, vra1;
  SCM lvra, *plvra;
  long *vinds;
  int k, kmax;
  switch (scm_ra_matchp (ra0, lra))
    {
    default:
    case 0:
      scm_wta (ra0, s_ra_mismatch, what);
    case 2:
    case 3:
    case 4:			/* Try unrolling arrays */
      kmax = (ARRAYP (ra0) ? ARRAY_NDIM (ra0) - 1 : 0);
      if (kmax < 0)
	goto gencase;
      vra0 = scm_array_contents (ra0, SCM_UNDEFINED);
      if IMP
	(vra0) goto gencase;
      if (!ARRAYP (vra0))
	{
	  vra1 = scm_make_ra (1);
	  ARRAY_BASE (vra1) = 0;
	  ARRAY_DIMS (vra1)->lbnd = 0;
	  ARRAY_DIMS (vra1)->ubnd = LENGTH (vra0) - 1;
	  ARRAY_DIMS (vra1)->inc = 1;
	  ARRAY_V (vra1) = vra0;
	  vra0 = vra1;
	}
      lvra = EOL;
      plvra = &lvra;
      for (z = lra; NIMP (z); z = CDR (z))
	{
	  ra1 = CAR (z);
	  vra1 = scm_make_ra (1);
	  ARRAY_DIMS (vra1)->lbnd = ARRAY_DIMS (vra0)->lbnd;
	  ARRAY_DIMS (vra1)->ubnd = ARRAY_DIMS (vra0)->ubnd;
	  if (!ARRAYP (ra1))
	    {
	      ARRAY_BASE (vra1) = 0;
	      ARRAY_DIMS (vra1)->inc = 1;
	      ARRAY_V (vra1) = ra1;
	    }
	  else if (!ARRAY_CONTP (ra1))
	    goto gencase;
	  else
	    {
	      ARRAY_BASE (vra1) = ARRAY_BASE (ra1);
	      ARRAY_DIMS (vra1)->inc = ARRAY_DIMS (ra1)[kmax].inc;
	      ARRAY_V (vra1) = ARRAY_V (ra1);
	    }
	  *plvra = scm_cons (vra1, EOL);
	  plvra = &CDR (*plvra);
	}
      return (UNBNDP (data) ? cproc (vra0, lvra) : cproc (vra0, data, lvra));
    case 1:
    gencase:			/* Have to loop over all dimensions. */
      vra0 = scm_make_ra (1);
      if ARRAYP
	(ra0)
	  {
	    kmax = ARRAY_NDIM (ra0) - 1;
	    if (kmax < 0)
	      {
		ARRAY_DIMS (vra0)->lbnd = 0;
		ARRAY_DIMS (vra0)->ubnd = 0;
		ARRAY_DIMS (vra0)->inc = 1;
	      }
	    else
	      {
		ARRAY_DIMS (vra0)->lbnd = ARRAY_DIMS (ra0)[kmax].lbnd;
		ARRAY_DIMS (vra0)->ubnd = ARRAY_DIMS (ra0)[kmax].ubnd;
		ARRAY_DIMS (vra0)->inc = ARRAY_DIMS (ra0)[kmax].inc;
	      }
	    ARRAY_BASE (vra0) = ARRAY_BASE (ra0);
	    ARRAY_V (vra0) = ARRAY_V (ra0);
	  }
      else
	{
	  kmax = 0;
	  ARRAY_DIMS (vra0)->lbnd = 0;
	  ARRAY_DIMS (vra0)->ubnd = LENGTH (ra0) - 1;
	  ARRAY_DIMS (vra0)->inc = 1;
	  ARRAY_BASE (vra0) = 0;
	  ARRAY_V (vra0) = ra0;
	  ra0 = vra0;
	}
      lvra = EOL;
      plvra = &lvra;
      for (z = lra; NIMP (z); z = CDR (z))
	{
	  ra1 = CAR (z);
	  vra1 = scm_make_ra (1);
	  ARRAY_DIMS (vra1)->lbnd = ARRAY_DIMS (vra0)->lbnd;
	  ARRAY_DIMS (vra1)->ubnd = ARRAY_DIMS (vra0)->ubnd;
	  if ARRAYP
	    (ra1)
	      {
		if (kmax >= 0)
		  ARRAY_DIMS (vra1)->inc = ARRAY_DIMS (ra1)[kmax].inc;
		ARRAY_V (vra1) = ARRAY_V (ra1);
	      }
	  else
	    {
	      ARRAY_DIMS (vra1)->inc = 1;
	      ARRAY_V (vra1) = ra1;
	    }
	  *plvra = scm_cons (vra1, EOL);
	  plvra = &CDR (*plvra);
	}
      inds = scm_make_uve (ARRAY_NDIM (ra0), MAKINUM (-1L));
      vinds = (long *) VELTS (inds);
      for (k = 0; k <= kmax; k++)
	vinds[k] = ARRAY_DIMS (ra0)[k].lbnd;
      k = kmax;
      do
	{
	  if (k == kmax)
	    {
	      SCM y = lra;
	      ARRAY_BASE (vra0) = cind (ra0, inds);
	      for (z = lvra; NIMP (z); z = CDR (z), y = CDR (y))
		ARRAY_BASE (CAR (z)) = cind (CAR (y), inds);
	      if (0 == (UNBNDP (data) ? cproc (vra0, lvra) : cproc (vra0, data, lvra)))
		return 0;
	      k--;
	      continue;
	    }
	  if (vinds[k] < ARRAY_DIMS (ra0)[k].ubnd)
	    {
	      vinds[k]++;
	      k++;
	      continue;
	    }
	  vinds[k] = ARRAY_DIMS (ra0)[k].lbnd - 1;
	  k--;
	}
      while (k >= 0);
      return 1;
    }
}

static char s_array_fill[] = "array-fill!";
int 
scm_rafill (ra, fill, ignore)
     SCM ra, fill, ignore;
{
  sizet i, n = ARRAY_DIMS (ra)->ubnd - ARRAY_DIMS (ra)->lbnd + 1;
  long inc = ARRAY_DIMS (ra)->inc;
  sizet base = ARRAY_BASE (ra);
  ra = ARRAY_V (ra);
  switch TYP7
    (ra)
      {
      default:
	for (i = base; n--; i += inc)
	  scm_aset (ra, fill, MAKINUM (i));
	break;
      case tc7_vector:
	for (i = base; n--; i += inc)
	  VELTS (ra)[i] = fill;
	break;
      case tc7_string:
	ASRTGO (ICHRP (fill), badarg2);
	for (i = base; n--; i += inc)
	  CHARS (ra)[i] = ICHR (fill);
	break;
      case tc7_bvect:
	{
	  long *ve = (long *) VELTS (ra);
	  if (1 == inc && (n >= LONG_BIT || n == LENGTH (ra)))
	    {
	      i = base / LONG_BIT;
	      if (BOOL_F == fill)
		{
		  if (base % LONG_BIT) /* leading partial word */
		    ve[i++] &= ~(~0L << (base % LONG_BIT));
		  for (; i < (base + n) / LONG_BIT; i++)
		    ve[i] = 0L;
		  if ((base + n) % LONG_BIT) /* trailing partial word */
		    ve[i] &= (~0L << ((base + n) % LONG_BIT));
		}
	      else if (BOOL_T == fill)
		{
		  if (base % LONG_BIT)
		    ve[i++] |= ~0L << (base % LONG_BIT);
		  for (; i < (base + n) / LONG_BIT; i++)
		    ve[i] = ~0L;
		  if ((base + n) % LONG_BIT)
		    ve[i] |= ~(~0L << ((base + n) % LONG_BIT));
		}
	      else
	      badarg2:scm_wta (fill, (char *) ARG2, s_array_fill);
	    }
	  else
	    {
	      if (BOOL_F == fill)
		for (i = base; n--; i += inc)
		  ve[i / LONG_BIT] &= ~(1L << (i % LONG_BIT));
	      else if (BOOL_T == fill)
		for (i = base; n--; i += inc)
		  ve[i / LONG_BIT] |= (1L << (i % LONG_BIT));
	      else
		goto badarg2;
	    }
	  break;
	}
      case tc7_uvect:
	ASRTGO (0 <= INUM (fill), badarg2);
      case tc7_ivect:
	ASRTGO (INUMP (fill), badarg2);
	{
	  long f = INUM (fill), *ve = (long *) VELTS (ra);
	  for (i = base; n--; i += inc)
	    ve[i] = f;
	  break;
	}
#ifdef FLOATS
#ifdef SINGLES
      case tc7_fvect:
	{
	  float f, *ve = (float *) VELTS (ra);
	  ASRTGO (NIMP (fill) && REALP (fill), badarg2);
	  f = REALPART (fill);
	  for (i = base; n--; i += inc)
	    ve[i] = f;
	  break;
	}
#endif /* SINGLES */
      case tc7_dvect:
	{
	  double f, *ve = (double *) VELTS (ra);
	  ASRTGO (NIMP (fill) && REALP (fill), badarg2);
	  f = REALPART (fill);
	  for (i = base; n--; i += inc)
	    ve[i] = f;
	  break;
	}
      case tc7_cvect:
	{
	  double fr, fi;
	  double (*ve)[2] = (double (*)[2]) VELTS (ra);
	  ASRTGO (NIMP (fill) && INEXP (fill), badarg2);
	  fr = REALPART (fill);
	  fi = (CPLXP (fill) ? IMAG (fill) : 0.0);
	  for (i = base; n--; i += inc)
	    {
	      ve[i][0] = fr;
	      ve[i][1] = fi;
	    }
	  break;
	}
#endif /* FLOATS */
      }
  return 1;
}
SCM 
scm_array_fill (ra, fill)
     SCM ra, fill;
{
  scm_ramapc (scm_rafill, fill, ra, EOL, s_array_fill);
  return UNSPECIFIED;
}

static char s_sarray_copy[] = "serial-array-copy!";
#define s_array_copy  (s_sarray_copy + 7)
static int 
racp (src, dst)
     SCM dst, src;
{
  long n = (ARRAY_DIMS (src)->ubnd - ARRAY_DIMS (src)->lbnd + 1);
  long inc_d, inc_s = ARRAY_DIMS (src)->inc;
  sizet i_d, i_s = ARRAY_BASE (src);
  dst = CAR (dst);
  inc_d = ARRAY_DIMS (dst)->inc;
  i_d = ARRAY_BASE (dst);
  src = ARRAY_V (src);
  dst = ARRAY_V (dst);
  switch TYP7
    (dst)
      {
      default:
      gencase: case tc7_vector:
	for (; n-- > 0; i_s += inc_s, i_d += inc_d)
	  scm_aset (dst, scm_cvref (src, i_s, SCM_UNDEFINED), MAKINUM (i_d));
	break;
      case tc7_string:
	if (tc7_string != TYP7 (dst))
	  goto gencase;
	for (; n-- > 0; i_s += inc_s, i_d += inc_d)
	  CHARS (dst)[i_d] = CHARS (src)[i_s];
	break;
      case tc7_bvect:
	if (tc7_bvect != TYP7 (dst))
	  goto gencase;
	if (1 == inc_d && 1 == inc_s && i_s % LONG_BIT == i_d % LONG_BIT && n >= LONG_BIT)
	  {
	    long *sv = (long *) VELTS (src);
	    long *dv = (long *) VELTS (dst);
	    sv += i_s / LONG_BIT;
	    dv += i_d / LONG_BIT;
	    if (i_s % LONG_BIT)
	      {			/* leading partial word */
		*dv = (*dv & ~(~0L << (i_s % LONG_BIT))) | (*sv & (~0L << (i_s % LONG_BIT)));
		dv++;
		sv++;
		n -= LONG_BIT - (i_s % LONG_BIT);
	      }
	    IVDEP (src != dst,
		   for (; n >= LONG_BIT; n -= LONG_BIT, sv++, dv++)
		   * dv = *sv;)
	      if (n)		/* trailing partial word */
		*dv = (*dv & (~0L << n)) | (*sv & ~(~0L << n));
	  }
	else
	  {
	    for (; n-- > 0; i_s += inc_s, i_d += inc_d)
	      if (VELTS (src)[i_s / LONG_BIT] & (1L << (i_s % LONG_BIT)))
		VELTS (dst)[i_d / LONG_BIT] |= (1L << (i_d % LONG_BIT));
	      else
		VELTS (dst)[i_d / LONG_BIT] &= ~(1L << (i_d % LONG_BIT));
	  }
	break;
      case tc7_uvect:
	if (tc7_uvect != TYP7 (src))
	  goto gencase;
	else
	  {
	    long *d = (long *) VELTS (dst), *s = (long *) VELTS (src);
	    IVDEP (src != dst,
		   for (; n-- > 0; i_s += inc_s, i_d += inc_d)
		   d[i_d] = s[i_s];)
	      break;
	  }
      case tc7_ivect:
	if (tc7_uvect != TYP7 (src) && tc7_ivect != TYP7 (src))
	  goto gencase;
	else
	  {
	    long *d = (long *) VELTS (dst), *s = (long *) VELTS (src);
	    IVDEP (src != dst,
		   for (; n-- > 0; i_s += inc_s, i_d += inc_d)
		   d[i_d] = s[i_s];)
	      break;
	  }
#ifdef FLOATS
#ifdef SINGLES
      case tc7_fvect:
	{
	  float *d = (float *) VELTS (dst);
	  float *s = (float *) VELTS (src);
	  switch TYP7
	    (src)
	      {
	      default:
		goto gencase;
	      case tc7_ivect:
	      case tc7_uvect:
		IVDEP (src != dst,
		       for (; n-- > 0; i_s += inc_s, i_d += inc_d)
		       d[i_d] = ((long *) s)[i_s];)
		  break;
	      case tc7_fvect:
		IVDEP (src != dst,
		       for (; n-- > 0; i_s += inc_s, i_d += inc_d)
		       d[i_d] = s[i_s];)
		  break;
	      case tc7_dvect:
		IVDEP (src != dst,
		       for (; n-- > 0; i_s += inc_s, i_d += inc_d)
		       d[i_d] = ((double *) s)[i_s];)
		  break;
	      }
	  break;
	}
#endif /* SINGLES */
      case tc7_dvect:
	{
	  double *d = (double *) VELTS (dst);
	  double *s = (double *) VELTS (src);
	  switch TYP7
	    (src)
	      {
	      default:
		goto gencase;
	      case tc7_ivect:
	      case tc7_uvect:
		IVDEP (src != dst,
		       for (; n-- > 0; i_s += inc_s, i_d += inc_d)
		       d[i_d] = ((long *) s)[i_s];)
		  break;
	      case tc7_fvect:
		IVDEP (src != dst,
		       for (; n-- > 0; i_s += inc_s, i_d += inc_d)
		       d[i_d] = ((float *) s)[i_s];)
		  break;
	      case tc7_dvect:
		IVDEP (src != dst,
		       for (; n-- > 0; i_s += inc_s, i_d += inc_d)
		       d[i_d] = s[i_s];)
		  break;
	      }
	  break;
	}
      case tc7_cvect:
	{
	  double (*d)[2] = (double (*)[2]) VELTS (dst);
	  double (*s)[2] = (double (*)[2]) VELTS (src);
	  switch TYP7
	    (src)
	      {
	      default:
		goto gencase;
	      case tc7_ivect:
	      case tc7_uvect:
		IVDEP (src != dst,
		       for (; n-- > 0; i_s += inc_s, i_d += inc_d)
		       {
			 d[i_d][0] = ((long *) s)[i_s];
			 d[i_d][1] = 0.0;
		       }
		       )
		  break;
	      case tc7_fvect:
		IVDEP (src != dst,
		       for (; n-- > 0; i_s += inc_s, i_d += inc_d)
		       {
			 d[i_d][0] = ((float *) s)[i_s];
			 d[i_d][1] = 0.0;
		       }
		       )
		  break;
	      case tc7_dvect:
		IVDEP (src != dst,
		       for (; n-- > 0; i_s += inc_s, i_d += inc_d)
		       {
			 d[i_d][0] = ((double *) s)[i_s];
			 d[i_d][1] = 0.0;
		       }
		       )
		  break;
	      case tc7_cvect:
		IVDEP (src != dst,
		       for (; n-- > 0; i_s += inc_s, i_d += inc_d)
		       {
			 d[i_d][0] = s[i_s][0];
			 d[i_d][1] = s[i_s][1];
		       }
		       )
	      }
	  break;
	}
      }
#endif /* FLOATS */
  return 1;
}
SCM scm_array_copy (src, dst)
     SCM src;
     SCM dst;
{
  scm_ramapc (racp, SCM_UNDEFINED, src, scm_cons (dst, EOL), s_array_copy);
  return SCM_UNDEFINED;
}

/* Functions callable by ARRAY-MAP! */
int scm_ra_eqp (ra0, ras)
     SCM ra0, ras;
{
  SCM ra1 = CAR (ras), ra2 = CAR (CDR (ras));
  long n = ARRAY_DIMS (ra0)->ubnd - ARRAY_DIMS (ra0)->lbnd + 1;
  sizet i0 = ARRAY_BASE (ra0), i1 = ARRAY_BASE (ra1), i2 = ARRAY_BASE (ra2);
  long inc0 = ARRAY_DIMS (ra0)->inc;
  long inc1 = ARRAY_DIMS (ra1)->inc;
  long inc2 = ARRAY_DIMS (ra1)->inc;
  ra0 = ARRAY_V (ra0);
  ra1 = ARRAY_V (ra1);
  ra2 = ARRAY_V (ra2);
  switch (TYP7 (ra1) == TYP7 (ra2) ? TYP7 (ra1) : 0)
    {
    default:
      {
	SCM e1 = SCM_UNDEFINED, e2 = SCM_UNDEFINED;
	for (; n-- > 0; i0 += inc0, i1 += inc1, i2 += inc2)
	  if BVE_REF
	    (ra0, i0)
	      if FALSEP
		(scm_eqp (RVREF (ra1, i1, e1), RVREF (ra2, i2, e2)))
		  BVE_CLR (ra0, i0);
	break;
      }
    case tc7_uvect:
    case tc7_ivect:
      for (; n-- > 0; i0 += inc0, i1 += inc1, i2 += inc2)
	if BVE_REF
	  (ra0, i0)
	    if (VELTS (ra1)[i1] != VELTS (ra2)[i2])
	      BVE_CLR (ra0, i0);
      break;
#ifdef FLOATS
#ifdef SINGLES
    case tc7_fvect:
      for (; n-- > 0; i0 += inc0, i1 += inc1, i2 += inc2)
	if BVE_REF
	  (ra0, i0)
	    if (((float *) VELTS (ra1))[i1] != ((float *) VELTS (ra2))[i2])
	      BVE_CLR (ra0, i0);
      break;
#endif /*SINGLES*/
    case tc7_dvect:
      for (; n-- > 0; i0 += inc0, i1 += inc1, i2 += inc2)
	if BVE_REF
	  (ra0, i0)
	    if (((double *) VELTS (ra1))[i1] != ((double *) VELTS (ra2))[i2])
	      BVE_CLR (ra0, i0);
      break;
    case tc7_cvect:
      for (; n-- > 0; i0 += inc0, i1 += inc1, i2 += inc2)
	if BVE_REF
	  (ra0, i0)
	    if (((double *) VELTS (ra1))[2 * i1] != ((double *) VELTS (ra2))[2 * i2] ||
		((double *) VELTS (ra1))[2 * i1 + 1] != ((double *) VELTS (ra2))[2 * i2 + 1])
	      BVE_CLR (ra0, i0);
      break;
#endif /*FLOATS*/
    }
  return 1;
}
/* opt 0 means <, nonzero means >= */
static int ra_compare (ra0, ra1, ra2, opt)
     SCM ra0, ra1, ra2;
     int opt;
{
  long n = ARRAY_DIMS (ra0)->ubnd - ARRAY_DIMS (ra0)->lbnd + 1;
  sizet i0 = ARRAY_BASE (ra0), i1 = ARRAY_BASE (ra1), i2 = ARRAY_BASE (ra2);
  long inc0 = ARRAY_DIMS (ra0)->inc;
  long inc1 = ARRAY_DIMS (ra1)->inc;
  long inc2 = ARRAY_DIMS (ra1)->inc;
  ra0 = ARRAY_V (ra0);
  ra1 = ARRAY_V (ra1);
  ra2 = ARRAY_V (ra2);
  switch (TYP7 (ra1) == TYP7 (ra2) ? TYP7 (ra1) : 0)
    {
    default:
      {
	SCM e1 = SCM_UNDEFINED, e2 = SCM_UNDEFINED;
	for (; n-- > 0; i0 += inc0, i1 += inc1, i2 += inc2)
	  if BVE_REF
	    (ra0, i0)
	      if (opt ?
		  NFALSEP (scm_lessp (RVREF (ra1, i1, e1), RVREF (ra2, i2, e2))) :
		  FALSEP (scm_lessp (RVREF (ra1, i1, e1), RVREF (ra2, i2, e2))))
		BVE_CLR (ra0, i0);
	break;
      }
    case tc7_uvect:
    case tc7_ivect:
      for (; n-- > 0; i0 += inc0, i1 += inc1, i2 += inc2)
	{
	  if BVE_REF
	    (ra0, i0)
	      if (opt ?
		  VELTS (ra1)[i1] < VELTS (ra2)[i2] :
		  VELTS (ra1)[i1] >= VELTS (ra2)[i2])
		BVE_CLR (ra0, i0);
	}
      break;
#ifdef FLOATS
#ifdef SINGLES
    case tc7_fvect:
      for (; n-- > 0; i0 += inc0, i1 += inc1, i2 += inc2)
	if BVE_REF
	  (ra0, i0)
	    if (opt ?
		((float *) VELTS (ra1))[i1] < ((float *) VELTS (ra2))[i2] :
		((float *) VELTS (ra1))[i1] >= ((float *) VELTS (ra2))[i2])
	      BVE_CLR (ra0, i0);
      break;
#endif /*SINGLES*/
    case tc7_dvect:
      for (; n-- > 0; i0 += inc0, i1 += inc1, i2 += inc2)
	if BVE_REF
	  (ra0, i0)
	    if (opt ?
		((double *) VELTS (ra1))[i1] < ((double *) VELTS (ra2))[i2] :
		((double *) VELTS (ra1))[i1] >= ((double *) VELTS (ra2))[i2])
	      BVE_CLR (ra0, i0);
      break;
#endif /*FLOATS*/
    }
  return 1;
}
int scm_ra_lessp (ra0, ras)
     SCM ra0, ras;
{
  return ra_compare (ra0, CAR (ras), CAR (CDR (ras)), 0);
}
int scm_ra_leqp (ra0, ras)
     SCM ra0, ras;
{
  return ra_compare (ra0, CAR (CDR (ras)), CAR (ras), 1);
}
int scm_ra_grp (ra0, ras)
     SCM ra0, ras;
{
  return ra_compare (ra0, CAR (CDR (ras)), CAR (ras), 0);
}
int scm_ra_greqp (ra0, ras)
     SCM ra0, ras;
{
  return ra_compare (ra0, CAR (ras), CAR (CDR (ras)), 1);
}

int scm_ra_sum (ra0, ras)
     SCM ra0, ras;
{
  long n = ARRAY_DIMS (ra0)->ubnd - ARRAY_DIMS (ra0)->lbnd + 1;
  sizet i0 = ARRAY_BASE (ra0);
  long inc0 = ARRAY_DIMS (ra0)->inc;
  ra0 = ARRAY_V (ra0);
  if NNULLP
    (ras)
      {
	SCM ra1 = CAR (ras);
	sizet i1 = ARRAY_BASE (ra1);
	long inc1 = ARRAY_DIMS (ra1)->inc;
	ra1 = ARRAY_V (ra1);
	switch (TYP7 (ra0) == TYP7 (ra1) ? TYP7 (ra0) : 0)
	  {
	  default:
	    {
	      SCM e0 = SCM_UNDEFINED, e1 = SCM_UNDEFINED;
	      for (; n-- > 0; i0 += inc0, i1 += inc1)
		scm_aset (ra0, scm_sum (RVREF (ra0, i0, e0), RVREF (ra1, i1, e1)),
		      MAKINUM (i0));
	      break;
	    }
	  case tc7_uvect:
	  case tc7_ivect:
	    {
	      long *v0 = VELTS (ra0);
	      long *v1 = VELTS (ra1);
	      IVDEP (ra0 != ra1,
		     for (; n-- > 0; i0 += inc0, i1 += inc1)
		     v0[i0] += v1[i1]);
	      break;
	    }
#ifdef FLOATS
#ifdef SINGLES
	  case tc7_fvect:
	    {
	      float *v0 = (float *) VELTS (ra0);
	      float *v1 = (float *) VELTS (ra1);
	      IVDEP (ra0 != ra1,
		     for (; n-- > 0; i0 += inc0, i1 += inc1)
		     v0[i0] += v1[i1]);
	      break;
	    }
#endif /* SINGLES */
	  case tc7_dvect:
	    {
	      double *v0 = (double *) VELTS (ra0);
	      double *v1 = (double *) VELTS (ra1);
	      IVDEP (ra0 != ra1,
		     for (; n-- > 0; i0 += inc0, i1 += inc1)
		     v0[i0] += v1[i1]);
	      break;
	    }
	  case tc7_cvect:
	    {
	      double (*v0)[2] = (double (*)[2]) VELTS (ra0);
	      double (*v1)[2] = (double (*)[2]) VELTS (ra1);
	      IVDEP (ra0 != ra1,
		     for (; n-- > 0; i0 += inc0, i1 += inc1)
		     {
		       v0[i0][0] += v1[i1][0];
		       v0[i0][1] += v1[i1][1];
		     }
		     );
	      break;
	    }
#endif /* FLOATS */
	  }
      }
  return 1;
}

int scm_ra_difference (ra0, ras)
     SCM ra0, ras;
{
  long n = ARRAY_DIMS (ra0)->ubnd - ARRAY_DIMS (ra0)->lbnd + 1;
  sizet i0 = ARRAY_BASE (ra0);
  long inc0 = ARRAY_DIMS (ra0)->inc;
  ra0 = ARRAY_V (ra0);
  if NULLP
    (ras)
      {
	switch TYP7
	  (ra0)
	    {
	    default:
	      {
		SCM e0 = SCM_UNDEFINED;
		for (; n-- > 0; i0 += inc0)
		  scm_aset (ra0, scm_difference (RVREF (ra0, i0, e0), SCM_UNDEFINED), MAKINUM (i0));
		break;
	      }
#ifdef FLOATS
#ifdef SINGLES
	    case tc7_fvect:
	      {
		float *v0 = (float *) VELTS (ra0);
		for (; n-- > 0; i0 += inc0)
		  v0[i0] = -v0[i0];
		break;
	      }
#endif /* SINGLES */
	    case tc7_dvect:
	      {
		double *v0 = (double *) VELTS (ra0);
		for (; n-- > 0; i0 += inc0)
		  v0[i0] = -v0[i0];
		break;
	      }
	    case tc7_cvect:
	      {
		double (*v0)[2] = (double (*)[2]) VELTS (ra0);
		for (; n-- > 0; i0 += inc0)
		  {
		    v0[i0][0] = -v0[i0][0];
		    v0[i0][1] = -v0[i0][1];
		  }
		break;
	      }
#endif /* FLOATS */
	    }
      }
  else
    {
      SCM ra1 = CAR (ras);
      sizet i1 = ARRAY_BASE (ra1);
      long inc1 = ARRAY_DIMS (ra1)->inc;
      ra1 = ARRAY_V (ra1);
      switch (TYP7 (ra0) == TYP7 (ra1) ? TYP7 (ra0) : 0)
	{
	default:
	  {
	    SCM e0 = SCM_UNDEFINED, e1 = SCM_UNDEFINED;
	    for (; n-- > 0; i0 += inc0, i1 += inc1)
	      scm_aset (ra0, scm_difference (RVREF (ra0, i0, e0), RVREF (ra1, i1, e1)), MAKINUM (i0));
	    break;
	  }
#ifdef FLOATS
#ifdef SINGLES
	case tc7_fvect:
	  {
	    float *v0 = (float *) VELTS (ra0);
	    float *v1 = (float *) VELTS (ra1);
	    IVDEP (ra0 != ra1,
		   for (; n-- > 0; i0 += inc0, i1 += inc1)
		   v0[i0] -= v1[i1]);
	    break;
	  }
#endif /* SINGLES */
	case tc7_dvect:
	  {
	    double *v0 = (double *) VELTS (ra0);
	    double *v1 = (double *) VELTS (ra1);
	    IVDEP (ra0 != ra1,
		   for (; n-- > 0; i0 += inc0, i1 += inc1)
		   v0[i0] -= v1[i1]);
	    break;
	  }
	case tc7_cvect:
	  {
	    double (*v0)[2] = (double (*)[2]) VELTS (ra0);
	    double (*v1)[2] = (double (*)[2]) VELTS (ra1);
	    IVDEP (ra0 != ra1,
		   for (; n-- > 0; i0 += inc0, i1 += inc1)
		   {
		     v0[i0][0] -= v1[i1][0];
		     v0[i0][1] -= v1[i1][1];
		   }
		   )
	      break;
	  }
#endif /* FLOATS */
	}
    }
  return 1;
}

int scm_ra_product (ra0, ras)
     SCM ra0, ras;
{
  long n = ARRAY_DIMS (ra0)->ubnd - ARRAY_DIMS (ra0)->lbnd + 1;
  sizet i0 = ARRAY_BASE (ra0);
  long inc0 = ARRAY_DIMS (ra0)->inc;
  ra0 = ARRAY_V (ra0);
  if NNULLP
    (ras)
      {
	SCM ra1 = CAR (ras);
	sizet i1 = ARRAY_BASE (ra1);
	long inc1 = ARRAY_DIMS (ra1)->inc;
	ra1 = ARRAY_V (ra1);
	switch (TYP7 (ra0) == TYP7 (ra1) ? TYP7 (ra0) : 0)
	  {
	  default:
	    {
	      SCM e0 = SCM_UNDEFINED, e1 = SCM_UNDEFINED;
	      for (; n-- > 0; i0 += inc0, i1 += inc1)
		scm_aset (ra0, scm_product (RVREF (ra0, i0, e0), RVREF (ra1, i1, e1)),
		      MAKINUM (i0));
	      break;
	    }
	  case tc7_uvect:
	  case tc7_ivect:
	    {
	      long *v0 = VELTS (ra0);
	      long *v1 = VELTS (ra1);
	      IVDEP (ra0 != ra1,
		     for (; n-- > 0; i0 += inc0, i1 += inc1)
		     v0[i0] *= v1[i1]);
	      break;
	    }
#ifdef FLOATS
#ifdef SINGLES
	  case tc7_fvect:
	    {
	      float *v0 = (float *) VELTS (ra0);
	      float *v1 = (float *) VELTS (ra1);
	      IVDEP (ra0 != ra1,
		     for (; n-- > 0; i0 += inc0, i1 += inc1)
		     v0[i0] *= v1[i1]);
	      break;
	    }
#endif /* SINGLES */
	  case tc7_dvect:
	    {
	      double *v0 = (double *) VELTS (ra0);
	      double *v1 = (double *) VELTS (ra1);
	      IVDEP (ra0 != ra1,
		     for (; n-- > 0; i0 += inc0, i1 += inc1)
		     v0[i0] *= v1[i1]);
	      break;
	    }
	  case tc7_cvect:
	    {
	      double (*v0)[2] = (double (*)[2]) VELTS (ra0);
	      register double r;
	      double (*v1)[2] = (double (*)[2]) VELTS (ra1);
	      IVDEP (ra0 != ra1,
		     for (; n-- > 0; i0 += inc0, i1 += inc1)
		     {
		       r = v0[i0][0] * v1[i1][0] - v0[i0][1] * v1[i1][1];
		       v0[i0][1] = v0[i0][0] * v1[i1][1] + v0[i0][1] * v1[i1][0];
		       v0[i0][0] = r;
		     }
		     );
	      break;
	    }
#endif /* FLOATS */
	  }
      }
  return 1;
}
int scm_ra_divide (ra0, ras)
     SCM ra0, ras;
{
  long n = ARRAY_DIMS (ra0)->ubnd - ARRAY_DIMS (ra0)->lbnd + 1;
  sizet i0 = ARRAY_BASE (ra0);
  long inc0 = ARRAY_DIMS (ra0)->inc;
  ra0 = ARRAY_V (ra0);
  if NULLP
    (ras)
      {
	switch TYP7
	  (ra0)
	    {
	    default:
	      {
		SCM e0 = SCM_UNDEFINED;
		for (; n-- > 0; i0 += inc0)
		  scm_aset (ra0, scm_divide (RVREF (ra0, i0, e0), SCM_UNDEFINED), MAKINUM (i0));
		break;
	      }
#ifdef FLOATS
#ifdef SINGLES
	    case tc7_fvect:
	      {
		float *v0 = (float *) VELTS (ra0);
		for (; n-- > 0; i0 += inc0)
		  v0[i0] = 1.0 / v0[i0];
		break;
	      }
#endif /* SINGLES */
	    case tc7_dvect:
	      {
		double *v0 = (double *) VELTS (ra0);
		for (; n-- > 0; i0 += inc0)
		  v0[i0] = 1.0 / v0[i0];
		break;
	      }
	    case tc7_cvect:
	      {
		register double d;
		double (*v0)[2] = (double (*)[2]) VELTS (ra0);
		for (; n-- > 0; i0 += inc0)
		  {
		    d = v0[i0][0] * v0[i0][0] + v0[i0][1] * v0[i0][1];
		    v0[i0][0] /= d;
		    v0[i0][1] /= -d;
		  }
		break;
	      }
#endif /* FLOATS */
	    }
      }
  else
    {
      SCM ra1 = CAR (ras);
      sizet i1 = ARRAY_BASE (ra1);
      long inc1 = ARRAY_DIMS (ra1)->inc;
      ra1 = ARRAY_V (ra1);
      switch (TYP7 (ra0) == TYP7 (ra1) ? TYP7 (ra0) : 0)
	{
	default:
	  {
	    SCM e0 = SCM_UNDEFINED, e1 = SCM_UNDEFINED;
	    for (; n-- > 0; i0 += inc0, i1 += inc1)
	      scm_aset (ra0, scm_divide (RVREF (ra0, i0, e0), RVREF (ra1, i1, e1)), MAKINUM (i0));
	    break;
	  }
#ifdef FLOATS
#ifdef SINGLES
	case tc7_fvect:
	  {
	    float *v0 = (float *) VELTS (ra0);
	    float *v1 = (float *) VELTS (ra1);
	    IVDEP (ra0 != ra1,
		   for (; n-- > 0; i0 += inc0, i1 += inc1)
		   v0[i0] /= v1[i1]);
	    break;
	  }
#endif /* SINGLES */
	case tc7_dvect:
	  {
	    double *v0 = (double *) VELTS (ra0);
	    double *v1 = (double *) VELTS (ra1);
	    IVDEP (ra0 != ra1,
		   for (; n-- > 0; i0 += inc0, i1 += inc1)
		   v0[i0] /= v1[i1]);
	    break;
	  }
	case tc7_cvect:
	  {
	    register double d, r;
	    double (*v0)[2] = (double (*)[2]) VELTS (ra0);
	    double (*v1)[2] = (double (*)[2]) VELTS (ra1);
	    IVDEP (ra0 != ra1,
		   for (; n-- > 0; i0 += inc0, i1 += inc1)
		   {
		     d = v1[i1][0] * v1[i1][0] + v1[i1][1] * v1[i1][1];
		     r = (v0[i0][0] * v1[i1][0] + v0[i0][1] * v1[i1][1]) / d;
		     v0[i0][1] = (v0[i0][1] * v1[i1][0] - v0[i0][0] * v1[i1][1]) / d;
		     v0[i0][0] = r;
		   }
		   )
	      break;
	  }
#endif /* FLOATS */
	}
    }
  return 1;
}
static int ra_identity (dst, src)
     SCM src, dst;
{
  return racp (CAR (src), scm_cons (dst, EOL));
}

static int ramap (ra0, proc, ras)
     SCM ra0, proc, ras;
{
  long i = ARRAY_DIMS (ra0)->lbnd;
  long inc = ARRAY_DIMS (ra0)->inc;
  long n = ARRAY_DIMS (ra0)->ubnd;
  long base = ARRAY_BASE (ra0) - i * inc;
  ra0 = ARRAY_V (ra0);
  if NULLP
    (ras)
      for (; i <= n; i++)
	scm_aset (ra0, scm_apply (proc, EOL, EOL), MAKINUM (i * inc + base));
  else
    {
      SCM ra1 = CAR (ras);
      SCM args, *ve = &ras;
      sizet k, i1 = ARRAY_BASE (ra1);
      long inc1 = ARRAY_DIMS (ra1)->inc;
      ra1 = ARRAY_V (ra1);
      ras = CDR (ras);
      if NULLP
	(ras)
	  ras = nullvect;
      else
	{
	  ras = scm_vector (ras);
	  ve = VELTS (ras);
	}
      for (; i <= n; i++, i1 += inc1)
	{
	  args = EOL;
	  for (k = LENGTH (ras); k--;)
	    args = scm_cons (scm_aref (ve[k], MAKINUM (i)), args);
	  args = scm_cons (scm_cvref (ra1, i1, SCM_UNDEFINED), args);
	  scm_aset (ra0, scm_apply (proc, args, EOL), MAKINUM (i * inc + base));
	}
    }
  return 1;
}
static int ramap_cxr (ra0, proc, ras)
     SCM ra0, proc, ras;
{
  SCM ra1 = CAR (ras);
  SCM e1 = SCM_UNDEFINED;
  sizet i0 = ARRAY_BASE (ra0), i1 = ARRAY_BASE (ra1);
  long inc0 = ARRAY_DIMS (ra0)->inc, inc1 = ARRAY_DIMS (ra1)->inc;
  long n = ARRAY_DIMS (ra0)->ubnd - ARRAY_DIMS (ra1)->lbnd + 1;
  ra0 = ARRAY_V (ra0);
  ra1 = ARRAY_V (ra1);
  switch TYP7
    (ra0)
      {
      default:
      gencase:
	for (; n-- > 0; i0 += inc0, i1 += inc1)
	  scm_aset (ra0, scm_apply (proc, RVREF (ra1, i1, e1), listofnull), MAKINUM (i0));
	break;
#ifdef FLOATS
#ifdef SINGLES
      case tc7_fvect:
	{
	  float *dst = (float *) VELTS (ra0);
	  switch TYP7
	    (ra1)
	      {
	      default:
		goto gencase;
	      case tc7_fvect:
		for (; n-- > 0; i0 += inc0, i1 += inc1)
		  dst[i0] = DSUBRF (proc) ((double) ((float *) VELTS (ra1))[i1]);
		break;
	      case tc7_uvect:
	      case tc7_ivect:
		for (; n-- > 0; i0 += inc0, i1 += inc1)
		  dst[i0] = DSUBRF (proc) ((double) VELTS (ra1)[i1]);
		break;
	      }
	  break;
	}
#endif /* SINGLES */
      case tc7_dvect:
	{
	  double *dst = (double *) VELTS (ra0);
	  switch TYP7
	    (ra1)
	      {
	      default:
		goto gencase;
	      case tc7_dvect:
		for (; n-- > 0; i0 += inc0, i1 += inc1)
		  dst[i0] = DSUBRF (proc) (((double *) VELTS (ra1))[i1]);
		break;
	      case tc7_uvect:
	      case tc7_ivect:
		for (; n-- > 0; i0 += inc0, i1 += inc1)
		  dst[i0] = DSUBRF (proc) ((double) VELTS (ra1)[i1]);
		break;
	      }
	  break;
	}
#endif /* FLOATS */
      }
  return 1;
}
static int ramap_rp (ra0, proc, ras)
     SCM ra0, proc, ras;
{
  SCM ra1 = CAR (ras), ra2 = CAR (CDR (ras));
  SCM e1 = SCM_UNDEFINED, e2 = SCM_UNDEFINED;
  long n = ARRAY_DIMS (ra0)->ubnd - ARRAY_DIMS (ra0)->lbnd + 1;
  sizet i0 = ARRAY_BASE (ra0), i1 = ARRAY_BASE (ra1), i2 = ARRAY_BASE (ra2);
  long inc0 = ARRAY_DIMS (ra0)->inc;
  long inc1 = ARRAY_DIMS (ra1)->inc;
  long inc2 = ARRAY_DIMS (ra1)->inc;
  ra0 = ARRAY_V (ra0);
  ra1 = ARRAY_V (ra1);
  ra2 = ARRAY_V (ra2);
  switch (TYP7 (ra1) == TYP7 (ra2) ? TYP7 (ra1) : 0)
    {
    default:
      for (; n-- > 0; i0 += inc0, i1 += inc1, i2 += inc2)
	if BVE_REF
	  (ra0, i0)
	    if FALSEP
	      (SUBRF (proc) (RVREF (ra1, i1, e1), RVREF (ra2, i2, e2)))
		BVE_CLR (ra0, i0);
      break;
    case tc7_uvect:
    case tc7_ivect:
      for (; n-- > 0; i0 += inc0, i1 += inc1, i2 += inc2)
	if BVE_REF
	  (ra0, i0)
	    {
	      if FALSEP
		(SUBRF (proc) (MAKINUM (VELTS (ra1)[i1]),
			       MAKINUM (VELTS (ra2)[i2])))
		  BVE_CLR (ra0, i0);
	    }
      break;
#ifdef FLOATS
#ifdef SINGLES
    case tc7_fvect:
      {
	SCM a1 = makflo (1.0), a2 = makflo (1.0);
	for (; n-- > 0; i0 += inc0, i1 += inc1, i2 += inc2)
	  if BVE_REF
	    (ra0, i0)
	      {
		FLO (a1) = ((float *) VELTS (ra1))[i1];
		FLO (a2) = ((float *) VELTS (ra2))[i2];
		if FALSEP
		  (SUBRF (proc) (a1, a2))
		    BVE_CLR (ra0, i0);
	      }
	break;
      }
#endif /*SINGLES*/
    case tc7_dvect:
      {
	SCM a1 = scm_makdbl (1.0 / 3.0, 0.0), a2 = scm_makdbl (1.0 / 3.0, 0.0);
	for (; n-- > 0; i0 += inc0, i1 += inc1, i2 += inc2)
	  if BVE_REF
	    (ra0, i0)
	      {
		REAL (a1) = ((double *) VELTS (ra1))[i1];
		REAL (a2) = ((double *) VELTS (ra2))[i2];
		if FALSEP
		  (SUBRF (proc) (a1, a2))
		    BVE_CLR (ra0, i0);
	      }
	break;
      }
    case tc7_cvect:
      {
	SCM a1 = scm_makdbl (1.0, 1.0), a2 = scm_makdbl (1.0, 1.0);
	for (; n-- > 0; i0 += inc0, i1 += inc1, i2 += inc2)
	  if BVE_REF
	    (ra0, i0)
	      {
		REAL (a1) = ((double *) VELTS (ra1))[2 * i1];
		IMAG (a1) = ((double *) VELTS (ra1))[2 * i1 + 1];
		REAL (a2) = ((double *) VELTS (ra2))[2 * i2];
		IMAG (a2) = ((double *) VELTS (ra2))[2 * i2 + 1];
		if FALSEP
		  (SUBRF (proc) (a1, a2))
		    BVE_CLR (ra0, i0);
	      }
	break;
      }
#endif /*FLOATS*/
    }
  return 1;
}
static int ramap_1 (ra0, proc, ras)
     SCM ra0, proc, ras;
{
  SCM ra1 = CAR (ras);
  SCM e1 = SCM_UNDEFINED;
  long n = ARRAY_DIMS (ra0)->ubnd - ARRAY_DIMS (ra0)->lbnd + 1;
  sizet i0 = ARRAY_BASE (ra0), i1 = ARRAY_BASE (ra1);
  long inc0 = ARRAY_DIMS (ra0)->inc, inc1 = ARRAY_DIMS (ra1)->inc;
  ra0 = ARRAY_V (ra0);
  ra1 = ARRAY_V (ra1);
  if (tc7_vector == TYP7 (ra0))
    for (; n-- > 0; i0 += inc0, i1 += inc1)
      scm_aset (ra0, SUBRF (proc) (scm_cvref (ra1, i1, SCM_UNDEFINED)), MAKINUM (i0));
  else
    for (; n-- > 0; i0 += inc0, i1 += inc1)
      scm_aset (ra0, SUBRF (proc) (RVREF (ra1, i1, e1)), MAKINUM (i0));
  return 1;
}
static int ramap_2o (ra0, proc, ras)
     SCM ra0, proc, ras;
{
  SCM ra1 = CAR (ras);
  SCM e1 = SCM_UNDEFINED;
  long n = ARRAY_DIMS (ra0)->ubnd - ARRAY_DIMS (ra0)->lbnd + 1;
  sizet i0 = ARRAY_BASE (ra0), i1 = ARRAY_BASE (ra1);
  long inc0 = ARRAY_DIMS (ra0)->inc, inc1 = ARRAY_DIMS (ra1)->inc;
  ra0 = ARRAY_V (ra0);
  ra1 = ARRAY_V (ra1);
  ras = CDR (ras);
  if NULLP
    (ras)
      {
	if (tc7_vector == TYP7 (ra0))
	  for (; n-- > 0; i0 += inc0, i1 += inc1)
	    scm_aset (ra0, SUBRF (proc) (scm_cvref (ra1, i1, SCM_UNDEFINED), SCM_UNDEFINED),
		  MAKINUM (i0));
	else
	  for (; n-- > 0; i0 += inc0, i1 += inc1)
	    scm_aset (ra0, SUBRF (proc) (RVREF (ra1, i1, e1), SCM_UNDEFINED),
		  MAKINUM (i0));
      }
  else
    {
      SCM ra2 = CAR (ras);
      SCM e2 = SCM_UNDEFINED;
      sizet i2 = ARRAY_BASE (ra2);
      long inc2 = ARRAY_DIMS (ra2)->inc;
      ra2 = ARRAY_V (ra2);
      if (tc7_vector == TYP7 (ra0))
	for (; n-- > 0; i0 += inc0, i1 += inc1, i2 += inc2)
	  scm_aset (ra0,
		SUBRF (proc) (scm_cvref (ra1, i1, SCM_UNDEFINED), scm_cvref (ra2, i2, SCM_UNDEFINED)),
		MAKINUM (i0));
      else
	for (; n-- > 0; i0 += inc0, i1 += inc1, i2 += inc2)
	  scm_aset (ra0,
		SUBRF (proc) (RVREF (ra1, i1, e1), RVREF (ra2, i2, e2)),
		MAKINUM (i0));
    }
  return 1;
}
static int ramap_a (ra0, proc, ras)
     SCM ra0, proc, ras;
{
  SCM e0 = SCM_UNDEFINED, e1 = SCM_UNDEFINED;
  long n = ARRAY_DIMS (ra0)->ubnd - ARRAY_DIMS (ra0)->lbnd + 1;
  sizet i0 = ARRAY_BASE (ra0);
  long inc0 = ARRAY_DIMS (ra0)->inc;
  ra0 = ARRAY_V (ra0);
  if NULLP
    (ras)
      for (; n-- > 0; i0 += inc0)
	scm_aset (ra0, SUBRF (proc) (RVREF (ra0, i0, e0), SCM_UNDEFINED), MAKINUM (i0));
  else
    {
      SCM ra1 = CAR (ras);
      sizet i1 = ARRAY_BASE (ra1);
      long inc1 = ARRAY_DIMS (ra1)->inc;
      ra1 = ARRAY_V (ra1);
      for (; n-- > 0; i0 += inc0, i1 += inc1)
	scm_aset (ra0, SUBRF (proc) (RVREF (ra0, i0, e0), RVREF (ra1, i1, e1)),
	      MAKINUM (i0));
    }
  return 1;
}

/* These tables are a kluge that will not scale well when more
   vectorized subrs are added.  It is tempting to steal some bits from
   the CAR of all subrs (like those selected by SMOBNUM) to hold an
   offset into a table of vectorized subrs.  */

static ra_iproc ra_rpsubrs[] =
{
  {"=", SCM_UNDEFINED, scm_ra_eqp},
  {"<", SCM_UNDEFINED, scm_ra_lessp},
  {"<=", SCM_UNDEFINED, scm_ra_leqp},
  {">", SCM_UNDEFINED, scm_ra_grp},
  {">=", SCM_UNDEFINED, scm_ra_greqp},
  {0, 0, 0}};
static ra_iproc ra_asubrs[] =
{
  {"+", SCM_UNDEFINED, scm_ra_sum},
  {"-", SCM_UNDEFINED, scm_ra_difference},
  {"*", SCM_UNDEFINED, scm_ra_product},
  {"/", SCM_UNDEFINED, scm_ra_divide},
  {0, 0, 0}};

static char s_sarray_map[] = "serial-scm_array-map!";
#define s_array_map  (s_sarray_map + 7)
SCM scm_array_map (ra0, proc, lra)
     SCM ra0, proc, lra;
{
  ASSERT (BOOL_T == scm_procedurep (proc), proc, ARG2, s_array_map);
  switch TYP7
    (proc)
      {
      default:
      gencase:
	scm_ramapc (ramap, proc, ra0, lra, s_array_map);
	return UNSPECIFIED;
      case tc7_subr_1:
	scm_ramapc (ramap_1, proc, ra0, lra, s_array_map);
	return UNSPECIFIED;
      case tc7_subr_2:
      case tc7_subr_2o:
	scm_ramapc (ramap_2o, proc, ra0, lra, s_array_map);
	return UNSPECIFIED;
      case tc7_cxr:
	if (!SUBRF (proc))
	  goto gencase;
	scm_ramapc (ramap_cxr, proc, ra0, lra, s_array_map);
	return UNSPECIFIED;
      case tc7_rpsubr:
	{
	  ra_iproc *p;
	  if (FALSEP (scm_arrayp (ra0, BOOL_T)))
	    goto gencase;
	  scm_array_fill (ra0, BOOL_T);
	  for (p = ra_rpsubrs; p->name; p++)
	    if (proc == p->sproc)
	      {
		while (NNULLP (lra) && NNULLP (CDR (lra)))
		  {
		    scm_ramapc (p->vproc, SCM_UNDEFINED, ra0, lra, s_array_map);
		    lra = CDR (lra);
		  }
		return UNSPECIFIED;
	      }
	  while (NNULLP (lra) && NNULLP (CDR (lra)))
	    {
	      scm_ramapc (ramap_rp, proc, ra0, lra, s_array_map);
	      lra = CDR (lra);
	    }
	  return UNSPECIFIED;
	}
      case tc7_asubr:
	if NULLP
	  (lra)
	    {
	      SCM prot, fill = SUBRF (proc) (SCM_UNDEFINED, SCM_UNDEFINED);
	      if INUMP
		(fill)
		  {
		    prot = scm_array_prot (ra0);
		    if (NIMP (prot) && INEXP (prot))
		      fill = scm_makdbl ((double) INUM (fill), 0.0);
		  }
	      scm_array_fill (ra0, fill);
	    }
	else
	  {
	    SCM tail, ra1 = CAR (lra);
	    SCM v0 = (NIMP (ra0) && ARRAYP (ra0) ? ARRAY_V (ra0) : ra0);
	    ra_iproc *p;
	    /* Check to see if order might matter.
	       This might be an argument for a separate
	       SERIAL-ARRAY-MAP! */
	    if (v0 == ra1 || (NIMP (ra1) && ARRAYP (ra1) && v0 == ARRAY_V (ra1)))
	      if (ra0 != ra1 || (ARRAYP(ra0) && !ARRAY_CONTP(ra0)))
		goto gencase;
	    for (tail = CDR (lra); NNULLP (tail); tail = CDR (tail))
	      {
		ra1 = CAR (tail);
		if (v0 == ra1 || (NIMP (ra1) && ARRAYP (ra1) && v0 == ARRAY_V (ra1)))
		  goto gencase;
	      }
	    for (p = ra_asubrs; p->name; p++)
	      if (proc == p->sproc)
		{
		  if (ra0 != CAR (lra))
		    scm_ramapc (ra_identity, SCM_UNDEFINED, ra0, scm_cons (CAR (lra), EOL), s_array_map);
		  lra = CDR (lra);
		  while (1)
		    {
		      scm_ramapc (p->vproc, SCM_UNDEFINED, ra0, lra, s_array_map);
		      if (IMP (lra) || IMP (CDR (lra)))
			return UNSPECIFIED;
		      lra = CDR (lra);
		    }
		}
	    scm_ramapc (ramap_2o, proc, ra0, lra, s_array_map);
	    lra = CDR (lra);
	    if NIMP
	      (lra)
		for (lra = CDR (lra); NIMP (lra); lra = CDR (lra))
		  scm_ramapc (ramap_a, proc, ra0, lra, s_array_map);
	  }
	return UNSPECIFIED;
      }
}

static int rafe (ra0, proc, ras)
     SCM ra0, proc, ras;
{
  long i = ARRAY_DIMS (ra0)->lbnd;
  sizet i0 = ARRAY_BASE (ra0);
  long inc0 = ARRAY_DIMS (ra0)->inc;
  long n = ARRAY_DIMS (ra0)->ubnd;
  ra0 = ARRAY_V (ra0);
  if NULLP
    (ras)
      for (; i <= n; i++, i0 += inc0)
	scm_apply (proc, scm_cvref (ra0, i0, SCM_UNDEFINED), listofnull);
  else
    {
      SCM ra1 = CAR (ras);
      SCM args, *ve = &ras;
      sizet k, i1 = ARRAY_BASE (ra1);
      long inc1 = ARRAY_DIMS (ra1)->inc;
      ra1 = ARRAY_V (ra1);
      ras = CDR (ras);
      if NULLP
	(ras)
	  ras = nullvect;
      else
	{
	  ras = scm_vector (ras);
	  ve = VELTS (ras);
	}
      for (; i <= n; i++, i0 += inc0, i1 += inc1)
	{
	  args = EOL;
	  for (k = LENGTH (ras); k--;)
	    args = scm_cons (scm_aref (ve[k], MAKINUM (i)), args);
	  args = scm_cons2 (scm_cvref (ra0, i0, SCM_UNDEFINED), scm_cvref (ra1, i1, SCM_UNDEFINED), args);
	  scm_apply (proc, args, EOL);
	}
    }
  return 1;
}
static char s_array_for_each[] = "array-for-each";
SCM scm_array_for_each (proc, ra0, lra)
     SCM proc, ra0, lra;
{
  ASSERT (BOOL_T == scm_procedurep (proc), proc, ARG1, s_array_for_each);
  scm_ramapc (rafe, proc, ra0, lra, s_array_for_each);
  return UNSPECIFIED;
}

static char s_array_imap[] = "array-index-map!";
SCM scm_array_imap (ra, proc)
     SCM ra, proc;
{
  sizet i;
  ASSERT (NIMP (ra), ra, ARG1, s_array_imap);
  ASSERT (BOOL_T == scm_procedurep (proc), proc, ARG2, s_array_imap);
  switch TYP7
    (ra)
      {
      default:
      badarg:scm_wta (ra, (char *) ARG1, s_array_imap);
      case tc7_vector:
	{
	  SCM *ve = VELTS (ra);
	  for (i = 0; i < LENGTH (ra); i++)
	    ve[i] = scm_apply (proc, MAKINUM (i), listofnull);
	  return UNSPECIFIED;
	}
      case tc7_string:
      case tc7_bvect:
      case tc7_uvect:
      case tc7_ivect:
      case tc7_fvect:
      case tc7_dvect:
      case tc7_cvect:
	for (i = 0; i < LENGTH (ra); i++)
	  scm_aset (ra, scm_apply (proc, MAKINUM (i), listofnull), MAKINUM (i));
	return UNSPECIFIED;
      case tc7_smob:
	ASRTGO (ARRAYP (ra), badarg);
	{
	  SCM args = EOL;
	  SCM inds = scm_make_uve (ARRAY_NDIM (ra), MAKINUM (-1L));
	  long *vinds = VELTS (inds);
	  int j, k, kmax = ARRAY_NDIM (ra) - 1;
	  for (k = 0; k <= kmax; k++)
	    vinds[k] = ARRAY_DIMS (ra)[k].lbnd;
	  k = kmax;
	  do
	    {
	      if (k == kmax)
		{
		  vinds[k] = ARRAY_DIMS (ra)[k].lbnd;
		  i = cind (ra, inds);
		  for (; vinds[k] <= ARRAY_DIMS (ra)[k].ubnd; vinds[k]++)
		    {
		      for (j = kmax + 1, args = EOL; j--;)
			args = scm_cons (MAKINUM (vinds[j]), args);
		      scm_aset (ARRAY_V (ra), scm_apply (proc, args, EOL), MAKINUM (i));
		      i += ARRAY_DIMS (ra)[k].inc;
		    }
		  k--;
		  continue;
		}
	      if (vinds[k] < ARRAY_DIMS (ra)[k].ubnd)
		{
		  vinds[k]++;
		  k++;
		  continue;
		}
	      vinds[k] = ARRAY_DIMS (ra)[k].lbnd - 1;
	      k--;
	    }
	  while (k >= 0);
	  return UNSPECIFIED;
	}
      }
}

SCM scm_array_equal P ((SCM ra0, SCM ra1));
static int raeql_1 (ra0, as_equal, ra1)
     SCM ra0, as_equal, ra1;
{
  SCM e0 = SCM_UNDEFINED, e1 = SCM_UNDEFINED;
  sizet i0 = 0, i1 = 0;
  long inc0 = 1, inc1 = 1;
  sizet n = LENGTH (ra0);
  ra1 = CAR (ra1);
  if ARRAYP
    (ra0)
      {
	n = ARRAY_DIMS (ra0)->ubnd - ARRAY_DIMS (ra0)->lbnd + 1;
	i0 = ARRAY_BASE (ra0);
	inc0 = ARRAY_DIMS (ra0)->inc;
	ra0 = ARRAY_V (ra0);
      }
  if ARRAYP
    (ra1)
      {
	i1 = ARRAY_BASE (ra1);
	inc1 = ARRAY_DIMS (ra1)->inc;
	ra1 = ARRAY_V (ra1);
      }
  switch TYP7
    (ra0)
      {
      case tc7_vector:
      default:
	for (; n--; i0 += inc0, i1 += inc1)
	  {
	    if FALSEP
	      (as_equal)
		{
		  if FALSEP
		    (scm_array_equal (RVREF (ra0, i0, e0), RVREF (ra1, i1, e1)))
		      return 0;
		}
	    else if FALSEP
	      (scm_equal (RVREF (ra0, i0, e0), RVREF (ra1, i1, e1)))
		return 0;
	  }
	return 1;
      case tc7_string:
	{
	  char *v0 = CHARS (ra0) + i0;
	  char *v1 = CHARS (ra1) + i1;
	  for (; n--; v0 += inc0, v1 += inc1)
	    if (*v0 != *v1)
	      return 0;
	  return 1;
	}
      case tc7_bvect:
	for (; n--; i0 += inc0, i1 += inc1)
	  if (BVE_REF (ra0, i0) != BVE_REF (ra1, i1))
	    return 0;
	return 1;
      case tc7_uvect:
      case tc7_ivect:
	{
	  long *v0 = (long *) VELTS (ra0) + i0;
	  long *v1 = (long *) VELTS (ra1) + i1;
	  for (; n--; v0 += inc0, v1 += inc1)
	    if (*v0 != *v1)
	      return 0;
	  return 1;
	}
#ifdef FLOATS
#ifdef SINGLES
      case tc7_fvect:
	{
	  float *v0 = (float *) VELTS (ra0) + i0;
	  float *v1 = (float *) VELTS (ra1) + i1;
	  for (; n--; v0 += inc0, v1 += inc1)
	    if (*v0 != *v1)
	      return 0;
	  return 1;
	}
#endif /* SINGLES */
      case tc7_dvect:
	{
	  double *v0 = (double *) VELTS (ra0) + i0;
	  double *v1 = (double *) VELTS (ra1) + i1;
	  for (; n--; v0 += inc0, v1 += inc1)
	    if (*v0 != *v1)
	      return 0;
	  return 1;
	}
      case tc7_cvect:
	{
	  double (*v0)[2] = (double (*)[2]) VELTS (ra0) + i0;
	  double (*v1)[2] = (double (*)[2]) VELTS (ra1) + i1;
	  for (; n--; v0 += inc0, v1 += inc1)
	    {
	      if ((*v0)[0] != (*v1)[0])
		return 0;
	      if ((*v0)[1] != (*v1)[1])
		return 0;
	    }
	  return 1;
	}
#endif /* FLOATS */
      }
}
static int raeql (ra0, as_equal, ra1)
     SCM ra0, as_equal, ra1;
{
  SCM v0 = ra0, v1 = ra1;
  scm_array_dim dim0, dim1;
  scm_array_dim *s0 = &dim0, *s1 = &dim1;
  sizet bas0 = 0, bas1 = 0;
  int k, unroll = 1, vlen = 1, ndim = 1;
  if ARRAYP
    (ra0)
      {
	ndim = ARRAY_NDIM (ra0);
	s0 = ARRAY_DIMS (ra0);
	bas0 = ARRAY_BASE (ra0);
	v0 = ARRAY_V (ra0);
      }
  else
    {
      s0->inc = 1;
      s0->lbnd = 0;
      s0->ubnd = LENGTH (v0) - 1;
      unroll = 0;
    }
  if ARRAYP
    (ra1)
      {
	if (ndim != ARRAY_NDIM (ra1))
	  return 0;
	s1 = ARRAY_DIMS (ra1);
	bas1 = ARRAY_BASE (ra1);
	v1 = ARRAY_V (ra1);
      }
  else
    {
      if (1 != ndim)
	return BOOL_F;
      s1->inc = 1;
      s1->lbnd = 0;
      s1->ubnd = LENGTH (v1) - 1;
      unroll = 0;
    }
  if (TYP7 (v0) != TYP7 (v1))
    return 0;
  for (k = ndim; k--;)
    {
      if (s0[k].lbnd != s1[k].lbnd || s0[k].ubnd != s1[k].ubnd)
	return 0;
      if (unroll)
	{
	  unroll = (s0[k].inc == s1[k].inc);
	  vlen *= s0[k].ubnd - s1[k].lbnd + 1;
	}
    }
  if (unroll && bas0 == bas1 && v0 == v1)
    return BOOL_T;
  return scm_ramapc (raeql_1, as_equal, ra0, scm_cons (ra1, EOL), "");
}

SCM scm_raequal (ra0, ra1)
     SCM ra0, ra1;
{
  return (raeql (ra0, BOOL_T, ra1) ? BOOL_T : BOOL_F);
}
static char s_array_equalp[] = "array-equal?";
SCM scm_array_equal (ra0, ra1)
     SCM ra0, ra1;
{
  if (IMP (ra0) || IMP (ra1))
  callequal:return scm_equal (ra0, ra1);
  switch TYP7
    (ra0)
      {
      default:
	goto callequal;
      case tc7_bvect:
      case tc7_string:
      case tc7_uvect:
      case tc7_ivect:
      case tc7_fvect:
      case tc7_dvect:
      case tc7_cvect:
      case tc7_vector:
	break;
      case tc7_smob:
	if (!ARRAYP (ra0))
	  goto callequal;
      }
  switch TYP7
    (ra1)
      {
      default:
	goto callequal;
      case tc7_bvect:
      case tc7_string:
      case tc7_uvect:
      case tc7_ivect:
      case tc7_fvect:
      case tc7_dvect:
      case tc7_cvect:
      case tc7_vector:
	break;
      case tc7_smob:
	if (!ARRAYP (ra1))
	  goto callequal;
      }
  return (raeql (ra0, BOOL_F, ra1) ? BOOL_T : BOOL_F);
}

static scm_iproc subr2s[] =
{
  {s_array_fill, scm_array_fill},
  {s_array_copy, scm_array_copy},
  {s_sarray_copy, scm_array_copy},
  {0, 0}};

static scm_iproc lsubr2s[] =
{
  {s_array_map, scm_array_map},
  {s_sarray_map, scm_array_map},
  {s_array_for_each, scm_array_for_each},
  {s_array_imap, scm_array_imap},
  {0, 0}};

static void init_raprocs (subra)
     ra_iproc *subra;
{
  for (; subra->name; subra++)
    subra->sproc = CDR (scm_intern (subra->name, strlen (subra->name)));
}

void init_ramap ()
{
  init_raprocs (ra_rpsubrs);
  init_raprocs (ra_asubrs);
  scm_init_iprocs (subr2s, tc7_subr_2);
  scm_init_iprocs (lsubr2s, tc7_lsubr_2);
  scm_make_subr (s_array_equalp, tc7_rpsubr, scm_array_equal);
  scm_smobs[0x0ff & (scm_tc16_array >> 8)].equalp = scm_raequal;
  scm_add_feature (s_array_for_each);
}

#endif /* ARRAYS */
