#
# Create Core.xs
# - needed since we allow bad pixel handling to be switched off
#

use strict;

use Config;
use File::Basename qw(&basename &dirname);

# check for bad value support
use vars qw( $bvalflag $usenan $bvalPerPdl );
require "badsupport.p";

# are we big or little endian?
require PDL::Core::Dev;
my $isbigendian = PDL::Core::Dev::isbigendian();

# This forces PL files to create target in same directory as PL file.
# This is so that make depend always knows where to find PL derivatives.
chdir(dirname($0));
my $file;
($file = basename($0)) =~ s/\.PL$//;
$file =~ s/\.pl$//
        if ($Config{'osname'} eq 'VMS' or
            $Config{'osname'} eq 'OS2');  # "case-forgiving"
 
if ( $bvalflag ) {
    print "Extracting $file (WITH bad value support)\n";
} else {		     
    print "Extracting $file (NO bad value support)\n";
}
open OUT,">$file" or die "Can't create $file: $!";
chmod 0644, $file;

print OUT <<"!WITH!SUBS!";

/*
   Core.xs
   - automatically generated by Core.xs.PL
   - bad value support = $bvalflag
*/
!WITH!SUBS!

print OUT <<'!NO!SUBS!';

#ifndef WIN32
#include <unistd.h>
#include <sys/mman.h>
#include <fcntl.h>
#define USE_MMAP
#endif

#include "EXTERN.h"   /* std perl include */
#include "perl.h"     /* std perl include */
#include "XSUB.h"     /* XSUB include */

#if defined(CONTEXT)
#undef CONTEXT
#endif

#define PDL_CORE      /* For certain ifdefs */
#include "pdl.h"      /* Data structure declarations */
#include "pdlcore.h"  /* Core declarations */

!NO!SUBS!

    if ( $bvalflag ) {
	print OUT "#include <float.h>\n" unless $usenan;
	print OUT "#include <limits.h>\n";

    } # if: $bvalflag

print OUT <<'!NO!SUBS!';

/* Return a integer or numeric scalar as approroate */

#define setflag(reg,flagval,val) (val?(reg |= flagval):(reg &= ~flagval))

#define SET_RETVAL_NV(x) x->datatype<PDL_F ? (RETVAL=newSViv( (IV)result )) : (RETVAL=newSVnv( result ))

Core PDL; /* Struct holding pointers to shared C routines */

#ifdef FOO
Core *pdl__Core_get_Core() /* INTERNAL TO CORE! DON'T CALL FROM OUTSIDE */
{
	return PDL;
}
#endif

int pdl_debugging=0;

#define CHECKP(p)    if ((p) == NULL) croak("Out of memory")

static int* pdl_packint( SV* sv, int *ndims ) {

   SV*  bar;
   AV*  array;
   int i;
   int *dims;

   if (!(SvROK(sv) && SvTYPE(SvRV(sv))==SVt_PVAV))  /* Test */
       return NULL;
   array = (AV *) SvRV(sv);   /* dereference */
     *ndims = (int) av_len(array) + 1;  /* Number of dimensions */
   /* Array space */
   dims = (int *) pdl_malloc( (*ndims) * sizeof(*dims) );
   CHECKP(dims);

   for(i=0; i<(*ndims); i++) {
      bar = *(av_fetch( array, i, 0 )); /* Fetch */
      dims[i] = (int) SvIV(bar);
   }
   return dims;
}

static SV* pdl_unpackint ( PDL_Long *dims, int ndims ) {

   AV*  array;
   int i;

   array = newAV();

   for(i=0; i<ndims; i++) /* if ndims == 0, nothing stored -> ok */
         av_store( array, i, newSViv( (IV)dims[i] ) );

   return (SV*) array;
}

/*
 * Free the data if possible; used by mmapper
 * Moved from pdlhash.c July 10 2006 DJB
 */
static void pdl_freedata (pdl *a) {
	if(a->datasv) {
		SvREFCNT_dec(a->datasv);
		a->datasv=0;
		a->data=0;
	} else if(a->data) {
		die("Trying to free data of untouchable (mmapped?) pdl");
	}
}

