FICO
FICO Xpress Optimization Examples Repository
FICO Optimization Community FICO Xpress Optimization Home
Back to examples browserPrevious exampleNext example

Definition of complex numbers and operators to work with them

Description
Language extensions provided by this module:
  • type: mathematical
  • operators: arithmetic, constructors, assignment, comparator
  • service: reset
Further explanation of this example: 'Mosel Native Interface User Guide', Chapter 5 Creating external types: second example


Source Files

Data Files





complex.c

/******************************************
  Mosel NI Examples
  =================
    
  File complex.c
  ``````````````
  Example module defining a new type
    complex
  with the corresponding arithmetic operators.

  (c) 2008 Fair Isaac Corporation
      author: Y. Colombani, 2002
*******************************************/

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include "xprm_ni.h"

#ifdef _WIN32
#define snprintf _snprintf
#endif

#define NCXL 20               /* Number of complex to allocate at once */

/**** Function prototypes ****/
static int cx_new0(XPRMcontext ctx,void *libctx);
static int cx_new1(XPRMcontext ctx,void *libctx);
static int cx_new2(XPRMcontext ctx,void *libctx);
static int cx_zero(XPRMcontext ctx,void *libctx);
static int cx_one(XPRMcontext ctx,void *libctx);
static int cx_asgn(XPRMcontext ctx,void *libctx);
static int cx_asgn_r(XPRMcontext ctx,void *libctx);
static int cx_pls(XPRMcontext ctx,void *libctx);
static int cx_pls_r(XPRMcontext ctx,void *libctx);
static int cx_neg(XPRMcontext ctx,void *libctx);
static int cx_mul(XPRMcontext ctx,void *libctx);
static int cx_mul_r(XPRMcontext ctx,void *libctx);
static int cx_div(XPRMcontext ctx,void *libctx);
static int cx_div_r1(XPRMcontext ctx,void *libctx);
static int cx_div_r2(XPRMcontext ctx,void *libctx);
static int cx_eql(XPRMcontext ctx,void *libctx);
static int cx_eql_r(XPRMcontext ctx,void *libctx);
static int cx_tostr(XPRMcontext ctx,void *,void *,char *,int,int);
static int cx_fromstr(XPRMcontext ctx,void *libctx,void *toinit,const char *str,int,const char **endptr);
static int cx_copy(XPRMcontext ctx,void *libctx,void *toinit,void *src,int typnum);
static int cx_compare(XPRMcontext ctx,void *libctx,void *c1,void *c2,int typnum);
static void *cx_create(XPRMcontext ctx,void *,void *,int);
static void cx_delete(XPRMcontext ctx,void *,void *,int);
static void *cx_reset(XPRMcontext ctx,void *libctx,int version);
static size_t cx_memuse(XPRMcontext ctx,void *cxctx,void *ref,int code);

/**** Structures for passing info to Mosel ****/
/* Subroutines */
static XPRMdsofct tabfct[]=
        {
         {"@&",1000,XPRM_TYP_EXTN,1,"complex:|complex|",cx_new0},
         {"@&",1001,XPRM_TYP_EXTN,1,"complex:r",cx_new1},
         {"@&",1002,XPRM_TYP_EXTN,2,"complex:rr",cx_new2},
         {"@0",1003,XPRM_TYP_EXTN,0,"complex:",cx_zero},
         {"@1",1004,XPRM_TYP_EXTN,0,"complex:",cx_one},
         {"@:",1005,XPRM_TYP_NOT,2,"|complex||complex|",cx_asgn},
         {"@:",1006,XPRM_TYP_NOT,2,"|complex|r",cx_asgn_r},
         {"@+",1007,XPRM_TYP_EXTN,2,"complex:|complex||complex|",cx_pls},
         {"@+",1008,XPRM_TYP_EXTN,2,"complex:|complex|r",cx_pls_r},
         {"@*",1009,XPRM_TYP_EXTN,2,"complex:|complex||complex|",cx_mul},
         {"@*",1010,XPRM_TYP_EXTN,2,"complex:|complex|r",cx_mul_r},
         {"@-",1011,XPRM_TYP_EXTN,1,"complex:|complex|",cx_neg},
         {"@/",1012,XPRM_TYP_EXTN,2,"complex:|complex||complex|",cx_div},
         {"@/",1013,XPRM_TYP_EXTN,2,"complex:|complex|r",cx_div_r1},
         {"@/",1014,XPRM_TYP_EXTN,2,"complex:r|complex|",cx_div_r2},
         {"@=",1015,XPRM_TYP_BOOL,2,"|complex||complex|",cx_eql},
         {"@=",1016,XPRM_TYP_BOOL,2,"|complex|r",cx_eql_r}
        };

