Annotation of OpenXM_contrib2/asir2000/engine/cplx.c, Revision 1.6
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.6 ! ohara 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/cplx.c,v 1.5 2003/02/14 22:29:08 ohara 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: #endif
55:
56: void toreim(a,rp,ip)
57: Num a;
58: Num *rp,*ip;
59: {
60: if ( !a )
61: *rp = *ip = 0;
1.4 saito 62: #if defined(INTERVAL)
63: else if ( NID(a) <= N_PRE_C ) {
64: #else
1.1 noro 65: else if ( NID(a) <= N_B ) {
1.4 saito 66: #endif
1.1 noro 67: *rp = a; *ip = 0;
68: } else {
69: *rp = ((C)a)->r; *ip = ((C)a)->i;
70: }
71: }
72:
73: void reimtocplx(r,i,cp)
74: Num r,i;
75: Num *cp;
76: {
77: C c;
78:
79: if ( !i )
80: *cp = r;
81: else {
82: NEWC(c); c->r = r; c->i = i; *cp = (Num)c;
83: }
84: }
85:
86: void addcplx(a,b,c)
87: Num a,b;
88: Num *c;
89: {
90: Num ar,ai,br,bi,cr,ci;
91:
92: if ( !a )
93: *c = b;
94: else if ( !b )
95: *c = a;
1.4 saito 96: #if defined(INTERVAL)
97: else if ( (NID(a) <= N_PRE_C) && (NID(b) <= N_PRE_C ) )
98: #else
1.1 noro 99: else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
1.4 saito 100: #endif
1.1 noro 101: addnum(0,a,b,c);
102: else {
103: toreim(a,&ar,&ai); toreim(b,&br,&bi);
104: addnum(0,ar,br,&cr); addnum(0,ai,bi,&ci);
105: reimtocplx(cr,ci,c);
106: }
107: }
108:
109: void subcplx(a,b,c)
110: Num a,b;
111: Num *c;
112: {
113: Num ar,ai,br,bi,cr,ci;
114:
115: if ( !a )
116: chsgnnum(b,c);
117: else if ( !b )
118: *c = a;
1.4 saito 119: #if defined(INTERVAL)
120: else if ( (NID(a) <= N_PRE_C) && (NID(b) <= N_PRE_C ) )
121: #else
1.1 noro 122: else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
1.4 saito 123: #endif
1.1 noro 124: subnum(0,a,b,c);
125: else {
126: toreim(a,&ar,&ai); toreim(b,&br,&bi);
127: subnum(0,ar,br,&cr); subnum(0,ai,bi,&ci);
128: reimtocplx(cr,ci,c);
129: }
130: }
131:
132: void mulcplx(a,b,c)
133: Num a,b;
134: Num *c;
135: {
136: Num ar,ai,br,bi,cr,ci,t,s;
137:
138: if ( !a || !b )
139: *c = 0;
1.4 saito 140: #if defined(INTERVAL)
141: else if ( (NID(a) <= N_PRE_C) && (NID(b) <= N_PRE_C ) )
142: #else
1.1 noro 143: else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
1.4 saito 144: #endif
1.1 noro 145: mulnum(0,a,b,c);
146: else {
147: toreim(a,&ar,&ai); toreim(b,&br,&bi);
148: mulnum(0,ar,br,&t); mulnum(0,ai,bi,&s); subnum(0,t,s,&cr);
149: mulnum(0,ar,bi,&t); mulnum(0,ai,br,&s); addnum(0,t,s,&ci);
150: reimtocplx(cr,ci,c);
151: }
152: }
153:
154: void divcplx(a,b,c)
155: Num a,b;
156: Num *c;
157: {
158: Num ar,ai,br,bi,cr,ci,t,s,u,w;
159:
160: if ( !b )
161: error("divcplx : division by 0");
162: else if ( !a )
163: *c = 0;
1.4 saito 164: #if defined(INTERVAL)
165: else if ( (NID(a) <= N_PRE_C) && (NID(b) <= N_PRE_C ) )
166: #else
1.1 noro 167: else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
1.4 saito 168: #endif
1.1 noro 169: divnum(0,a,b,c);
170: else {
171: toreim(a,&ar,&ai); toreim(b,&br,&bi);
172: mulnum(0,br,br,&t); mulnum(0,bi,bi,&s); addnum(0,t,s,&u);
173: mulnum(0,ar,br,&t); mulnum(0,ai,bi,&s);
174: addnum(0,t,s,&w); divnum(0,w,u,&cr);
175: mulnum(0,ar,bi,&t); mulnum(0,ai,br,&s);
176: subnum(0,s,t,&w); divnum(0,w,u,&ci);
177: reimtocplx(cr,ci,c);
178: }
179: }
180:
181: void pwrcplx(a,e,c)
182: Num a,e;
183: Num *c;
184: {
185: int ei;
186: Num t;
187:
188: if ( !e )
189: *c = (Num)ONE;
190: else if ( !a )
191: *c = 0;
192: else if ( !INT(e) ) {
1.5 ohara 193: #if defined(PARI)
1.6 ! ohara 194: gpui_ri((Obj)a,(Obj)c,(Obj *)c);
1.1 noro 195: #else
196: error("pwrcplx : can't calculate a fractional power");
197: #endif
198: } else {
199: ei = SGN((Q)e)*QTOS((Q)e);
200: pwrcplx0(a,ei,&t);
201: if ( SGN((Q)e) < 0 )
202: divnum(0,(Num)ONE,t,c);
203: else
204: *c = t;
205: }
206: }
207:
208: void pwrcplx0(a,e,c)
209: Num a;
210: int e;
211: Num *c;
212: {
213: Num t,s;
214:
215: if ( e == 1 )
216: *c = a;
217: else {
218: pwrcplx0(a,e/2,&t);
219: if ( !(e%2) )
220: mulnum(0,t,t,c);
221: else {
222: mulnum(0,t,t,&s); mulnum(0,s,a,c);
223: }
224: }
225: }
226:
227: void chsgncplx(a,c)
228: Num a,*c;
229: {
230: Num r,i;
231:
232: if ( !a )
233: *c = 0;
1.4 saito 234: #if defined(INTERVAL)
235: else if ( NID(a) <= N_PRE_C )
236: #else
1.1 noro 237: else if ( NID(a) <= N_B )
1.4 saito 238: #endif
1.1 noro 239: chsgnnum(a,c);
240: else {
241: chsgnnum(((C)a)->r,&r); chsgnnum(((C)a)->i,&i);
242: reimtocplx(r,i,c);
243: }
244: }
245:
246: int cmpcplx(a,b)
247: Num a,b;
248: {
249: Num ar,ai,br,bi;
250: int s;
251:
252: if ( !a ) {
1.4 saito 253: #if defined(INTERVAL)
254: if ( !b || (NID(b)<=N_PRE_C) )
255: #else
1.1 noro 256: if ( !b || (NID(b)<=N_B) )
1.4 saito 257: #endif
1.1 noro 258: return compnum(0,a,b);
259: else
260: return -1;
261: } else if ( !b ) {
1.4 saito 262: #if defined(INTERVAL)
263: if ( !a || (NID(a)<=N_PRE_C) )
264: #else
1.1 noro 265: if ( !a || (NID(a)<=N_B) )
1.4 saito 266: #endif
1.1 noro 267: return compnum(0,a,b);
268: else
269: return 1;
270: } else {
271: toreim(a,&ar,&ai); toreim(b,&br,&bi);
272: s = compnum(0,ar,br);
273: if ( s )
274: return s;
275: else
276: return compnum(0,ai,bi);
277: }
278: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>