!NO!SUBS!

    if ( $bvalflag ) {
	print OUT <<'!NO!SUBS!';

#ifdef FOOFOO_PROPOGATE_BADFLAG

/*
 * this seems to cause an infinite loop in between tests 42 & 43 of
 * t/bad.t - ie
 *
 * $a = sequence( byte, 2, 3 );
 * $b = $a->slice("(1),:");
 * my $mask = sequence( byte, 2, 3 );
 * $mask = $mask->setbadif( ($mask % 3) == 2 );
 * print "a,b == ", $a->badflag, ",", $b->badflag, "\n";
 * $a->inplace->copybad( $mask );                          <-- think this is the call
 * print "a,b == ", $a->badflag, ",", $b->badflag, "\n";
 * print "$a $b\n";
 * ok( $b->badflag, 1 );
 * 
 */

/* used by propogate_badflag() */

void propogate_badflag_children( pdl *it, int newval ) {
    PDL_DECL_CHILDLOOP(it)
    PDL_START_CHILDLOOP(it)
    {
	pdl_trans *trans = PDL_CHILDLOOP_THISCHILD(it);
	int i;

	for( i = trans->vtable->nparents;
	     i < trans->vtable->npdls; 
	     i++ ) {
	    
	    pdl *child = trans->pdls[i];

	    if ( newval ) child->state |=  PDL_BADVAL;
            else          child->state &= ~PDL_BADVAL;

	    /* make sure we propogate to grandchildren, etc */
	    propogate_badflag_children( child, newval );

        } /* for: i */
    }
    PDL_END_CHILDLOOP(it)
} /* propogate_badflag_children */

/* used by propogate_badflag() */

void propogate_badflag_parents( pdl *it ) {
    PDL_DECL_CHILDLOOP(it)
    PDL_START_CHILDLOOP(it)
    {
	pdl_trans *trans = PDL_CHILDLOOP_THISCHILD(it);
	int i;

	for( i = 0;
	     i < trans->vtable->nparents; 
	     i++ ) {
	    
	    pdl *parent = trans->pdls[i];

	    /* only sets allowed here */
	    parent->state |= PDL_BADVAL;

	    /* make sure we propogate to grandparents, etc */
	    propogate_badflag_parents( parent );

        } /* for: i */
    }
    PDL_END_CHILDLOOP(it)
} /* propogate_badflag_parents */

/*
 * we want to change the bad flag of the children
 * (newval = 1 means set flag, 0 means clear it).
 * If newval == 1, then we also loop through the
 * parents, setting their bad flag
 *
 * thanks to Christian Soeller for this 
 */

void propogate_badflag( pdl *it, int newval ) {
   /* only do anything if the flag has changed - do we need this check ? */
   if ( newval ) {
      if ( (it->state & PDL_BADVAL) == 0 ) {
         propogate_badflag_parents( it );
         propogate_badflag_children( it, newval );
      }
   } else {
      if ( (it->state & PDL_BADVAL) > 0 ) {
         propogate_badflag_children( it, newval );
      }

   }

} /* propogate_badflag */

#else        /* FOOFOO_PROPOGATE_BADFLAG */

/* newval = 1 means set flag, 0 means clear it */
/* thanks to Christian Soeller for this */

void propogate_badflag( pdl *it, int newval ) {
    PDL_DECL_CHILDLOOP(it)
    PDL_START_CHILDLOOP(it)
    {
	pdl_trans *trans = PDL_CHILDLOOP_THISCHILD(it);
	int i;

	for( i = trans->vtable->nparents;
	     i < trans->vtable->npdls; i++ ) {
	    
	    pdl *child = trans->pdls[i];

	    if ( newval ) child->state |=  PDL_BADVAL;
            else          child->state &= ~PDL_BADVAL;

	    /* make sure we propogate to grandchildren, etc */
	    propogate_badflag( child, newval );

        } /* for: i */
    }
    PDL_END_CHILDLOOP(it)
} /* propogate_badflag */

#endif    /* FOOFOO_PROPOGATE_BADFLAG */

void propogate_badvalue( pdl *it ) {
    PDL_DECL_CHILDLOOP(it)
    PDL_START_CHILDLOOP(it)
    {
	pdl_trans *trans = PDL_CHILDLOOP_THISCHILD(it);
	int i;

	for( i = trans->vtable->nparents;
	     i < trans->vtable->npdls; i++ ) {
	    
	    pdl *child = trans->pdls[i];

            child->has_badvalue = 1;
            child->badvalue = it->badvalue;

	    /* make sure we propogate to grandchildren, etc */
	    propogate_badvalue( child );

        } /* for: i */
    }
    PDL_END_CHILDLOOP(it)
} /* propogate_badvalue */


/* this is horrible - the routines from bad should perhaps be here instead ? */
double pdl_get_badvalue( int datatype ) {
    double retval;
    switch ( datatype ) {
!NO!SUBS!


use PDL::Types;
my $ntypes = $#PDL::Types::names;

my $str;
foreach my $i ( 0 .. $ntypes ) {
    my $type = PDL::Type->new( $i );
    my $typesym = $type->symbol;

    my $cname = $type->ctype;
    $cname =~ s/^PDL_//;
    my $storage = "PDL.bvals.$cname";

    print OUT "\tcase $typesym: retval = $storage; break;\n";
}

print OUT <<'!NO!SUBS!';

      default:
	croak("Unknown type sent to pdl_get_badvalue\n");
    }
    return retval;
} /* pdl_get_badvalue() */


double pdl_get_pdl_badvalue( pdl *it ) {
    double retval;
    int datatype;

!NO!SUBS!

	if ($bvalPerPdl) {
	    print OUT '
    if (it->has_badvalue) {
        retval = it->badvalue;
    } else {
        datatype = it->datatype;
        retval = pdl_get_badvalue( datatype );
    }
';
	} else {
    print OUT '
    datatype = it->datatype;
    retval = pdl_get_badvalue( datatype );
';
        }
print OUT <<'!NO!SUBS!';
    return retval;
} /* pdl_get_pdl_badvalue() */


!NO!SUBS!

} # if: $bvalflag

print OUT <<'!NO!SUBS!';

MODULE = PDL::Core     PACKAGE = PDL


# Destroy a PDL - note if a hash do nothing, the $$x{PDL} component
# will be destroyed anyway on a separate call

void
DESTROY(sv)
  SV *	sv;
  CODE:
    pdl *self;
    if (SvROK(sv) && SvTYPE(SvRV(sv)) == SVt_PVHV)
       1; /* Do nothing */
    else {
       self = SvPDLV(sv);
       PDLDEBUG_f(printf("DESTROYING %d\n",self);)
       if (self != NULL)
          pdl_destroy(self);
    }

# Return the transformation object or an undef otherwise.

SV *
get_trans(self)
	pdl *self;
	CODE:
	ST(0) = sv_newmortal();
	if(self->trans)  {
		sv_setref_pv(ST(0), "PDL::Trans", (void*)(self->trans));
	} else {
               ST(0) = &PL_sv_undef;
	}

# This will change in the future, as can be seen from the name ;)
# the argument passing is a real quick hack: you can pass 3 integers
# and nothing else.

MODULE = PDL::Core	PACKAGE = PDL::Trans
void
call_trans_foomethod(trans,i1,i2,i3)
	pdl_trans *trans
	int i1
	int i2
	int i3
	CODE:
	PDL_TR_CHKMAGIC(trans);
	pdl_trans_changesoon(trans,PDL_PARENTDIMSCHANGED|PDL_PARENTDATACHANGED);
	if(trans->vtable->foomethod == NULL) {
		croak("This transformation doesn't have a foomethod!");
	}
	(trans->vtable->foomethod)(trans,i1,i2,i3);
	pdl_trans_changed(trans,PDL_PARENTDIMSCHANGED|PDL_PARENTDATACHANGED);

MODULE = PDL::Core	PACKAGE = PDL

int
iscontig(x)
   pdl*	x
   CODE:
      RETVAL = 1;
      pdl_make_physvaffine( x );
	if PDL_VAFFOK(x) {
	   int i, inc=1;
	   printf("vaff check...\n");
  	   for (i=0;i<x->ndims;i++) {
     	      if (PDL_REPRINC(x,i) != inc) {
		     RETVAL = 0;
		     break;
              }
     	      inc *= x->dims[i];
  	   }
        }
  OUTPUT:
    RETVAL

!NO!SUBS!