/* Types */
static XPRMdsotyp tabtyp[]=
        {
         {"complex",1,XPRM_DTYP_PNCTX|XPRM_DTYP_RFCNT|XPRM_DTYP_APPND,cx_create,cx_delete,cx_tostr,cx_fromstr,cx_copy,cx_compare}
        };

/* Services */
static XPRMdsoserv tabserv[]=
        {
         {XPRM_SRV_RESET,(void *)cx_reset},
         {XPRM_SRV_MEMUSE,(void *)cx_memuse}
        };

/* Interface structure */
static XPRMdsointer dsointer= 
        { 
         0,NULL,
         sizeof(tabfct)/sizeof(XPRMdsofct),tabfct,
         sizeof(tabtyp)/sizeof(XPRMdsotyp),tabtyp,
         sizeof(tabserv)/sizeof(XPRMdsoserv),tabserv
        };

/**** Structures used by this module ****/
static XPRMnifct mm;             /* For storing Mosel NI function table */

typedef struct                  /* Where we store a complex number */
        {
         int refcnt;            /* For reference count */
         double re,im;
        } s_complex;

typedef union Freelist          /* List of allocated but not used complex */
        {
         s_complex cx;
         union Freelist *next;
        } u_freelist;
        
typedef struct Nmlist           /* A block of memory */
        {
         s_complex list[NCXL];
         int nextfree;
         struct Nmlist *next;
        } s_nmlist;

typedef struct                  /* A context for this module */
        {
         u_freelist *freelist;
         s_nmlist *nmlist;
        } s_cxctx;

/************************************************/
/* Initialize the library just after loading it */
/************************************************/
DSO_INIT complex_init(XPRMnifct nifct, int *interver,int *libver, XPRMdsointer **interf)
{
 mm=nifct;                      /* Save the table of Mosel NI functions */
 *interver=XPRM_NIVERS;         /* The interface version we are using */
 *libver=XPRM_MKVER(0,0,1);     /* The version of the module: 0.0.1 */
 *interf=&dsointer;             /* Our interface */

 return 0;
}

/******** Functions implementing the operators ********/

/*******************/
/* Clone a complex */
/*******************/
static int cx_new0(XPRMcontext ctx,void *libctx)
{
 s_complex *complex,*new_complex;

 complex=XPRM_POP_REF(ctx);
 if(complex!=NULL)
 {
  new_complex=cx_create(ctx,libctx,NULL,0);
  new_complex->im=complex->im;
  new_complex->re=complex->re;
  XPRM_PUSH_REF(ctx,new_complex);
 }
 else
  XPRM_PUSH_REF(ctx,NULL);
 return XPRM_RT_OK;
}

/********************************/
/* Create a complex from a real */
/********************************/
static int cx_new1(XPRMcontext ctx,void *libctx)
{
 s_complex *complex;

 complex=cx_create(ctx,libctx,NULL,0);
 complex->re=XPRM_POP_REAL(ctx);
 XPRM_PUSH_REF(ctx,complex);
 return XPRM_RT_OK;
}

/***********************************/
/* Create a complex from two reals */
/***********************************/
static int cx_new2(XPRMcontext ctx,void *libctx)
{
 s_complex *complex;

 complex=cx_create(ctx,libctx,NULL,0);
 complex->re=XPRM_POP_REAL(ctx);
 complex->im=XPRM_POP_REAL(ctx);
 XPRM_PUSH_REF(ctx,complex);
 return XPRM_RT_OK;
}

/**********************************************/
/* Zero in complex (used to initialise `sum') */
/**********************************************/
static int cx_zero(XPRMcontext ctx,void *libctx)
{
 XPRM_PUSH_REF(ctx,cx_create(ctx,libctx,NULL,0));
 return XPRM_RT_OK;
}

/**********************************************/
/* One in complex (used to initialise `prod') */
/**********************************************/
static int cx_one(XPRMcontext ctx,void *libctx)
{
 s_complex *complex;

 complex=cx_create(ctx,libctx,NULL,0);
 complex->re=1;
 XPRM_PUSH_REF(ctx,complex);
 return XPRM_RT_OK;
}

