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

version 1.2, 2000/08/21 08:31:27 version 1.9, 2019/06/04 07:11:22
Line 23 
Line 23 
  * shall be made on your publication or presentation in any form of the   * shall be made on your publication or presentation in any form of the
  * results obtained by use of the SOFTWARE.   * results obtained by use of the SOFTWARE.
  * (4) In the event that you modify the SOFTWARE, you shall notify FLL by   * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
  * e-mail at risa-admin@flab.fujitsu.co.jp of the detailed specification   * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
  * for such modification or the source code of the modified part of the   * for such modification or the source code of the modified part of the
  * SOFTWARE.   * SOFTWARE.
  *   *
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.1.1.1 1999/12/03 07:39:08 noro 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 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;
         else if ( NID(a) <= N_B ) {  #if defined(INTERVAL)
                 *rp = a; *ip = 0;    else if ( NID(a) < N_C ) {
         } else {  #else
                 *rp = ((C)a)->r; *ip = ((C)a)->i;    else if ( NID(a) <= N_B ) {
         }  #endif
       *rp = a; *ip = 0;
     } else {
       *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;
         else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )  #if defined(INTERVAL)
                 addnum(0,a,b,c);    else if ( (NID(a) < N_C) && (NID(b) < N_C ) )
         else {  #else
                 toreim(a,&ar,&ai); toreim(b,&br,&bi);    else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
                 addnum(0,ar,br,&cr); addnum(0,ai,bi,&ci);  #endif
                 reimtocplx(cr,ci,c);      addnum(0,a,b,c);
         }    else {
       toreim(a,&ar,&ai); toreim(b,&br,&bi);
       addnum(0,ar,br,&cr); addnum(0,ai,bi,&ci);
       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;
         else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )  #if defined(INTERVAL)
                 subnum(0,a,b,c);    else if ( (NID(a) < N_C) && (NID(b) <= N_C ) )
         else {  #else
                 toreim(a,&ar,&ai); toreim(b,&br,&bi);    else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
                 subnum(0,ar,br,&cr); subnum(0,ai,bi,&ci);  #endif
                 reimtocplx(cr,ci,c);      subnum(0,a,b,c);
         }    else {
       toreim(a,&ar,&ai); toreim(b,&br,&bi);
       subnum(0,ar,br,&cr); subnum(0,ai,bi,&ci);
       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;
         else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )  #if defined(INTERVAL)
                 mulnum(0,a,b,c);    else if ( (NID(a) < N_C) && (NID(b) < N_C ) )
         else {  #else
                 toreim(a,&ar,&ai); toreim(b,&br,&bi);    else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
                 mulnum(0,ar,br,&t); mulnum(0,ai,bi,&s); subnum(0,t,s,&cr);  #endif
                 mulnum(0,ar,bi,&t); mulnum(0,ai,br,&s); addnum(0,t,s,&ci);      mulnum(0,a,b,c);
                 reimtocplx(cr,ci,c);    else {
         }      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,bi,&t); mulnum(0,ai,br,&s); addnum(0,t,s,&ci);
       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;
         else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )  #if defined(INTERVAL)
                 divnum(0,a,b,c);    else if ( (NID(a) < N_C) && (NID(b) < N_C ) )
         else {  #else
                 toreim(a,&ar,&ai); toreim(b,&br,&bi);    else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
                 mulnum(0,br,br,&t); mulnum(0,bi,bi,&s); addnum(0,t,s,&u);  #endif
                 mulnum(0,ar,br,&t); mulnum(0,ai,bi,&s);      divnum(0,a,b,c);
                 addnum(0,t,s,&w); divnum(0,w,u,&cr);    else {
                 mulnum(0,ar,bi,&t); mulnum(0,ai,br,&s);      toreim(a,&ar,&ai); toreim(b,&br,&bi);
                 subnum(0,s,t,&w); divnum(0,w,u,&ci);      mulnum(0,br,br,&t); mulnum(0,bi,bi,&s); addnum(0,t,s,&u);
                 reimtocplx(cr,ci,c);      mulnum(0,ar,br,&t); mulnum(0,ai,bi,&s);
         }      addnum(0,t,s,&w); divnum(0,w,u,&cr);
       mulnum(0,ar,bi,&t); mulnum(0,ai,br,&s);
       subnum(0,s,t,&w); divnum(0,w,u,&ci);
       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 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 200  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;
         else if ( NID(a) <= N_B )  #if defined(INTERVAL)
                 chsgnnum(a,c);    else if ( NID(a) < N_C )
         else {  #else
                 chsgnnum(((C)a)->r,&r); chsgnnum(((C)a)->i,&i);    else if ( NID(a) <= N_B )
                 reimtocplx(r,i,c);  #endif
         }      chsgnnum(a,c);
     else {
       chsgnnum(((C)a)->r,&r); chsgnnum(((C)a)->i,&i);
       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 ( !b || (NID(b)<=N_B) )  #if defined(INTERVAL)
                         return compnum(0,a,b);      if ( !b || (NID(b)<N_C) )
                 else  #else
                         return -1;      if ( !b || (NID(b)<=N_B) )
         } else if ( !b ) {  #endif
                 if ( !a || (NID(a)<=N_B) )        return compnum(0,a,b);
                         return compnum(0,a,b);      else
                 else        return -1;
                         return 1;    } else if ( !b ) {
         } else {  #if defined(INTERVAL)
                 toreim(a,&ar,&ai); toreim(b,&br,&bi);      if ( !a || (NID(a)<N_C) )
                 s = compnum(0,ar,br);  #else
                 if ( s )      if ( !a || (NID(a)<=N_B) )
                         return s;  #endif
                 else        return compnum(0,a,b);
                         return compnum(0,ai,bi);      else
         }        return 1;
     } else {
       toreim(a,&ar,&ai); toreim(b,&br,&bi);
       s = compnum(0,ar,br);
       if ( s )
         return s;
       else
         return compnum(0,ai,bi);
     }
 }  }
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.9
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>