# access (read, if set is true then write as well; if postset true then
#         read first and write new value after that)
# to piddle's state
#

my %flags = 
    ( 
      hdrcpy => { set => 1 },
      fflows => { FLAG => "DATAFLOW_F" },
      bflows => { FLAG => "DATAFLOW_B" },
      is_inplace => { FLAG => "INPLACE", postset => 1 },
      donttouch => { FLAG => "DONTTOUCHDATA" },
      allocated => { },
      vaffine => { FLAG => "OPT_VAFFTRANSOK" },
      anychgd => { FLAG => "ANYCHANGED" },
      dimschgd => { FLAG => "PARENTDIMSCHANGED" },
      tracedebug => { FLAG => "TRACEDEBUG", set => 1},
     );

#if ( $bvalflag ) { $flags{baddata} = { set => 1, FLAG => "BADVAL" }; }

foreach my $name ( keys %flags ) {
    my $flag = "PDL_" . ($flags{$name}{FLAG} || uc($name));
    if ( $flags{$name}{set} ) {
	print OUT <<"!WITH!SUBS!";
int
$name(x,mode=0)
	pdl *x
	int mode
	CODE:
	if (items>1) 
	   { setflag(x->state,$flag,mode); }
	RETVAL = ((x->state & $flag) > 0);
	OUTPUT:
	RETVAL

!WITH!SUBS!

} elsif ($flags{$name}{postset}) {

	print OUT <<"!WITH!SUBS!";
int
$name(x,mode=0)
	pdl *x
	int mode
	CODE:
	RETVAL = ((x->state & $flag) > 0);
	if (items>1) 
	   { setflag(x->state,$flag,mode); }
	OUTPUT:
	RETVAL

!WITH!SUBS!

} else {
	print OUT <<"!WITH!SUBS!";
int
$name(self)
	pdl *self
	CODE:
	RETVAL = ((self->state & $flag) > 0);
	OUTPUT:
	RETVAL

!WITH!SUBS!

} 

} # foreach: keys %flags 

=begin windows_mmap

I found this at http://mollyrocket.com/forums/viewtopic.php?p=2529&sid=973b8e0a1e639e3008d7ef05f686c6fa
and thougt we might consider using it to make windows mmapping possible.

-David Mertens

 /*
 This code was placed in the public domain by the author,
 Sean Barrett, in November 2007. Do with it as you will.
 (Seee the page for stb_vorbis or the mollyrocket source
 page for a longer description of the public domain non-license).
 */
 
 #define WIN32_LEAN_AND_MEAN
 #include <windows.h>
 
 typedef struct
 {
    HANDLE f;
    HANDLE m;
    void *p;
 } SIMPLE_UNMMAP;
 
 // map 'filename' and return a pointer to it. fill out *length and *un if not-NULL
 void *simple_mmap(const char *filename, int *length, SIMPLE_UNMMAP *un)
 {
    HANDLE f = CreateFile(filename, GENERIC_READ, FILE_SHARE_READ,  NULL, OPEN_EXISTING, FILE_ATTRIBUTE_NORMAL, NULL);
    HANDLE m;
    void *p;
    if (!f) return NULL;
    m = CreateFileMapping(f, NULL, PAGE_READONLY, 0,0, NULL);
    if (!m) { CloseHandle(f); return NULL; }
    p = MapViewOfFile(m, FILE_MAP_READ, 0,0,0);
    if (!p) { CloseHandle(m); CloseHandle(f); return NULL; }
    if (n) *n = GetFileSize(f, NULL);
    if (un) {
       un->f = f;
       un->m = m;
       un->p = p;
    }
    return p;
 }
 
 void simple_unmmap(SIMPLE_UNMMAP *un)
 {
    UnmapViewOfFile(un->p);
    CloseHandle(un->m);
    CloseHandle(un->f);
 } 
 
=end windows_mmap

=cut

print OUT <<'!NO!SUBS!';

void
set_inplace(self,val)
  pdl *self;
  int val;
  CODE:
    setflag(self->state,PDL_INPLACE,val);

int
address(self)
  pdl *self;
  CODE:
    RETVAL = (int) self;
  OUTPUT:
    RETVAL

pdl *
pdl_hard_copy(src)
	pdl *src;

pdl *
sever(src)
	pdl *src;
	CODE:
		if(src->trans) {
			pdl_make_physvaffine(src);
			pdl_destroytransform(src->trans,1);
		}
		RETVAL=src;
	OUTPUT:
		RETVAL

int
set_data_by_mmap(it,fname,len,writable,shared,creat,mode,trunc)
	pdl *it
	char *fname
	int len
	int writable
	int shared
	int creat
	int mode
	int trunc
	CODE:
#ifdef USE_MMAP
       int fd;
       pdl_freedata(it);
       fd = open(fname,(writable && shared ? O_RDWR : O_RDONLY)|
               (creat ? O_CREAT : 0),mode);
       if(fd < 0) {
               croak("Error opening file");
       }
       if(trunc) {
               ftruncate(fd,0);   /* Clear all previous data */
               ftruncate(fd,len); /* And make it long enough */
       }
       if(len) {
		it->data = mmap(0,len,PROT_READ | (writable ?
					PROT_WRITE : 0),
				(shared ? MAP_SHARED : MAP_PRIVATE),
				fd,0);
		if(!it->data)
			croak("Error mmapping!");
       } else {
               /* Special case: zero-length file */
               it->data = NULL;
       }
       PDLDEBUG_f(printf("PDL::MMap: mapped to %d\n",it->data);)
       it->state |= PDL_DONTTOUCHDATA | PDL_ALLOCATED;
       pdl_add_deletedata_magic(it, pdl_delete_mmapped_data, len);
       close(fd);
#else
	croak("mmap not supported on this architecture");
#endif
       RETVAL = 1;
OUTPUT:
       RETVAL


int
set_data_by_offset(it,orig,offset)
      pdl *it
      pdl *orig
      int offset
      CODE:
              pdl_freedata(it);
              it->data = ((char *) orig->data) + offset;
              it->state |= PDL_DONTTOUCHDATA | PDL_ALLOCATED;
              RETVAL = 1;
      OUTPUT:
              RETVAL

int
nelem(x)
	pdl *x
	CODE:
		pdl_make_physdims(x);
		RETVAL = x->nvals;
	OUTPUT:
		RETVAL

