Annotation of OpenXM_contrib2/asir2000/engine/cplx.c, Revision 1.5
1.2 noro 1: /*
2: * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
3: * All rights reserved.
4: *
5: * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
6: * non-exclusive and royalty-free license to use, copy, modify and
7: * redistribute, solely for non-commercial and non-profit purposes, the
8: * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
9: * conditions of this Agreement. For the avoidance of doubt, you acquire
10: * only a limited right to use the SOFTWARE hereunder, and FLL or any
11: * third party developer retains all rights, including but not limited to
12: * copyrights, in and to the SOFTWARE.
13: *
14: * (1) FLL does not grant you a license in any way for commercial
15: * purposes. You may use the SOFTWARE only for non-commercial and
16: * non-profit purposes only, such as academic, research and internal
17: * business use.
18: * (2) The SOFTWARE is protected by the Copyright Law of Japan and
19: * international copyright treaties. If you make copies of the SOFTWARE,
20: * with or without modification, as permitted hereunder, you shall affix
21: * to all such copies of the SOFTWARE the above copyright notice.
22: * (3) An explicit reference to this SOFTWARE and its copyright owner
23: * shall be made on your publication or presentation in any form of the
24: * results obtained by use of the SOFTWARE.
25: * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
1.3 noro 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.2 noro 27: * for such modification or the source code of the modified part of the
28: * SOFTWARE.
29: *
30: * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
31: * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
32: * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
33: * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
34: * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
35: * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
36: * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
37: * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
38: * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
39: * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
40: * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
41: * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
42: * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
43: * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
44: * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
45: * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
46: * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
47: *
1.5 ! ohara 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/cplx.c,v 1.4 2000/12/22 10:03:28 saito Exp $
1.2 noro 49: */
1.1 noro 50: #include "ca.h"
51: #include "base.h"
1.5 ! ohara 52: #if defined(PARI)
1.1 noro 53: #include "genpari.h"
54: void patori(GEN,Obj *);
55: void patori_i(GEN,N *);
56: void ritopa(Obj,GEN *);
57: void ritopa_i(N,int,GEN *);
58: #endif
59:
60: void toreim(a,rp,ip)
61: Num a;
62: Num *rp,*ip;
63: {
64: if ( !a )
65: *rp = *ip = 0;
1.4 saito 66: #if defined(INTERVAL)
67: else if ( NID(a) <= N_PRE_C ) {
68: #else
1.1 noro 69: else if ( NID(a) <= N_B ) {
1.4 saito 70: #endif
1.1 noro 71: *rp = a; *ip = 0;
72: } else {
73: *rp = ((C)a)->r; *ip = ((C)a)->i;
74: }
75: }
76:
77: void reimtocplx(r,i,cp)
78: Num r,i;
79: Num *cp;
80: {
81: C c;
82:
83: if ( !i )
84: *cp = r;
85: else {
86: NEWC(c); c->r = r; c->i = i; *cp = (Num)c;
87: }
88: }
89:
90: void addcplx(a,b,c)
91: Num a,b;
92: Num *c;
93: {
94: Num ar,ai,br,bi,cr,ci;
95:
96: if ( !a )
97: *c = b;
98: else if ( !b )
99: *c = a;
1.4 saito 100: #if defined(INTERVAL)
101: else if ( (NID(a) <= N_PRE_C) && (NID(b) <= N_PRE_C ) )
102: #else
1.1 noro 103: else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
1.4 saito 104: #endif
1.1 noro 105: addnum(0,a,b,c);
106: else {
107: toreim(a,&ar,&ai); toreim(b,&br,&bi);
108: addnum(0,ar,br,&cr); addnum(0,ai,bi,&ci);
109: reimtocplx(cr,ci,c);
110: }
111: }
112:
113: void subcplx(a,b,c)
114: Num a,b;
115: Num *c;
116: {
117: Num ar,ai,br,bi,cr,ci;
118:
119: if ( !a )
120: chsgnnum(b,c);
121: else if ( !b )
122: *c = a;
1.4 saito 123: #if defined(INTERVAL)
124: else if ( (NID(a) <= N_PRE_C) && (NID(b) <= N_PRE_C ) )
125: #else
1.1 noro 126: else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
1.4 saito 127: #endif
1.1 noro 128: subnum(0,a,b,c);
129: else {
130: toreim(a,&ar,&ai); toreim(b,&br,&bi);
131: subnum(0,ar,br,&cr); subnum(0,ai,bi,&ci);
132: reimtocplx(cr,ci,c);
133: }
134: }
135:
136: void mulcplx(a,b,c)
137: Num a,b;
138: Num *c;
139: {
140: Num ar,ai,br,bi,cr,ci,t,s;
141:
142: if ( !a || !b )
143: *c = 0;
1.4 saito 144: #if defined(INTERVAL)
145: else if ( (NID(a) <= N_PRE_C) && (NID(b) <= N_PRE_C ) )
146: #else
1.1 noro 147: else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
1.4 saito 148: #endif
1.1 noro 149: mulnum(0,a,b,c);
150: else {
151: toreim(a,&ar,&ai); toreim(b,&br,&bi);
152: mulnum(0,ar,br,&t); mulnum(0,ai,bi,&s); subnum(0,t,s,&cr);
153: mulnum(0,ar,bi,&t); mulnum(0,ai,br,&s); addnum(0,t,s,&ci);
154: reimtocplx(cr,ci,c);
155: }
156: }
157:
158: void divcplx(a,b,c)
159: Num a,b;
160: Num *c;
161: {
162: Num ar,ai,br,bi,cr,ci,t,s,u,w;
163:
164: if ( !b )
165: error("divcplx : division by 0");
166: else if ( !a )
167: *c = 0;
1.4 saito 168: #if defined(INTERVAL)
169: else if ( (NID(a) <= N_PRE_C) && (NID(b) <= N_PRE_C ) )
170: #else
1.1 noro 171: else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
1.4 saito 172: #endif
1.1 noro 173: divnum(0,a,b,c);
174: else {
175: toreim(a,&ar,&ai); toreim(b,&br,&bi);
176: mulnum(0,br,br,&t); mulnum(0,bi,bi,&s); addnum(0,t,s,&u);
177: mulnum(0,ar,br,&t); mulnum(0,ai,bi,&s);
178: addnum(0,t,s,&w); divnum(0,w,u,&cr);
179: mulnum(0,ar,bi,&t); mulnum(0,ai,br,&s);
180: subnum(0,s,t,&w); divnum(0,w,u,&ci);
181: reimtocplx(cr,ci,c);
182: }
183: }
184:
185: void pwrcplx(a,e,c)
186: Num a,e;
187: Num *c;
188: {
189: int ei;
190: Num t;
191: extern long prec;
192:
193: if ( !e )
194: *c = (Num)ONE;
195: else if ( !a )
196: *c = 0;
197: else if ( !INT(e) ) {
1.5 ! ohara 198: #if defined(PARI)
1.1 noro 199: GEN pa,pe,z;
200: int ltop,lbot;
201:
202: ltop = avma; ritopa((Obj)a,&pa); ritopa((Obj)e,&pe); lbot = avma;
203: z = gerepile(ltop,lbot,gpui(pa,pe,prec));
204: patori(z,(Obj *)c); cgiv(z);
205: #else
206: error("pwrcplx : can't calculate a fractional power");
207: #endif
208: } else {
209: ei = SGN((Q)e)*QTOS((Q)e);
210: pwrcplx0(a,ei,&t);
211: if ( SGN((Q)e) < 0 )
212: divnum(0,(Num)ONE,t,c);
213: else
214: *c = t;
215: }
216: }
217:
218: void pwrcplx0(a,e,c)
219: Num a;
220: int e;
221: Num *c;
222: {
223: Num t,s;
224:
225: if ( e == 1 )
226: *c = a;
227: else {
228: pwrcplx0(a,e/2,&t);
229: if ( !(e%2) )
230: mulnum(0,t,t,c);
231: else {
232: mulnum(0,t,t,&s); mulnum(0,s,a,c);
233: }
234: }
235: }
236:
237: void chsgncplx(a,c)
238: Num a,*c;
239: {
240: Num r,i;
241:
242: if ( !a )
243: *c = 0;
1.4 saito 244: #if defined(INTERVAL)
245: else if ( NID(a) <= N_PRE_C )
246: #else
1.1 noro 247: else if ( NID(a) <= N_B )
1.4 saito 248: #endif
1.1 noro 249: chsgnnum(a,c);
250: else {
251: chsgnnum(((C)a)->r,&r); chsgnnum(((C)a)->i,&i);
252: reimtocplx(r,i,c);
253: }
254: }
255:
256: int cmpcplx(a,b)
257: Num a,b;
258: {
259: Num ar,ai,br,bi;
260: int s;
261:
262: if ( !a ) {
1.4 saito 263: #if defined(INTERVAL)
264: if ( !b || (NID(b)<=N_PRE_C) )
265: #else
1.1 noro 266: if ( !b || (NID(b)<=N_B) )
1.4 saito 267: #endif
1.1 noro 268: return compnum(0,a,b);
269: else
270: return -1;
271: } else if ( !b ) {
1.4 saito 272: #if defined(INTERVAL)
273: if ( !a || (NID(a)<=N_PRE_C) )
274: #else
1.1 noro 275: if ( !a || (NID(a)<=N_B) )
1.4 saito 276: #endif
1.1 noro 277: return compnum(0,a,b);
278: else
279: return 1;
280: } else {
281: toreim(a,&ar,&ai); toreim(b,&br,&bi);
282: s = compnum(0,ar,br);
283: if ( s )
284: return s;
285: else
286: return compnum(0,ai,bi);
287: }
288: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to cplx.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / engine