/* data.cc */

#include "data.h"
#include "data_classes.h"
#include <math.h>
#include "units.h"

inline Real sqr( Real r ) { return r * r; }

dc_data::dc_data() {
  valid = false;
}

ostream &dc_data::display( ostream &stream = cout ) {
  return stream << full_type();
}

/* returns true if both have units and units are consistent.
   prints warning message and returns true if only one has units or if both do 
   and units are inconsistent.
   returns true without warning if neither have units */
bool check_units( const dc_data *d1, const dc_data *d2,
		  const char *func_name = "check_units" ) {
  if( d1->has_units() ) {
    if( d2->has_units() ) {
      if( ( ( dc_udata * )d1 )->get_units() ==
	  ( ( dc_udata * )d2 )->get_units() ) return true;

      dc_trace( TRACE_WARNING ) {
	cerr << func_name << " -- " << ( ( dc_udata * )d1 )->get_units()
	     << " " << d1->full_type()  << " and "
	     << ( ( dc_udata * )d2 )->get_units() << " " << d2->full_type()
	     << " have different units\n";
      }
    } else {
      dc_trace( TRACE_WARNING ) {
	cerr << func_name << " -- " << ( ( dc_udata * )d1 )->get_units()
	     << " " << d1->full_type()  << " and " << d2->full_type()
	     << " have different units\n";
      }
    }
  } else if( d2->has_units() ) {
    dc_trace( TRACE_WARNING ) {
      cerr << func_name << " -- " << " " << d1->full_type()  << " and "
	   << ( ( dc_udata * )d2 )->get_units() << " " << d2->full_type()
	   << " have different units\n";
    }
  } else return true;

  return false;
}

/* returns true if both have quivalent units or no units, false otherwise */
bool equiv_units( const dc_data *d1, const dc_data *d2 ) {
  if( d1->has_units() ) {
    if( d2->has_units() && ( ( dc_udata * )d1 )->get_units() ==
	( ( dc_udata * )d2 )->get_units() ) return true; 
  } else if( !d2->has_units() ) return true;

  return false;
}

#define clean2 if(is_temporary()) delete(this); if(D.is_temporary()) delete(&D);
#define opclean2 if(D1.is_temporary()) delete(&D1); if(D2.is_temporary()) delete(&D2);
#define clean  if(is_temporary()) delete(this);

bool set_data( dc_data &dest, dc_data &src ) {
  switch( dest.type() ) {
  case Boolean_t :
    if( src.type() != Boolean_t ) return true;
    ( ( dc_boolean * )&dest )->set( ( ( dc_boolean * ) &src )->get() );
    break;
  case Int_t :
    if( src.type() == Int_t ) {
      ( ( dc_int * )&dest )->set( ( ( dc_int * ) &src )->get() );
    } else if( src.sub_type() == Real_t ) {
      ( ( dc_int * )&dest )->set( ( int )( ( dc_real * ) &src )->get() );
    } else return true;
    break;
  case Pointer_t :
    if( src.type() != Pointer_t ) return true;
    *( ( dc_pointer * )&dest ) = *( ( dc_pointer * )&src );
    break;
  case Real_t :
    switch( dest.sub_type() ) {
    case Real_t :
      if( src.type() == Int_t ) {
	( ( dc_real * )&dest )->set( ( double )( ( dc_int * ) &src )->get() );
      } else if( src.type() == Real_t ) {
	( ( dc_real * )&dest )->set( ( ( dc_real * )&src )->get() );
      } else return true;
      break;
    case Distrib_t :
      if( src.sub_type() != Distrib_t ) return true;
      *( ( dc_distrib * )&dest ) = src;
      break;
    default : return true;
    }
    break;
  case Set_t :
    if( src.type() != Set_t ) return true;
    *( ( dc_set * )&dest ) = *( ( dc_set * )&src );
    break;
  case Shape_t :
    if( src.type() != Shape_t ) return true;
    ( ( dc_shape * )&dest )->set( *( ( dc_shape * )&src ) );
    break;
  case String_t :
    if( src.type() != String_t ) return true;
    ( ( dc_string * )&dest )->set( *( ( dc_string * )&src ) );
    break;
  case Symbol_t :
    if( src.type() != Symbol_t ) return true;
    ( ( dc_symbol * )&dest )->set( *( ( dc_symbol * )&src ) );
    break;
  case Matrix_t :
    switch( src.sub_type() ) {
    case Rect_Matrix_t :
      if( src.sub_type() != Rect_Matrix_t ) return true;
      ( ( dc_rect_matrix *)&dest )->set( ( ( dc_rect_matrix * ) &src )->get() );
      break;
    case Vector_t :
      if( src.sub_type() != Vector_t ) return true;
      ( ( dc_vector * )&dest )->set( ( ( dc_vector * ) &src )->get() );
      break;
    case Triple_t :
      if( src.sub_type() != Triple_t ) return true;
      ( ( dc_triple * )&dest )->set( ( ( dc_triple * ) &src )->get() );
      break;
    default : return true;
    }
    break;
  default : return true;
  }

  if( ( ( dc_udata * )&dest )->has_units() )
    ( ( dc_udata * )&dest )->set_units( ( ( dc_udata * )&src )->get_units() );
  return false;
}

dc_data &dc_data::operator=( dc_data &D ) {
  if( set_data( *this, D ) ) {
    dc_trace( TRACE_ERROR ) {
      cerr << "dc_data::operator= -- cannot assign " << D.full_type() << " to "
	   << full_type() << "\n";
      exit( 1 );
    }
  }
  if( D.is_temporary() ) { delete( &D ); }
  return *this;
}

dc_data &assign( dc_data *&dest, dc_data &source ) {
  if( source.is_temporary() ) {
    if( dest != nil ) {
      if( dest->is_temporary() ) {
	delete( dest );
	dest = &source;
      } else {
	string label = dest->local_label();
	dc_node *parent = dest->get_parent();
	bool visible = dest->is_visible();
	delete( dest );
	dest = &source;
	dest->set_both( label, parent );
	dest->set_visible( visible );
      }
    } else {
      dest = &source;
    }
  } else {
    if( dest == nil ) {
      dest = ( dc_data * )allocate( source.sub_type() );
      *dest = source;
    } else {
      if( dest->is_temporary() ) {
	if( dest->sub_type() != source.sub_type() ) {
	  delete( dest );
	  dest = ( dc_data * )allocate( source.sub_type() );
	  *dest = source;
	} else {
	  *dest = source;
	}
      } else {
	if( dest->sub_type() != source.sub_type() ) {
	  string label = dest->local_label();
	  dc_node *parent = dest->get_parent();
	  bool visible = dest->is_visible();
	  delete( dest );
	  dest = ( dc_data * )allocate( source.sub_type(), label, parent );
	  *dest = source;
	  dest->set_visible( visible );
	}else {
	  *dest = source;
	}
      }
    }
  }
  return *dest;
}