# Convert PDL to new datatype (called by float(), int() etc.)

# SV *
# convert(a,datatype)
#    pdl*	a
#    int	datatype
#    CODE:
#     pdl* b;
#     pdl_make_physical(a);
#     RETVAL = pdl_copy(a,""); /* Init value to return */
#     b = SvPDLV(RETVAL);      /* Map */
#     pdl_converttype( &b, datatype, PDL_PERM );
#     PDLDEBUG_f(printf("converted %d, %d, %d, %d\n",a, b, a->datatype, b->datatype));

#     OUTPUT:
#      RETVAL


# Call my howbig function

int
howbig_c(datatype)
   int	datatype
   CODE:
     RETVAL = pdl_howbig(datatype);
   OUTPUT:
     RETVAL

MODULE = PDL::Core     PACKAGE = PDL::Core

unsigned int
is_scalar_SvPOK(arg)
	SV* arg;
	CODE:
	RETVAL = SvPOK(arg);
	OUTPUT:
	RETVAL


int
set_debugging(i)
	int i;
	CODE:
	RETVAL = pdl_debugging;
	pdl_debugging = i;
	OUTPUT:
	RETVAL


SV *
sclr_c(it)
   pdl* it
   CODE:
	PDL_Long nullp = 0;
	PDL_Long dummyd = 1;
	PDL_Long dummyi = 1;
	double result;

        /* get the first element of a piddle and return as
         * Perl double scalar (NV)
         */
        pdl_make_physvaffine( it );
	if (it->nvals < 1)
           croak("piddle must have at least one element");
	/* offs = PDL_REPROFFS(it); */
        /* result = pdl_get_offs(PDL_REPRP(it),offs); */
        result=pdl_at(PDL_REPRP(it), it->datatype, &nullp, &dummyd,
        &dummyi, PDL_REPROFFS(it),1);
        SET_RETVAL_NV(it) ;

    OUTPUT:
        RETVAL


SV *
at_c(x,position)
   pdl*	x
   PDL_Long *	pos = NO_INIT
   CODE:
    int npos, ipos;
    double result;

      pdl_make_physvaffine( x );

    pos = pdl_packdims( ST(1), &npos);
    
    if (pos == NULL || npos < x->ndims)
       croak("Invalid position");

    /*  allow additional trailing indices
     *  which must be all zero, i.e. a
     *  [3,1,5] piddle is treated as an [3,1,5,1,1,1,....]
     *  infinite dim piddle
     */
    for (ipos=x->ndims; ipos<npos; ipos++)
      if (pos[ipos] != 0)
         croak("Invalid position");

    result=pdl_at(PDL_REPRP(x), x->datatype, pos, x->dims,
        (PDL_VAFFOK(x) ? x->vafftrans->incs : x->dimincs), PDL_REPROFFS(x),
	x->ndims);

    SET_RETVAL_NV(x) ;

    OUTPUT:
     RETVAL

SV *
at_bad_c(x,position)
   pdl*	x
   PDL_Long *	pos = NO_INIT
   CODE:
    int npos, ipos, badflag;
    double result;

      pdl_make_physvaffine( x );

    pos = pdl_packdims( ST(1), &npos);
    
    if (pos == NULL || npos < x->ndims)
       croak("Invalid position");

    /*  allow additional trailing indices
     *  which must be all zero, i.e. a
     *  [3,1,5] piddle is treated as an [3,1,5,1,1,1,....]
     *  infinite dim piddle
     */
    for (ipos=x->ndims; ipos<npos; ipos++)
      if (pos[ipos] != 0)
         croak("Invalid position");

    result=pdl_at(PDL_REPRP(x), x->datatype, pos, x->dims,
        (PDL_VAFFOK(x) ? x->vafftrans->incs : x->dimincs), PDL_REPROFFS(x),
	x->ndims);
!NO!SUBS!

if ( $bvalflag ) {
    print OUT 
'
   badflag = (x->state & PDL_BADVAL) > 0;

';

    # do we have to bother about NaN's?
    if ( $usenan ) { 
       print OUT
'
   if ( badflag &&
        ( ( x->datatype < 4 && ( result == pdl_get_badvalue( x->datatype ) ) ) ||
          ( x->datatype >= 4 && ( finite(result) == 0 ) )
        )
      ) {
	 RETVAL = newSVpvn( "BAD", 3 );
   } else
';
   } else {
       print OUT
'
   if ( badflag && 
        ( pdl_get_badvalue( x->datatype ) == result )
      ) {
	 RETVAL = newSVpvn( "BAD", 3 );
   } else
';

   } # if: $usenan
} # if: $bvalflag

print OUT <<'!NO!SUBS!';
    SET_RETVAL_NV(x) ;

    OUTPUT:
     RETVAL


void
list_c(x)
	pdl *x
	PPCODE:
	PDL_Long *inds,*incs,offs;
	void *data;
	int ind;
	int stop = 0;
        pdl_make_physvaffine( x );
	inds = pdl_malloc(sizeof(PDL_Long) * x->ndims); /* GCC -> on stack :( */

	data = PDL_REPRP(x);
	incs = (PDL_VAFFOK(x) ? x->vafftrans->incs : x->dimincs);
	offs = PDL_REPROFFS(x);
	EXTEND(sp,x->nvals);
	for(ind=0; ind < x->ndims; ind++) inds[ind] = 0;
	while(!stop) {
		PUSHs(sv_2mortal(newSVnv(pdl_at( data, x->datatype,
			inds, x->dims, incs, offs, x->ndims))));
		stop = 1;
		for(ind = 0; ind < x->ndims; ind++)
			if(++(inds[ind]) >= x->dims[ind])
				inds[ind] = 0;
			else
				{stop = 0; break;}
	}

# returns the string 'BAD' if an element is bad
#

SV *
listref_c(x)
   pdl *x
  CODE:
   PDL_Long *inds,*incs,offs;
   void *data;
   int ind, lind;
   int stop = 0;
   AV *av;

!NO!SUBS!

