[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.1.1.1 and 1.8

version 1.1.1.1, 1999/12/03 07:39:08 version 1.8, 2018/03/29 01:32:51
Line 1 
Line 1 
 /* $OpenXM: OpenXM/src/asir99/engine/cplx.c,v 1.1.1.1 1999/11/10 08:12:26 noro Exp $ */  /*
    * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
    * All rights reserved.
    *
    * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
    * non-exclusive and royalty-free license to use, copy, modify and
    * redistribute, solely for non-commercial and non-profit purposes, the
    * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
    * conditions of this Agreement. For the avoidance of doubt, you acquire
    * only a limited right to use the SOFTWARE hereunder, and FLL or any
    * third party developer retains all rights, including but not limited to
    * copyrights, in and to the SOFTWARE.
    *
    * (1) FLL does not grant you a license in any way for commercial
    * purposes. You may use the SOFTWARE only for non-commercial and
    * non-profit purposes only, such as academic, research and internal
    * business use.
    * (2) The SOFTWARE is protected by the Copyright Law of Japan and
    * international copyright treaties. If you make copies of the SOFTWARE,
    * with or without modification, as permitted hereunder, you shall affix
    * to all such copies of the SOFTWARE the above copyright notice.
    * (3) An explicit reference to this SOFTWARE and its copyright owner
    * shall be made on your publication or presentation in any form of the
    * results obtained by use of the SOFTWARE.
    * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
    * 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
    * SOFTWARE.
    *
    * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
    * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
    * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
    * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
    * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
    * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
    * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
    * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
    * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
    * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
    * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
    * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
    * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
    * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
    * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
    * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
    * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
    *
    * $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 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_PRE_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_PRE_C) && (NID(b) <= N_PRE_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_PRE_C) && (NID(b) <= N_PRE_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_PRE_C) && (NID(b) <= N_PRE_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_PRE_C) && (NID(b) <= N_PRE_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 152  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_PRE_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_PRE_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_PRE_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.1.1.1  
changed lines
  Added in v.1.8
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>