/* arithmetic operations */
/* operations performed based on types of operands */
dc_data &dc_data::operator+( dc_data &D ) { /* addition */
  dc_data *d = op_add( *this, D );
  if( d == nil ) {
    clean2;
    exit( 1 );
  }
  return *d;
}

dc_data *op_add( dc_data &D1, dc_data &D2 ) {
  dc_data *rVal;

  double d;
  switch( D1.type() ) { 
  case Int_t :
    switch( D2.type() ) {
    case Int_t :
      d = ( ( dc_int * )&D1 )->get() + ( ( dc_int * )&D2 )->get();
      rVal = new dc_int( ( long int )d );
      break;
    case Real_t :
      d = ( ( dc_int * )&D1 )->get_real() + ( ( dc_real * )&D2 )->get();
      rVal = new dc_real( d );
      break;
    default :
      cerr << "op_add -- operation not defined on " 
	   << D1.type_string() << " and " << D2.type_string() << "\n";
      return nil;
    } 
    break;
  case Real_t :
    switch( D2.type() ) {
    case Int_t :
      return op_add( D2, D1 );
    case Real_t :
      d = ( ( dc_real * )&D1 )->get() + ( ( dc_real * )&D2 )->get();
      rVal = new dc_real( d );
      break;
    default :
      cerr << "op_add -- operation not defined on " 
	   << D1.type_string() << " and " << D2.type_string() << "\n";
      return nil;
    } 
    break;
  case Matrix_t :
    switch( D1.sub_type() ) {
    case Rect_Matrix_t :
      if( ( D2.sub_type() != Rect_Matrix_t ) || 
	  ( ( ( dc_rect_matrix * )&D1 )->Nrows() != 
	    ( ( dc_rect_matrix * )&D2 )->Nrows() ) ||
	  ( ( ( dc_rect_matrix * )&D1 )->Ncols() != 
	    ( ( dc_rect_matrix * )&D2 )->Ncols() ) ) {
	cerr << "op_add -- operation not defined on " 
	     << D1.type_string() << " and " << D2.type_string() << "\n";
	return nil;
      }
      rVal = new dc_rect_matrix( ( ( dc_rect_matrix * )&D1 )->get() + 
				 ( ( dc_rect_matrix * )&D2 )->get() );
      break;
    case Triple_t :
      if( ( D2.sub_type() != Triple_t ) || 
	  ( ( ( dc_triple * )&D1 )->dim() != 
	    ( ( dc_triple * )&D2 )->dim() ) ) {
	cerr << "op_add -- operation not defined on " 
	     << D1.type_string() << " and " << D2.type_string() << "\n";
	return nil;
      }
      rVal = new dc_triple( ( ( dc_triple * )&D1 )->get() + 
			    ( ( dc_triple * )&D2 )->get() );
      break;
    case Vector_t :
      if( ( D2.sub_type() != Vector_t ) || 
	  ( ( ( dc_vector * )&D1 )->dim() != 
	    ( ( dc_vector * )&D2 )->dim() ) ) {
	cerr << "op_add -- operation not defined on " 
	     << D1.type_string() << " and " << D2.type_string() << "\n";
	return nil;
      }
      rVal = new dc_vector( ( ( dc_vector * )&D1 )->get() + 
			    ( ( dc_vector * )&D2 )->get() );
      break;
    default :
      cerr << "op_add -- operation not defined on " 
	   << D1.type_string() << " and " << D2.type_string() << "\n";
      return nil;
    }
    break;
  case String_t :
    if( D2.type() == String_t ) {
      rVal = ( dc_data * )&( ( *( ( dc_string * )&D1 ) ) + 
			     ( *( ( dc_string * )&D2 ) ) );
      break;
    }
  default :
    cerr << "op_add -- operation not defined on " 
	 << D1.type_string() << " and " << D2.type_string() << "\n";
    return nil;
  }
  if( rVal->has_units() ) {
    if( check_units( &D1, &D2, "op_add" ) )
      ( ( dc_udata * )rVal )->set_units( ( ( dc_udata * )&D1 )->get_units() );
    else {
      delete( rVal );
      return nil;
    }
  }

  /* clean after error checking */
  opclean2;
  return rVal;
}

dc_data &dc_data::operator-( dc_data &D ) { /* subtraction */
  dc_data *d = op_subtract( *this, D );
  if( d == nil ) {
    clean2;
    exit( 1 );
  }
  return *d;
}