if ( $bvalflag ) {
    # note: 
    #  the badvalue is stored in a double, but that's what pdl_at()
    #  returns

    print OUT 
'
   SV *sv;
   double pdl_val, pdl_badval;
   int badflag = (x->state & PDL_BADVAL) > 0;

';

    # do we have to bother about NaN's?
    if ( $usenan ) { 
       print OUT
'
   if ( badflag && x->datatype < 4 ) {
      pdl_badval = pdl_get_pdl_badvalue( x ); 
   }
';
   } else {
       print OUT
'
   if ( badflag ) {
      pdl_badval = pdl_get_pdl_badvalue( x ); 
   }
';

   } # if: $usenan
} # if: $bvalflag

print OUT <<'!NO!SUBS!';

   pdl_make_physvaffine( x );
   inds = pdl_malloc(sizeof(PDL_Long) * x->ndims); /* GCC -> on stack :( */
   data = PDL_REPRP(x);
   incs = (PDL_VAFFOK(x) ? x->vafftrans->incs : x->dimincs);
   offs = PDL_REPROFFS(x);
   av = newAV();
   av_extend(av,x->nvals);
   lind=0;
   for(ind=0; ind < x->ndims; ind++) inds[ind] = 0;
   while(!stop) {

!NO!SUBS!

    if ( $bvalflag ) {

	my $condition;
	if ( $usenan ) {
	    $condition = '( (x->datatype < 4 && pdl_val == pdl_badval) ||
			(x->datatype >= 4 && finite(pdl_val) == 0) )';
	} else {
	    $condition = 'pdl_val == pdl_badval';
	}

	print OUT <<"!WITH!SUBS!";
      pdl_val = pdl_at( data, x->datatype, inds, x->dims, incs, offs, x->ndims );
      if ( badflag && $condition ) {
	 sv = newSVpvn( "BAD", 3 );
      } else {
	 sv = newSVnv( pdl_val );
      }
      av_store( av, lind, sv );
!WITH!SUBS!

    } else {

	print OUT <<'!NO!SUBS!';
      av_store(av,lind,
               newSVnv( pdl_at( data, x->datatype,
	       inds, x->dims, incs, offs, x->ndims ) )
               );
!NO!SUBS!

} # bvalflag

print OUT <<'!NO!SUBS!';

      lind++;
      stop = 1;
      for(ind = 0; ind < x->ndims; ind++) {
	 if(++(inds[ind]) >= x->dims[ind]) {
       	    inds[ind] = 0;
         } else {
       	    stop = 0; break;
         }
      }
   }
   RETVAL = newRV_noinc((SV *)av);
  OUTPUT:
   RETVAL

void
set_c(x,position,value)
    pdl*	x
    PDL_Long *	pos = NO_INIT
    double	value
   CODE:
    int npos,ipos;

    pdl_make_physvaffine( x );

    pos = pdl_packdims( ST(1), &npos);
    if (pos == NULL || npos < x->ndims)
       croak("Invalid position");

    /*  allow additional trailing indices
     *  which must be all zero, i.e. a
     *  [3,1,5] piddle is treated as an [3,1,5,1,1,1,....]
     *  infinite dim piddle
     */
    for (ipos=x->ndims; ipos<npos; ipos++)
      if (pos[ipos] != 0)
         croak("Invalid position");

    pdl_children_changesoon( x , PDL_PARENTDATACHANGED );
    pdl_set(PDL_REPRP(x), x->datatype, pos, x->dims,
        (PDL_VAFFOK(x) ? x->vafftrans->incs : x->dimincs), PDL_REPROFFS(x),
	x->ndims,value);
    if (PDL_VAFFOK(x))
       pdl_vaffinechanged(x, PDL_PARENTDATACHANGED);
    else
       pdl_changed( x , PDL_PARENTDATACHANGED , 0 );

BOOT:

   /* Initialize structure of pointers to core C routines */

   PDL.Version     = PDL_CORE_VERSION;
   PDL.SvPDLV      = SvPDLV;
   PDL.SetSV_PDL   = SetSV_PDL;
   PDL.create      = pdl_create;
   PDL.pdlnew      = pdl_external_new;
   PDL.tmp         = pdl_external_tmp;
   PDL.destroy     = pdl_destroy;
   PDL.null        = pdl_null;
   PDL.copy        = pdl_copy;
   PDL.hard_copy   = pdl_hard_copy;
   PDL.converttype = pdl_converttype;
   PDL.twod        = pdl_twod;
   PDL.smalloc     = pdl_malloc;
   PDL.howbig      = pdl_howbig;
   PDL.packdims    = pdl_packdims;
   PDL.unpackdims  = pdl_unpackdims;
   PDL.setdims     = pdl_setdims;
   PDL.grow        = pdl_grow;
   PDL.flushcache  = NULL;
   PDL.reallocdims = pdl_reallocdims;
   PDL.reallocthreadids = pdl_reallocthreadids;
   PDL.resize_defaultincs = pdl_resize_defaultincs;
   PDL.get_threadoffsp = pdl_get_threadoffsp;
   PDL.thread_copy = pdl_thread_copy;
   PDL.clearthreadstruct = pdl_clearthreadstruct;
   PDL.initthreadstruct = pdl_initthreadstruct;
   PDL.startthreadloop = pdl_startthreadloop;
   PDL.iterthreadloop = pdl_iterthreadloop;
   PDL.freethreadloop = pdl_freethreadloop;
   PDL.thread_create_parameter = pdl_thread_create_parameter;
   PDL.add_deletedata_magic = pdl_add_deletedata_magic;

   PDL.setdims_careful = pdl_setdims_careful;
   PDL.put_offs = pdl_put_offs;
   PDL.get_offs = pdl_get_offs;
   PDL.get = pdl_get;
   PDL.set_trans_childtrans = pdl_set_trans_childtrans;
   PDL.set_trans_parenttrans = pdl_set_trans_parenttrans;
   PDL.make_now = pdl_make_now;

   PDL.get_convertedpdl = pdl_get_convertedpdl;

   PDL.make_trans_mutual = pdl_make_trans_mutual;
   PDL.trans_mallocfreeproc = pdl_trans_mallocfreeproc;
   PDL.make_physical = pdl_make_physical;
   PDL.make_physdims = pdl_make_physdims;
   PDL.make_physvaffine = pdl_make_physvaffine;
   PDL.pdl_barf      = pdl_barf;
   PDL.allocdata     = pdl_allocdata;
   PDL.safe_indterm  = pdl_safe_indterm;
   PDL.children_changesoon = pdl_children_changesoon;
   PDL.changed       = pdl_changed;
   PDL.vaffinechanged = pdl_vaffinechanged;

   PDL.NaN_float  = union_nan_float.f;
   PDL.NaN_double = union_nan_double.d;

