[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.5 and 1.9

version 1.5, 2003/02/14 22:29:08 version 1.9, 2019/06/04 07:11:22
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.4 2000/12/22 10:03:28 saito Exp $   * $OpenXM: OpenXM_contrib2/asir2000/engine/cplx.c,v 1.8 2018/03/29 01:32:51 noro Exp $
 */  */
 #include "ca.h"  #include "ca.h"
 #include "base.h"  #include "base.h"
 #if defined(PARI)  
 #include "genpari.h"  
 void patori(GEN,Obj *);  
 void patori_i(GEN,N *);  
 void ritopa(Obj,GEN *);  
 void ritopa_i(N,int,GEN *);  
 #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_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_C) && (NID(b) < N_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_C) && (NID(b) <= N_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_C) && (NID(b) < N_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_C) && (NID(b) < N_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;
         extern long prec;  
   
         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()).");
                 GEN pa,pe,z;    else {
                 int ltop,lbot;      ei = SGN((Q)e)*QTOS((Q)e);
       pwrcplx0(a,ei,&t);
                 ltop = avma; ritopa((Obj)a,&pa); ritopa((Obj)e,&pe); lbot = avma;      if ( SGN((Q)e) < 0 )
                 z = gerepile(ltop,lbot,gpui(pa,pe,prec));        divnum(0,(Num)ONE,t,c);
                 patori(z,(Obj *)c); cgiv(z);      else
 #else        *c = t;
                 error("pwrcplx : can't calculate a fractional power");    }
 #endif  
         } else {  
                 ei = SGN((Q)e)*QTOS((Q)e);  
                 pwrcplx0(a,ei,&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 220  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_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_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_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.5  
changed lines
  Added in v.1.9
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>