home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / mitsch75.zip / scheme-7_5_17-src.zip / scheme-7.5.17 / src / microcode / array.c < prev    next >
C/C++ Source or Header  |  1999-01-02  |  63KB  |  2,250 lines
  1. /* -*-C-*-
  2.  
  3. $Id: array.c,v 9.46 1999/01/02 06:11:34 cph Exp $
  4.  
  5. Copyright (c) 1987-1999 Massachusetts Institute of Technology
  6.  
  7. This program is free software; you can redistribute it and/or modify
  8. it under the terms of the GNU General Public License as published by
  9. the Free Software Foundation; either version 2 of the License, or (at
  10. your option) any later version.
  11.  
  12. This program is distributed in the hope that it will be useful, but
  13. WITHOUT ANY WARRANTY; without even the implied warranty of
  14. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  15. General Public License for more details.
  16.  
  17. You should have received a copy of the GNU General Public License
  18. along with this program; if not, write to the Free Software
  19. Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  20. */
  21.  
  22. #include "scheme.h"
  23. #include "prims.h"
  24. #include "array.h"
  25. #include <math.h>
  26. #include <values.h>
  27. /* <values.h> contains some math constants */
  28.  
  29. /* ARRAY (as a scheme object)
  30.    is a usual array (in C) containing REAL numbers (float/double)
  31.    and tagged as a non-marked vector.
  32.  
  33.    Basic contents:
  34.    constructors, selectors, arithmetic operations,
  35.    conversion routines between C_Array, and Scheme_Vector
  36.  
  37.    see array.h for macros, NM_VECTOR, and extern */
  38.  
  39. /* mathematical constants */
  40. #ifdef PI
  41. #undef PI
  42. #endif
  43. #define PI           3.141592653589793238462643
  44. #define PI_OVER_2    1.570796326794896619231322
  45. #define TWOPI        6.283185307179586476925287
  46. #define SQRT_2          1.4142135623730950488
  47. #define ONE_OVER_SQRT_2  .7071067811865475244
  48. /* Abramowitz and Stegun p.3 */
  49.  
  50. REAL
  51. flonum_to_real (argument, arg_number)
  52.      fast SCHEME_OBJECT argument;
  53.      int arg_number;
  54. {
  55.   switch (OBJECT_TYPE (argument))
  56.     {
  57.     case TC_FIXNUM:
  58.       return ((REAL) (FIXNUM_TO_DOUBLE (argument)));
  59.  
  60.     case TC_BIG_FIXNUM:
  61.       if (! (BIGNUM_TO_DOUBLE_P (argument)))
  62.       error_bad_range_arg (arg_number);
  63.       return ((REAL) (bignum_to_double (argument)));
  64.  
  65.     case TC_BIG_FLONUM:
  66.       return ((REAL) (FLONUM_TO_DOUBLE (argument)));
  67.  
  68.     default:
  69.       error_wrong_type_arg (arg_number);
  70.       /* NOTREACHED */
  71.     }
  72. }
  73.  
  74. SCHEME_OBJECT
  75. allocate_array (length)
  76.      long length;
  77. {
  78. #if (REAL_IS_DEFINED_DOUBLE == 0)
  79.  
  80.   fast SCHEME_OBJECT result =
  81.     (allocate_non_marked_vector
  82.      (TC_NON_MARKED_VECTOR, ((length * REAL_SIZE) + 1), true));
  83.   FAST_MEMORY_SET (result, 1, length);
  84.   return (result);
  85.  
  86. #else /* (REAL_IS_DEFINED_DOUBLE != 0) */
  87.   
  88.   long n_words = (length * DOUBLE_SIZE);
  89.   ALIGN_FLOAT (Free);
  90.   Primitive_GC_If_Needed (n_words + 1);
  91.   {
  92.     SCHEME_OBJECT result = (MAKE_POINTER_OBJECT (TC_BIG_FLONUM, (Free)));
  93.     (*Free++) = (MAKE_OBJECT (TC_MANIFEST_NM_VECTOR, n_words));
  94.     Free += n_words;
  95.     return (result);
  96.   }
  97.  
  98. #endif /* (REAL_IS_DEFINED_DOUBLE != 0) */
  99. }
  100.  
  101. DEFINE_PRIMITIVE ("VECTOR->ARRAY", Prim_vector_to_array, 1, 1, 0)
  102. {
  103.   PRIMITIVE_HEADER (1);
  104.   CHECK_ARG (1, VECTOR_P);
  105.   {
  106.     SCHEME_OBJECT vector = (ARG_REF (1));
  107.     long length = (VECTOR_LENGTH (vector));
  108.     SCHEME_OBJECT result = (allocate_array (length));
  109.     fast SCHEME_OBJECT * scan_source = (& (VECTOR_REF (vector, 0)));
  110.     fast SCHEME_OBJECT * end_source = (scan_source + length);
  111.     fast REAL * scan_target = (ARRAY_CONTENTS (result));
  112.     while (scan_source < end_source)
  113.       (*scan_target++) = (flonum_to_real ((*scan_source++), 1));
  114.     PRIMITIVE_RETURN (result);
  115.   }
  116. }
  117.  
  118. DEFINE_PRIMITIVE ("ARRAY->VECTOR", Prim_array_to_vector, 1, 1, 0)
  119. {
  120.   PRIMITIVE_HEADER (1);
  121.   CHECK_ARG (1, ARRAY_P);
  122.   {
  123.     SCHEME_OBJECT array = (ARG_REF (1));
  124.     long length = (ARRAY_LENGTH (array));
  125.     fast REAL * scan_source = (ARRAY_CONTENTS (array));
  126.     fast REAL * end_source = (scan_source + length);
  127.     SCHEME_OBJECT result = (allocate_marked_vector (TC_VECTOR, length, true));
  128.     fast SCHEME_OBJECT * scan_result = (MEMORY_LOC (result, 1));
  129.     while (scan_source < end_source)
  130.       (*scan_result++) = (double_to_flonum ((double) (*scan_source++)));
  131.     PRIMITIVE_RETURN (result);
  132.   }
  133. }
  134.  
  135. DEFINE_PRIMITIVE ("ARRAY-ALLOCATE", Prim_array_allocate, 1,1, 0)
  136. {
  137.   fast REAL * scan;
  138.   long length;
  139.   SCHEME_OBJECT result;
  140.   PRIMITIVE_HEADER (1);
  141.  
  142.   length = (arg_nonnegative_integer (1));
  143.   result = (allocate_array (length));
  144.   for (scan = (ARRAY_CONTENTS (result)); --length >= 0; )
  145.     *scan++ = ((REAL) 0.0);
  146.   PRIMITIVE_RETURN (result);
  147. }
  148.  
  149. DEFINE_PRIMITIVE ("ARRAY-CONS-REALS", Prim_array_cons_reals, 3, 3, 0)
  150. {
  151.   PRIMITIVE_HEADER (3);
  152.   {
  153.     fast double from = (arg_real_number (1));
  154.     fast double dt = (arg_real_number (2));
  155.     long length = (arg_nonnegative_integer (3));
  156.     SCHEME_OBJECT result = (allocate_array (length));
  157.     fast REAL * scan_result = (ARRAY_CONTENTS (result));
  158.     fast int i;
  159.     for (i = 0; (i < length); i += 1)
  160.       {
  161.     (*scan_result++) = ((REAL) from);
  162.     from += dt;
  163.       }
  164.     PRIMITIVE_RETURN (result);
  165.   }
  166. }
  167.  
  168. DEFINE_PRIMITIVE ("ARRAY-LENGTH", Prim_array_length, 1, 1, 0)
  169. {
  170.   PRIMITIVE_HEADER (1);
  171.   CHECK_ARG (1, ARRAY_P);
  172.   PRIMITIVE_RETURN (LONG_TO_UNSIGNED_FIXNUM (ARRAY_LENGTH (ARG_REF (1))));
  173. }
  174.  
  175. DEFINE_PRIMITIVE ("ARRAY-REF", Prim_array_ref, 2, 2, 0)
  176. {
  177.   SCHEME_OBJECT array;
  178.   PRIMITIVE_HEADER (2);
  179.   CHECK_ARG (1, ARRAY_P);
  180.   array = (ARG_REF (1));
  181.   PRIMITIVE_RETURN
  182.     (double_to_flonum
  183.      ((double)
  184.       ((ARRAY_CONTENTS (array))
  185.        [arg_index_integer (2, (ARRAY_LENGTH (array)))])));
  186. }
  187.  
  188. DEFINE_PRIMITIVE ("ARRAY-SET!", Prim_array_set, 3, 3, 0)
  189. {
  190.   SCHEME_OBJECT array;
  191.   REAL * array_ptr;
  192.   double old_value, new_value;
  193.   PRIMITIVE_HEADER (3);
  194.   CHECK_ARG (1, ARRAY_P);
  195.   array = (ARG_REF (1));
  196.   array_ptr =
  197.     (& ((ARRAY_CONTENTS (array))
  198.     [arg_index_integer (2, (ARRAY_LENGTH (array)))]));
  199.   old_value = (*array_ptr);
  200.   new_value = (arg_real_number (3));
  201. #if (REAL_IS_DEFINED_DOUBLE == 0)
  202.   if ((new_value >= 0.0)
  203.       ? (new_value < ((double) FLT_MIN))
  204.       : (new_value > (0.0 - ((double) FLT_MIN))))
  205.     new_value = ((REAL) 0.0);
  206. #endif
  207.   (*array_ptr) = ((REAL) new_value);
  208.   PRIMITIVE_RETURN (double_to_flonum (old_value));
  209. }
  210.  
  211. /*____________________ file readers ___________
  212.   ascii and 2bint formats 
  213.   ______________________________________________*/
  214.  
  215. /* Reading data from files 
  216.    To read REAL numbers, use "lf" for double, "%f" for float 
  217.    */
  218. #if (REAL_IS_DEFINED_DOUBLE == 1)
  219. #define REALREAD  "%lf"
  220. #define REALREAD2 "%lf %lf"
  221. #else
  222. #define REALREAD  "%f"
  223. #define REALREAD2 "%f %f"
  224. #endif
  225.  
  226. static void
  227. C_Array_Read_Ascii_File (a, N, fp)          /* 16 ascii decimal digits */
  228.      REAL * a;
  229.      long N;
  230.      FILE * fp;
  231. { 
  232.   fast long i;
  233.   for (i = 0; (i < N); i += 1)
  234.     {
  235.       if ((fscanf (fp, REALREAD, (&(a[i])))) != 1)
  236.     { printf("Not enough values read ---\n last value a[%d] = % .16e \n", (i-1), a[i-1]);
  237.       error_external_return (); }
  238.     }
  239.   return;
  240. }
  241.  
  242. /* 2BINT FORMAT = integer stored in 2 consecutive bytes.
  243.    On many machines, "putw" and "getw" use 4 byte integers (C int)
  244.    so use "putc" "getc" as shown below.
  245.    */
  246.  
  247. static void
  248. C_Array_Read_2bint_File (a, N, fp)
  249.      REAL * a;
  250.      long N;
  251.      FILE * fp;
  252. {
  253.   fast long i;
  254.   fast int msd;
  255.   for (i = 0; (i < N); i += 1)
  256.     {
  257.       if (feof (fp))
  258.     error_external_return ();
  259.       msd = (getc (fp));
  260.       (a [i]) = ((REAL) ((msd << 8) | (getc (fp))));
  261.     }
  262.   return;
  263. }
  264.  
  265. DEFINE_PRIMITIVE ("ARRAY-READ-FROM-FILE", Prim_array_read_from_file, 3,3, 0)
  266. {
  267.   PRIMITIVE_HEADER (3);
  268.   CHECK_ARG (1, STRING_P);    /* 1 = filename */
  269.   /*                               2 = length of data */
  270.   CHECK_ARG (3, FIXNUM_P);    /* 3 = format of data   0=ascii 1=2bint  */
  271.   {
  272.     fast long length = (arg_nonnegative_integer (2));
  273.     fast SCHEME_OBJECT result = (allocate_array (length));
  274.     int format;
  275.     fast FILE * fp;
  276.     if ( (fp = fopen((STRING_ARG (1)), "r")) == NULL)
  277.       error_bad_range_arg (1);
  278.     
  279.     format = arg_nonnegative_integer(3);
  280.     if (format==0)
  281.       C_Array_Read_Ascii_File ((ARRAY_CONTENTS (result)), length, fp);
  282.     else if (format==1)
  283.       C_Array_Read_2bint_File ((ARRAY_CONTENTS (result)), length, fp);
  284.     else
  285.       error_bad_range_arg(3);    /* illegal format code */
  286.     
  287.     if ((fclose (fp)) != 0)
  288.       error_external_return ();
  289.     PRIMITIVE_RETURN (result);
  290.   }
  291. }
  292.  
  293. static void
  294. C_Array_Write_Ascii_File (a, N, fp)           /* 16 ascii decimal digits */
  295.      REAL * a;
  296.      long N;
  297.      FILE * fp;
  298. {
  299.   fast long i;
  300.   for (i = 0; (i < N); i += 1)
  301.     {
  302.       if (feof (fp))
  303.     error_external_return ();
  304.       fprintf (fp, "% .16e \n", a[i]);
  305.     }
  306.   return;
  307. }
  308.  
  309. DEFINE_PRIMITIVE ("ARRAY-WRITE-ASCII-FILE", Prim_array_write_ascii_file, 2, 2, 0)
  310. {
  311.   PRIMITIVE_HEADER (2);
  312.   CHECK_ARG (1, ARRAY_P);
  313.   CHECK_ARG (2, STRING_P);
  314.   {
  315.     fast SCHEME_OBJECT array = (ARG_REF (1));
  316.     fast FILE * fp = (fopen((STRING_ARG (2)), "w"));
  317.     if (fp == ((FILE *) 0))
  318.       error_bad_range_arg (2);
  319.     C_Array_Write_Ascii_File
  320.       ((ARRAY_CONTENTS (array)),
  321.        (ARRAY_LENGTH (array)),
  322.        fp);
  323.     if ((fclose (fp)) != 0)
  324.       error_external_return ();
  325.   }
  326.   PRIMITIVE_RETURN (UNSPECIFIC);
  327. }
  328.  
  329.  
  330.  
  331.  
  332. DEFINE_PRIMITIVE ("SUBARRAY-COPY!", Prim_subarray_copy, 5, 5, 0)
  333. {
  334.   PRIMITIVE_HEADER (5);
  335.   CHECK_ARG (1, ARRAY_P);    /* source array */
  336.   CHECK_ARG (2, ARRAY_P);    /* destination array */
  337.   {
  338.     REAL * source = (ARRAY_CONTENTS (ARG_REF (1)));
  339.     REAL * target = (ARRAY_CONTENTS (ARG_REF (2)));
  340.     long start_source = (arg_nonnegative_integer (3));
  341.     long start_target = (arg_nonnegative_integer (4));
  342.     long n_elements = (arg_nonnegative_integer (5));
  343.     if ((start_source + n_elements) > (ARRAY_LENGTH (ARG_REF (1))))
  344.       error_bad_range_arg (3);
  345.     if ((start_target + n_elements) > (ARRAY_LENGTH (ARG_REF (2))))
  346.       error_bad_range_arg (4);
  347.     {
  348.       fast REAL * scan_source = (source + start_source);
  349.       fast REAL * end_source = (scan_source + n_elements);
  350.       fast REAL * scan_target = (target + start_target);
  351.       while (scan_source < end_source)
  352.     (*scan_target++) = (*scan_source++);
  353.     }
  354.   }
  355.   PRIMITIVE_RETURN (UNSPECIFIC);
  356. }
  357.  
  358. DEFINE_PRIMITIVE ("ARRAY-REVERSE!", Prim_array_reverse, 1, 1, 0)
  359. {
  360.   PRIMITIVE_HEADER (1);
  361.   CHECK_ARG (1, ARRAY_P);
  362.   {
  363.     SCHEME_OBJECT array = (ARG_REF (1));
  364.     long length = (ARRAY_LENGTH (array));
  365.     long half_length = (length / 2);
  366.     fast REAL * array_ptr = (ARRAY_CONTENTS (array));
  367.     fast long i;
  368.     fast long j;
  369.     for (i = 0, j = (length - 1); (i < half_length); i += 1, j -= 1)
  370.       {
  371.     fast REAL Temp = (array_ptr [j]);
  372.     (array_ptr [j]) = (array_ptr [i]);
  373.     (array_ptr [i]) = Temp;
  374.       }
  375.   }
  376.   PRIMITIVE_RETURN (UNSPECIFIC);
  377. }
  378.  
  379. DEFINE_PRIMITIVE ("ARRAY-TIME-REVERSE!", Prim_array_time_reverse, 1, 1, 0)
  380. {
  381.   void C_Array_Time_Reverse ();
  382.   PRIMITIVE_HEADER (1);
  383.   CHECK_ARG (1, ARRAY_P);
  384.   C_Array_Time_Reverse
  385.     ((ARRAY_CONTENTS (ARG_REF (1))), (ARRAY_LENGTH (ARG_REF (1))));
  386.   PRIMITIVE_RETURN (UNSPECIFIC);
  387. }
  388.  
  389. /* time-reverse
  390.    x[0] remains fixed. (time=0)
  391.    x[i] swapped with x[n-i]. (mirror image around x[0]) */
  392.  
  393. void
  394. C_Array_Time_Reverse (x,n)
  395.      REAL *x;
  396.      long n;
  397. { long i, ni, n2;
  398.   REAL xt;
  399.   if ((n % 2) == 0)        /* even length */
  400.   { n2 = (n/2);
  401.     for (i=1; i<n2; i++)    /* i=1,2,..,n/2-1 */
  402.     {  ni = n-i;
  403.        xt    = x[i];
  404.        x[i]  = x[ni];
  405.        x[ni] = xt; }}
  406.   else                /* odd length */
  407.   { n2 = (n+1)/2;        /* (n+1)/2 = (n-1)/2 + 1 */
  408.     for (i=1; i<n2; i++)    /* i=1,2,..,(n-1)/2 */
  409.     {  ni = n-i;
  410.        xt   = x[i];
  411.        x[i] = x[ni];
  412.        x[ni] = xt; }}
  413. }
  414.  
  415. /* The following is smart
  416.    and avoids computation when offset or scale are degenerate 0,1 */
  417.  
  418. DEFINE_PRIMITIVE ("SUBARRAY-OFFSET-SCALE!", Prim_subarray_offset_scale, 5, 5, 0)
  419. {
  420.   long i, at, m,mplus;
  421.   REAL *a, offset,scale;
  422.   PRIMITIVE_HEADER (5);
  423.   CHECK_ARG (1, ARRAY_P);
  424.   CHECK_ARG (2, FIXNUM_P);
  425.   CHECK_ARG (3, FIXNUM_P);
  426.   a = ARRAY_CONTENTS(ARG_REF(1));
  427.   at = arg_nonnegative_integer(2); /*       at = starting index             */
  428.   m  = arg_nonnegative_integer(3); /*       m  = number of points to change */
  429.   mplus = at + m;
  430.   if (mplus > (ARRAY_LENGTH(ARG_REF(1)))) error_bad_range_arg(3);
  431.   offset = (arg_real (4));
  432.   scale = (arg_real (5));
  433.   if ((offset == 0.0) && (scale == 1.0))
  434.     ;                /* be smart */
  435.   else if (scale == 0.0)
  436.     for (i=at; i<mplus; i++)  a[i] = offset;
  437.   else if (offset == 0.0)
  438.     for (i=at; i<mplus; i++)  a[i] = scale * a[i];
  439.   else if (scale == 1.0)
  440.     for (i=at; i<mplus; i++)  a[i] = offset + a[i];
  441.   else
  442.     for (i=at; i<mplus; i++)  a[i] = offset + scale * a[i];
  443.   PRIMITIVE_RETURN (UNSPECIFIC);
  444. }
  445.  
  446. DEFINE_PRIMITIVE ("COMPLEX-SUBARRAY-COMPLEX-SCALE!", Prim_complex_subarray_complex_scale, 6,6, 0)
  447. { long i, at,m,mplus;
  448.   REAL *a,*b;            /* (a,b) = (real,imag) arrays */
  449.   double temp, minus_y,  x, y;    /* (x,y) = (real,imag) scale */
  450.   PRIMITIVE_HEADER (6);
  451.   CHECK_ARG (1, ARRAY_P);
  452.   CHECK_ARG (2, ARRAY_P);
  453.   at = (arg_nonnegative_integer (3)); /* starting index */
  454.   m  = (arg_nonnegative_integer (4)); /* number of points to change */
  455.   mplus = at + m;
  456.   if (mplus > (ARRAY_LENGTH(ARG_REF(1)))) error_bad_range_arg(4);
  457.   x = (arg_real_number (5));
  458.   y = (arg_real_number (6));
  459.   a = ARRAY_CONTENTS(ARG_REF(1));
  460.   b = ARRAY_CONTENTS(ARG_REF(2));
  461.   if ((ARRAY_LENGTH(ARG_REF(1))) != (ARRAY_LENGTH(ARG_REF(2))))
  462.     error_bad_range_arg(2);
  463.   if (x==0.0)            /* imaginary only */
  464.     if       (y==0.0)
  465.       for (i=at; i<mplus; i++)
  466.       { a[i] = 0.0;
  467.     b[i] = 0.0; }
  468.     else if  (y==1.0)
  469.       for (i=at; i<mplus; i++)
  470.       { temp = b[i];
  471.     b[i] = a[i];
  472.     a[i] = (-temp); }
  473.     else if  (y==-1.0)
  474.       for (i=at; i<mplus; i++)
  475.       { temp = b[i];
  476.     b[i] = (-a[i]);
  477.     a[i] = temp; }
  478.     else
  479.     { minus_y = (-y);
  480.       for (i=at; i<mplus; i++)
  481.       { temp =               y * ((double) a[i]);
  482.     a[i] = (REAL) (minus_y * ((double) b[i]));
  483.     b[i] = (REAL) temp; }}
  484.   else if (y==0.0)        /* real only */
  485.     if (x==1.0) ;
  486.     else for (i=at; i<mplus; i++)
  487.     { a[i] = (REAL) (x * ((double) a[i]));
  488.       b[i] = (REAL) (x * ((double) b[i])); }
  489.   else                /* full complex scale */
  490.     for (i=at; i<mplus; i++)
  491.     { temp =         ((double) a[i])*x - ((double) b[i])*y;
  492.       b[i] = (REAL) (((double) b[i])*x + ((double) a[i])*y);
  493.       a[i] = (REAL) temp; }
  494.   PRIMITIVE_RETURN (UNSPECIFIC);
  495. }
  496.  
  497. /* Accumulate
  498.    using combinators              *
  499.    corresponding type codes       1 */
  500.  
  501. DEFINE_PRIMITIVE ("COMPLEX-SUBARRAY-ACCUMULATE!", Prim_complex_subarray_accumulate, 6,6, 0)
  502. { long  at,m,mplus, tc, i;
  503.   REAL *a,*b;            /* (a,b) = (real,imag) input arrays */
  504.   REAL *c;    /* result = output array of length 2, holds a complex number */
  505.   double x, y, temp;
  506.   PRIMITIVE_HEADER (6);
  507.   CHECK_ARG (1, ARRAY_P);    /* a = input array (real) */
  508.   CHECK_ARG (2, ARRAY_P);    /* b = input array (imag) */
  509.   a = ARRAY_CONTENTS(ARG_REF(1));
  510.   b = ARRAY_CONTENTS(ARG_REF(2));
  511.   if ((ARRAY_LENGTH(ARG_REF(1))) != (ARRAY_LENGTH(ARG_REF(2))))
  512.     error_bad_range_arg(2);
  513.   tc = arg_nonnegative_integer(3); /*       tc = type code 0 or 1            */
  514.   at = arg_nonnegative_integer(4); /*       at = starting index              */
  515.   m  = arg_nonnegative_integer(5); /*       m  = number of points to process */
  516.   CHECK_ARG (6, ARRAY_P);    /* c = output array of length 2 */
  517.   c = ARRAY_CONTENTS(ARG_REF(6));
  518.   if ((ARRAY_LENGTH(ARG_REF(6))) != 2) error_bad_range_arg(6);
  519.   mplus = at + m;
  520.   if (mplus > (ARRAY_LENGTH(ARG_REF(1)))) error_bad_range_arg(5);
  521.   if (tc==1)
  522.   { x = 1.0;            /* real part of accumulator */
  523.     y = 0.0;            /* imag part of accumulator */
  524.     for (i=at;i<mplus;i++)
  525.     { temp = ((double) a[i])*x - ((double) b[i])*y;
  526.       y    = ((double) b[i])*x + ((double) a[i])*y;
  527.       x    = temp; }
  528.   }
  529.   else
  530.     error_bad_range_arg(3);
  531.   c[0] = ((REAL) x);        /* mechanism for returning complex number */
  532.   c[1] = ((REAL) y);        /* do not use lists, avoid heap pointer   */
  533.   PRIMITIVE_RETURN (UNSPECIFIC);
  534. }
  535.  
  536. DEFINE_PRIMITIVE ("CS-ARRAY-TO-COMPLEX-ARRAY!", Prim_cs_array_to_complex_array, 3, 3, 0)
  537. { long n,n2,n2_1, i;
  538.   REAL *a, *b,*c;
  539.   PRIMITIVE_HEADER (3);
  540.   CHECK_ARG (1, ARRAY_P);
  541.   CHECK_ARG (2, ARRAY_P);
  542.   CHECK_ARG (3, ARRAY_P);
  543.   a = ARRAY_CONTENTS(ARG_REF(1));
  544.   n = ARRAY_LENGTH(ARG_REF(1));
  545.   b = ARRAY_CONTENTS(ARG_REF(2));
  546.   c = ARRAY_CONTENTS(ARG_REF(3));
  547.   if (n!=(ARRAY_LENGTH(ARG_REF(2)))) error_bad_range_arg(2);
  548.   if (n!=(ARRAY_LENGTH(ARG_REF(3)))) error_bad_range_arg(3);
  549.   b[0] = a[0];
  550.   c[0] = 0.0;
  551.   n2 = n/2;            /* integer division truncates down */
  552.   n2_1 = n2+1;
  553.   if (2*n2 == n)        /* even length, n2 is only real */
  554.   { b[n2]  = a[n2];  c[n2]  = 0.0; }
  555.   else /* odd length, make the loop include the n2 index */
  556.   { n2   = n2+1;
  557.     n2_1 = n2; }
  558.   for (i=1; i<n2; i++)   { b[i] = a[i];
  559.                c[i] = a[n-i]; }
  560.   for (i=n2_1; i<n; i++) { b[i] =  a[n-i];
  561.                c[i] = (-a[i]); }
  562.   PRIMITIVE_RETURN (UNSPECIFIC);
  563. }
  564.  
  565. DEFINE_PRIMITIVE ("CS-ARRAY-MULTIPLY-INTO-SECOND-ONE!", Prim_cs_array_multiply_into_second_one, 2, 2, 0)
  566. { long n,n2;
  567.   REAL *a, *b;
  568.   void cs_array_multiply_into_second_one();
  569.   PRIMITIVE_HEADER (2);
  570.   CHECK_ARG (1, ARRAY_P);
  571.   CHECK_ARG (2, ARRAY_P);
  572.   a = ARRAY_CONTENTS(ARG_REF(1));
  573.   n = ARRAY_LENGTH(ARG_REF(1));
  574.   b = ARRAY_CONTENTS(ARG_REF(2));
  575.   if (n!=(ARRAY_LENGTH(ARG_REF(2)))) error_bad_range_arg(2);
  576.   n2 = n/2;            /* integer division truncates down */
  577.   cs_array_multiply_into_second_one(a,b, n,n2);
  578.   PRIMITIVE_RETURN (UNSPECIFIC);
  579. }
  580.  
  581. void
  582. cs_array_multiply_into_second_one (a,b, n,n2)
  583.      REAL *a, *b;
  584.      long n,n2;
  585. {
  586.   REAL temp;
  587.   long i,ni;
  588.   b[0]   = a[0]  * b[0];
  589.   if (2*n2 == n)        /* even length, n2 is only real */
  590.     b[n2]  = a[n2] * b[n2];
  591.   else
  592.     n2 = n2+1; /* odd length, make the loop include the n2 index */
  593.   for (i=1; i<n2; i++)
  594.     {
  595.       ni = n-i;
  596.       temp   = a[i]*b[i]   -  a[ni]*b[ni]; /* real part */
  597.       b[ni]  = a[i]*b[ni]  +  a[ni]*b[i]; /*  imag part */
  598.       b[i]   = temp;
  599.     }
  600. }
  601.  
  602. DEFINE_PRIMITIVE ("CS-ARRAY-DIVIDE-INTO-XXX!", Prim_cs_array_divide_into_xxx, 4, 4, 0)
  603. {
  604.   long n,n2, one_or_two;
  605.   REAL *a, *b, inf;
  606.   void cs_array_divide_into_z();
  607.   PRIMITIVE_HEADER (4);
  608.   CHECK_ARG (1, ARRAY_P);
  609.   CHECK_ARG (2, ARRAY_P);
  610.   inf = (arg_real (3));
  611.   /* where to store result of division */
  612.   one_or_two = (arg_nonnegative_integer (4));
  613.   a = ARRAY_CONTENTS(ARG_REF(1));
  614.   b = ARRAY_CONTENTS(ARG_REF(2));
  615.   n = ARRAY_LENGTH(ARG_REF(1));
  616.   if (n!=(ARRAY_LENGTH(ARG_REF(2)))) error_bad_range_arg(2);
  617.   n2 = n/2;            /* integer division truncates down */
  618.   if (one_or_two == 1)
  619.     cs_array_divide_into_z(a,b, a,  n,n2, inf);
  620.   else if (one_or_two == 2)
  621.     cs_array_divide_into_z(a,b, b,  n,n2, inf);
  622.   else
  623.     error_bad_range_arg(4);
  624.   PRIMITIVE_RETURN (UNSPECIFIC);
  625. }
  626.  
  627. void
  628. cs_array_divide_into_second_one (a,b, n,n2,inf)    /* used in image.c */
  629.      REAL *a,*b, inf;
  630.      long n,n2;
  631. {
  632.   void cs_array_divide_into_z ();
  633.   cs_array_divide_into_z (a,b, b, n,n2,inf);
  634. }
  635.  
  636. void
  637. cs_array_divide_into_z (a,b, z, n,n2, inf) /* z can be either a or b */
  638.      REAL *a,*b,*z, inf;
  639.      long n,n2;
  640. { long i,ni;
  641.   REAL temp, radius;
  642.  
  643.   if (b[0] == 0.0)
  644.     if (a[0] == 0.0) z[0] = 1.0;
  645.     else             z[0] = a[0] * inf;
  646.   else               z[0] = a[0] / b[0];
  647.  
  648.   if (2*n2 == n)        /* even length, n2 is only real */
  649.     if (b[n2] == 0.0)
  650.       if (a[n2] == 0.0) z[n2] = 1.0;
  651.       else              z[n2] = a[n2] * inf;
  652.     else                z[n2] = a[n2] / b[n2];
  653.   else
  654.     n2 = n2+1; /* odd length, make the loop include the n2 index */
  655.  
  656.   for (i=1; i<n2; i++)
  657.   { ni = n-i;
  658.     radius = b[i]*b[i] + b[ni]*b[ni]; /* b^2 denominator = real^2 + imag^2 */
  659.  
  660.     if (radius == 0.0) {
  661.       if (a[i]  == 0.0) z[i]  = 1.0;
  662.       else              z[i]  = a[i] * inf;
  663.       if (a[ni] == 0.0) z[ni] = 1.0;
  664.       else              z[ni] = a[ni] * inf; }
  665.     else {
  666.       temp  = a[i]*b[i]    +  a[ni]*b[ni];
  667.       z[ni] = (a[ni]*b[i]  -  a[i]*b[ni]) / radius; /* imag part */
  668.       z[i]  = temp                        / radius; /* real part */
  669.     }}
  670. }
  671.  
  672. /* ARRAY-UNARY-FUNCTION!
  673.    apply unary-function elementwise on array
  674.    Available functions : */
  675.  
  676. void
  677. REALabs (a,b)
  678.      REAL *a,*b;
  679. {
  680.   (*b) = ( (REAL) fabs( (double) (*a)) );
  681. }
  682.  
  683. void
  684. REALexp (a,b)
  685.      REAL *a,*b;
  686. {
  687.   fast double y;
  688.   if ((y = exp((double) (*a))) == HUGE)
  689.     error_bad_range_arg (1);    /* OVERFLOW */
  690.   (*b) = ((REAL) y);
  691. }
  692.  
  693. void
  694. REALlog (a,b)
  695.      REAL *a,*b;
  696. {
  697.   if ((*a) < 0.0)
  698.     error_bad_range_arg(1);    /* log(negative) */
  699.   (*b) = ( (REAL) log( (double) (*a)) );
  700. }
  701.  
  702. void
  703. REALtruncate (a,b)
  704.      REAL *a,*b;        /* towards zero */
  705. {
  706.   double integral_part, modf();
  707.   modf( ((double) (*a)), &integral_part);
  708.   (*b) = ( (REAL) integral_part);
  709. }
  710.  
  711. void
  712. REALround (a,b)
  713.      REAL *a,*b;        /* towards nearest integer */
  714. {
  715.   double integral_part, modf();
  716.   if ((*a) >= 0.0)        /* It may be faster to look at the sign
  717.                    of mantissa, and dispatch */
  718.     modf( ((double) ((*a)+0.5)), &integral_part);
  719.   else
  720.     modf( ((double) ((*a)-0.5)), &integral_part);
  721.   (*b) = ( (REAL) integral_part);
  722. }
  723.  
  724. void
  725. REALsquare (a,b)
  726.      REAL *a,*b;
  727. {
  728.   (*b) = ( (REAL) ((*a) * (*a)) );
  729. }
  730.  
  731. void
  732. REALsqrt (a,b)
  733.      REAL *a,*b;
  734. {
  735.   if ((*a) < 0.0)
  736.     error_bad_range_arg(1);    /* sqrt(negative) */
  737.   (*b) = ( (REAL) sqrt( (double) (*a)) );
  738. }
  739.  
  740. void
  741. REALsin (a,b)
  742.      REAL *a,*b;
  743. {
  744.   (*b) = ( (REAL) sin( (double) (*a)) );
  745. }
  746.  
  747. void
  748. REALcos (a,b)
  749.      REAL *a,*b;
  750. {
  751.   (*b) = ( (REAL) cos( (double) (*a)) );
  752. }
  753.  
  754. void
  755. REALtan (a,b)
  756.      REAL *a,*b;
  757. {
  758.   (*b) = ( (REAL) tan( (double) (*a)) );
  759. }
  760.  
  761. void
  762. REALasin (a,b)
  763.      REAL *a,*b;
  764. {
  765.   (*b) = ( (REAL) asin( (double) (*a)) );
  766. }
  767.  
  768. void
  769. REALacos (a,b)
  770.      REAL *a,*b;
  771. {
  772.   (*b) = ( (REAL) acos( (double) (*a)) );
  773. }
  774.  
  775. void
  776. REALatan (a,b)
  777.      REAL *a,*b;
  778. {
  779.   (*b) = ( (REAL) atan( (double) (*a)) );
  780. }
  781.  
  782. void
  783. REALgamma (a,b)
  784.      REAL *a,*b;
  785. {
  786.   fast double y;
  787.   if ((y = gamma(((double) (*a)))) > LN_MAXDOUBLE)
  788.     error_bad_range_arg(1);    /* gamma( non-positive integer ) */
  789.   (*b) = ((REAL) (signgam * exp(y))); /* see HPUX Section 3 */
  790. }
  791.  
  792. void
  793. REALerf (a,b)
  794.      REAL *a,*b;
  795. {
  796.   (*b) = ( (REAL) erf((double) (*a)) );
  797. }
  798.  
  799. void
  800. REALerfc (a,b)
  801.      REAL *a,*b;
  802. {
  803.   (*b) = ( (REAL) erfc((double) (*a)) );
  804. }
  805.  
  806. void
  807. REALbessel1 (order,a,b)        /* Bessel of first kind */
  808.      long order;
  809.      REAL *a,*b;
  810. {
  811.   if (order == 0)
  812.     (*b) = ( (REAL) j0((double) (*a)) );
  813.   if (order == 1)
  814.     (*b) = ( (REAL) j1((double) (*a)) );
  815.   else
  816.     (*b) = ( (REAL) jn(((int) order), ((double) (*a))) );
  817. }
  818.  
  819. void
  820. REALbessel2 (order,a,b)        /* Bessel of second kind */
  821.      long order;
  822.      REAL *a,*b;
  823. {
  824.   if ((*a) <= 0.0)
  825.     error_bad_range_arg(1);    /* Blows Up */
  826.   if (order == 0)
  827.     (*b) = ( (REAL) y0((double) (*a)) );
  828.   if (order == 1)
  829.     (*b) = ( (REAL) y1((double) (*a)) );
  830.   else
  831.     (*b) = ( (REAL) yn(((int) order), ((double) (*a))) );
  832. }
  833.  
  834. /* Table to store the available unary-functions.
  835.    Also some binary functions at the end -- not available right now.
  836.    The (1 and 2)s denote the numofargs (1 for unary 2 for binary) */
  837.  
  838. struct array_func_table
  839. {
  840.   long numofargs;
  841.   void (*func)();
  842. }
  843. Array_Function_Table [] =
  844. {
  845.   1, REALabs,            /*0*/
  846.   1, REALexp,            /*1*/
  847.   1, REALlog,            /*2*/
  848.   1, REALtruncate,        /*3*/
  849.   1, REALround,            /*4*/
  850.   1, REALsquare,        /*5*/
  851.   1, REALsqrt,            /*6*/
  852.   1, REALsin,            /*7*/
  853.   1, REALcos,            /*8*/
  854.   1, REALtan,            /*9*/
  855.   1, REALasin,            /*10*/
  856.   1, REALacos,            /*11*/
  857.   1, REALatan,            /*12*/
  858.   1, REALgamma,            /*13*/
  859.   1, REALerf,            /*14*/
  860.   1, REALerfc,            /*15*/
  861.   2, REALbessel1,        /*16*/
  862.   2, REALbessel2        /*17*/
  863.   };
  864.  
  865. #define MAX_ARRAY_FUNCTC 17
  866.  
  867. /* array-unary-function!    could be called        array-operation-1!
  868.    following the naming convention for other similar procedures
  869.    but it is specialized to mappings only, so we have special name. */
  870.  
  871. DEFINE_PRIMITIVE ("ARRAY-UNARY-FUNCTION!", Prim_array_unary_function, 2,2, 0)
  872. {
  873.   long n, i;
  874.   REAL *a,*b;
  875.   long tc;
  876.   void (*f)();
  877.   PRIMITIVE_HEADER (2);
  878.   CHECK_ARG (1, ARRAY_P);    /*      a  = input (and output) array    */
  879.   CHECK_ARG (2, FIXNUM_P);    /*      tc = type code                   */
  880.   tc = arg_nonnegative_integer(2);
  881.   if (tc > MAX_ARRAY_FUNCTC) error_bad_range_arg(2);
  882.   f = ((Array_Function_Table[tc]).func);
  883.   if (1 != (Array_Function_Table[tc]).numofargs) error_wrong_type_arg(2);
  884.   /* check it is a unary function */
  885.   a = ARRAY_CONTENTS(ARG_REF(1));
  886.   b = a;
  887.   n = ARRAY_LENGTH(ARG_REF(1));
  888.   for (i=0; i<n; i++)
  889.     (*f) ( &(a[i]), &(b[i]) );    /* a into b */
  890.   PRIMITIVE_RETURN (UNSPECIFIC);
  891. }
  892.  
  893. /* Accumulate
  894.    using combinators              +  or  *
  895.    corresponding type codes       0      1 */
  896.  
  897. DEFINE_PRIMITIVE ("SUBARRAY-ACCUMULATE", Prim_subarray_accumulate, 4,4, 0)
  898. {
  899.   long at,m,mplus, tc, i;
  900.   REAL *a;
  901.   double result;
  902.   PRIMITIVE_HEADER (4);
  903.   CHECK_ARG (1, ARRAY_P);    /*           a = input array                 */
  904.   a  = ARRAY_CONTENTS(ARG_REF(1));
  905.   tc = arg_nonnegative_integer(2); /*       tc = type code 0 or 1            */
  906.   at = arg_nonnegative_integer(3); /*       at = starting index              */
  907.   m  = arg_nonnegative_integer(4); /*       m  = number of points to process */
  908.   mplus = at + m;
  909.   if (mplus > (ARRAY_LENGTH(ARG_REF(1)))) error_bad_range_arg(4);
  910.   if (tc == 0)
  911.     {
  912.       result = 0.0;
  913.       for (i=at;i<mplus;i++)
  914.     result = result + ((double) a[i]);
  915.     }
  916.   else if (tc == 1)
  917.     {
  918.       result = 1.0;
  919.       for (i=at;i<mplus;i++)
  920.     result = result * ((double) a[i]);
  921.     }
  922.   else
  923.     error_bad_range_arg (2);
  924.   PRIMITIVE_RETURN (double_to_flonum (result));
  925. }
  926.  
  927. /* The following searches for value within tolerance
  928.    starting from index=from in array.
  929.    Returns first index where match occurs (useful for finding zeros). */
  930.  
  931. DEFINE_PRIMITIVE ("ARRAY-SEARCH-VALUE-TOLERANCE-FROM", Prim_array_search_value_tolerance_from, 4, 4, 0)
  932. {
  933.   SCHEME_OBJECT array;
  934.   fast long length;
  935.   fast REAL * a;
  936.   fast REAL value;        /* value to search for */
  937.   fast double tolerance;    /* tolerance allowed */
  938.   PRIMITIVE_HEADER (4);
  939.   CHECK_ARG (1, ARRAY_P);
  940.   array = (ARG_REF (1));
  941.   length = (ARRAY_LENGTH (array));
  942.   a = (ARRAY_CONTENTS (array));
  943.   value = (arg_real (2));
  944.   tolerance = (arg_real (3));
  945.   {
  946.     fast long i;
  947.     for (i = (arg_index_integer (4, length)); i<length; i+=1)
  948.       if (tolerance >= (fabs ((double) ((a [i]) - value))))
  949.     PRIMITIVE_RETURN (LONG_TO_UNSIGNED_FIXNUM (i));
  950.   }
  951.   PRIMITIVE_RETURN (SHARP_F);
  952. }
  953.  
  954. DEFINE_PRIMITIVE ("SUBARRAY-MIN-MAX-INDEX", Prim_subarray_min_max_index, 3, 3, 0)
  955. {
  956.   PRIMITIVE_HEADER (3);
  957.   CHECK_ARG (1, ARRAY_P);
  958.   {
  959.     REAL * a = (ARRAY_CONTENTS (ARG_REF (1)));
  960.     long at = (arg_nonnegative_integer (2)); /* starting index */
  961.     long m  = (arg_nonnegative_integer (3)); /* number of points to process */
  962.     long mplus = (at + m);
  963.     long nmin;
  964.     long nmax;
  965.     if (mplus > (ARRAY_LENGTH (ARG_REF (1))))
  966.       error_bad_range_arg (3);
  967.     
  968.     C_Array_Find_Min_Max ((& (a [at])), m, (&nmin), (&nmax));
  969.     nmin = nmin + at;        /* offset appropriately */
  970.     nmax = nmax + at;
  971.     
  972.     PRIMITIVE_RETURN
  973.       (cons ((LONG_TO_FIXNUM (nmin)),
  974.          (cons ((LONG_TO_FIXNUM (nmax)),
  975.             EMPTY_LIST))));
  976.   }
  977. }
  978.  
  979. void
  980. C_Array_Find_Min_Max (x, n, nmin, nmax)
  981.      fast REAL * x;
  982.      fast long n;
  983.      long * nmin;
  984.      long * nmax;
  985. { REAL *xold = x;
  986.   register REAL xmin, xmax;
  987.   register long nnmin, nnmax;
  988.   register long count;
  989.  
  990.   nnmin = nnmax = 0;
  991.   xmin = xmax = *x++;
  992.   n--;
  993.   count = 1;
  994.   if(n>0)
  995.   {
  996.     do {
  997.       if(*x < xmin) {
  998.     nnmin = count++ ;
  999.     xmin = *x++ ;
  1000.       } else if(*x > xmax) {
  1001.     nnmax = count++ ;
  1002.     xmax = *x++ ;
  1003.       } else {
  1004.     count++ ;
  1005.     x++ ;
  1006.       }
  1007.     } while( --n > 0 ) ;
  1008.   }
  1009.   *nmin = nnmin ;
  1010.   *nmax = nnmax ;
  1011. }
  1012.  
  1013.  
  1014. /* array-average
  1015.    can be done with (array-reduce +) and division by array-length.
  1016.    But there is also this C primitive.
  1017.    Keep it around, may be useful someday. */
  1018.  
  1019. /* Computes the average in pieces, so as to reduce
  1020.    roundoff smearing in cumulative sum.
  1021.    example= first huge positive numbers, then small nums,
  1022.    then huge negative numbers. */
  1023.  
  1024. static void
  1025. C_Array_Find_Average (Array, Length, pAverage)
  1026.      long Length;
  1027.      REAL * Array;
  1028.      REAL * pAverage;
  1029. {
  1030.   long i;
  1031.   long array_index;
  1032.   REAL average_n, sum;
  1033.  
  1034.   average_n = 0.0;
  1035.   array_index = 0;
  1036.   while (array_index < Length)
  1037.     {
  1038.       sum = 0.0;
  1039.       for (i=0;((array_index<Length) && (i<2000));i++) {
  1040.     sum += Array[array_index];
  1041.     array_index++;
  1042.       }
  1043.       average_n += (sum / ((REAL) Length));
  1044.     }
  1045.   (*pAverage) = average_n;
  1046.   return;
  1047. }
  1048.  
  1049. DEFINE_PRIMITIVE ("ARRAY-AVERAGE", Prim_array_find_average, 1, 1, 0)
  1050. {
  1051.   SCHEME_OBJECT array;
  1052.   REAL average;
  1053.   PRIMITIVE_HEADER (1);
  1054.   CHECK_ARG (1, ARRAY_P);
  1055.   array = (ARG_REF (1));
  1056.   C_Array_Find_Average
  1057.     ((ARRAY_CONTENTS (array)),
  1058.      (ARRAY_LENGTH (array)),
  1059.      (&average));
  1060.   PRIMITIVE_RETURN (double_to_flonum ((double) average));
  1061. }
  1062.  
  1063. void
  1064. C_Array_Make_Histogram (Array, Length, Histogram, npoints)
  1065.      REAL Array[], Histogram[];
  1066.      long Length, npoints;
  1067. {
  1068.   REAL Max, Min, Offset, Scale;
  1069.   long i, nmin, nmax, index;
  1070.   C_Array_Find_Min_Max (Array, Length, (&nmin), (&nmax));
  1071.   Min = (Array [nmin]);
  1072.   Max = (Array [nmax]);
  1073.   Find_Offset_Scale_For_Linear_Map
  1074.     (Min, Max, 0.0, ((REAL) npoints), (&Offset), (&Scale));
  1075.   for (i = 0; (i < npoints); i += 1)
  1076.     (Histogram [i]) = 0.0;
  1077.   for (i = 0; (i < Length); i += 1)
  1078.     {
  1079.       /* Everything from 0 to 1 maps to bin 0, and so on */
  1080.       index = ((long) (floor ((double) ((Scale * (Array [i])) + Offset))));
  1081.       /* max that won't floor to legal array index */
  1082.       if (index == npoints)
  1083.     index = (index - 1);
  1084.       (Histogram [index]) += 1.0;
  1085.     }
  1086.   return;
  1087. }
  1088.  
  1089. DEFINE_PRIMITIVE ("ARRAY-MAKE-HISTOGRAM", Prim_array_make_histogram, 2, 2, 0)
  1090. {
  1091.   PRIMITIVE_HEADER (2);
  1092.   CHECK_ARG (1, ARRAY_P);
  1093.   {
  1094.     fast SCHEME_OBJECT array = (ARG_REF (1));
  1095.     long length = (ARRAY_LENGTH (array));
  1096.     long npoints = (arg_integer_in_range (2, 1, ((2 * length) + 1)));
  1097.     SCHEME_OBJECT result = (allocate_array (npoints));
  1098.     C_Array_Make_Histogram
  1099.       ((ARRAY_CONTENTS (array)),
  1100.        length,
  1101.        (ARRAY_CONTENTS (result)),
  1102.        npoints);
  1103.     PRIMITIVE_RETURN (result);
  1104.   }
  1105. }
  1106.  
  1107. /* The following geometrical map is slightly tricky. */
  1108. void
  1109. Find_Offset_Scale_For_Linear_Map (Min, Max, New_Min, New_Max, Offset, Scale)
  1110.      REAL Min, Max, New_Min, New_Max, *Offset, *Scale;
  1111. {
  1112.   if (Min != Max)
  1113.     {
  1114.       (*Scale)  = ((New_Max - New_Min) / (Max - Min));
  1115.       (*Offset) = (New_Min - ((*Scale) * Min));
  1116.     }
  1117.   else
  1118.     {
  1119.       (*Scale) =
  1120.     ((Max == 0.0)
  1121.      ? 0.0
  1122.      : (0.25 * (mabs ((New_Max - New_Min) / Max))));
  1123.       (*Offset) = ((New_Max + New_Min) / 2.0);
  1124.     }
  1125.   return;
  1126. }
  1127.  
  1128.  
  1129. DEFINE_PRIMITIVE ("ARRAY-CLIP-MIN-MAX!", Prim_array_clip_min_max, 3, 3, 0)
  1130. {
  1131.   PRIMITIVE_HEADER (3);
  1132.   CHECK_ARG (1, ARRAY_P);
  1133.   {
  1134.     SCHEME_OBJECT array = (ARG_REF (1));
  1135.     REAL xmin = (arg_real (2));
  1136.     REAL xmax = (arg_real (3));
  1137.     long Length = (ARRAY_LENGTH (array));
  1138.     REAL * From_Here = (ARRAY_CONTENTS (array));
  1139.     REAL * To_Here = From_Here;
  1140.     long i;
  1141.     if (xmin > xmax)
  1142.       error_bad_range_arg (3);
  1143.     for (i = 0; (i < Length); i += 1)
  1144.       {
  1145.     if ((*From_Here) < xmin)
  1146.       (*To_Here++) = xmin;
  1147.     else if ((*From_Here) > xmax)
  1148.       (*To_Here++) = xmax;
  1149.     else
  1150.       (*To_Here++) = (*From_Here);
  1151.     From_Here += 1;
  1152.       }
  1153.   }
  1154.   PRIMITIVE_RETURN (UNSPECIFIC);
  1155. }
  1156.  
  1157. /* complex-array-operation-1!
  1158.    groups together procedures that use 1 complex-array
  1159.    and store the result in place */
  1160.  
  1161. DEFINE_PRIMITIVE ("COMPLEX-ARRAY-OPERATION-1!", Prim_complex_array_operation_1, 3,3, 0)
  1162. { long n, i, opcode;
  1163.   REAL *a, *b;
  1164.   void complex_array_to_polar(), complex_array_exp(), complex_array_sqrt();
  1165.   void complex_array_sin(), complex_array_cos();
  1166.   void complex_array_asin(), complex_array_acos(), complex_array_atan();
  1167.   PRIMITIVE_HEADER (3);
  1168.   CHECK_ARG (1, FIXNUM_P);    /* operation opcode */
  1169.   CHECK_ARG (2, ARRAY_P);    /* input array -- n      real part         */
  1170.   CHECK_ARG (3, ARRAY_P);    /* input array -- n      imag part         */
  1171.   n = ARRAY_LENGTH(ARG_REF(2));
  1172.   if (n != ARRAY_LENGTH(ARG_REF(3))) error_bad_range_arg(3);
  1173.   a  = ARRAY_CONTENTS(ARG_REF(2)); /*  real part */
  1174.   b  = ARRAY_CONTENTS(ARG_REF(3)); /*  imag part */
  1175.   opcode = arg_nonnegative_integer(1);
  1176.   if (opcode==1)
  1177.     complex_array_to_polar(a,b,n);
  1178.   else if (opcode==2)
  1179.     complex_array_exp(a,b,n);
  1180.   else if (opcode==3)
  1181.     complex_array_sqrt(a,b,n);
  1182.  
  1183.   else if (opcode==4)
  1184.     complex_array_sin(a,b,n);
  1185.   else if (opcode==5)
  1186.     complex_array_cos(a,b,n);
  1187.   /* for tan(z) use sin(z)/cos(z) */
  1188.   
  1189.   else if (opcode==6)
  1190.     complex_array_asin(a,b,n);
  1191.   else if (opcode==7)
  1192.     complex_array_acos(a,b,n);
  1193.   else if (opcode==8)
  1194.     complex_array_atan(a,b,n);
  1195.   
  1196.   else
  1197.     error_bad_range_arg(1);    /* illegal opcode */
  1198.   PRIMITIVE_RETURN (UNSPECIFIC);
  1199. }
  1200.  
  1201. void
  1202. complex_array_to_polar (a,b,n)
  1203.      REAL *a,*b;
  1204.      long n;
  1205. {
  1206.   long i;
  1207.   double x,y, temp;
  1208.   for (i=0; i<n; i++)
  1209.     {
  1210.       x = (double) a[i];
  1211.       y = (double) b[i];
  1212.       temp = sqrt(x*x + y*y);
  1213.       if (temp == 0.0)
  1214.     b[i] = 0.0;        /* choose angle = 0    for x,y=0,0 */
  1215.       else
  1216.     b[i] = (REAL) atan2(y,x);
  1217.       a[i]   = (REAL) temp;
  1218.     }
  1219. }
  1220.  
  1221. void
  1222. complex_array_exp (a,b,n)
  1223.      REAL *a,*b;
  1224.      long n;
  1225. {
  1226.   long i;
  1227.   double x,y, temp;
  1228.  
  1229.   for (i=0; i<n; i++)
  1230.     {
  1231.       x = (double) a[i];
  1232.       y = (double) b[i];
  1233.       if ((temp = exp(x)) == HUGE) error_bad_range_arg(2); /* overflow */
  1234.       a[i] = (REAL) (temp*cos(y));
  1235.       b[i] = (REAL) (temp*sin(y));
  1236.     }
  1237. }
  1238.  
  1239. void
  1240. complex_array_sqrt (a,b,n)
  1241.      REAL *a,*b;
  1242.      long n;
  1243. {
  1244.   long i;
  1245.   double x,y, r;
  1246.  
  1247.   for (i=0; i<n; i++)
  1248.     {
  1249.       x = (double) a[i];
  1250.       y = (double) b[i];
  1251.       r = sqrt( x*x + y*y);
  1252.       a[i] = sqrt((r+x)/2.0);
  1253.       if (y>=0.0)
  1254.     b[i] =  sqrt((r-x)/2.0); /* choose principal root */
  1255.       else            /* see Abramowitz (p.17 3.7.27) */
  1256.     b[i] = -sqrt((r-x)/2.0);
  1257.     }
  1258. }
  1259.  
  1260. void
  1261. complex_array_sin (a,b,n)
  1262.      REAL *a,*b;
  1263.      long n;
  1264. {
  1265.   long i;
  1266.   double x, ey,fy;
  1267.   REAL temp;
  1268.  
  1269.   for (i=0; i<n; i++)
  1270.     {
  1271.       x = (double) a[i];
  1272.       ey = exp((double) b[i]);    /* radius should be small to avoid overflow */
  1273.       fy = 1.0/ey;        /* exp(-y) */
  1274.       /* expanded (e(iz)-e(-iz))*(-.5i) formula */
  1275.       temp = (REAL) (sin(x) * (ey + fy) * 0.5);
  1276.       /* see my notes in Abram.p.71 */
  1277.       b[i] = (REAL) (cos(x) * (ey - fy) * 0.5);
  1278.       a[i] = temp;
  1279.     }
  1280. }
  1281.  
  1282. void
  1283. complex_array_cos (a,b,n)
  1284.      REAL *a,*b;
  1285.      long n;
  1286. {
  1287.   long i;
  1288.   double x, ey,fy;
  1289.   REAL temp;
  1290.  
  1291.   for (i=0; i<n; i++)
  1292.     {
  1293.       x = (double) a[i];
  1294.       ey = exp((double) b[i]);    /* radius should be small to avoid overflow */
  1295.       fy = 1.0/ey;        /* exp(-y) */
  1296.       /* expanded (e(iz)+e(-iz))*.5 formula */
  1297.       temp = (REAL) (cos(x) * (ey + fy) * 0.5);
  1298.       b[i] = (REAL) (sin(x) * (fy - ey) * 0.5); /* see my notes in Abram.p.71*/
  1299.       a[i] = temp;
  1300.     }
  1301. }
  1302.  
  1303.  
  1304. void
  1305. complex_array_asin (a,b,n)
  1306.      REAL *a,*b;
  1307.      long n;
  1308. { /* logarithmic formula as in R3.99, about 21ops plus log,atan - see my notes */
  1309.   long i; 
  1310.   double oldx,oldy, x,y, real,imag, r;
  1311.   
  1312.   for (i=0; i<n; i++)
  1313.   {
  1314.     oldx = (double) a[i];
  1315.     oldy = (double) b[i];
  1316.     
  1317.     x = 1.0 - oldx*oldx + oldy*oldy; /* 1 - z*z */
  1318.     y = -2.0 * oldx * oldy;
  1319.     
  1320.     r = sqrt(x*x + y*y);    /* sqrt(1-z*z)  */
  1321.     real = sqrt((r+x)/2.0);
  1322.     if (y>=0.0)
  1323.       imag =  sqrt((r-x)/2.0);    /* choose principal root */
  1324.     else            /* see Abramowitz (p.17 3.7.27) */
  1325.       imag = -sqrt((r-x)/2.0);
  1326.     
  1327.     real = real - oldy;        /* i*z + sqrt(...) */
  1328.     imag = imag + oldx;
  1329.     
  1330.     b[i] = (REAL) (- log (sqrt (real*real + imag*imag))); /* -i*log(...) */
  1331.     a[i] = (REAL) atan2( imag, real); /* chosen angle is okay 
  1332.                      Also 0/0 doesnot occur */
  1333.   }
  1334. }
  1335.  
  1336. void
  1337. complex_array_acos (a,b,n)
  1338.      REAL *a,*b;
  1339.      long n;
  1340. {
  1341.   long i;
  1342.  
  1343.   complex_array_asin (a,b,n);
  1344.   
  1345.   for (i=0; i<n; i++)
  1346.     {
  1347.       a[i] = PI_OVER_2 - a[i];
  1348.       b[i] =           - b[i];
  1349.     }
  1350. }
  1351.   
  1352.  
  1353. void
  1354. complex_array_atan (a,b,n)
  1355.      REAL *a,*b;
  1356.      long n;
  1357. { /* logarithmic formula, expanded, simplified - see my notes */
  1358.   long i; 
  1359.   double x,y, xx, real,imag, d;
  1360.   
  1361.   for (i=0; i<n; i++)
  1362.   {
  1363.     x = (double) a[i];
  1364.     y = (double) b[i];
  1365.     
  1366.     xx = x*x;
  1367.     imag = 1.0 + y;        /* temp var */
  1368.     d  = xx + imag*imag;
  1369.     
  1370.     real = (1 - y*y - xx) / d;
  1371.     imag = (2.0 * x)  / d;
  1372.     
  1373.     b[i] = (REAL) ((log (sqrt (real*real + imag*imag))) / -2.0);
  1374.     a[i] = (atan2 (imag,real)) / 2.0;
  1375.   }
  1376. }
  1377.  
  1378.  
  1379. /* complex-array-operation-1b!
  1380.    groups together procedures that use 1 complex-array & 1 number
  1381.    and store the result in place (e.g. invert 1/x) */
  1382.  
  1383. DEFINE_PRIMITIVE ("COMPLEX-ARRAY-OPERATION-1B!", Prim_complex_array_operation_1b, 4,4, 0)
  1384. {
  1385.   long n, i, opcode;
  1386.   REAL *a, *b, inf;
  1387.   void complex_array_invert ();
  1388.   PRIMITIVE_HEADER (4);
  1389.   CHECK_ARG (1, FIXNUM_P);    /* operation opcode */
  1390.   CHECK_ARG (2, ARRAY_P);    /* input array -- n      real part         */
  1391.   CHECK_ARG (3, ARRAY_P);    /* input array -- n      imag part         */
  1392.   inf = (arg_real (4));        /* User-Provided Infinity */
  1393.   n = ARRAY_LENGTH(ARG_REF(2));
  1394.   if (n != ARRAY_LENGTH(ARG_REF(3))) error_bad_range_arg(3);
  1395.   a  = ARRAY_CONTENTS(ARG_REF(2)); /*  real part */
  1396.   b  = ARRAY_CONTENTS(ARG_REF(3)); /*  imag part */
  1397.   opcode = arg_nonnegative_integer(1);
  1398.   if (opcode==1)
  1399.     complex_array_invert(a,b, n, inf);  /* performs 1/x */
  1400.   else if (opcode==2)
  1401.     error_bad_range_arg(1);    /* illegal opcode */
  1402.   else
  1403.     error_bad_range_arg(1);    /* illegal opcode */
  1404.   PRIMITIVE_RETURN (UNSPECIFIC);
  1405. }
  1406.  
  1407. void
  1408. complex_array_invert (a,b, n, inf)
  1409.      REAL *a,*b, inf;
  1410.      long n;
  1411. {
  1412.   long i;
  1413.   double x,y, r;
  1414.   for (i=0; i<n; i++)
  1415.     {
  1416.       x = (double) a[i];
  1417.       y = (double) b[i];
  1418.       r = (x*x + y*y);
  1419.       if (r==0.0)
  1420.     {
  1421.       a[i] = inf;
  1422.       b[i] = inf;
  1423.     }
  1424.       else
  1425.     {
  1426.       a[i] = (REAL)  x/r;
  1427.       b[i] = (REAL) -y/r;
  1428.     }
  1429.     }
  1430. }
  1431.  
  1432. /* complex-array-operation-1a
  1433.    groups together procedures that use 1 complex-array
  1434.    and store result in a 3rd real array. */
  1435.  
  1436. DEFINE_PRIMITIVE ("COMPLEX-ARRAY-OPERATION-1A", Prim_complex_array_operation_1a, 4,4, 0)
  1437. {
  1438.   long n, i, opcode;
  1439.   REAL *a, *b, *c;
  1440.   void complex_array_magnitude(), complex_array_angle();
  1441.   PRIMITIVE_HEADER (4);
  1442.   CHECK_ARG (1, FIXNUM_P);    /* operation opcode */
  1443.   CHECK_ARG (2, ARRAY_P);    /* input array -- n      real part         */
  1444.   CHECK_ARG (3, ARRAY_P);    /* input array -- n      imag part         */
  1445.   CHECK_ARG (4, ARRAY_P);    /* output array -- n                       */
  1446.   n = ARRAY_LENGTH(ARG_REF(2));
  1447.   if (n != ARRAY_LENGTH(ARG_REF(3))) error_bad_range_arg(3);
  1448.   if (n != ARRAY_LENGTH(ARG_REF(4))) error_bad_range_arg(4);
  1449.   a  = ARRAY_CONTENTS(ARG_REF(2)); /*  real part */
  1450.   b  = ARRAY_CONTENTS(ARG_REF(3)); /*  imag part */
  1451.   c  = ARRAY_CONTENTS(ARG_REF(4)); /*  output    */
  1452.   opcode = arg_nonnegative_integer(1);
  1453.   if (opcode==1)
  1454.     complex_array_magnitude(a,b,c,n);
  1455.   else if (opcode==2)
  1456.     complex_array_angle(a,b,c,n);
  1457.   else
  1458.     error_bad_range_arg(1);    /* illegal opcode */
  1459.   PRIMITIVE_RETURN (UNSPECIFIC);
  1460. }
  1461.  
  1462. void
  1463. complex_array_magnitude (a,b,c,n)
  1464.      REAL *a,*b,*c;
  1465.      long n;
  1466. {
  1467.   long i;
  1468.   for (i=0; i<n; i++)
  1469.     c[i] = (REAL) sqrt( (double) a[i]*a[i] + b[i]*b[i] );
  1470. }
  1471.  
  1472. void
  1473. complex_array_angle (a,b,c,n)
  1474.      REAL *a,*b,*c;
  1475.      long n;
  1476. {
  1477.   long i;
  1478.   for (i=0; i<n; i++)
  1479.   {
  1480.     if ((a[i] == 0.0) && (b[i]==0.0))
  1481.       c[i] = 0.0;        /* choose angle=0 for point (0,0) */
  1482.     else
  1483.       c[i] = (REAL) atan2( (double) b[i], (double) a[i]);
  1484.     /* angle ==   -pi (exclusive) to +pi (inclusive) */
  1485.   }
  1486. }
  1487.  
  1488. DEFINE_PRIMITIVE ("CS-ARRAY-MAGNITUDE!", Prim_cs_array_magnitude, 1, 1, 0)
  1489. {
  1490.   long n, i;
  1491.   REAL *a;
  1492.   void cs_array_magnitude();
  1493.   PRIMITIVE_HEADER (1);
  1494.   CHECK_ARG (1, ARRAY_P);
  1495.   a = ARRAY_CONTENTS(ARG_REF(1)); /* input cs-array */
  1496.   n = ARRAY_LENGTH(ARG_REF(1));    /* becomes a standard array on return  */
  1497.   cs_array_magnitude(a,n);
  1498.   PRIMITIVE_RETURN (UNSPECIFIC);
  1499. }
  1500.  
  1501. /* result is a standard array (even signal, real data) */
  1502. void
  1503. cs_array_magnitude (a,n)
  1504.      REAL *a;
  1505.      long n;
  1506. {
  1507.   long i, n2, ni;
  1508.   n2 = n/2;            /* integer division truncates down */
  1509.   a[0]  = (REAL) fabs((double) a[0]); /*   imag=0 */
  1510.   if (2*n2 == n)        /* even length, n2 is only real */
  1511.     a[n2] = (REAL) fabs((double) a[n2]); /*  imag=0 */
  1512.   else
  1513.     /* odd length, make the loop include the n2 index */
  1514.     n2 = n2+1;
  1515.   for (i=1; i<n2; i++)
  1516.     {
  1517.       ni = n-i;
  1518.       a[i]   = (REAL)  sqrt( (double) a[i]*a[i] + (double) a[ni]*a[ni] );
  1519.       a[ni]  = a[i];        /* even signal */
  1520.     }
  1521. }
  1522.  
  1523. /* Rectangular and Polar
  1524.    A cs-array has even magnitude and odd angle (almost)
  1525.    hence a polar cs-array stores magnitude in the first half (real part)
  1526.    and angle in the second half (imag part)
  1527.    except for a[0] real-only and a[n2] (n even)
  1528.    The angle of a[0] is either 0 (pos. sign) or pi (neg. sign),
  1529.    but there is no place in an n-point cs-array to store this, so
  1530.    a[0] and a[n2] when n even are left unchanged  when going polar.
  1531.    as opposed to taking their absolute values, for magnitude. */
  1532.  
  1533. DEFINE_PRIMITIVE ("CS-ARRAY-TO-POLAR!", Prim_cs_array_to_polar, 1,1, 0)
  1534. {
  1535.   long n, i;
  1536.   REAL *a;
  1537.   void cs_array_to_polar();
  1538.   PRIMITIVE_HEADER (1);
  1539.   CHECK_ARG (1, ARRAY_P);
  1540.   a  = ARRAY_CONTENTS(ARG_REF(1));
  1541.   n = ARRAY_LENGTH(ARG_REF(1));
  1542.   cs_array_to_polar(a,n);
  1543.   PRIMITIVE_RETURN (UNSPECIFIC);
  1544. }
  1545.  
  1546. void
  1547. cs_array_to_polar (a,n)
  1548.      REAL *a;
  1549.      long n;
  1550. {
  1551.   long i, n2;
  1552.   double real, imag;        /* temporary variables */
  1553.   n2 = n/2;            /* integer division truncates down */
  1554.   /* a[0] stores both magnitude and angle
  1555.      (pos. sign angle=0 , neg. sign angle=pi) */
  1556.   if (2*n2 == n)        /* even length, n2 is only real */
  1557.     ;                /* a[n2] stores sign information like a[0] */
  1558.   else
  1559.     /* odd length, make the loop include the n2 index */
  1560.     n2 = n2+1;
  1561.   for (i=1; i<n2; i++)
  1562.     {
  1563.       real = (double) a[i];
  1564.       imag = (double) a[n-i];
  1565.       a[i]   = (REAL)  sqrt( real*real + imag*imag );
  1566.       if (a[i] == 0.0)
  1567.     a[n-i] = 0.0;
  1568.       else
  1569.     a[n-i] = (REAL) atan2( imag, real );
  1570.     }
  1571. }
  1572.  
  1573. DEFINE_PRIMITIVE ("CS-ARRAY-TO-RECTANGULAR!", Prim_cs_array_to_rectangular, 1,1, 0)
  1574. {
  1575.   long n,n2, i;
  1576.   double magn,angl;        /* temporary variables */
  1577.   REAL *a;
  1578.   PRIMITIVE_HEADER (1);
  1579.   CHECK_ARG (1, ARRAY_P);
  1580.   a  = ARRAY_CONTENTS(ARG_REF(1));
  1581.   n = ARRAY_LENGTH(ARG_REF(1));
  1582.   n2 = n/2;            /* integer division truncates down */
  1583.   ;                /* a[0] is okay */
  1584.   if (2*n2 == n)        /* even length, n2 is real only */
  1585.     ;                /* a[n2] is okay */
  1586.   else
  1587.     n2 = n2+1; /* odd length, make the loop include the n2 index */
  1588.   for (i=1; i<n2; i++)
  1589.   { magn = (double) a[i];
  1590.     angl = (double) a[n-i];
  1591.     a[i]   = (REAL)  magn * cos(angl);
  1592.     a[n-i] = (REAL)  magn * sin(angl); }
  1593.   PRIMITIVE_RETURN (UNSPECIFIC);
  1594. }
  1595.  
  1596. /* Convolution in the Time-Domain */
  1597. /* In the following macro
  1598.    To1 and To2 should be (Length1-1) and (Length2-1) respectively. */
  1599.  
  1600. #define C_Convolution_Point_Macro(X, Y, To1, To2, N, Result)        \
  1601. {                                    \
  1602.   long Min_of_N_To1 = (min ((N), (To1)));                \
  1603.   long mi, N_minus_mi;                            \
  1604.   REAL Sum = 0.0;                            \
  1605.   for (mi = (max (0, ((N) - (To2)))), N_minus_mi = ((N) - mi);        \
  1606.        (mi <= Min_of_N_To1);                        \
  1607.        mi += 1, N_minus_mi -= 1)                    \
  1608.     Sum += ((X [mi]) * (Y [N_minus_mi]));                \
  1609.   (Result) = Sum;                            \
  1610. }
  1611.  
  1612. DEFINE_PRIMITIVE ("CONVOLUTION-POINT", Prim_convolution_point, 3, 3, 0)
  1613. {
  1614.   long Length1, Length2, N;
  1615.   REAL *Array1, *Array2;
  1616.   REAL C_Result;
  1617.   PRIMITIVE_HEADER (3);
  1618.   CHECK_ARG (1, ARRAY_P);
  1619.   CHECK_ARG (2, ARRAY_P);
  1620.   Length1 = ARRAY_LENGTH (ARG_REF (1));
  1621.   Length2 = ARRAY_LENGTH (ARG_REF (2));
  1622.   N = (arg_nonnegative_integer (3));
  1623.   Array1 = ARRAY_CONTENTS (ARG_REF (1));
  1624.   Array2 = ARRAY_CONTENTS (ARG_REF (2));
  1625.   C_Convolution_Point_Macro(Array1, Array2, Length1-1, Length2-1, N, C_Result);
  1626.   PRIMITIVE_RETURN (double_to_flonum ((double) C_Result));
  1627. }
  1628.  
  1629. DEFINE_PRIMITIVE ("ARRAY-CONVOLUTION-IN-TIME!", Prim_array_convolution_in_time, 3, 3, 0)
  1630. {
  1631.   long n,m,l, n_1,m_1, i;
  1632.   REAL *a,*b,*c;
  1633.   PRIMITIVE_HEADER (3);
  1634.   CHECK_ARG (1, ARRAY_P);
  1635.   CHECK_ARG (2, ARRAY_P);
  1636.   CHECK_ARG (3, ARRAY_P);
  1637.   a = ARRAY_CONTENTS(ARG_REF(1));
  1638.   b = ARRAY_CONTENTS(ARG_REF(2));
  1639.   c = ARRAY_CONTENTS(ARG_REF(3));
  1640.   n = ARRAY_LENGTH(ARG_REF(1));
  1641.   m = ARRAY_LENGTH(ARG_REF(2));
  1642.   l = n+m-1;            /* resulting length */
  1643.   if (l != ARRAY_LENGTH(ARG_REF(3))) error_bad_range_arg(3);
  1644.   n_1 = n-1; m_1 = m-1;
  1645.   for (i=0; i<l; i++)
  1646.   { C_Convolution_Point_Macro(a, b, n_1, m_1, i, c[i]); }
  1647.   PRIMITIVE_RETURN (UNSPECIFIC);
  1648. }
  1649.  
  1650. DEFINE_PRIMITIVE ("ARRAY-MULTIPLY-INTO-SECOND-ONE!", Prim_array_multiply_into_second_one, 2, 2, 0)
  1651. {
  1652.   long Length, i;
  1653.   REAL *To_Here;
  1654.   REAL *From_Here_1, *From_Here_2;
  1655.   PRIMITIVE_HEADER (2);
  1656.   CHECK_ARG (1, ARRAY_P);
  1657.   CHECK_ARG (2, ARRAY_P);
  1658.   Length = ARRAY_LENGTH (ARG_REF (1));
  1659.   if (Length != ARRAY_LENGTH (ARG_REF (2))) error_bad_range_arg (2);
  1660.   From_Here_1 = ARRAY_CONTENTS (ARG_REF (1));
  1661.   From_Here_2 = ARRAY_CONTENTS (ARG_REF (2));
  1662.   To_Here = From_Here_2;
  1663.   for (i=0; i < Length; i++)
  1664.     {
  1665.       *To_Here++ = (*From_Here_1) * (*From_Here_2);
  1666.       From_Here_1++ ;
  1667.       From_Here_2++ ;
  1668.     }
  1669.   PRIMITIVE_RETURN (UNSPECIFIC);
  1670. }
  1671.  
  1672. /* complex-array-operation-2!
  1673.    groups together procedures that use 2 complex-arrays
  1674.    and store result in either 1st or 2nd */
  1675.  
  1676. DEFINE_PRIMITIVE ("COMPLEX-ARRAY-OPERATION-2!", Prim_complex_array_operation_2, 5,5, 0)
  1677. {
  1678.   long n, opcode;
  1679.   REAL *ax,*ay, *bx,*by;
  1680.   void complex_array_multiply_into_second_one();
  1681.   PRIMITIVE_HEADER (5);
  1682.   CHECK_ARG (1, FIXNUM_P);    /* operation opcode */
  1683.   CHECK_ARG (2, ARRAY_P);    /* ax array -- n      real         */
  1684.   CHECK_ARG (3, ARRAY_P);    /* ay array -- n      imag         */
  1685.   CHECK_ARG (4, ARRAY_P);    /* bx array -- n      real         */
  1686.   CHECK_ARG (5, ARRAY_P);    /* by array -- n      imag         */
  1687.   n = ARRAY_LENGTH(ARG_REF(2));
  1688.   if (n != ARRAY_LENGTH(ARG_REF(3))) error_bad_range_arg(3);
  1689.   if (n != ARRAY_LENGTH(ARG_REF(4))) error_bad_range_arg(4);
  1690.   if (n != ARRAY_LENGTH(ARG_REF(4))) error_bad_range_arg(5);
  1691.   ax  = ARRAY_CONTENTS(ARG_REF(2)); /*  real */
  1692.   ay  = ARRAY_CONTENTS(ARG_REF(3)); /*  imag */
  1693.   bx  = ARRAY_CONTENTS(ARG_REF(4)); /*  real */
  1694.   by  = ARRAY_CONTENTS(ARG_REF(5)); /*  imag */
  1695.   opcode = arg_nonnegative_integer(1);
  1696.   if (opcode==1)
  1697.     complex_array_multiply_into_second_one(ax,ay,bx,by, n);
  1698.   else if (opcode==2)
  1699.     error_bad_range_arg(1);    /* illegal opcode */
  1700.   else
  1701.     error_bad_range_arg(1);    /* illegal opcode */
  1702.   PRIMITIVE_RETURN (UNSPECIFIC);
  1703. }
  1704.  
  1705. void
  1706. complex_array_multiply_into_second_one (ax,ay,bx,by, n)
  1707.      REAL *ax,*ay,*bx,*by;
  1708.      long n;
  1709. {
  1710.   long i;
  1711.   REAL temp;
  1712.   for (i=0;i<n;i++)
  1713.     {
  1714.       temp   = ax[i]*bx[i]  -  ay[i]*by[i]; /*  real part */
  1715.       by[i]  = ax[i]*by[i]  +  ay[i]*bx[i]; /*  imag part */
  1716.       bx[i]  = temp;
  1717.     }
  1718. }
  1719.  
  1720. void
  1721. C_Array_Complex_Multiply_Into_First_One (a,b,c,d, length) /* used in fft.c */
  1722.      REAL *a,*b,*c,*d;
  1723.      long length;
  1724. {
  1725.   long i;
  1726.   REAL temp;
  1727.   for (i=0;i<length;i++)
  1728.     {
  1729.       temp = a[i]*c[i] - b[i]*d[i];
  1730.       b[i] = a[i]*d[i] + b[i]*c[i];
  1731.       a[i] = temp;
  1732.     }
  1733. }
  1734.  
  1735. DEFINE_PRIMITIVE ("ARRAY-DIVIDE-INTO-XXX!", Prim_array_divide_into_xxx, 4,4, 0)
  1736. {
  1737.   long n, i, one_or_two;
  1738.   REAL *x,*y,*z, inf;
  1739.   void array_divide_into_z();
  1740.   PRIMITIVE_HEADER (4);
  1741.   CHECK_ARG (1, ARRAY_P);
  1742.   CHECK_ARG (2, ARRAY_P);
  1743.   inf = (arg_real (3));
  1744.   /* where to store result of division */
  1745.   one_or_two = (arg_nonnegative_integer (4));
  1746.   x = ARRAY_CONTENTS(ARG_REF(1));
  1747.   y = ARRAY_CONTENTS(ARG_REF(2));
  1748.   n = ARRAY_LENGTH(ARG_REF(1));
  1749.   if (n!=(ARRAY_LENGTH(ARG_REF(2)))) error_bad_range_arg(2);
  1750.   if (one_or_two == 1)
  1751.     array_divide_into_z( x,y, x,  n, inf);
  1752.   else if (one_or_two == 2)
  1753.     array_divide_into_z( x,y, y,  n, inf);
  1754.   else
  1755.     error_bad_range_arg(4);
  1756.   PRIMITIVE_RETURN (UNSPECIFIC);
  1757. }
  1758.  
  1759. void
  1760. array_divide_into_z (x,y, z, n, inf) /* z can either x or y */
  1761.      REAL *x,*y,*z, inf;
  1762.      long n;
  1763. {
  1764.   long i;
  1765.   for (i=0; i<n; i++)
  1766.     {
  1767.       if (y[i] == 0.0)
  1768.     {
  1769.       if (x[i] == 0.0)
  1770.         z[i] = 1.0;
  1771.       else
  1772.         z[i] = inf * x[i];
  1773.     }
  1774.       else
  1775.     z[i] = x[i] / y[i];
  1776.     }
  1777. }
  1778.  
  1779. /* complex-array-operation-2b!
  1780.    groups together procedures that use 2 complex-arrays
  1781.    & 1 additional real number
  1782.    and store result in either 1st or 2nd (e.g. division) */
  1783.  
  1784. DEFINE_PRIMITIVE ("COMPLEX-ARRAY-OPERATION-2B!", Prim_complex_array_operation_2b, 6,6, 0)
  1785. { long n, opcode;
  1786.   REAL *ax,*ay, *bx,*by,  inf;
  1787.   void complex_array_divide_into_z();
  1788.   PRIMITIVE_HEADER (6);
  1789.   CHECK_ARG (1, FIXNUM_P);    /* operation opcode */
  1790.   CHECK_ARG (2, ARRAY_P);    /* ax array -- n      real         */
  1791.   CHECK_ARG (3, ARRAY_P);    /* ay array -- n      imag         */
  1792.   CHECK_ARG (4, ARRAY_P);    /* bx array -- n      real         */
  1793.   CHECK_ARG (5, ARRAY_P);    /* by array -- n      imag         */
  1794.   inf = (arg_real (6));        /* User-Provided Infinity */
  1795.   n = ARRAY_LENGTH(ARG_REF(2));
  1796.   if (n != ARRAY_LENGTH(ARG_REF(3))) error_bad_range_arg(3);
  1797.   if (n != ARRAY_LENGTH(ARG_REF(4))) error_bad_range_arg(4);
  1798.   if (n != ARRAY_LENGTH(ARG_REF(4))) error_bad_range_arg(5);
  1799.   ax  = ARRAY_CONTENTS(ARG_REF(2)); /*  real */
  1800.   ay  = ARRAY_CONTENTS(ARG_REF(3)); /*  imag */
  1801.   bx  = ARRAY_CONTENTS(ARG_REF(4)); /*  real */
  1802.   by  = ARRAY_CONTENTS(ARG_REF(5)); /*  imag */
  1803.   opcode = arg_nonnegative_integer(1);
  1804.   if (opcode==1)
  1805.     complex_array_divide_into_z (ax,ay,bx,by, ax,ay,  n, inf);
  1806.   else if (opcode==2)
  1807.     complex_array_divide_into_z (ax,ay,bx,by, bx,by,  n, inf);
  1808.   else
  1809.     error_bad_range_arg(1);    /* illegal opcode */
  1810.   PRIMITIVE_RETURN (UNSPECIFIC);
  1811. }
  1812.  
  1813. void
  1814. complex_array_divide_into_z (xr,xi, yr,yi, zr,zi, n, inf)
  1815.      REAL *xr,*xi, *yr,*yi, *zr,*zi, inf;
  1816.      long n;
  1817. {
  1818.   long i;
  1819.   fast double temp, radius;
  1820.   for (i=0; i<n; i++)
  1821.     {
  1822.       radius = (double) (yr[i] * yr[i]) + (yi[i] * yi[i]); /* denominator */
  1823.       if (radius == 0.0)
  1824.     {
  1825.       if (xr[i] == 0.0)
  1826.         zr[i] = 1.0;
  1827.       else
  1828.         zr[i] = inf * xr[i];
  1829.       if (xi[i] == 0.0)
  1830.         zi[i] = 1.0;
  1831.       else
  1832.         zi[i] = inf * xi[i];
  1833.     }
  1834.       else
  1835.     {
  1836.       temp =  (double) (xr[i] * yr[i]  +  xi[i] * yi[i]);
  1837.       zi[i] = (REAL) (xi[i] * yr[i]  -  xr[i] * yi[i]) / radius;
  1838.       zr[i] = (REAL) temp                              / radius;
  1839.     }
  1840.     }
  1841. }
  1842.  
  1843. DEFINE_PRIMITIVE ("ARRAY-LINEAR-SUPERPOSITION-INTO-SECOND-ONE!", Prim_array_linear_superposition_into_second_one, 4, 4, 0)
  1844. {
  1845.   long n, i;
  1846.   REAL *To_Here, Coeff1, Coeff2;
  1847.   REAL *From_Here_1, *From_Here_2;
  1848.   PRIMITIVE_HEADER (4);
  1849.   Coeff1 = (arg_real (1));
  1850.   CHECK_ARG (2, ARRAY_P);
  1851.   Coeff2 = (arg_real (3));
  1852.   CHECK_ARG (4, ARRAY_P);
  1853.   n = (ARRAY_LENGTH (ARG_REF (2)));
  1854.   if (n != (ARRAY_LENGTH (ARG_REF (4))))
  1855.     error_bad_range_arg (4);
  1856.   From_Here_1 = (ARRAY_CONTENTS (ARG_REF (2)));
  1857.   From_Here_2 = (ARRAY_CONTENTS (ARG_REF (4)));
  1858.   To_Here = From_Here_2;
  1859.   for (i=0; i < n; i++)
  1860.     {
  1861.       *To_Here++ = (Coeff1 * (*From_Here_1)) + (Coeff2 * (*From_Here_2));
  1862.       From_Here_1++ ;
  1863.       From_Here_2++ ;
  1864.     }
  1865.   PRIMITIVE_RETURN (UNSPECIFIC);
  1866. }
  1867.  
  1868. /*  m_pi = 3.14159265358979323846264338327950288419716939937510; */
  1869.  
  1870. DEFINE_PRIMITIVE ("SAMPLE-PERIODIC-FUNCTION", Prim_sample_periodic_function, 4, 4, 0)
  1871. {
  1872.   long N, i, Function_Number;
  1873.   double Signal_Frequency, Sampling_Frequency, DT, DTi;
  1874.   double twopi = 6.28318530717958;
  1875.   SCHEME_OBJECT Result, Pfunction_number, Psignal_frequency;
  1876.   SCHEME_OBJECT Pfunction_Number;
  1877.   REAL *To_Here;
  1878.   double unit_square_wave(), unit_triangle_wave();
  1879.   PRIMITIVE_HEADER (4);
  1880.   Function_Number = (arg_index_integer (1, 11));
  1881.   Signal_Frequency = (arg_real_number (2));
  1882.   if (Signal_Frequency == 0)
  1883.     error_bad_range_arg (2);
  1884.   Sampling_Frequency = (arg_real_number (3));
  1885.   if (Sampling_Frequency == 0)
  1886.     error_bad_range_arg (3);
  1887.   N = (arg_nonnegative_integer (4));
  1888.   Result = (allocate_array (N));
  1889.   To_Here = ARRAY_CONTENTS(Result);
  1890.   DT = (double) (twopi * Signal_Frequency * (1 / Sampling_Frequency));
  1891.   if (Function_Number == 0)
  1892.     for (i=0, DTi=0.0; i < N; i++, DTi += DT)
  1893.       *To_Here++ = (REAL) cos(DTi);
  1894.   else if (Function_Number == 1)
  1895.     for (i=0, DTi=0.0; i < N; i++, DTi += DT)
  1896.       *To_Here++ = (REAL) sin(DTi);
  1897.   else if (Function_Number == 2)
  1898.     for (i=0, DTi=0.0; i < N; i++, DTi += DT)
  1899.       *To_Here++ = (REAL) unit_square_wave(DTi);
  1900.   else if (Function_Number == 3)
  1901.     for (i=0, DTi=0.0; i < N; i++, DTi += DT)
  1902.       *To_Here++ = (REAL) unit_triangle_wave(DTi);
  1903.   else
  1904.     error_bad_range_arg (1);
  1905.   PRIMITIVE_RETURN (Result);
  1906. }
  1907.  
  1908. double
  1909. hamming (t, length)
  1910.      double t, length;
  1911. {
  1912.   double twopi = 6.28318530717958;
  1913.   double pi = twopi/2.;
  1914.   double t_bar = cos(twopi * (t / length));
  1915.   if ((t<length) && (t>0.0)) return(.08 + .46 * (1 - t_bar));
  1916.   else return (0);
  1917. }
  1918.  
  1919. double
  1920. unit_square_wave (t)
  1921.      double t;
  1922. {
  1923.   double twopi = 6.28318530717958;
  1924.   double fmod(), fabs();
  1925.   double pi = twopi/2.;
  1926.   double t_bar = ((REAL) fabs(fmod( ((double) t), twopi)));
  1927.   if (t_bar < pi) return(1);
  1928.   else return(-1);
  1929. }
  1930.  
  1931. double
  1932. unit_triangle_wave (t)
  1933.      double t;
  1934. {
  1935.   double twopi = 6.28318530717958;
  1936.   double pi = twopi/2.;
  1937.   double pi_half = pi/2.;
  1938.   double three_pi_half = pi+pi_half;
  1939.   double t_bar = ((double) fabs(fmod( ((double) t), twopi)));
  1940.   if (t_bar<pi_half)
  1941.     return (-(t_bar/pi));
  1942.   else if (t_bar<pi)
  1943.     return (t_bar/pi);
  1944.   else if (t_bar<three_pi_half)
  1945.     return ((twopi-t_bar)/pi);
  1946.   else
  1947.     return (-((twopi-t_bar)/pi));
  1948. }
  1949.  
  1950. DEFINE_PRIMITIVE ("ARRAY-HANNING!", Prim_array_hanning, 2,2, 0)
  1951. {
  1952.   long n, hanning_power;
  1953.   REAL *a;
  1954.   void C_Array_Make_Hanning();
  1955.   PRIMITIVE_HEADER (2);
  1956.   CHECK_ARG (1, ARRAY_P);    /* input array -- n */
  1957.   CHECK_ARG (2, FIXNUM_P);    /* hanning power */
  1958.   a  = ARRAY_CONTENTS(ARG_REF(1));
  1959.   n = ARRAY_LENGTH(ARG_REF(1));
  1960.   hanning_power = arg_nonnegative_integer(2);
  1961.   C_Array_Make_Hanning (a, n, hanning_power);
  1962.   PRIMITIVE_RETURN (UNSPECIFIC);
  1963. }
  1964.  
  1965. void
  1966. C_Array_Make_Hanning (f1, length, power)
  1967.      REAL f1[];
  1968.      long length, power;
  1969. {
  1970.   double window_length;
  1971.   long i;
  1972.   double integer_power(), hanning();
  1973.   window_length = ((double) length);
  1974.   for (i=0;i<length;i++)
  1975.     {
  1976.       f1[i] = ((REAL) hanning(((double) i), window_length));
  1977.       f1[i] = (REAL) integer_power(((double) f1[i]), power);
  1978.     }
  1979. }
  1980.  
  1981. double
  1982. hanning (t, length)
  1983.      double t, length;
  1984. {
  1985.   double twopi = 6.283185307179586476925287;
  1986.   double t_bar;
  1987.   t_bar = cos(twopi * (t / length));
  1988.   if ((t<length) && (t>0.0))
  1989.     return(.5 * (1 - t_bar));
  1990.   else
  1991.     return (0.0);
  1992. }
  1993.  
  1994. double
  1995. integer_power (a, n)
  1996.      double a;
  1997.      long n;
  1998. {
  1999.   double b;
  2000.   double integer_power();
  2001.  
  2002.   if (n<0) exit(-1);
  2003.   else if (n==0) return(1.0);
  2004.   else if (n==1) return(a);
  2005.   else if ((n%2) == 0)
  2006.     {
  2007.       b = integer_power(a, n/2);
  2008.       return(b*b);
  2009.     }
  2010.   else
  2011.     return(a * integer_power(a, (n-1)));
  2012. }
  2013.  
  2014. /* array-operation-1!
  2015.    groups together procedures that use 1 array
  2016.    and store the result in place (e.g. random) */
  2017.  
  2018. DEFINE_PRIMITIVE ("ARRAY-OPERATION-1!", Prim_array_operation_1, 2,2, 0)
  2019. {
  2020.   long n, opcode;
  2021.   REAL *a;
  2022.   void array_random();
  2023.   PRIMITIVE_HEADER (2);
  2024.   CHECK_ARG (1, FIXNUM_P);    /* operation opcode */
  2025.   CHECK_ARG (2, ARRAY_P);    /* input array -- n */
  2026.   n = ARRAY_LENGTH(ARG_REF(2));
  2027.   a  = ARRAY_CONTENTS(ARG_REF(2));
  2028.   opcode = arg_nonnegative_integer(1);
  2029.   if (opcode==1)
  2030.     array_random(a,n);
  2031.   else if (opcode==2)
  2032.     error_bad_range_arg(1);    /* illegal opcode */
  2033.   else
  2034.     error_bad_range_arg(1);    /* illegal opcode */
  2035.   PRIMITIVE_RETURN (UNSPECIFIC);
  2036. }
  2037.  
  2038. void
  2039. array_random (a,n)
  2040.      REAL *a;
  2041.      long n;
  2042. {
  2043.   long i;
  2044.   /* HPUX 3: Rand uses a multiplicative congruential random-number generator
  2045.      with period 2^32 that returns successive pseudo-random numbers in the
  2046.      range from 0 to 2^15-1 */
  2047.   for (i=0;i<n;i++)
  2048.     a[i] = ((REAL) rand()) * (3.0517578125e-5);
  2049.   /* 3.051xxx = 2^(-15) makes the range from 0 to 1 */
  2050. }
  2051.  
  2052. /* The following should go away.
  2053.    superceded by ARRAY-CONS-INTEGERS, ARRAY-UNARY-FUNCTION and array-random */
  2054. DEFINE_PRIMITIVE ("SAMPLE-APERIODIC-FUNCTION", Prim_sample_aperiodic_function, 3, 3, 0)
  2055. {
  2056.   long N, i, Function_Number;
  2057.   double Sampling_Frequency, DT, DTi;
  2058.   double twopi = 6.28318530717958;
  2059.   SCHEME_OBJECT Result;
  2060.   REAL *To_Here, twopi_dt;
  2061.   PRIMITIVE_HEADER (3);
  2062.   Function_Number = (arg_index_integer (1, 7));
  2063.   Sampling_Frequency = (arg_real_number (2));
  2064.   if (Sampling_Frequency == 0)
  2065.     error_bad_range_arg (2);
  2066.   N = (arg_nonnegative_integer (3));
  2067.   Result = (allocate_array (N));
  2068.   To_Here = ARRAY_CONTENTS(Result);
  2069.   DT = (twopi * (1 / Sampling_Frequency));
  2070.   if      (Function_Number == 0)
  2071.     /* HPUX 3: Rand uses a multiplicative congruential random-number generator
  2072.        with period 2^32 that returns successive pseudo-random numbers in the
  2073.        range from 0 to 2^15-1 */
  2074.     for (i=0; i<N; i++)
  2075.       /* 2^(-15) makes range from 0 to 1 */
  2076.       *To_Here++ = 3.0517578125e-5 * ((REAL) rand());
  2077.   else if (Function_Number == 1)
  2078.   { double length=DT*N;
  2079.     for (i=0, DTi=0.0; i < N; i++, DTi += DT)
  2080.       *To_Here++ = (REAL) hanning(DTi, length);
  2081.   }
  2082.   else if (Function_Number == 2)
  2083.   { double length=DT*N;
  2084.     for (i=0, DTi=0.0; i < N; i++, DTi += DT)
  2085.       *To_Here++ = (REAL) hamming(DTi, length);
  2086.   }
  2087.   else if (Function_Number == 3)
  2088.     for (i=0, DTi=0.0; i < N; i++, DTi += DT)
  2089.       *To_Here++ = (REAL) sqrt(DTi);
  2090.   else if (Function_Number == 4)
  2091.     for (i=0, DTi=0.0; i < N; i++, DTi += DT)
  2092.       *To_Here++ = (REAL) log(DTi);
  2093.   else if (Function_Number == 5)
  2094.     for (i=0, DTi=0.0; i < N; i++, DTi += DT)
  2095.       *To_Here++ = (REAL) exp(DTi);
  2096.   else
  2097.     error_bad_range_arg (1);
  2098.   PRIMITIVE_RETURN (Result);
  2099. }
  2100.  
  2101. DEFINE_PRIMITIVE ("ARRAY-PERIODIC-DOWNSAMPLE", Prim_array_periodic_downsample, 2, 2, 0)
  2102. {
  2103.   long Length, Pseudo_Length, Sampling_Ratio;
  2104.   REAL *Array, *To_Here;
  2105.   SCHEME_OBJECT Result;
  2106.   long i, array_index;
  2107.   PRIMITIVE_HEADER (2);
  2108.   CHECK_ARG (1, ARRAY_P);
  2109.   Length = ARRAY_LENGTH(ARG_REF (1));
  2110.   Sampling_Ratio = ((arg_integer (2)) % Length);
  2111.   if (Sampling_Ratio < 1)
  2112.     error_bad_range_arg (2);
  2113.   Array = ARRAY_CONTENTS(ARG_REF (1));
  2114.   Result = (allocate_array (Length));
  2115.   To_Here = ARRAY_CONTENTS(Result);
  2116.   Pseudo_Length = Length * Sampling_Ratio;
  2117.   /* new Array has the same Length by assuming periodicity */
  2118.   for (i=0; i<Pseudo_Length; i += Sampling_Ratio)
  2119.   { array_index = i % Length;
  2120.     *To_Here++ = Array[array_index]; }
  2121.   PRIMITIVE_RETURN (Result);
  2122. }
  2123.  
  2124. /* Shift is not done in place (no side-effects). */
  2125. DEFINE_PRIMITIVE ("ARRAY-PERIODIC-SHIFT", Prim_array_periodic_shift, 2, 2, 0)
  2126. {
  2127.   long Length, Shift;
  2128.   REAL *Array, *To_Here;
  2129.   SCHEME_OBJECT Result;
  2130.   long i, array_index;
  2131.   PRIMITIVE_HEADER (2);
  2132.   CHECK_ARG (1, ARRAY_P);
  2133.   Length = ARRAY_LENGTH(ARG_REF (1));
  2134.   /* periodic waveform, same sign as dividend */
  2135.   Shift = ((arg_integer (2)) % Length);
  2136.   Array = ARRAY_CONTENTS(ARG_REF (1));
  2137.   Result = (allocate_array (Length));
  2138.   To_Here = ARRAY_CONTENTS(Result);
  2139.   /* new Array has the same Length by assuming periodicity */
  2140.   for (i=0; i<Length; i++) {
  2141.     array_index = (i+Shift) % Length;
  2142.     if (array_index<0) array_index = Length + array_index; /* wrap around */
  2143.     *To_Here++ = Array[array_index]; }
  2144.   PRIMITIVE_RETURN (Result);
  2145. }
  2146.  
  2147. /* This is done here because array-map is very slow */
  2148. DEFINE_PRIMITIVE ("ARRAY-APERIODIC-DOWNSAMPLE", Prim_array_aperiodic_downsample, 2, 2, 0)
  2149. {
  2150.   long Length, New_Length, Sampling_Ratio;
  2151.   REAL *Array, *To_Here;
  2152.   SCHEME_OBJECT Result;
  2153.   long i, array_index;
  2154.   PRIMITIVE_HEADER (2);
  2155.   CHECK_ARG (1, ARRAY_P);
  2156.   Array = ARRAY_CONTENTS(ARG_REF (1));
  2157.   Length = ARRAY_LENGTH(ARG_REF (1));
  2158.   Sampling_Ratio = (arg_integer_in_range (2, 1, (Length + 1)));
  2159.   if (Length < 1) error_bad_range_arg (1);
  2160.   /* 1 for first one and then the rest --- integer division chops */
  2161.   New_Length = 1 + ((Length-1)/Sampling_Ratio);
  2162.   Result = (allocate_array (New_Length));
  2163.   To_Here = ARRAY_CONTENTS(Result);
  2164.   for (i=0; i<Length; i += Sampling_Ratio)
  2165.     *To_Here++ = Array[i];
  2166.   PRIMITIVE_RETURN (Result);
  2167. }
  2168.  
  2169. /* one more hack for speed */
  2170.  
  2171. /* (SOLVE-SYSTEM A B N)
  2172.     Solves the system of equations Ax = b.  A and B are
  2173.     arrays and b is the order of the system.  Returns x.
  2174.     From the Fortran procedure in Strang. */
  2175.  
  2176. DEFINE_PRIMITIVE ("SOLVE-SYSTEM", Prim_gaussian_elimination, 2, 2, 0)
  2177. {
  2178.   REAL *A, *B, *X;
  2179.   long Length;
  2180.   SCHEME_OBJECT Result;
  2181.   PRIMITIVE_HEADER (2);
  2182.   CHECK_ARG (1, ARRAY_P);
  2183.   CHECK_ARG (2, ARRAY_P);
  2184.   Length  = ARRAY_LENGTH(ARG_REF (2));
  2185.   if ((Length*Length) != ARRAY_LENGTH(ARG_REF (1))) error_bad_range_arg (2);
  2186.   A = ARRAY_CONTENTS(ARG_REF (1));
  2187.   B = ARRAY_CONTENTS(ARG_REF (2));
  2188.   Result = (allocate_array (Length));
  2189.   X = ARRAY_CONTENTS(Result);
  2190.   C_Array_Copy(B, X, Length);
  2191.   C_Gaussian_Elimination(A, X, Length);
  2192.   PRIMITIVE_RETURN (Result);
  2193. }
  2194.  
  2195. /* C routine side-effects b. */
  2196. C_Gaussian_Elimination (a, b, n)
  2197.      REAL *a, *b;
  2198.      long n;
  2199. {
  2200.   long *pvt;
  2201.   REAL p, t;
  2202.   long i, j, k, m;
  2203.   Primitive_GC_If_Needed (n);
  2204.   pvt = ((long *) Free);
  2205.   *(pvt+n-1) = 1;
  2206.   if (n != 1) {
  2207.     for (k=1; k<n; k++) {
  2208.       m = k;
  2209.       for (i=k+1; i<=n; i++)
  2210.     if (fabs(*(a+i+(k-1)*n-1)) > fabs(*(a+m+(k-1)*n-1)))
  2211.       m = i;
  2212.       *(pvt+k-1) = m;
  2213.       if (m != k)
  2214.     *(pvt+n-1) = - *(pvt+n-1);
  2215.       p = *(a+m+(k-1)*n-1);
  2216.       *(a+m+(k-1)*n-1) = *(a+k+(k-1)*n-1);
  2217.       *(a+k+(k-1)*n-1) = p;
  2218.       if (p != 0.0) {
  2219.     for (i=k+1; i<=n; i++)
  2220.       *(a+i+(k-1)*n-1) = - *(a+i+(k-1)*n-1) / p;
  2221.     for (j=k+1; j<=n; j++) {
  2222.       t = *(a+m+(j-1)*n-1);
  2223.       *(a+m+(j-1)*n-1) = *(a+k+(j-1)*n-1);
  2224.       *(a+k+(j-1)*n-1) = t;
  2225.       if (t != 0.0)
  2226.         for (i=k+1; i<=n; i++)
  2227.           *(a+i+(j-1)*n-1) = *(a+i+(j-1)*n-1) + *(a+i+(k-1)*n-1) * t;
  2228.     }
  2229.       }
  2230.     }
  2231.     for (k=1; k<n; k++) {
  2232.       m = *(pvt+k-1);
  2233.       t = *(b+m-1);
  2234.       *(b+m-1) = *(b+k-1);
  2235.       *(b+k-1) = t;
  2236.       for (i=k+1; i<=n; i++)
  2237.     *(b+i-1) = *(b+i-1) + *(a+i+(k-1)*n-1) * t;
  2238.     }
  2239.     for (j=1; j<n; j++) {
  2240.       k = n - j + 1;
  2241.       *(b+k-1) = *(b+k-1) / *(a+k+(k-1)*n-1);
  2242.       t = - *(b+k-1);
  2243.       for (i=1; i <= n-j; i++)
  2244.     *(b+i-1) = *(b+i-1) + *(a+i+(k-1)*n-1) * t;
  2245.     }
  2246.   }
  2247.   *b = *b / *a;
  2248.   return;
  2249. }
  2250.