!NO!SUBS!

if ( $bvalflag ) {
    print OUT <<'!NO!SUBS!';

    PDL.propogate_badflag = propogate_badflag;
    PDL.propogate_badvalue = propogate_badvalue;
    PDL.get_pdl_badvalue = pdl_get_pdl_badvalue;
!NO!SUBS!

  for my $type (PDL::Types::types) {
    my $typename = $type->ctype;
    $typename =~ s/^PDL_//;
    my $bval = $type->defbval;
    my $ctype = $type->ctype;
    if ($usenan && $type->usenan) {
      # note: no defaults if usenan
      print OUT "\tPDL.bvals.$typename = PDL.NaN_$type;\n"; #Core NaN value
    } else {
      print OUT
	"\tPDL.bvals.$typename = PDL.bvals.default_$typename = $bval;\n";
    }
  }
#    PDL.bvals.Byte   = PDL.bvals.default_Byte   = UCHAR_MAX;
#    PDL.bvals.Short  = PDL.bvals.default_Short  = SHRT_MIN;
#    PDL.bvals.Ushort = PDL.bvals.default_Ushort = USHRT_MAX;
#    PDL.bvals.Long   = PDL.bvals.default_Long   = INT_MIN;

} # if: $bvalflag

print OUT <<'!NO!SUBS!';
   /*
      "Publish" pointer to this structure in perl variable for use
       by other modules
   */

   sv_setiv(get_sv("PDL::SHARE",TRUE|GV_ADDMULTI), PTR2IV(&PDL));

# version of eval() which propogates errors encountered in
# any internal eval(). Must be passed a code reference - could
# be use perl_eval_sv() but that is still buggy. This subroutine is
# primarily for the perlDL shell to use.
#
# Thanks to Sarathy (gsar@engin.umich.edu) for suggesting this, though
# it needs to be wrapped up in the stack stuff to avoid certain SEGVs!

void
myeval(code)
  SV *	code;
  PROTOTYPE: $
  CODE:
   PUSHMARK(sp) ;
   perl_call_sv(code, G_EVAL|G_KEEPERR|GIMME_V);


# make piddle belonging to 'class' and of type 'type'
# from avref 'array_ref' which is checked for being
# rectangular first

SV*
pdl_avref(array_ref, class, type)
     SV* array_ref
     char* class
     int type
  CODE:
     /* make a piddle from a Perl array ref */
     AV *dims, *av;
     int i, depth;
     int datalevel = -1;
     SV* psv;
     pdl* p;

     if (!SvROK(array_ref))
       croak("pdl_avref: not a reference");
     if (SvTYPE(SvRV(array_ref)) != SVt_PVAV)
       croak("pdl_avref: not an array reference");
     av = (AV *) SvRV(array_ref);
     dims = (AV *) sv_2mortal( (SV *) newAV());

     av_store(dims,0,newSViv((IV) av_len(av)+1));
     /* even if we contain nothing depth is one */
     depth = 1 + av_ndcheck(av,dims,0,&datalevel);

     /* printf("will make type %s\n",class); */
     /* 
	at this stage start making a piddle and populate it with
	values from the array (which has already been checked in av_check)
     */
     if (strcmp(class,"PDL") == 0) {
        p = pdl_from_array(av,dims,type,NULL); /* populate with data */
        ST(0) = sv_newmortal();
        SetSV_PDL(ST(0),p);
     } else {
       /* call class->initialize method */
       PUSHMARK(SP);
       XPUSHs(sv_2mortal(newSVpv(class, 0)));
       PUTBACK;
       perl_call_method("initialize", G_SCALAR);
       SPAGAIN;
       psv = POPs;
       PUTBACK;
       p = SvPDLV(psv); /* and get piddle from returned object */
       ST(0) = psv;
       pdl_from_array(av,dims,type,p); /* populate ;) */
     }

MODULE = PDL::Core	PACKAGE = PDL

# pdl_null is created/imported with no PREFIX  as pdl_null.
#  'null' is supplied in Core.pm that calls 'initialize' which calls
#   the pdl_null here

pdl *
pdl_null(...)


MODULE = PDL::Core     PACKAGE = PDL::Core     PREFIX = pdl_

int
pdl_pthreads_enabled()

MODULE = PDL::Core	PACKAGE = PDL	PREFIX = pdl_

int
isnull(self)
	pdl *self;
	CODE:
		RETVAL= !!(self->state & PDL_NOMYDIMS);
	OUTPUT:
		RETVAL

pdl *
make_physical(self)
	pdl *self;
	CODE:
		pdl_make_physical(self);
		RETVAL = self;
	OUTPUT:
		RETVAL

pdl *
make_physvaffine(self)
	pdl *self;
	CODE:
		pdl_make_physvaffine(self);
		RETVAL = self;
	OUTPUT:
		RETVAL


pdl *
make_physdims(self)
	pdl *self;
	CODE:
		pdl_make_physdims(self);
		RETVAL = self;
	OUTPUT:
		RETVAL

void
pdl_dump(x)
  pdl *x;

void
pdl_add_threading_magic(it,nthdim,nthreads)
	pdl *it
	int nthdim
	int nthreads

void
pdl_remove_threading_magic(it)
	pdl *it
	CODE:
		pdl_add_threading_magic(it,-1,-1);

MODULE = PDL::Core	PACKAGE = PDL	

SV *
initialize(class)
	SV *class

        PPCODE:
	HV *bless_stash;

        if (SvROK(class)) { /* a reference to a class */
	  bless_stash = SvSTASH(SvRV(class));
        } else {            /* a class name */
          bless_stash = gv_stashsv(class, 0);
        }
        ST(0) = sv_newmortal();
        SetSV_PDL(ST(0),pdl_null());   /* set a null PDL to this SV * */
        ST(0) = sv_bless(ST(0), bless_stash); /* bless appropriately  */
	XSRETURN(1);