/*******************************/
/* Assignment complex:=complex */
/*******************************/
static int cx_asgn(XPRMcontext ctx,void *libctx)
{
 s_complex *c1,*c2;

 c1=XPRM_POP_REF(ctx);
 c2=XPRM_POP_REF(ctx);
 c1->im=c2->im;
 c1->re=c2->re;
 cx_delete(ctx,libctx,c2,0);
 return XPRM_RT_OK;
}

/****************************/
/* Assignment complex:=real */
/****************************/
static int cx_asgn_r(XPRMcontext ctx,void *libctx)
{
 s_complex *c1;

 c1=XPRM_POP_REF(ctx);
 c1->re=XPRM_POP_REAL(ctx);
 c1->im=0;
 return XPRM_RT_OK;
}

/***************************************/
/* Addition complex+complex -> complex */
/***************************************/
static int cx_pls(XPRMcontext ctx,void *libctx)
{
 s_complex *c1,*c2;

 c1=XPRM_POP_REF(ctx);
 c2=XPRM_POP_REF(ctx);
 if(c1!=NULL)
 {
  if(c2!=NULL)
  {
   c1->re+=c2->re;
   c1->im+=c2->im;
   cx_delete(ctx,libctx,c2,0);
  }
  XPRM_PUSH_REF(ctx,c1);
 }
 else
  XPRM_PUSH_REF(ctx,c2);
 return XPRM_RT_OK;
}

/************************************/
/* Addition complex+real -> complex */
/************************************/
static int cx_pls_r(XPRMcontext ctx,void *libctx)
{
 s_complex *c1;

 c1=XPRM_POP_REF(ctx);
 if(c1!=NULL)
  c1->re+=XPRM_POP_REAL(ctx);
 XPRM_PUSH_REF(ctx,c1);
 return XPRM_RT_OK;
}

/**************************************/
/* Change of sign complex -> -complex */
/**************************************/
static int cx_neg(XPRMcontext ctx,void *libctx)
{
 s_complex *c1;

 c1=XPRM_POP_REF(ctx);
 if(c1!=NULL)
 {
  c1->re=-c1->re;
  c1->im=-c1->im;
 }
 XPRM_PUSH_REF(ctx,c1);
 return XPRM_RT_OK;
}

/**************************************/
/* Product complex*complex -> complex */
/**************************************/
static int cx_mul(XPRMcontext ctx,void *libctx)
{
 s_complex *c1,*c2;
 double re,im;

 c1=XPRM_POP_REF(ctx);
 c2=XPRM_POP_REF(ctx);
 if(c1!=NULL)
 {
  if(c2!=NULL)
  {
   re=c1->re*c2->re-c1->im*c2->im;
   im=c1->re*c2->im+c1->im*c2->re;
   c1->re=re;
   c1->im=im;
  }
  else
   c1->re=c2->im=0;
 }
 cx_delete(ctx,libctx,c2,0);
 XPRM_PUSH_REF(ctx,c1);
 return XPRM_RT_OK;
}

/***********************************/
/* Product complex*real -> complex */
/***********************************/
static int cx_mul_r(XPRMcontext ctx,void *libctx)
{
 s_complex *c1;
 double r;

 c1=XPRM_POP_REF(ctx);
 r=XPRM_POP_REAL(ctx);
 if(c1!=NULL)
 {
  c1->re*=r;
  c1->im*=r;
 }
 XPRM_PUSH_REF(ctx,c1);
 return XPRM_RT_OK;
}

/***************************************/
/* Division complex/complex -> complex */
/***************************************/
static int cx_div(XPRMcontext ctx,void *libctx)
{
 s_complex *c1,*c2;
 double re,im;

 c1=XPRM_POP_REF(ctx);
 c2=XPRM_POP_REF(ctx);
 if((c2==NULL)||((c2->re==0)&&(c2->im==0)))
 {
  mm->dispmsg(ctx,"Complex: Division by 0.");
  return XPRM_RT_ERROR;
 }
 else
 {
  if(c1!=NULL)
  {                             /* Compute 1/c2 then c1* 1/c2 */
   re=c2->re/(c2->re*c2->re+c2->im*c2->im);
   im=-c2->im/(c2->re*c2->re+c2->im*c2->im);
   c2->re=re;
   c2->im=im;
   XPRM_PUSH_REF(ctx,c2);
   XPRM_PUSH_REF(ctx,c1);
   return cx_mul(ctx,libctx);
  }
  else
  {
   cx_delete(ctx,libctx,c2,0);
   XPRM_PUSH_REF(ctx,c1);
   return XPRM_RT_OK;
  }
 }
}