dc_data *op_subtract( dc_data &D1, dc_data &D2 ) {
  dc_data *rVal;

  double d;
  switch( D1.type() ) { 
  case Int_t :
    switch( D2.type() ) {
    case Int_t :
      d = ( ( dc_int * )&D1 )->get() - ( ( dc_int * )&D2 )->get_real();
      rVal = new dc_int( ( long int )d );
      break;
    case Real_t :
      d = ( ( dc_int * )&D1 )->get_real() - ( ( dc_real * )&D2 )->get();
      rVal = new dc_real( d );
      break;
    default :
      cerr << "op_subtract -- operation not defined on " 
	   << D1.type_string() << " and " << D2.type_string() << "\n";
      return nil;
    } 
    break;
  case Real_t :
    switch( D2.type() ) {
    case Int_t :
      d = ( ( dc_real * )&D1 )->get() - ( ( dc_int * )&D2 )->get_real();
      rVal = new dc_real( d );
      break;
    case Real_t :
      d = ( ( dc_real * )&D1 )->get() - ( ( dc_real * )&D2 )->get();
      rVal = new dc_real( d );
      break;
    default :
      cerr << "op_subtract -- operation not defined on " 
	   << D1.type_string() << " and " << D2.type_string() << "\n";
      return nil;
    } 
    break;
  case Matrix_t :
    switch( D1.sub_type() ) {
    case Rect_Matrix_t :
      if( ( D2.sub_type() != Rect_Matrix_t ) || 
	  ( ( ( dc_rect_matrix * )&D1 )->Nrows() != 
	    ( ( dc_rect_matrix * )&D2 )->Nrows() ) ||
	  ( ( ( dc_rect_matrix * )&D1 )->Ncols() != 
	    ( ( dc_rect_matrix * )&D2 )->Ncols() ) ) {
	cerr << "op_subtract -- operation not defined on " 
	     << D1.type_string() << " and " << D2.type_string() << "\n";
	return nil;
      }
      rVal = new dc_rect_matrix( ( ( dc_rect_matrix * )&D1 )->get() - 
				 ( ( dc_rect_matrix * )&D2 )->get() );
      break;
    case Triple_t :
      if( ( D2.sub_type() != Triple_t ) || 
	  ( ( ( dc_triple * )&D1 )->dim() != 
	    ( ( dc_triple * )&D2 )->dim() ) ) {
	cerr << "op_subtract -- operation not defined on " 
	     << D1.type_string() << " and " << D2.type_string() << "\n";
	return nil;
      }
      rVal = new dc_triple( ( ( dc_triple * )&D1 )->get() - 
			    ( ( dc_triple * )&D2 )->get() );
      break;
    case Vector_t :
      if( ( D2.sub_type() != Vector_t ) || 
	  ( ( ( dc_vector * )&D1 )->dim() != 
	    ( ( dc_vector * )&D2 )->dim() ) ) {
	cerr << "op_subtract -- operation not defined on " 
	     << D1.type_string() << " and " << D2.type_string() << "\n";
	return nil;
      }
      rVal = new dc_vector( ( ( dc_vector * )&D1 )->get() - 
			    ( ( dc_vector * )&D2 )->get() );
      break;
    default :
      cerr << "op_subtract -- operation not defined on " 
	   << D1.type_string() << " and " << D2.type_string() << "\n";
      return nil;
    }
    break;
  default :
    cerr << "op_subtract -- operation not defined on " 
	 << D1.type_string() << " and " << D2.type_string() << "\n";
    return nil;
  }
  if( rVal->has_units() ) {
    if( check_units( &D1, &D2, "op_subtract" ) )
      ( ( dc_udata * )rVal )->set_units( ( ( dc_udata * )&D1 )->get_units() );
    else {
      delete( rVal );
      return nil;
    }
  }

  /* clean after error checking */
  opclean2;
  return rVal;
}

dc_data &dc_data::operator*( dc_data &D ) { /* multiplication */
  dc_data *rVal;

  switch( type() ) { 
  case Int_t :
    switch( D.type() ) {
    case Int_t :
      rVal = new dc_int( ( ( dc_int * )this )->get() * 
			 ( ( dc_int * )&D )->get() );
      break;
    case Matrix_t :
      switch( D.sub_type() ) {
      case Rect_Matrix_t :
	if( !D.is_valid() ){
	  cerr << "operator* -- operation not defined on " 
	       << type_string() << " and " << D.type_string() << "\n";
	  exit( 1 );
	}
	rVal = new dc_rect_matrix( ( ( dc_rect_matrix * )&D )->get() *
				   ( ( dc_int * )this )->get() );
	break; 
      case Triple_t :
	rVal = new dc_triple( ( ( dc_triple * )&D )->get() * 
			      ( ( dc_int * )this )->get() );
	break;
      case Vector_t :
	rVal = new dc_vector( ( ( dc_vector * )&D )->get() * 
			      ( ( dc_int * )this )->get() );
	break;
      default :
	cerr << "operator* -- operation not defined on " 
	     << type_string() << " and " << D.type_string() << "\n";
	exit( 1 );
      } 
      break;
    case Real_t :
      rVal = new dc_real( ( ( dc_int * )this )->get() * 
			  ( ( dc_real * )&D )->get() );
      break;
    default :
      cerr << "operator* -- operation not defined on " 
	   << type_string() << " and " << D.type_string() << "\n";
      exit( 1 );
    } 
    break;
  case Real_t :
    switch( D.type() ) {
    case Int_t :
      return( D * *this );
    case Matrix_t :
      switch( D.sub_type() ) {
      case Rect_Matrix_t :
	if( !D.is_valid() ){
	  cerr << "operator* -- operation not defined on " 
	       << type_string() << " and " << D.type_string() << "\n";
	  exit( 1 );
	}
	rVal = new dc_rect_matrix( ( ( dc_rect_matrix * )&D )->get() *
				   ( ( dc_real * )this )->get() );
	break;
      case Triple_t :
	rVal = new dc_triple( ( ( dc_triple * )&D )->get() * 
			      ( ( dc_real * )this )->get() );
	break;
      case Vector_t :
	rVal = new dc_vector( ( ( dc_vector * )&D )->get() * 
			      ( ( dc_real * )this )->get() );
	break;
      default :
	cerr << "operator* -- operation not defined on " 
	     << type_string() << " and " << D.type_string() << "\n";
	exit( 1 );
      } 
      break;
    case Real_t :
      rVal = new dc_real( ( ( dc_real * )this )->get() * 
			  ( ( dc_real * )&D )->get() );
      break;
    default :
      cerr << "operator* -- operation not defined on " 
	   << type_string() << " and " << D.type_string() << "\n";
      exit( 1 );
    } 
    break;
  case Matrix_t :
    switch( sub_type() ) {
    case Rect_Matrix_t :
      switch( D.type() ) {
      case Int_t :
	return ( D * *this );
	break;
      case Matrix_t :
	switch( D.sub_type() ) {
	case Rect_Matrix_t :
	  if( ( ( ( dc_rect_matrix * )this )->Ncols() != 
		( ( dc_rect_matrix * )&D )->Nrows() ) ) {
	    cerr << "operator* -- operation not defined on " 
		 << type_string() << " and " << D.type_string() << "\n";
	    exit( 1 );
	  }
	  rVal = new dc_rect_matrix( ( ( dc_rect_matrix * )this )->get() *
				     ( ( dc_rect_matrix * )&D )->get() );
	  break;
	case Vector_t :
	  if( ( ( ( dc_rect_matrix * )this )->Ncols() != 
		( ( dc_vector * )&D )->dim() ) ) {
	    cerr << "operator* -- operation not defined on " 
		 << type_string() << " and " << D.type_string() << "\n";
	    exit( 1 );
	  }
	  rVal = new dc_vector( ( ( dc_rect_matrix * )this )->get() *
				( ( dc_vector * )&D )->get() );
	  break;
	default :
	  cerr << "operator* -- operation not defined on " 
	       << type_string() << " and " << D.type_string() << "\n";
	  exit( 1 );
	}
	break;
      case Real_t :
	return ( D * *this );
	break;
      default :
	cerr << "operator* -- operation not defined on " 
	     << type_string() << " and " << D.type_string() << "\n";
	exit( 1 );
      }
      break;
    case Vector_t :
      switch( D.type() ) {
      case Int_t :
	return ( D * *this );
	break;
      case Real_t :
	return ( D * *this );
	break;
      case Matrix_t :
	switch( D.sub_type() ) {
	case Rect_Matrix_t :
	  rVal = new dc_vector( ( ( dc_vector * )this )->get_as_row() *
				( ( dc_rect_matrix * )&D )->get() );
	  break;
	case Vector_t : /* cross multiplication */
	  rVal = new dc_rect_matrix( ( ( ( dc_vector * )this )->get() * 
				       ( ( dc_vector * )&D )->get_as_row()) );
	  break;
	default :
	  cerr << "operator* -- operation not defined on " 
	       << type_string() << " and " << D.type_string() << "\n";
	  exit( 1 );
	}
	break;
      default :
	cerr << "operator* -- operation not defined on " 
	     << type_string() << " and " << D.type_string() << "\n";
	exit( 1 );
      }
      break;
    case Triple_t :
      switch( D.type() ) {
      case Int_t :
	return ( D * *this );
	break;
      case Real_t :
	return ( D * *this );
	break;
      case Matrix_t :
	switch( D.sub_type() ) {
	case Rect_Matrix_t :
	  rVal = new dc_triple( ( ( dc_triple * )this )->get_as_row() *
				( ( dc_rect_matrix * )&D )->get() );
	  break;
	case Triple_t : /* cross multiplication */
	  rVal = new dc_rect_matrix( ( ( ( dc_triple * )this )->get() * 
				       ( ( dc_triple * )&D )->get_as_row()) );
	  break;
	case Vector_t : /* cross multiplication */
	  rVal = new dc_rect_matrix( ( ( ( dc_triple * )this )->get() * 
				       ( ( dc_vector * )&D )->get_as_row()) );
	  break;
	default :
	  cerr << "operator* -- operation not defined on " 
	       << type_string() << " and " << D.type_string() << "\n";
	  exit( 1 );
	}
	break;
      default :
	cerr << "operator* -- operation not defined on " 
	     << type_string() << " and " << D.type_string() << "\n";
	exit( 1 );
      }
      break;
    default :
      cerr << "operator* -- operation not defined on " 
	   << type_string() << " and " << D.type_string() << "\n";
      exit( 1 );
    }
    break;
  default :
    cerr << "operator* -- operation not defined on " 
	 << type_string() << " and " << D.type_string() << "\n";
    exit( 1 );
  }
  ( ( dc_udata * )rVal )->set_units( ( ( dc_udata * )this )->get_units() );
  ( ( dc_udata * )rVal )->get_units() *= ( ( dc_udata * )&D )->get_units();

  clean2;
  return *rVal;
}

