/*  pack								*/
/* & pjt								*/
/* : low-level-i/o							*/
/* +									*/
/*									*/
/*		This converts data between disk and internal		*/
/*		format. Disk format is IEEE reals and 16 or 32		*/
/*		bit integers (most significant byte first).		*/
/*									*/
/*		This assumes that these are the local machine format	*/
/*		(float == IEEE real, int == 32 bit integer,		*/
/*		short int == 16 bit integer).				*/
/*									*/
/*		packx_c, unpackx_c, pack32_c and unpack32_c are		*/
/*		implemented as macros (calling bcopy) in the		*/
/*		system dependent include file.				*/
/*--									*/
/*									*/
/*  History:								*/
/*    rjs  Dark-ages Original version.					*/
/*    bs   ?????89   Improved efficiency using "register" declarations.	*/
/*    rjs   1nov89   Incoporated Brian's changes.			*/
/*    mjs  28feb91   Merge Sun and Cray versions.			*/
/*    mjs  18mar91   Added convex definition.                           */
/*    mjs  19feb93   Added mips definition.                             */
/*    pjt  25jan95   linux kludge to include packALPHA.c                */
/*    pjt  14jun01   packALPHA.c now included in this source code       */
/*                   and using the standard WORDS_BIGENDIAN macro       */
/*    pjt  21jun02   MIR4 prototyping                                   */
/************************************************************************/

#include <Python.h>
#include "sysdep.h"
#include "miriad.h"

#if defined(WORDS_BIGENDIAN)

/* static int words_bigendian = 1; never used actually, but handy symbol to find via nm(1) */

void	pack16_c(register int *from,char *to,int n)
{
	register short int	*tto;
	register int		i;

	tto = (short int *)to;
	for (i=0; i < n; i++)	*tto++ = *from++;
}
void	unpack16_c(char *from,register int *to,int n)
{
	register short int	*ffrom;
	register int		i;

	ffrom = (short int *)from;
	for (i=0; i < n; i++)	*to++ = *ffrom++;
}

void	pack64_c(register int *from,char *to,int n)
{
	register short int	*tto;
	register int		i;

	tto = (short int *)to;
	for (i=0; i < n; i++)	*tto++ = *from++;
}
void	unpack64_c(char *from,register int *to,int n)
{
	register short int	*ffrom;
	register int		i;

	ffrom = (short int *)from;
	for (i=0; i < n; i++)	*to++ = *ffrom++;
}

#endif