/************************************/
/* Division complex/real -> complex */
/************************************/
static int cx_div_r1(XPRMcontext ctx,void *libctx)
{
 s_complex *c1;
 double r;

 c1=XPRM_POP_REF(ctx);
 r=XPRM_POP_REAL(ctx);
 if(r==0)
 {
  mm->dispmsg(ctx,"Complex: Division by 0.");
  return XPRM_RT_ERROR;
 }
 else
 {
  if(c1!=NULL)
  {
   c1->re/=r;
   c1->im/=r;
  }
  XPRM_PUSH_REF(ctx,c1);
  return XPRM_RT_OK;
 }
}

/************************************/
/* Division real/complex -> complex */
/************************************/
static int cx_div_r2(XPRMcontext ctx,void *libctx)
{
 s_complex *c1;
 double r,re,im;

 r=XPRM_POP_REAL(ctx);
 c1=XPRM_POP_REF(ctx);
 if((c1==NULL)||((c1->re==0)&&(c1->im==0)))
 {
  mm->dispmsg(ctx,"Complex: Division by 0.");
  return XPRM_RT_ERROR;
 }
 else
 {
  re=(c1->re*r)/(c1->re*c1->re+c1->im*c1->im);
  im=-(c1->im*r)/(c1->re*c1->re+c1->im*c1->im);
  c1->re=re;
  c1->im=im;
  XPRM_PUSH_REF(ctx,c1);
  return XPRM_RT_OK;
 }
}

/******************************/
/* Comparison complex=complex */
/******************************/
static int cx_eql(XPRMcontext ctx,void *libctx)
{
 s_complex *c1,*c2;
 int b;

 c1=XPRM_POP_REF(ctx);
 c2=XPRM_POP_REF(ctx);
 if(c1!=NULL)
 {
  if(c2!=NULL)
   b=(c1->re==c2->re)&&(c1->im==c2->im);
  else
   b=0;
 }
 else
  b=(c2==NULL);
 XPRM_PUSH_INT(ctx,b);
 return XPRM_RT_OK;
}

/***************************/
/* Comparison complex=real */
/***************************/
static int cx_eql_r(XPRMcontext ctx,void *libctx)
{
 s_complex *c1;
 double r;
 int b;

 c1=XPRM_POP_REF(ctx);
 r=XPRM_POP_REAL(ctx);
 if(c1!=NULL)
  b=(c1->im==0)&&(c1->re==r);
 else
  b=(r==0);
 XPRM_PUSH_INT(ctx,b);
 return XPRM_RT_OK;
}

/**************** Type-related functions ****************/

/*****************************/
/* Allocate a complex number */
/*****************************/
static void *cx_create(XPRMcontext ctx,void *libctx,void *todup,int typnum)
{
 s_cxctx *cxctx;
 s_complex *complex;
 s_nmlist *nmlist;

 if(todup!=NULL)
 {
  ((s_complex *)todup)->refcnt++;
  return todup;
 }
 else
 {
  cxctx=libctx;
  if(cxctx->freelist!=NULL)              /* We have got some free complex */
  {
   complex=&(cxctx->freelist->cx);
   cxctx->freelist=cxctx->freelist->next;
  }
  else                                   /* We must allocate a new block */
   if((cxctx->nmlist==NULL)||(cxctx->nmlist->nextfree>=NCXL))
   {
    nmlist=malloc(sizeof(s_nmlist));
    nmlist->next=cxctx->nmlist;
    cxctx->nmlist=nmlist;
    nmlist->nextfree=1;
    complex=nmlist->list;
   }
   else                                  /* We can take one from the block */
    complex=&(cxctx->nmlist->list[cxctx->nmlist->nextfree++]);
  complex->re=complex->im=0;
  complex->refcnt=1;
  return complex;
 }
}

/*******************************/
/* Deallocate a complex number */
/*******************************/
static void cx_delete(XPRMcontext ctx,void *libctx,void *todel,int typnum)
{
 s_cxctx *cxctx;
 u_freelist *freelist;

 if((todel!=NULL)&&((--((s_complex *)todel)->refcnt)<1))
 {
  cxctx=libctx;
  freelist=todel;
  freelist->next=cxctx->freelist;
  cxctx->freelist=freelist;
 }
}

