[BACK]Return to cplx.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib2 / asir2000 / engine

Diff for /OpenXM_contrib2/asir2000/engine/cplx.c between version 1.6 and 1.8

version 1.6, 2009/03/27 14:42:29 version 1.8, 2018/03/29 01:32:51
Line 45 
Line 45 
  * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,   * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
  * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.   * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
  *   *
  * $OpenXM: OpenXM_contrib2/asir2000/engine/cplx.c,v 1.5 2003/02/14 22:29:08 ohara Exp $   * $OpenXM: OpenXM_contrib2/asir2000/engine/cplx.c,v 1.7 2015/08/20 08:42:07 noro Exp $
 */  */
 #include "ca.h"  #include "ca.h"
 #include "base.h"  #include "base.h"
 #if defined(PARI)  
 #include "genpari.h"  
 #endif  
   
 void toreim(a,rp,ip)  void toreim(a,rp,ip)
 Num a;  Num a;
 Num *rp,*ip;  Num *rp,*ip;
 {  {
         if ( !a )    if ( !a )
                 *rp = *ip = 0;      *rp = *ip = 0;
 #if defined(INTERVAL)  #if defined(INTERVAL)
         else if ( NID(a) <= N_PRE_C ) {    else if ( NID(a) <= N_PRE_C ) {
 #else  #else
         else if ( NID(a) <= N_B ) {    else if ( NID(a) <= N_B ) {
 #endif  #endif
                 *rp = a; *ip = 0;      *rp = a; *ip = 0;
         } else {    } else {
                 *rp = ((C)a)->r; *ip = ((C)a)->i;      *rp = ((C)a)->r; *ip = ((C)a)->i;
         }    }
 }  }
   
 void reimtocplx(r,i,cp)  void reimtocplx(r,i,cp)
 Num r,i;  Num r,i;
 Num *cp;  Num *cp;
 {  {
         C c;    C c;
   
         if ( !i )    if ( !i )
                 *cp = r;      *cp = r;
         else {    else {
                 NEWC(c); c->r = r; c->i = i; *cp = (Num)c;      NEWC(c); c->r = r; c->i = i; *cp = (Num)c;
         }    }
 }  }
   
 void addcplx(a,b,c)  void addcplx(a,b,c)
 Num a,b;  Num a,b;
 Num *c;  Num *c;
 {  {
         Num ar,ai,br,bi,cr,ci;    Num ar,ai,br,bi,cr,ci;
   
         if ( !a )    if ( !a )
                 *c = b;      *c = b;
         else if ( !b )    else if ( !b )
                 *c = a;      *c = a;
 #if defined(INTERVAL)  #if defined(INTERVAL)
         else if ( (NID(a) <= N_PRE_C) && (NID(b) <= N_PRE_C ) )    else if ( (NID(a) <= N_PRE_C) && (NID(b) <= N_PRE_C ) )
 #else  #else
         else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )    else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
 #endif  #endif
                 addnum(0,a,b,c);      addnum(0,a,b,c);
         else {    else {
                 toreim(a,&ar,&ai); toreim(b,&br,&bi);      toreim(a,&ar,&ai); toreim(b,&br,&bi);
                 addnum(0,ar,br,&cr); addnum(0,ai,bi,&ci);      addnum(0,ar,br,&cr); addnum(0,ai,bi,&ci);
                 reimtocplx(cr,ci,c);      reimtocplx(cr,ci,c);
         }    }
 }  }
   
 void subcplx(a,b,c)  void subcplx(a,b,c)
 Num a,b;  Num a,b;
 Num *c;  Num *c;
 {  {
         Num ar,ai,br,bi,cr,ci;    Num ar,ai,br,bi,cr,ci;
   
         if ( !a )    if ( !a )
                 chsgnnum(b,c);      chsgnnum(b,c);
         else if ( !b )    else if ( !b )
                 *c = a;      *c = a;
 #if defined(INTERVAL)  #if defined(INTERVAL)
         else if ( (NID(a) <= N_PRE_C) && (NID(b) <= N_PRE_C ) )    else if ( (NID(a) <= N_PRE_C) && (NID(b) <= N_PRE_C ) )
 #else  #else
         else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )    else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
 #endif  #endif
                 subnum(0,a,b,c);      subnum(0,a,b,c);
         else {    else {
                 toreim(a,&ar,&ai); toreim(b,&br,&bi);      toreim(a,&ar,&ai); toreim(b,&br,&bi);
                 subnum(0,ar,br,&cr); subnum(0,ai,bi,&ci);      subnum(0,ar,br,&cr); subnum(0,ai,bi,&ci);
                 reimtocplx(cr,ci,c);      reimtocplx(cr,ci,c);
         }    }
 }  }
   
 void mulcplx(a,b,c)  void mulcplx(a,b,c)
 Num a,b;  Num a,b;
 Num *c;  Num *c;
 {  {
         Num ar,ai,br,bi,cr,ci,t,s;    Num ar,ai,br,bi,cr,ci,t,s;
   
         if ( !a || !b )    if ( !a || !b )
                 *c = 0;      *c = 0;
 #if defined(INTERVAL)  #if defined(INTERVAL)
         else if ( (NID(a) <= N_PRE_C) && (NID(b) <= N_PRE_C ) )    else if ( (NID(a) <= N_PRE_C) && (NID(b) <= N_PRE_C ) )
 #else  #else
         else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )    else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
 #endif  #endif
                 mulnum(0,a,b,c);      mulnum(0,a,b,c);
         else {    else {
                 toreim(a,&ar,&ai); toreim(b,&br,&bi);      toreim(a,&ar,&ai); toreim(b,&br,&bi);
                 mulnum(0,ar,br,&t); mulnum(0,ai,bi,&s); subnum(0,t,s,&cr);      mulnum(0,ar,br,&t); mulnum(0,ai,bi,&s); subnum(0,t,s,&cr);
                 mulnum(0,ar,bi,&t); mulnum(0,ai,br,&s); addnum(0,t,s,&ci);      mulnum(0,ar,bi,&t); mulnum(0,ai,br,&s); addnum(0,t,s,&ci);
                 reimtocplx(cr,ci,c);      reimtocplx(cr,ci,c);
         }    }
 }  }
   
 void divcplx(a,b,c)  void divcplx(a,b,c)
 Num a,b;  Num a,b;
 Num *c;  Num *c;
 {  {
         Num ar,ai,br,bi,cr,ci,t,s,u,w;    Num ar,ai,br,bi,cr,ci,t,s,u,w;
   
         if ( !b )    if ( !b )
                 error("divcplx : division by 0");      error("divcplx : division by 0");
         else if ( !a )    else if ( !a )
                 *c = 0;      *c = 0;
 #if defined(INTERVAL)  #if defined(INTERVAL)
         else if ( (NID(a) <= N_PRE_C) && (NID(b) <= N_PRE_C ) )    else if ( (NID(a) <= N_PRE_C) && (NID(b) <= N_PRE_C ) )
 #else  #else
         else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )    else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
 #endif  #endif
                 divnum(0,a,b,c);      divnum(0,a,b,c);
         else {    else {
                 toreim(a,&ar,&ai); toreim(b,&br,&bi);      toreim(a,&ar,&ai); toreim(b,&br,&bi);
                 mulnum(0,br,br,&t); mulnum(0,bi,bi,&s); addnum(0,t,s,&u);      mulnum(0,br,br,&t); mulnum(0,bi,bi,&s); addnum(0,t,s,&u);
                 mulnum(0,ar,br,&t); mulnum(0,ai,bi,&s);      mulnum(0,ar,br,&t); mulnum(0,ai,bi,&s);
                 addnum(0,t,s,&w); divnum(0,w,u,&cr);      addnum(0,t,s,&w); divnum(0,w,u,&cr);
                 mulnum(0,ar,bi,&t); mulnum(0,ai,br,&s);      mulnum(0,ar,bi,&t); mulnum(0,ai,br,&s);
                 subnum(0,s,t,&w); divnum(0,w,u,&ci);      subnum(0,s,t,&w); divnum(0,w,u,&ci);
                 reimtocplx(cr,ci,c);      reimtocplx(cr,ci,c);
         }    }
 }  }
   
 void pwrcplx(a,e,c)  void pwrcplx(a,e,c)
 Num a,e;  Num a,e;
 Num *c;  Num *c;
 {  {
         int ei;    int ei;
         Num t;    Num t;
   
         if ( !e )    if ( !e )
                 *c = (Num)ONE;      *c = (Num)ONE;
         else if ( !a )    else if ( !a )
                 *c = 0;      *c = 0;
         else if ( !INT(e) ) {    else if ( !INT(e) )
 #if defined(PARI)      error("pwrcplx : not implemented (use eval()).");
                 gpui_ri((Obj)a,(Obj)c,(Obj *)c);    else {
 #else      ei = SGN((Q)e)*QTOS((Q)e);
                 error("pwrcplx : can't calculate a fractional power");      pwrcplx0(a,ei,&t);
 #endif      if ( SGN((Q)e) < 0 )
         } else {        divnum(0,(Num)ONE,t,c);
                 ei = SGN((Q)e)*QTOS((Q)e);      else
                 pwrcplx0(a,ei,&t);        *c = t;
                 if ( SGN((Q)e) < 0 )    }
                         divnum(0,(Num)ONE,t,c);  
                 else  
                         *c = t;  
         }  
 }  }
   
 void pwrcplx0(a,e,c)  void pwrcplx0(a,e,c)