#ifndef WORDS_BIGENDIAN 
#ifndef unicos
/* static int words_littleendian = 1; never used actually, but handy symbol to find via nm(1) */
/************************************************************************/
/*									*/
/*  The pack routines -- these convert between the host format and	*/
/*  the disk format. Disk format is IEEE 32 and 64 bit reals, and 2's	*/
/*  complement integers. Byte order is the FITS byte order (most	*/
/*  significant bytes first).						*/
/*									*/
/*  This version is for a machine which uses IEEE internally, but which	*/
/*  uses least significant bytes first (little endian), e.g. PCs and	*/
/*  Alphas.								*/
/*									*/
/*  History:								*/
/*    rjs  21nov94 Original version.					*/
/************************************************************************/
void pack16_c(int *in,char *out,int n)
/*
  Pack an integer array into 16 bit integers.
------------------------------------------------------------------------*/
{
  int i;
  char *s;

  s = (char *)in;
  for(i=0; i < n; i++){
    *out++ = *(s+1);
    *out++ = *s;
    s += sizeof(int);
  }
}
/************************************************************************/
void unpack16_c(char *in,int *out,int n)
/*
  Unpack an array of 16 bit integers into integers.
------------------------------------------------------------------------*/
{
  int i;
  char *s;

  s = (char *)out;
  for(i=0; i < n; i++){
    *s++ = *(in+1);
    *s++ = *in;
    if(0x80 & *in){
      *s++ = 0xFF;
      *s++ = 0xFF;
    } else {
      *s++ = 0;
      *s++ = 0;
    }
    in += 2;
  }
}
/************************************************************************/
void pack32_c(int *in,char *out,int n)
/*
  Pack an array of integers into 32 bit integers.
------------------------------------------------------------------------*/
{
  int i;
  char *s;

  s = (char *)in;
  for(i = 0; i < n; i++){
    *out++ = *(s+3);
    *out++ = *(s+2);
    *out++ = *(s+1);
    *out++ = *s;
    s += 4;
  }
}
/************************************************************************/
void unpack32_c(char *in,int *out,int n)
/*
  Unpack an array of 32 bit integers into integers.
------------------------------------------------------------------------*/
{
  int i;
  char *s;

  s = (char *)out;
  for(i = 0; i < n; i++){
    *s++ = *(in+3);
    *s++ = *(in+2);
    *s++ = *(in+1);
    *s++ = *in;
    in += 4;
  }
}
/************************************************************************/
void pack64_c(int8 *in,char *out,int n)
/*
  Pack an integer array into 64 bit integers.
------------------------------------------------------------------------*/
{
  int i;
  char *s;

  s = (char *)in;
  for(i=0; i < n; i++){
    *out++ = *(s+7);
    *out++ = *(s+6);
    *out++ = *(s+5);
    *out++ = *(s+4);
    *out++ = *(s+3);
    *out++ = *(s+2);
    *out++ = *(s+1);
    *out++ = *s;
    s += 8;
  }
}
/************************************************************************/
void unpack64_c(char *in,int8 *out,int n)
/*
  Unpack an array of 64 bit integers into integers.
------------------------------------------------------------------------*/
{
  int i;
  char *s;

  s = (char *)out;
  for(i=0; i < n; i++){
    *s++ = *(in+7);
    *s++ = *(in+6);
    *s++ = *(in+5);
    *s++ = *(in+4);
    *s++ = *(in+3);
    *s++ = *(in+2);
    *s++ = *(in+1);
    *s++ = *in;
    in += 8;
  }
}
/************************************************************************/
void packr_c(float *in,char *out,int n)
/*
  Pack an array of reals into IEEE reals -- just do byte reversal.
------------------------------------------------------------------------*/
{
  int i;
  char *s;

  s = (char *)in;
  for(i = 0; i < n; i++){
    *out++ = *(s+3);
    *out++ = *(s+2);
    *out++ = *(s+1);
    *out++ = *s;
    s += 4;
  }
}
/************************************************************************/
void unpackr_c(char *in,float *out,int n)
/*
  Unpack an array of IEEE reals into reals -- just do byte reversal.
------------------------------------------------------------------------*/
{
  int i;
  char *s;

  s = (char *)out;
  for(i = 0; i < n; i++){
    *s++ = *(in+3);
    *s++ = *(in+2);
    *s++ = *(in+1);
    *s++ = *in;
    in += 4;
  }
}
/************************************************************************/
void packd_c(double *in,char *out,int n)
/*
  Pack an array of doubles -- this involves simply performing byte
  reversal.
------------------------------------------------------------------------*/
{
  int i;
  char *s;

  s = (char *)in;
  for(i = 0; i < n; i++){
    *out++ = *(s+7);
    *out++ = *(s+6);
    *out++ = *(s+5);
    *out++ = *(s+4);
    *out++ = *(s+3);
    *out++ = *(s+2);
    *out++ = *(s+1);
    *out++ = *s;
    s += 8;
  }
}
/************************************************************************/
void unpackd_c(char *in,double *out,int n)
/*
  Unpack an array of doubles -- this involves simply performing byte
  reversal.
------------------------------------------------------------------------*/
{
  int i;
  char *s;

  s = (char *)out;
  for(i = 0; i < n; i++){
    *s++ = *(in+7);
    *s++ = *(in+6);
    *s++ = *(in+5);
    *s++ = *(in+4);
    *s++ = *(in+3);
    *s++ = *(in+2);
    *s++ = *(in+1);
    *s++ = *in;
    in += 8;
  }
}
#endif
#endif


#if defined(unicos)
static int words_unicos = 1;
#define TWO15  0x8000
#define TWO16  0x10000
#define TWO31  0x80000000
#define TWO32  0x100000000
#define HILONG 0xFFFFFFFF00000000
#define LOLONG 0x00000000FFFFFFFF
#define WORD0  0x000000000000FFFF
#define WORD1  0x00000000FFFF0000
#define WORD2  0x0000FFFF00000000
#define WORD3  0xFFFF000000000000

/* Masks for IEEE floating format (both hi and lo types). */