SV *
get_dataref(self)
	pdl *self
	CODE:
	if(self->state & PDL_DONTTOUCHDATA) {
		croak("Trying to get dataref to magical (mmaped?) pdl");
	}
	pdl_make_physical(self); /* XXX IS THIS MEMLEAK WITHOUT MORTAL? */
	RETVAL = (newRV(self->datasv));
	OUTPUT:
	RETVAL

int
get_datatype(self)
	pdl *self
	CODE:
	RETVAL = self->datatype;
	OUTPUT:
	RETVAL

int
upd_data(self)
	pdl *self
	CODE:
       STRLEN n_a;
	if(self->state & PDL_DONTTOUCHDATA) {
		croak("Trying to touch dataref of magical (mmaped?) pdl");
	}
       self->data = SvPV((SV*)self->datasv,n_a);
	XSRETURN(0);

void
set_dataflow_f(self,value)
	pdl *self;
	int value;
	CODE:
	if(value)
		self->state |= PDL_DATAFLOW_F;
	else
		self->state &= ~PDL_DATAFLOW_F;

void
set_dataflow_b(self,value)
	pdl *self;
	int value;
	CODE:
	if(value)
		self->state |= PDL_DATAFLOW_B;
	else
		self->state &= ~PDL_DATAFLOW_B;

int
getndims(x)
	pdl *x
	ALIAS:
	     PDL::ndims = 1
	CODE:
		pdl_make_physdims(x);
		RETVAL = x->ndims;
	OUTPUT:
		RETVAL

int
getdim(x,y)
	pdl *x
	int y
	ALIAS:
	     PDL::dim = 1
	CODE:
		pdl_make_physdims(x);
		if (y < 0) y = x->ndims + y;
		if (y < 0) croak("negative dim index too large");
		if (y < x->ndims)
                   RETVAL = x->dims[y];
                else
		   RETVAL = 1; /* return size 1 for all other dims */
	OUTPUT:
		RETVAL

int
getnthreadids(x)
	pdl *x
	CODE:
		pdl_make_physdims(x);
		RETVAL = x->nthreadids;
	OUTPUT:
		RETVAL

int
getthreadid(x,y)
	pdl *x
	int y
	CODE:
		RETVAL = x->threadids[y];
	OUTPUT:
		RETVAL

void
setdims(x,dims)
	pdl *x
	PDL_Long *dims = NO_INIT
	CODE:
	{
		int ndims; int i;
		pdl_children_changesoon(x,PDL_PARENTDIMSCHANGED|PDL_PARENTDATACHANGED);
		dims = pdl_packdims(ST(1),&ndims);
		pdl_reallocdims(x,ndims);
		for(i=0; i<ndims; i++) x->dims[i] = dims[i];
		pdl_resize_defaultincs(x);
		x->threadids[0] = ndims;
 /* make null != dims = [0] */
#ifndef ELIFJELFIJSEJIF
		x->state &= ~PDL_NOMYDIMS;
#else
		   if(ndims == 1 && dims[0] == 0) {
			x->state |= PDL_NOMYDIMS;
		   } else {
			x->state &= ~PDL_NOMYDIMS;
		   }
#endif
		pdl_changed(x,PDL_PARENTDIMSCHANGED|PDL_PARENTDATACHANGED,0);
	}

void
dowhenidle()
	CODE:
		pdl_run_delayed_magic();
		XSRETURN(0);

void
bind(p,c)
	pdl *p
	SV *c
	PROTOTYPE: $&
	CODE:
		pdl_add_svmagic(p,c);
		XSRETURN(0);

void
sethdr(p,h)
	pdl *p
	SV *h
	CODE:
	HV* hash;
		if(p->hdrsv == NULL) {             
		      p->hdrsv =  &PL_sv_undef; /*(void*) newSViv(0);*/
		} 

		/* Throw an error if we're not either undef or hash */
                if ( (h != &PL_sv_undef && h != NULL) && 
		     ( !SvROK(h) || SvTYPE(SvRV(h)) != SVt_PVHV )
		   ) 
		      croak("Not a HASH reference");		

		/* Clear the old header */
		SvREFCNT_dec(p->hdrsv);

		/* Put the new header (or undef) in place */
		if(h == &PL_sv_undef || h == NULL)
		   p->hdrsv = NULL;
		else  
		   p->hdrsv = (void*) newRV( (SV*) SvRV(h) );

SV *
hdr(p)
	pdl *p
	CODE:
		pdl_make_physdims(p);

                /* Make sure that in the undef case we return not */
                /* undef but an empty hash ref. */

                if((p->hdrsv==NULL) || (p->hdrsv == &PL_sv_undef)) {
	            p->hdrsv = (void*) newRV_noinc( (SV*)newHV() );
                }

		RETVAL = newRV( (SV*) SvRV((SV*)p->hdrsv) );

	OUTPUT:
	 RETVAL

# fhdr(p) is implemented in perl; see Core.pm.PL if you're looking for it 
#   --CED 9-Feb-2003  
#

SV *
gethdr(p)
	pdl *p
	CODE:
		pdl_make_physdims(p);

                if((p->hdrsv==NULL) || (p->hdrsv == &PL_sv_undef)) {
	            RETVAL = &PL_sv_undef;
                } else {
		    RETVAL = newRV( (SV*) SvRV((SV*)p->hdrsv) );
                }

	OUTPUT:
	 RETVAL

void
set_datatype(a,datatype)
   pdl *a
   int datatype
   CODE:
    pdl_make_physical(a);
    if(a->trans)
	    pdl_destroytransform(a->trans,1);
/*     if(! (a->state && PDL_NOMYDIMS)) { */
    pdl_converttype( &a, datatype, PDL_PERM );
/*     } */