dc_data &dc_data::dot( dc_data &D ) {
  dc_data *rVal = nil;
  switch( type() ) {
  case Int_t : case Real_t :
    return ( *this * D );
    break;
  case Matrix_t :
    switch( sub_type() ) {
    case Rect_Matrix_t :
      return ( *this * D );
      break;
    case Triple_t :
      if( D.sub_type() == Triple_t ) {
	rVal = new dc_real( ( ( *( dc_triple * )this )( 1 ) * 
			      ( *( dc_triple * )&D )( 1 ) ) + 
			    ( ( *( dc_triple * )this )( 2 ) * 
			      ( *( dc_triple * )&D )( 2 ) ) + 
			    ( ( *( dc_triple * )this )( 3 ) * 
			      ( *( dc_triple * )&D )( 3 ) ) );
      } else if( ( D.sub_type() == Vector_t ) && 
		 ( ( ( dc_vector * )&D )->dim() == 3 ) ) {
	rVal = new dc_real( ( ( *( dc_triple * )this )( 1 ) * 
			      ( *( dc_vector * )&D )( 1 ) ) + 
			    ( ( *( dc_triple * )this )( 2 ) * 
			      ( *( dc_vector * )&D )( 2 ) ) + 
			    ( ( *( dc_triple * )this )( 3 ) * 
			      ( *( dc_vector * )&D )( 3 ) ) );
      } else {
	cerr << "dot -- operation not defined on " 
	     << type_string() << " and " << D.type_string() << "\n";
	exit( 1 );
      }
      break;
    case Vector_t :
      if( D.sub_type() == Triple_t && ( ( dc_vector * )this )->dim() == 3 ) {
	rVal = new dc_real( ( ( *( dc_vector * )this )( 1 ) * 
			      ( *( dc_triple * )&D )( 1 ) ) + 
			    ( ( *( dc_vector * )this )( 2 ) * 
			      ( *( dc_triple * )&D )( 2 ) ) + 
			    ( ( *( dc_vector * )this )( 3 ) * 
			      ( *( dc_triple * )&D )( 3 ) ) );
      } else if( ( D.sub_type() == Vector_t ) && 
		 ( ( ( dc_vector * )&D )->dim() == 
		 ( ( dc_vector * )this )->dim() ) ) {
	Real r = 0;
	for( int i = 1 ; i <= ( ( dc_vector * )this )->dim() ; i++ ) {
	  r +=  ( *( dc_vector * )this )( i ) * ( *( dc_vector * )&D )( i );
	}
	rVal = new dc_real( r );
	break;
      } else {
	cerr << "dot -- operation not defined on " 
	     << type_string() << " and " << D.type_string() << "\n";
	exit( 1 );
      }
      break;
    default :
      cerr << "dot -- operation not defined on " 
	   << type_string() << " and " << D.type_string() << "\n";
      exit( 1 );
    }
    break;
  default :
    cerr << "dot -- operation not defined on " 
	 << type_string() << " and " << D.type_string() << "\n";
    exit( 1 );
  }

  ( ( dc_udata * )rVal )->set_units( ( ( dc_udata * )this )->get_units() );
  ( ( dc_udata * )rVal )->get_units() *= ( ( dc_udata * )&D )->get_units();
  clean2;
  return *rVal;
}