Line 210  Num a;
Line 203  Num a;
 int e;  int e;
 Num *c;  Num *c;
 {  {
         Num t,s;    Num t,s;
   
         if ( e == 1 )    if ( e == 1 )
                 *c = a;      *c = a;
         else {    else {
                 pwrcplx0(a,e/2,&t);      pwrcplx0(a,e/2,&t);
                 if ( !(e%2) )      if ( !(e%2) )
                         mulnum(0,t,t,c);        mulnum(0,t,t,c);
                 else {      else {
                         mulnum(0,t,t,&s); mulnum(0,s,a,c);        mulnum(0,t,t,&s); mulnum(0,s,a,c);
                 }      }
         }    }
 }  }
   
 void chsgncplx(a,c)  void chsgncplx(a,c)
 Num a,*c;  Num a,*c;
 {  {
         Num r,i;    Num r,i;
   
         if ( !a )    if ( !a )
                 *c = 0;      *c = 0;
 #if defined(INTERVAL)  #if defined(INTERVAL)
         else if ( NID(a) <= N_PRE_C )    else if ( NID(a) <= N_PRE_C )
 #else  #else
         else if ( NID(a) <= N_B )    else if ( NID(a) <= N_B )
 #endif  #endif
                 chsgnnum(a,c);      chsgnnum(a,c);
         else {    else {
                 chsgnnum(((C)a)->r,&r); chsgnnum(((C)a)->i,&i);      chsgnnum(((C)a)->r,&r); chsgnnum(((C)a)->i,&i);
                 reimtocplx(r,i,c);      reimtocplx(r,i,c);
         }    }
 }  }
   
 int cmpcplx(a,b)  int cmpcplx(a,b)
 Num a,b;  Num a,b;
 {  {
         Num ar,ai,br,bi;    Num ar,ai,br,bi;
         int s;    int s;
   
         if ( !a ) {    if ( !a ) {
 #if defined(INTERVAL)  #if defined(INTERVAL)
                 if ( !b || (NID(b)<=N_PRE_C) )      if ( !b || (NID(b)<=N_PRE_C) )
 #else  #else
                 if ( !b || (NID(b)<=N_B) )      if ( !b || (NID(b)<=N_B) )
 #endif  #endif
                         return compnum(0,a,b);        return compnum(0,a,b);
                 else      else
                         return -1;        return -1;
         } else if ( !b ) {    } else if ( !b ) {
 #if defined(INTERVAL)  #if defined(INTERVAL)
                 if ( !a || (NID(a)<=N_PRE_C) )      if ( !a || (NID(a)<=N_PRE_C) )
 #else  #else
                 if ( !a || (NID(a)<=N_B) )      if ( !a || (NID(a)<=N_B) )
 #endif  #endif
                         return compnum(0,a,b);        return compnum(0,a,b);
                 else      else
                         return 1;        return 1;
         } else {    } else {
                 toreim(a,&ar,&ai); toreim(b,&br,&bi);      toreim(a,&ar,&ai); toreim(b,&br,&bi);
                 s = compnum(0,ar,br);      s = compnum(0,ar,br);
                 if ( s )      if ( s )
                         return s;        return s;
                 else      else
                         return compnum(0,ai,bi);        return compnum(0,ai,bi);
         }    }
 }  }
Legend:
Removed from v.1.6  
changed lines
  Added in v.1.8
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>