void
threadover_n(...)
   CODE:
   {
    int npdls = items - 1;
    if(npdls <= 0)
    	croak("Usage: threadover_n(pdl[,pdl...],sub)");
    {
	    int i,sd;
	    pdl **pdls = malloc(sizeof(pdl *) * npdls);
	    int *realdims = malloc(sizeof(int) * npdls);
	    pdl_thread pdl_thr;
	    SV *code = ST(items-1);
	    for(i=0; i<npdls; i++) {
		pdls[i] = SvPDLV(ST(i));
		/* XXXXXXXX Bad */
		pdl_make_physical(pdls[i]);
		realdims[i] = 0;
	    }
	    PDL_THR_CLRMAGIC(&pdl_thr);
	    pdl_initthreadstruct(0,pdls,realdims,realdims,npdls,NULL,&pdl_thr,NULL);
	    pdl_startthreadloop(&pdl_thr,NULL,NULL);
	    sd = pdl_thr.ndims;
	    do {
	    	dSP;
		PUSHMARK(sp);
		EXTEND(sp,items);
		PUSHs(sv_2mortal(newSViv((sd-1))));
		for(i=0; i<npdls; i++) {
			PUSHs(sv_2mortal(newSVnv(
				pdl_get_offs(pdls[i],pdl_thr.offs[i]))));
		}
	    	PUTBACK;
		perl_call_sv(code,G_DISCARD);
	    } while(sd = pdl_iterthreadloop(&pdl_thr,0));
	    pdl_freethreadloop(&pdl_thr);
	    free(pdls);
	    free(realdims);
    }
   }

void
threadover(...)
   CODE:
   {
    int npdls, nothers = -1;
    int targs = items - 4;
    if (items > 0) nothers = SvIV(ST(0));
    if(targs <= 0 || nothers < 0 || nothers >= targs)
    	croak("Usage: threadover(nothers,pdl[,pdl...][,otherpars..],realdims,creating,sub)");
    npdls = targs-nothers;
    {
	    int i,j,nd1,nd2,dtype=0,nc=npdls;
	    SV* rdimslist = ST(items-3);
	    SV* cdimslist = ST(items-2);
	    SV *code = ST(items-1);
	    pdl_thread pdl_thr;
	    pdl **pdls = malloc(sizeof(pdl *) * npdls);
	    pdl **child = malloc(sizeof(pdl *) * npdls);
	    SV **csv = malloc(sizeof(SV *) * npdls);
	    SV **dims = malloc(sizeof(SV *) * npdls);
	    SV **incs = malloc(sizeof(SV *) * npdls);
	    SV **others = malloc(sizeof(SV *) * nothers);
	    int *creating = pdl_packint(cdimslist,&nd2);
	    int *realdims = pdl_packint(rdimslist,&nd1);
	    CHECKP(pdls); CHECKP(child); CHECKP(dims);
	    CHECKP(incs); CHECKP(csv);

	    if (nd1 != npdls || nd2 < npdls)
		croak("threadover: need one realdim and creating flag "
		      "per pdl!");
	    for(i=0; i<npdls; i++) {
		pdls[i] = SvPDLV(ST(i+1));
		if (creating[i])
		  nc += realdims[i];
		else {
		  pdl_make_physical(pdls[i]); /* is this what we want?XXX */
		  dtype = PDLMAX(dtype,pdls[i]->datatype);
		}
	    }
	    for (i=npdls+1; i<=targs; i++)
		others[i-npdls-1] = ST(i);
	    if (nd2 < nc)
		croak("Not enough dimension info to create pdls");
#ifdef DEBUG_PTHREAD
		for (i=0;i<npdls;i++) { /* just for debugging purposes */
		printf("pdl %d Dims: [",i);
		for (j=0;j<realdims[i];j++)
			printf("%d ",pdls[i]->dims[j]);
		printf("] Incs: [");
		for (j=0;j<realdims[i];j++)
			printf("%d ",PDL_REPRINC(pdls[i],j));
		printf("]\n");
	        }
#endif
	    PDL_THR_CLRMAGIC(&pdl_thr);
	    pdl_initthreadstruct(0,pdls,realdims,creating,npdls,
				NULL,&pdl_thr,NULL);
	    for(i=0, nc=npdls; i<npdls; i++)  /* create as necessary */
              if (creating[i]) {
		int *cp = creating+nc;
		pdls[i]->datatype = dtype;
		pdl_thread_create_parameter(&pdl_thr,i,cp,0);
		nc += realdims[i];
		pdl_make_physical(pdls[i]);
		PDLDEBUG_f(pdl_dump(pdls[i]));
		/* And make it nonnull, now that we've created it */
		pdls[i]->state &= (~PDL_NOMYDIMS);
	      }
	    pdl_startthreadloop(&pdl_thr,NULL,NULL);
	    for(i=0; i<npdls; i++) { /* will the SV*'s be properly freed? */
		dims[i] = newRV(pdl_unpackint(pdls[i]->dims,realdims[i]));
		incs[i] = newRV(pdl_unpackint(PDL_VAFFOK(pdls[i]) ?
		pdls[i]->vafftrans->incs: pdls[i]->dimincs,realdims[i]));
		/* need to make sure we get the vaffine (grand)parent */
		if (PDL_VAFFOK(pdls[i]))
		   pdls[i] = pdls[i]->vafftrans->from;
		child[i]=pdl_null();
		/*  instead of pdls[i] its vaffine parent !!!XXX */
		PDL.affine_new(pdls[i],child[i],pdl_thr.offs[i],dims[i],
						incs[i]);
		pdl_make_physical(child[i]); /* make sure we can get at
						the vafftrans          */
		csv[i] = sv_newmortal();
		SetSV_PDL(csv[i], child[i]); /* pdl* into SV* */
	    }
	    do {  /* the actual threadloop */
		pdl_trans_affine *traff;
	    	dSP;
		PUSHMARK(sp);
		EXTEND(sp,npdls);
		for(i=0; i<npdls; i++) {
		   /* just twiddle the offset - quick and dirty */
		   /* we must twiddle both !! */
		   traff = (pdl_trans_affine *) child[i]->trans;
		   traff->offs = pdl_thr.offs[i];
		   child[i]->vafftrans->offs = pdl_thr.offs[i];
		   child[i]->state |= PDL_PARENTDATACHANGED;
		   PUSHs(csv[i]);
		}
		for (i=0; i<nothers; i++)
		  PUSHs(others[i]);   /* pass the OtherArgs onto the stack */
	    	PUTBACK;
		perl_call_sv(code,G_DISCARD);
	    } while (pdl_iterthreadloop(&pdl_thr,0));
	    pdl_freethreadloop(&pdl_thr);
	    free(pdls);  /* should all these be done with pdl_malloc */
	    free(dims);  /* in case the sub barfs ? XXXX            */
	    free(child);
	    free(csv);
	    free(incs);
	    free(others);
    }
   }

!NO!SUBS!