dc_data &dc_data::operator/( dc_data &D ) { /* division */
  dc_data *rVal;
  switch( type() ) { 
  case Int_t :
    switch( D.type() ) {
    case Int_t :
      rVal = new dc_int( ( ( dc_int * )this )->get() / 
			 ( ( dc_int * )&D )->get() );
      ( ( dc_int * )rVal )->set_units( ( ( dc_int * )this )->get_units() );
      ( ( dc_int * )rVal )->get_units() /= ( ( dc_int * )&D )->get_units();
      break;
    case Matrix_t :
      switch( D.sub_type() ) {
      case Rect_Matrix_t :
	if( !D.is_valid() ){
	  cerr << "operator/ -- operation not defined on " 
	       << type_string() << " and " << D.type_string() << "\n";
	  exit( 1 );
	}
	rVal = new dc_rect_matrix( ( ( dc_rect_matrix * )&D )->get() * 
				   ( 1.0 / ( double )( ( ( dc_int * )this )
						       ->get() ) ) ); 
	break;
      case Vector_t :
	if( !D.is_valid() ){
	  cerr << "operator/ -- operation not defined on " 
	       << type_string() << " and " << D.type_string() << "\n";
	  exit( 1 );
	}
	rVal = new dc_rect_matrix( ( ( dc_vector * )&D )->get() * 
				   ( 1.0 / ( double )( ( ( dc_int * )this )
						       ->get() ) ) ); 
	break;
      default :
	cerr << "operator/ -- operation not defined on " 
	     << type_string() << " and " << D.type_string() << "\n";
	exit( 1 );
      }
      break;
    case Real_t :
      rVal = new dc_real( ( ( dc_int * )this )->get() / 
			  ( ( dc_real * )&D )->get() );
      ( ( dc_real * )rVal )->set_units( ( ( dc_int * )this )->get_units() );
      ( ( dc_real * )rVal )->get_units() /= ( ( dc_real * )&D )->get_units();
      break;
    default :
      cerr << "operator/ -- operation not defined on " 
	   << type_string() << " and " << D.type_string() << "\n";
      exit( 1 );
    } 
    break;
  case Real_t :
    switch( D.type() ) {
    case Int_t :
      rVal = new dc_real( ( ( dc_real * )this )->get() / 
			  ( ( dc_int * )&D )->get() );
      ( ( dc_real * )rVal )->set_units( ( ( dc_real * )this )->get_units() );
      ( ( dc_real * )rVal )->get_units() /= ( ( dc_real * )&D )->get_units();
      break;
    case Matrix_t :
      switch( D.sub_type() ) {
      case Rect_Matrix_t :
	if( !D.is_valid() ){
	  cerr << "operator/ -- operation not defined on " 
	       << type_string() << " and " << D.type_string() << "\n";
	  exit( 1 );
	}
	rVal = new dc_rect_matrix( ( ( dc_rect_matrix * )&D )->get() *
				   ( 1.0 / ( ( dc_real * )this )->get() ) );
	break;
      case Vector_t :
	if( !D.is_valid() ){
	  cerr << "operator/ -- operation not defined on " 
	       << type_string() << " and " << D.type_string() << "\n";
	  exit( 1 );
	}
	rVal = new dc_rect_matrix( ( ( dc_vector * )&D )->get() * 
				   ( 1.0 / ( ( ( dc_real * )this )->get() ) ) ); 
	break;
      default :
	cerr << "operator/ -- operation not defined on " 
	     << type_string() << " and " << D.type_string() << "\n";
	exit( 1 );
      } 
      break;
    case Real_t :
      rVal = new dc_real( ( ( dc_real * )this )->get() / 
			  ( ( dc_real * )&D )->get() );
      ( ( dc_real * )this )->get_units() /= ( ( dc_real * )&D )->get_units();
      break;
    default :
      cerr << "operator/ -- operation not defined on " 
	   << type_string() << " and " << D.type_string() << "\n";
      exit( 1 );
    } 
    break;
  default :
    cerr << "operator/ -- operation not defined on " << type_string() 
	 << " and " << D.type_string() << "\n";
    exit( 1 );
  }
  
  clean2;
  return *rVal;
}

dc_data &dc_data::operator%( dc_data &D ) { /* modulus / dot product */
  dc_data *rVal;
  if( ( type() != Int_t ) || ( D.type() != Int_t ) ) {
    cerr << "operator% -- operation not defined on " << type_string() 
	 << " and " << D.type_string() << "\n";
    exit( 1 );
  }
  rVal = new dc_int( ( ( dc_int * )this )->get() % 
			 ( ( ( dc_int * )&D )->get() ) );
  ( ( dc_int * )rVal )->set_units( ( ( dc_int * )this )->get_units() );
  clean2;
  return *rVal;
}

/* arithmetic assignment operations. only work with non-temporaries */
dc_data &dc_data::operator+=( dc_data &D ) { /* add and assign */
  return ( *this = *this + D );
}

dc_data &dc_data::operator-=( dc_data &D ) { /* subtract and assign */ 
  return ( *this = *this - D );
}

dc_data &dc_data::operator*=( dc_data &D ) { /* multiply and assign */
  return ( *this = *this * D );
}

dc_data &dc_data::operator/=( dc_data &D ) { /* divide and assign */
  return ( *this = *this / D );
}

/* logical operations */
dc_data &dc_data::operator!( void ) { /* logical not */
  dc_data *rVal;
  
  if( type() != Boolean_t ) {
    cerr << "operator! -- operation not defined on " 
	 << type_string() << "\n";
    exit( 1 );
  } else {
    rVal = new dc_boolean( !( ( dc_boolean * )this )->get() );
  }

  clean;
  return *rVal;
}