/*********************/
/* Complex -> String */
/*********************/
static int cx_tostr(XPRMcontext ctx,void *libctx,void *toprt,char *str,int len,int typnum)
{
 s_complex *c;

 if(toprt==NULL)
 {
  strncpy(str,"0+0i",len);
  return 4;
 }
 else
 {
  c=toprt;
  return snprintf(str,len,"%g%+gi",c->re,c->im);
 }
}

/*********************/
/* String -> Complex */
/*********************/
static int cx_fromstr(XPRMcontext ctx,void *libctx,void *toinit,const char *str,int typnum,const char **endptr)
{
 double re,im;
 s_complex *c;
 int len;
 struct Info
 {
  char dummy[4];
  s_complex *c;
 } *ref;

 c=toinit;
 if((str[0]=='r') && (str[1]=='a') && (str[2]=='w') && (str[3]=='\0'))
 {
  if(endptr!=NULL) *endptr=NULL;
  ref=(struct Info *)str;
  if(ref->c==NULL)
   c->re=c->im=0;
  else
  {
   c->re=ref->c->re;
   c->im=ref->c->im;
  }
  return XPRM_RT_OK;
 }
 else
  if(sscanf(str,"%lf%lf%ni%n",&re,&im,&len,&len)<2)
  {
   if(endptr!=NULL) *endptr=str;
   return XPRM_RT_ERROR;
  }
  else
  {
   if(endptr!=NULL) *endptr=str+len;
   c->re=re;
   c->im=im;
   return XPRM_RT_OK;
  }
}

/******************/
/* Copy a complex */
/******************/
static int cx_copy(XPRMcontext ctx,void *libctx,void *toinit,void *src,int typnum)
{
 s_complex *c_dst;

 c_dst=(s_complex *)toinit;
 if(XPRM_CPY(typnum)==XPRM_CPY_APPEND)
 {
  if(src!=NULL)
  {
   c_dst->re+=((s_complex *)src)->re;
   c_dst->im+=((s_complex *)src)->im;
  }
 }
 else
  if(src!=NULL)
  {
   c_dst->re=((s_complex *)src)->re;
   c_dst->im=((s_complex *)src)->im;
  }
  else
  {
   c_dst->re=0;
   c_dst->im=0;
  }
 return 0;
}

/****************************/
/* Compare 2 complex values */
/****************************/
static int cx_compare(XPRMcontext ctx,void *libctx,void *c1,void *c2,int typnum)
{
 int b;

 if(c1!=NULL)
 {
  if(c2!=NULL)
   b=(((s_complex *)c1)->re==((s_complex *)c2)->re)&&
     (((s_complex *)c1)->im==((s_complex *)c2)->im);
  else
   b=((((s_complex *)c1)->re==0)&&(((s_complex *)c1)->im==0));
 }
 else
  b=(c2==NULL)||((((s_complex *)c2)->re==0)&&(((s_complex *)c2)->im==0));
 
 switch(XPRM_COMPARE(typnum))
 {
  case MM_COMPARE_EQ:
    return b;
  case MM_COMPARE_NEQ:
    return !b;
  default:
    return XPRM_COMPARE_ERROR;
 }
}

/******************** Services ********************/

/**************************************/
/* Reset the Complex module for a run */
/**************************************/
static void *cx_reset(XPRMcontext ctx,void *libctx,int version)
{
 s_cxctx *cxctx;
 s_nmlist *nmlist;

 if(libctx==NULL)               /* libctx==NULL => initialisation */
 {
  cxctx=malloc(sizeof(s_cxctx));
  memset(cxctx,0,sizeof(s_cxctx));
  return cxctx;
 }
 else                           /* otherwise release the resources we use */
 {
  cxctx=libctx;
  while(cxctx->nmlist!=NULL)
  {
   nmlist=cxctx->nmlist;
   cxctx->nmlist=nmlist->next;
   free(nmlist);
  }
  free(cxctx);
  return NULL;
 }
}

/*****************************/
/* Memory used by the module */
/*****************************/
static size_t cx_memuse(XPRMcontext ctx,void *cxctx,void *ref,int code)
{
 switch(code)
 {
  case 0:
    {
     size_t s;
     s_nmlist *nmlist;
     
     s=sizeof(s_cxctx);
     nmlist=((s_cxctx*)cxctx)->nmlist;
     while(nmlist!=NULL)
     {
      s+=sizeof(s_nmlist);
      nmlist=nmlist->next;
     }
     return s;
    }
  case 1:
    return sizeof(s_complex);
  default:
    return -1;
 }
}

Back to examples browserPrevious exampleNext example