=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/cplx.c,v retrieving revision 1.1.1.1 retrieving revision 1.9 diff -u -p -r1.1.1.1 -r1.9 --- OpenXM_contrib2/asir2000/engine/cplx.c 1999/12/03 07:39:08 1.1.1.1 +++ OpenXM_contrib2/asir2000/engine/cplx.c 2019/06/04 07:11:22 1.9 @@ -1,150 +1,201 @@ -/* $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.8 2018/03/29 01:32:51 noro Exp $ +*/ #include "ca.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) Num a; Num *rp,*ip; { - if ( !a ) - *rp = *ip = 0; - else if ( NID(a) <= N_B ) { - *rp = a; *ip = 0; - } else { - *rp = ((C)a)->r; *ip = ((C)a)->i; - } + if ( !a ) + *rp = *ip = 0; +#if defined(INTERVAL) + else if ( NID(a) < N_C ) { +#else + 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) Num r,i; Num *cp; { - C c; + C c; - if ( !i ) - *cp = r; - else { - NEWC(c); c->r = r; c->i = i; *cp = (Num)c; - } + if ( !i ) + *cp = r; + else { + NEWC(c); c->r = r; c->i = i; *cp = (Num)c; + } } void addcplx(a,b,c) Num a,b; Num *c; { - Num ar,ai,br,bi,cr,ci; + Num ar,ai,br,bi,cr,ci; - if ( !a ) - *c = b; - else if ( !b ) - *c = a; - else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) ) - 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); - } + if ( !a ) + *c = b; + else if ( !b ) + *c = a; +#if defined(INTERVAL) + else if ( (NID(a) < N_C) && (NID(b) < N_C ) ) +#else + else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) ) +#endif + 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) Num a,b; Num *c; { - Num ar,ai,br,bi,cr,ci; + Num ar,ai,br,bi,cr,ci; - if ( !a ) - chsgnnum(b,c); - else if ( !b ) - *c = a; - else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) ) - 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); - } + if ( !a ) + chsgnnum(b,c); + else if ( !b ) + *c = a; +#if defined(INTERVAL) + else if ( (NID(a) < N_C) && (NID(b) <= N_C ) ) +#else + else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) ) +#endif + 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) Num a,b; Num *c; { - Num ar,ai,br,bi,cr,ci,t,s; + Num ar,ai,br,bi,cr,ci,t,s; - if ( !a || !b ) - *c = 0; - else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) ) - mulnum(0,a,b,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); - } + if ( !a || !b ) + *c = 0; +#if defined(INTERVAL) + else if ( (NID(a) < N_C) && (NID(b) < N_C ) ) +#else + else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) ) +#endif + mulnum(0,a,b,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) Num a,b; 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 ) - error("divcplx : division by 0"); - else if ( !a ) - *c = 0; - else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) ) - divnum(0,a,b,c); - else { - 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,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); - } + if ( !b ) + error("divcplx : division by 0"); + else if ( !a ) + *c = 0; +#if defined(INTERVAL) + else if ( (NID(a) < N_C) && (NID(b) < N_C ) ) +#else + else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) ) +#endif + divnum(0,a,b,c); + else { + 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,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) Num a,e; Num *c; { - int ei; - Num t; - extern long prec; + int ei; + Num t; - if ( !e ) - *c = (Num)ONE; - else if ( !a ) - *c = 0; - else if ( !INT(e) ) { -#if PARI - GEN pa,pe,z; - int ltop,lbot; - - ltop = avma; ritopa((Obj)a,&pa); ritopa((Obj)e,&pe); lbot = avma; - z = gerepile(ltop,lbot,gpui(pa,pe,prec)); - patori(z,(Obj *)c); cgiv(z); -#else - 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; - } + if ( !e ) + *c = (Num)ONE; + else if ( !a ) + *c = 0; + else if ( !INT(e) ) + error("pwrcplx : not implemented (use eval())."); + 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) @@ -152,57 +203,69 @@ Num a; int e; Num *c; { - Num t,s; + Num t,s; - if ( e == 1 ) - *c = a; - else { - pwrcplx0(a,e/2,&t); - if ( !(e%2) ) - mulnum(0,t,t,c); - else { - mulnum(0,t,t,&s); mulnum(0,s,a,c); - } - } + if ( e == 1 ) + *c = a; + else { + pwrcplx0(a,e/2,&t); + if ( !(e%2) ) + mulnum(0,t,t,c); + else { + mulnum(0,t,t,&s); mulnum(0,s,a,c); + } + } } void chsgncplx(a,c) Num a,*c; { - Num r,i; + Num r,i; - if ( !a ) - *c = 0; - else if ( NID(a) <= N_B ) - chsgnnum(a,c); - else { - chsgnnum(((C)a)->r,&r); chsgnnum(((C)a)->i,&i); - reimtocplx(r,i,c); - } + if ( !a ) + *c = 0; +#if defined(INTERVAL) + else if ( NID(a) < N_C ) +#else + else if ( NID(a) <= N_B ) +#endif + chsgnnum(a,c); + else { + chsgnnum(((C)a)->r,&r); chsgnnum(((C)a)->i,&i); + reimtocplx(r,i,c); + } } int cmpcplx(a,b) Num a,b; { - Num ar,ai,br,bi; - int s; + Num ar,ai,br,bi; + int s; - if ( !a ) { - if ( !b || (NID(b)<=N_B) ) - return compnum(0,a,b); - else - return -1; - } else if ( !b ) { - if ( !a || (NID(a)<=N_B) ) - return compnum(0,a,b); - 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); - } + if ( !a ) { +#if defined(INTERVAL) + if ( !b || (NID(b)