dc_data &dc_data::operator&&( dc_data &D ) { /* logical and */
  dc_data *rVal;

  if( ( type() != Boolean_t ) || ( D.type() != Boolean_t ) ) {
    cerr << "operator&& -- operation not defined on " 
	 << type_string() << " and " << D.type_string() << "\n";
    exit( 1 );
  }
  rVal = new dc_boolean( ( ( dc_boolean * )this )->get() && 
			   ( ( dc_boolean * )&D )->get() );
  
  clean2;
  return *rVal;
}

dc_data &dc_data::operator||( dc_data &D ) { /* logical or */
  dc_data *rVal;

  if( ( type() != Boolean_t ) || ( D.type() != Boolean_t ) ) {
    cerr << "operator|| -- operation not defined on " 
	 << type_string() << " and " << D.type_string() << "\n";
    exit( 1 );
  }
  rVal = new dc_boolean( ( ( dc_boolean * )this )->get() 
			   || ( ( dc_boolean * )&D )->get() );

  clean2;
  return *rVal;
}

dc_data &dc_data::operator^( dc_data &D ) { /* logical xor */
  dc_data *rVal;

  if( ( type() != Boolean_t ) || ( D.type() != Boolean_t ) ) {
    cerr << "operator^ -- operation not defined on " 
	 << type_string() << " and " << D.type_string() << "\n";
    exit( 1 );
  }
  rVal = new dc_boolean( ( ( dc_boolean * )this )->get() != 
			   ( ( dc_boolean * )&D )->get() );

  clean2;
  return *rVal;
}

/* equality operations */
dc_data &dc_data::operator==( dc_data &D ) { /* equality */
  dc_data *rVal;

  if( !equiv_units( this, &D ) ) {
    rVal = new dc_boolean( false );
  } else {
    switch( type() ) {
    case Boolean_t :
      if( D.type() != Boolean_t ) {
	cerr << "operator== -- operation not defined on " 
	     << type_string() << " and " << D.type_string() << "\n";
	exit( 1 );
      }
      rVal = new dc_boolean( ( ( dc_boolean * )this )->get() == 
			     ( ( dc_boolean * )&D )->get() );
      break;
    case Int_t :
      switch( D.type() ) {
      case Int_t :
	rVal = new dc_boolean( ( ( double )( ( dc_int * )this )->get() ) == 
			       ( ( dc_int * )&D )->get() );
	break;
      case Real_t :
	rVal = new dc_boolean( ( ( dc_int * )this )->get() == 
			       ( ( dc_real * )&D )->get() );
	break;
      default :
	cerr << "operator== -- operation not defined on " 
	     << type_string() << " and " << D.type_string() << "\n";
	exit( 1 );
      }
      break;
    case Matrix_t :
      switch( sub_type() ) {
      case Rect_Matrix_t :
	if( D.sub_type() != Rect_Matrix_t ) {
	  cerr << "operator== -- operation not defined on " 
	       << type_string() << " and " << D.type_string() << "\n";
	  exit( 1 );
	}
	
	rVal = new dc_boolean( ( ( ( dc_rect_matrix * )this )->get() ) == 
			       ( ( ( dc_rect_matrix * )this )->get() ) ); 
	break;
      case Vector_t :
	if( D.sub_type() != Vector_t || 
	    !( is_valid() && ( D.is_valid() ) ) ){
	  cerr << "operator== -- operation not defined on " 
	       << type_string() << " and " << D.type_string() << "\n";
	  exit( 1 );
	}
	rVal = new dc_boolean( ( ( dc_vector * )this )->get() ==
			       ( ( dc_vector * )&D )->get() );
	break;
      default :
	cerr << "operator== -- operation not defined on " 
	     << type_string() << " and " << D.type_string() << "\n";
	exit( 1 );
	break;
      }
    case Pointer_t :
      if( D.type() != Pointer_t ) {
	cerr << "operator== -- operation not defined on " 
	     << type_string() << " and " << D.type_string() << "\n";
	exit( 1 );
      }
      rVal = new dc_boolean( ( ( dc_pointer * )this )->get() ==
			     ( ( dc_pointer * )&D )->get() );
      break;
    case Real_t :
      if( sub_type() == Distrib_t && D.sub_type() == Distrib_t ) {
	dc_distrib *d1 = ( dc_distrib * )this, *d2 = ( dc_distrib * )&D;
	rVal = new dc_boolean( d1->get_mean() == d2->get_mean() && 
			       d1->get_min() == d2->get_min() && 
			       d1->get_max() == d2->get_max() );
	break;
      }
      switch( D.type() ) {
      case Int_t :
	return ( D == *this );
	break;
      case Real_t :
	rVal = new dc_boolean( ( ( dc_real * )this )->get() == 
			       ( ( dc_real * )&D )->get() );
	break;
      default :
	cerr << "operator== -- operation not defined on " 
	     << type_string() << " and " << D.type_string() << "\n";
	exit( 1 );
      }
      break;
    case String_t :
      if( D.type() == String_t ) {
	rVal = new dc_boolean( ( *( ( dc_string * )this ) ) == 
			       ( *( ( dc_string * )&D ) ) );
	break;
      }
    case Symbol_t :
      if( D.type() == Symbol_t ) {
	rVal = new dc_boolean( ( ( dc_symbol * )this )->get() == 
			       ( ( dc_symbol * )&D )->get() );
	break;
      }
    default :
      cerr << "operator== -- operation not defined on " 
	   << type_string() << " and " << D.type_string() << "\n";
      exit( 1 );
      break;
    }
  }
  clean2;
  return *rVal;
}

dc_data &dc_data::operator!=( dc_data &D ) { /* inequality */
  dc_boolean *rVal = ( dc_boolean * )( &( *this == D ) );
  if( rVal ) {
    rVal->set( !rVal->get() );
  }
  return *rVal;
}