#define IEEE_HISIGN     0x8000000000000000
#define IEEE_HIEXPO     0x7F80000000000000
#define IEEE_HIMANT     0x007FFFFF00000000
#define IEEE_LOSIGN     0x0000000080000000
#define IEEE_LOEXPO     0x000000007F800000
#define IEEE_LOMANT     0x00000000007FFFFF
#define IEEE_DMANT      0x000FFFFFFFFFFFF0
#define IEEE_DEXPO      0x7FF0000000000000

/* Masks for Cray floating format. */

#define CRAY_MANT       0x0000FFFFFF000000      /* Including unhidden bit. */
#define CRAY_MANT1      0x00007FFFFF000000      /* No unhidden bit. */
#define CRAY_DMANT      0x0000FFFFFFFFFFFF
#define CRAY_DMANT1     0x00007FFFFFFFFFFF
#define CRAY_EXPO       0x7FFF000000000000
#define SIGN            0x8000000000000000

/* Mask of a pointer to char giving the character offset in a Cray word. */

#define CHAR_OFFSET 0xE000000000000000

/************************************************************************/
void pack16_c(int *in,char *out,int n)
/*
  Pack an integer array into 16 bit integers for unicos
------------------------------------------------------------------------*/
{
  int temp,offset,*outd,in1,in2,in3,in4,i;

  if(n <= 0)return;                             /* Return if nothing to do. */
  temp = (int)out;
  offset = ( temp & CHAR_OFFSET ) >> 62;        /* Get offset of first word. */
  outd = (int *)(temp & ~CHAR_OFFSET);  /* Get address of words. */

/* Handle the first few which are not aligned on a Cray word. */

  switch(offset){
    case 1: *outd = (*outd & ~WORD2) | ((*in++ << 32) & WORD2);
            if(--n == 0)break;
    case 2: *outd = (*outd & ~WORD1) | ((*in++ << 16) & WORD1);
            if(--n == 0)break;
    case 3: *outd = (*outd & ~WORD0) | ((*in++    ) & WORD0);
            outd++;
  }

/* Handle the ones which are aligned on a Cray word. */

  for(i=0; i < n-3; i=i+4){
    in1 = *in++ << 48;
    in2 = *in++ << 32;
    in3 = *in++ << 16;
    in4 = *in++;
    *outd++ = (in1 & WORD3) | (in2 & WORD2) | (in3 & WORD1) | (in4 & WORD0);
  }
  n -= i;

/* Handle the last ones which are not aligned on a Cray word. */

  if(n-- > 0){
    *outd = (*outd & ~WORD3) | ((*in++ << 48) & WORD3);
    if(n-- > 0){
      *outd = (*outd & ~WORD2) | ((*in++ << 32) & WORD2);
      if(n-- > 0){
        *outd = (*outd & ~WORD1) | ((*in++ << 16) & WORD1);
      }
    }
  }
}
/************************************************************************/
void unpack16_c(char *in,int *out,int n)
/*
  Unpack an array of 16 bit integers into integers for unicos
------------------------------------------------------------------------*/
{
  int temp,offset,*ind,i;

  if(n <= 0)return;                             /* Return if nothing to do. */
  temp = (int)in;
  offset = ( temp & CHAR_OFFSET ) >> 62;        /* Get offset of first word. */
  ind = (int *)(temp & ~CHAR_OFFSET);   /* Get address of words. */

/* Handle the first few which are not word aligned. */

  switch(offset){
    case 1:  temp = (*ind >> 32) & WORD0;
             *out++ = (temp < TWO15 ? temp : temp - TWO16);
             if(--n == 0) break;
    case 2:  temp = (*ind >> 16) & WORD0;
             *out++ = (temp < TWO15 ? temp : temp - TWO16);
             if(--n == 0) break;
    case 3:  temp = (*ind++    ) & WORD0;
             *out++ = (temp < TWO15 ? temp : temp - TWO16);
             if(--n == 0) break;
  }

/* Handle those that are Cray-word-aligned. */

  for(i=0; i < n-3; i=i+4){
    temp = (*ind >> 48) & WORD0;
    *out++ = (temp < TWO15 ? temp : temp - TWO16);
    temp = (*ind >> 32) & WORD0;
    *out++ = (temp < TWO15 ? temp : temp - TWO16);
    temp = (*ind >> 16) & WORD0;
    *out++ = (temp < TWO15 ? temp : temp - TWO16);
    temp = (*ind++    ) & WORD0;
    *out++ = (temp < TWO15 ? temp : temp - TWO16);
  }
  n -= i;

/* Handle the last few which are not Cray-word-aligned. */

  if(n-- > 0){
    temp = (*ind >> 48) & WORD0;
    *out++ = (temp < TWO15 ? temp : temp - TWO16);
    if(n-- > 0){
      temp = (*ind >> 32) & WORD0;
      *out++ = (temp < TWO15 ? temp : temp - TWO16);
      if(n-- > 0){
        temp = (*ind >> 16) & WORD0;
        *out++ = (temp < TWO15 ? temp : temp - TWO16);
      }
    }
  }
}
/************************************************************************/
void pack32_c(int *in,char *out,int n)
/*
  Pack an array of integers into 32 bit integers for unicos
------------------------------------------------------------------------*/
{
  int temp,offset,*outd,i,in1,in2;

  if(n <= 0)return;                             /* Return if nothing to do. */
  temp = (int)out;
  offset = ( temp & CHAR_OFFSET ) >> 63;        /* Get offset of first long. */
  outd = (int *)(temp & ~CHAR_OFFSET);  /* Get address of words. */

/* Do the first one, if it is not aligned on a Cray word. */

  if(offset==1){
    *outd = (*outd & ~LOLONG) | (*in++ & LOLONG);
    outd++;
  }
  n -= offset;

/* Do those which are Cray word aligned. */

  for(i=0; i < n-1; i=i+2){
    in1 = *in++ << 32;
    in2 = *in++;
    *outd++ = (in1 & HILONG) | (in2 & LOLONG);
  }
  n -= i;

/* Handle the last one, if there is one. */

  if(n==1)*outd =  (*outd & ~HILONG) | ((*in++ << 32) & HILONG);
}
/************************************************************************/
void unpack32_c(char *in,int *out,int n)
/*
  Unpack an array of 32 bit integers into integers for unicos
------------------------------------------------------------------------*/
{
  int temp,offset,*ind,i;

  if(n <= 0)return;                             /* Return if nothing to do. */
  temp = (int)in;
  offset = ( temp & CHAR_OFFSET ) >> 63;        /* Get offset of first word. */
  ind = (int *)(temp & ~CHAR_OFFSET);   /* Get address of words. */

/* Handle one which is not Cray word aligned. */

  if(offset==1){
    temp = (*ind++ & LOLONG);
    *out++ = (temp < TWO31 ? temp : temp - TWO32);
  }
  n -= offset;

/* Handle those which are Cray word aligned. */

  for(i=0; i < n-1; i=i+2){
    temp = (*ind >> 32) & LOLONG;
    *out++ = (temp < TWO31 ? temp : temp - TWO32);
    temp = (*ind++    ) & LOLONG;
    *out++ = (temp < TWO31 ? temp : temp - TWO32);
  }
  n -= i;

/* Possibly handle a last one which is not Cray word aligned. */

  if(n==1){
    temp = (*ind >> 32) & LOLONG;
    *out++ = (temp < TWO31 ? temp : temp - TWO32);
  }
}
/************************************************************************/
void packr_c(float *in,char *out,int n)
/*
  Pack an array of Cray reals into IEEE reals.
------------------------------------------------------------------------*/
{
  int temp,offset,*outd,bias,*ind,tin,tout,i;

  if(n <= 0)return;                             /* Return if nothing to do. */
  temp = (int)out;
  offset = ( temp & CHAR_OFFSET ) >> 63;        /* Get offset of first long. */
  outd = (int *)(temp & ~CHAR_OFFSET);  /* Get address of words. */
  bias = (16384 - 126) << 48;
  ind = (int *)in;

/* Do the first one, if it is not aligned on a Cray word. */

  if(offset==1){
    tin     = *ind++;
    *outd    =   (*outd & ~LOLONG) |
                 (tin & CRAY_MANT ? (((tin & CRAY_EXPO)-bias) >> 25) |
                ((tin & CRAY_MANT1) >> 24) | ((tin & SIGN) >> 32) : 0);
    outd++;
  }
  n -= offset;

/* Do those which are Cray word aligned. */

  for(i=0; i < n-1; i=i+2){
    tin = *ind++;
    tout =      (tin & CRAY_MANT ? (((tin & CRAY_EXPO)-bias) << 7) |
                ((tin & CRAY_MANT1) << 8) |   (tin & SIGN)        : 0);
    tin = *ind++;
    *outd++ =   tout |
                (tin & CRAY_MANT ? (((tin & CRAY_EXPO)-bias) >> 25) |
                ((tin & CRAY_MANT1) >> 24) | ((tin & SIGN) >> 32) : 0);
  }
  n -= i;

/* Handle the last one, if there is one. */

  if(n==1){
    tin = *ind;
    *outd = (*outd & ~HILONG) |
            (tin & CRAY_MANT ? (((tin & CRAY_EXPO)-bias) << 7) |
            ((tin & CRAY_MANT1) << 8) |   (tin & SIGN)        : 0);
  }
}
/************************************************************************/
void unpackr_c(char *in,float *out,int n)
/*
  Unpack an array of IEEE reals into Cray reals.
------------------------------------------------------------------------*/
{
  int temp,tin,*ind,*outd,offset,i,bias;

  if(n <= 0)return;                             /* Return if nothing to do. */
  temp = (int)in;
  offset = ( temp & CHAR_OFFSET ) >> 63;        /* Get offset of first word. */
  ind = (int *)(temp & ~CHAR_OFFSET);   /* Get address of words. */
  outd = (int *)out;
  bias = ((16384-126) <<48) + (1 << 47);

/* Handle the first one if it is not aligned on a Cray word. */

  if(offset==1){
    tin = *ind++;
    *outd++ = (tin & IEEE_LOEXPO ? (((tin & IEEE_LOEXPO) << 25)+bias) |
        ((tin & IEEE_LOMANT) << 24) | ((tin & IEEE_LOSIGN) << 32) : 0);
  }
  n -= offset;

/* Handle the bulk of them that are aligned on Cray words. */

  for(i=0; i < n-1; i=i+2){
    tin = *ind++;
    *outd++ = (tin & IEEE_HIEXPO ? (((tin & IEEE_HIEXPO) >> 7)+bias) |
        ((tin & IEEE_HIMANT) >> 8 ) |  (tin & IEEE_HISIGN)        : 0);
    *outd++ = (tin & IEEE_LOEXPO ? (((tin & IEEE_LOEXPO) << 25)+bias) |
        ((tin & IEEE_LOMANT) << 24) | ((tin & IEEE_LOSIGN) << 32) : 0);
  }
  n -= i;

/* Handle the last one, if needed, which is not aligned on a Cray word. */

  if(n==1){
    tin = *ind;
    *outd++ = (tin & IEEE_HIEXPO ? (((tin & IEEE_HIEXPO) >> 7)+bias) |
        ((tin & IEEE_HIMANT) >> 8 ) |  (tin & IEEE_HISIGN)        : 0);
  }
}
/************************************************************************/
void packd_c(double *in,char *out,int n)
/*
  Pack an array of Cray reals into IEEE double precision. This assumes
  that a "double" and a "float" are identical.
------------------------------------------------------------------------*/
{
  int *ind,*outd,bias,i,tin;

  ind = (int *)in;
  outd = (int *)out;
  bias = (16384 - 1022) << 48;

  for(i=0; i < n; i++){
    tin = *ind++;
    *outd++ = (tin & CRAY_DMANT ? (tin & SIGN) |
      (((tin & CRAY_EXPO)-bias) << 4) | ((tin & CRAY_DMANT1) << 5) : 0 );
  }
}
/************************************************************************/
void unpackd_c(char *in,double *out,int n)
/*
  Unpack an array of IEEE double precision numbers into Cray reals. This
  assumes that a "double" and a "float" are identical.
------------------------------------------------------------------------*/
{
  int *ind,*outd,bias,i,tin;

  ind = (int *)in;
  outd = (int *)out;
  bias = ((16384 - 1022) << 48) | (1 << 47);

  for(i=0; i < n; i++){
    tin = *ind++;
    *outd++ = (tin & IEEE_DEXPO ? (tin & SIGN) |
      (((tin & IEEE_DEXPO) >> 4) + bias) | ((tin & IEEE_DMANT) >> 5) : 0 );
  }
}

#endif