bool dc_data::lessthan( dc_data &D ) {
  bool b;

  switch( type() ) {
  case Int_t :
    switch( D.type() ) {
    case Int_t :
      b = ( ( dc_int * )this )->get() < ( ( dc_int * )&D )->get();
      break;
    case Real_t :
      b = ( ( dc_int * )this )->get() < ( ( dc_real * )&D )->get();
      break;
    default :
      cerr << "lessthan -- operation not defined on " 
	   << type_string() << " and " << D.type_string() << "\n";
      exit( 1 );
    }
    break;
  case Real_t :
    switch( D.type() ) {
    case Int_t :
      b = ( ( dc_real * )this )->get() < ( ( dc_int * )&D )->get();
      break;
    case Real_t :
      b = ( ( dc_real * )this )->get() < ( ( dc_real * )&D )->get();
      break;
    default :
      cerr << "lessthan -- operation not defined on " 
	   << type_string() << " and " << D.type_string() << "\n";
      exit( 1 );
    }
    break;
  case String_t :
    if( D.type() != String_t ) {
      cerr << "lessthan -- operation not defined on " 
	   << type_string() << " and " << D.type_string() << "\n";
      exit( 1 );
    }
    b = ( *( ( dc_string * )this ) ) < ( *( ( dc_string * )&D ) );
    break;
  case Symbol_t :
    if( D.type() != Symbol_t || 
	( ( dc_symbol * )this )->get_type() == undef_svec || 
	( ( dc_symbol * )this )->get_type() != 
	( ( dc_symbol * )&D )->get_type() ) {
      cerr << "lessthan -- operation not defined on " 
	   << type_string() << " and " << D.type_string() << "\n";
      exit( 1 );
    }
    b = ( ( dc_symbol * )this )->get_rank() < ( ( dc_symbol * )&D )->get_rank();
    break;
  default :
    cerr << "lessthan -- operation not defined on " 
	 << type_string() << " and " << D.type_string() << "\n";
    exit( 1 );
  }

  clean2;
  return b;
}

dc_data &dc_data::operator<( dc_data &D ) { /* less than */
  if( !check_units( this, &D, "lessthan" ) ) {
    exit( 1 );
  }

  return *( new dc_boolean( lessthan( D ) ) );
}

dc_data &dc_data::operator<=( dc_data &D ) { /* less than or equal to */
  if( !check_units( this, &D, "lessthan" ) ) {
    exit( 1 );
  }

  return *( new dc_boolean( !D.lessthan( *this ) ) );
}

dc_data &dc_data::operator>( dc_data &D ) { /* greater than */
  if( !check_units( this, &D, "lessthan" ) ) {
    exit( 1 );
  }

  return *( new dc_boolean( D.lessthan( *this ) ) );
}

dc_data &dc_data::operator>=( dc_data &D ) { /* greater than or equal to */
  if( !check_units( this, &D, "lessthan" ) ) {
    exit( 1 );
  }

  return *( new dc_boolean( !lessthan( D ) ) );
}

/* conditional operator since ?: cannot be overloaded */
/* does not require that type T = type E 
   perform cleanup on argument not used since no temps allocated */
dc_data &dc_data::If( dc_data &T, dc_data &E ) { /* conditional */
  if( type() != Boolean_t ) {
    cerr << "If -- First argument must be Bool\n\tFirst argument was " 
	 << full_type() << "\n";
    exit( 1 );
  }
  if( T.sub_type() != E.sub_type() ) {
    cerr << "If -- Warning : options not of same type\n\t"
	 << T.full_type() << " and " << E.full_type() << "\n";
  }
  bool i = ( ( dc_boolean * )this )->get();

  if( is_temporary() ) /* cleanup switch argument */
    delete( this );

  if( i ) {
    if( E.is_temporary() ) /* cleanup non-returned argument */
      delete( &E );
    return T;
  } else {
    if( T.is_temporary() ) /* cleanup non-returned argument */
      delete( &T );
    return E;
  }
}

/* magnitude operator */
dc_data &dc_data::mag( void ) { /* magnitude */
  dc_data *rVal = nil;
  int n; double d, e;
  switch( type() ){
  case Int_t :
    n = ( ( dc_int * )this )->get();
    rVal = new dc_int( n >=0 ? n : -n );
    break;
  case Matrix_t :
    switch( sub_type() ) {
    case Rect_Matrix_t :
      /* determinant of matrix */
      cerr << "DC_DATA::MAG -- NOT DEFINED YET. EXITING\n";
      exit( 1 );
      break;
    case Triple_t :
      rVal = new dc_real( sqrt( sqr( ( *( dc_triple * )this )( 1 ) ) + 
				sqr( ( * ( dc_triple * )this )( 2 ) ) +
				sqr( ( * ( dc_triple * )this )( 3 ) ) ) );
      break;
    case Vector_t :
      d = 0;
      for( int i = 0 ; i < ( ( dc_vector * )this )->dim() ; i++ ) {
	e = ( *( dc_vector * )this )(i);
	d += e * e;
      }
      rVal = new dc_real( sqrt( d ) );
      break;
    default :
      cerr << "mag -- operation not defined on " << type_string() << "\n";
      exit( 1 );
    }
    break;
  case Real_t :
    d = ( ( dc_real * )this )->get();
    rVal = new dc_real( d >=0 ? d : -d );
    break;
  case String_t :
    rVal = new dc_int( ( ( dc_string * )this )->length() );
    break;
  default :
    cerr << "mag -- operation not defined on " << type_string() << "\n";
    exit( 1 );
  }
  if( this->has_units() ) {
    ( ( dc_udata * )rVal )->set_units( ( ( dc_udata * )this )->get_units() );
  }
  clean;
  return *rVal;
}

dc_data &dc_data::operator-( void ) { /* negation */
  dc_data *rVal = nil;
  int n; double d;
  switch( type() ){
  case Int_t :
    n = ( ( dc_int * )this )->get();
    rVal = new dc_int( -n );
    ( ( dc_int * )rVal )->set_units( ( ( dc_int * )this )->get_units() );
     break;
  case Matrix_t :
    return( *this * *( new dc_real( -1 ) ) );
    break;
  case Real_t :
    d = ( ( dc_real * )this )->get();
    rVal = new dc_real( -d );
    ( ( dc_real * )rVal )->set_units( ( ( dc_real * )this )->get_units() );
    break;
  default :
    cerr << "operator- -- operation not defined on " << type_string() << "\n";
    exit( 1 );
  }
  ( ( dc_udata * )rVal )->set_units( ( ( dc_udata * )this )->get_units() );
  clean;
  return *rVal;
}

dc_data &dc_data::solve( dc_data &D ) {
  if( type() != Matrix_t || D.type() != Matrix_t ) {
    cerr << "solve -- operation not defined on " << type_string() << " and "
	 << D.type_string() << "\n";
    exit( 1 );
  }
  dc_data *rVal = nil;
  dc_matrix *m1 = ( dc_matrix * )this;
  dc_matrix *m2 = ( dc_matrix * )&D;
  rVal = new dc_rect_matrix( ( m1->get() ).i() * m2->get() );
  ( ( dc_udata * )rVal )->set_units( ( ( dc_udata * )&D )->get_units() );
  clean2;
  return *rVal;
}

/* returns a temporary and deletes d if it is temporary */
dc_data *cast( dc_data &d, const dc_type to ) {
  dc_type d_t = d.sub_type();
  if( to == d_t || to == Data_t ) return &d;

  dc_data *casted = nil;
  switch( d_t ) {
  case Boolean_t :
    break;
  case Distrib_t :
    if( to == Real_t ) {
      casted = new dc_real( ( ( dc_distrib * )&d )->get() );
    } else if( to == Int_t ) {
      casted = new dc_int( ( long int )( ( dc_distrib * )&d )->get() );
    }
    break;
  case Int_t :
    if( to == Real_t ) {
      casted = new dc_real( ( double )( ( dc_int * )&d )->get() );
    }
    break;
  case Real_t :
    if( to == Int_t ) {
      casted = new dc_int( ( long int )( ( dc_real * )&d )->get() );
    }
    break;
  case Rect_Matrix_t :
    if( ( ( dc_rect_matrix * )&d )->Nrows() == 1 ) {
      if( to == Vector_t ) {
	casted = new dc_vector( ( ( ( dc_rect_matrix * )&d )->get() )
				.AsColumn() );
      } else if( to == Triple_t && ( ( dc_rect_matrix * )&d )->Ncols() == 3 ) {
	casted = new dc_triple( ( ( ( dc_rect_matrix * )&d )->get() )
				.AsColumn() );
      }
    } else if( ( ( dc_rect_matrix * )&d )->Ncols() == 1 ) {
      if( to == Vector_t ) {
	casted = new dc_vector( ( ( ( dc_rect_matrix * )&d )->get_T() )
				.AsColumn() );
      } else if( to == Triple_t && ( ( dc_rect_matrix * )&d )->Nrows() == 3 ) {
	casted = new dc_triple( ( ( ( dc_rect_matrix * )&d )->get_T() )
				.AsColumn() );
      }
    }
    break;
  case Triple_t :
    if( to == Vector_t ) {
      casted = new dc_vector( ( ( dc_triple * )&d )->get() );
    } else if( to == Matrix_t || to == Rect_Matrix_t ) {
      casted = new dc_rect_matrix( ( ( dc_triple * )&d )->get() );
    }
    break;
  case Vector_t :
    if( to == Triple_t && ( ( dc_vector * )&d )->dim() == 3 ) {
      casted = new dc_triple( ( ( dc_vector * )&d )->get() );
    } else if( to == Matrix_t || to == Rect_Matrix_t ) {
      casted = new dc_rect_matrix( ( ( dc_vector * )&d )->get() );
    }
    break;
  default :;
  }
  if( casted == nil ) {
    dc_trace( TRACE_ERROR ) {
      cerr << "cast -- error trying to cast " << d.type_string() << " to "
	   << dc_type_string[ to ] << ".\n";
    }
    return nil;
  }

  if( d.has_units() )
    ( ( dc_udata * )casted )->set_units( ( ( dc_udata * )&d )->get_units() );
  if( d.is_temporary() )
    delete( &d );

  return casted;
}

bool castable( dc_type from, dc_type to ) {
  if( to == Data_t ) return true;
  switch( from ) {
  case Boolean_t :
    return to == Boolean_t;
  case Distrib_t :
    return to == Distrib_t || to == Int_t || to == Real_t;
  case Int_t : case Real_t : 
    return to == Int_t || to == Real_t;
  case Triple_t : case Vector_t :
    return to == Triple_t || to == Vector_t || to == Matrix_t || 
      to == Rect_Matrix_t;
  case Matrix_t : case Rect_Matrix_t :
    return to == Matrix_t || to == Rect_Matrix_t;
  default : /* Shape_t, String_t and Symbol_t not castable */
    return false;
  }
}

void dc_data::get_xtype( xtype &xt ) {
  xt.T = type();
  xt.subT = sub_type();
}

void dc_udata::get_xtype( xtype &xt ) {
  xt.T = type();
  xt.subT = sub_type();
  xt.uv = units;
}

#include "root.h"

bool set_default_units( dc_udata &d, cstring lbl, dc_node *origin ) {
  if( d.has_units() && !d.get_units() ) {
    dc_default *def = get_default( lbl, origin );
    dc_data *val;
    if( ( def != nil ) && ( ( val = def->get() ) != nil ) && val->has_units() )
    {
      d.set_units( ( ( dc_udata * )val )->get_units() );
      return true;
    }
  }
  return false;
}

dc_label *allocate( dc_type T, cstring label, dc_node *parent = nil ) {
  dc_label *new_data = allocate( T );
  if( new_data == nil ) {
    dc_trace( TRACE_ERROR ) {
      cerr << "allocate -- trying to allocate unknown type : " << ( int )T
	   << " : " << dc_type_string[T] << "\n";
    }
    exit( 1 );
  }
  new_data->set_both( label, parent );
  return new_data;
}

dc_label *allocate( dc_type T ) {
  dc_label *new_data;
  switch( T ){
  case Boolean_t :
    new_data = new dc_boolean();
    break;
  case Distrib_t :
    new_data = new dc_distrib();
    break;
  case Int_t :
    new_data = new dc_int();
    break;
  case Pointer_t :
    new_data = new dc_pointer();
    break;
  case Real_t :
    new_data = new dc_real();
    break;
  case Rect_Matrix_t :
    new_data = new dc_rect_matrix();
    break;
  case Set_t :
    new_data = new dc_set();
    break;
  case Shape_t :
    new_data = new dc_shape();
    break;
  case String_t :
    new_data = new dc_string();
    break;
  case Symbol_t :
    new_data = new dc_symbol();
    break;
  case Triple_t :
    new_data = new dc_triple();
    break;
  case Vector_t :
    new_data = new dc_vector();
    break;
  default :
    dc_trace( TRACE_ERROR ) {
      cerr << "allocate -- trying to allocate temporary of unknown type : " 
    	   << ( int )T << " : " << dc_type_string[T] << "\n";
    }
    exit( 1 );
  }
  return new_data;
}
