![[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|
Return to cplx.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.1 noro 1: /* $OpenXM: OpenXM/src/asir99/engine/cplx.c,v 1.1.1.1 1999/11/10 08:12:26 noro Exp $ */
2: #include "ca.h"
3: #include "base.h"
4: #if PARI
5: #include "genpari.h"
6: void patori(GEN,Obj *);
7: void patori_i(GEN,N *);
8: void ritopa(Obj,GEN *);
9: void ritopa_i(N,int,GEN *);
10: #endif
11:
12: void toreim(a,rp,ip)
13: Num a;
14: Num *rp,*ip;
15: {
16: if ( !a )
17: *rp = *ip = 0;
18: else if ( NID(a) <= N_B ) {
19: *rp = a; *ip = 0;
20: } else {
21: *rp = ((C)a)->r; *ip = ((C)a)->i;
22: }
23: }
24:
25: void reimtocplx(r,i,cp)
26: Num r,i;
27: Num *cp;
28: {
29: C c;
30:
31: if ( !i )
32: *cp = r;
33: else {
34: NEWC(c); c->r = r; c->i = i; *cp = (Num)c;
35: }
36: }
37:
38: void addcplx(a,b,c)
39: Num a,b;
40: Num *c;
41: {
42: Num ar,ai,br,bi,cr,ci;
43:
44: if ( !a )
45: *c = b;
46: else if ( !b )
47: *c = a;
48: else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
49: addnum(0,a,b,c);
50: else {
51: toreim(a,&ar,&ai); toreim(b,&br,&bi);
52: addnum(0,ar,br,&cr); addnum(0,ai,bi,&ci);
53: reimtocplx(cr,ci,c);
54: }
55: }
56:
57: void subcplx(a,b,c)
58: Num a,b;
59: Num *c;
60: {
61: Num ar,ai,br,bi,cr,ci;
62:
63: if ( !a )
64: chsgnnum(b,c);
65: else if ( !b )
66: *c = a;
67: else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
68: subnum(0,a,b,c);
69: else {
70: toreim(a,&ar,&ai); toreim(b,&br,&bi);
71: subnum(0,ar,br,&cr); subnum(0,ai,bi,&ci);
72: reimtocplx(cr,ci,c);
73: }
74: }
75:
76: void mulcplx(a,b,c)
77: Num a,b;
78: Num *c;
79: {
80: Num ar,ai,br,bi,cr,ci,t,s;
81:
82: if ( !a || !b )
83: *c = 0;
84: else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
85: mulnum(0,a,b,c);
86: else {
87: toreim(a,&ar,&ai); toreim(b,&br,&bi);
88: mulnum(0,ar,br,&t); mulnum(0,ai,bi,&s); subnum(0,t,s,&cr);
89: mulnum(0,ar,bi,&t); mulnum(0,ai,br,&s); addnum(0,t,s,&ci);
90: reimtocplx(cr,ci,c);
91: }
92: }
93:
94: void divcplx(a,b,c)
95: Num a,b;
96: Num *c;
97: {
98: Num ar,ai,br,bi,cr,ci,t,s,u,w;
99:
100: if ( !b )
101: error("divcplx : division by 0");
102: else if ( !a )
103: *c = 0;
104: else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
105: divnum(0,a,b,c);
106: else {
107: toreim(a,&ar,&ai); toreim(b,&br,&bi);
108: mulnum(0,br,br,&t); mulnum(0,bi,bi,&s); addnum(0,t,s,&u);
109: mulnum(0,ar,br,&t); mulnum(0,ai,bi,&s);
110: addnum(0,t,s,&w); divnum(0,w,u,&cr);
111: mulnum(0,ar,bi,&t); mulnum(0,ai,br,&s);
112: subnum(0,s,t,&w); divnum(0,w,u,&ci);
113: reimtocplx(cr,ci,c);
114: }
115: }
116:
117: void pwrcplx(a,e,c)
118: Num a,e;
119: Num *c;
120: {
121: int ei;
122: Num t;
123: extern long prec;
124:
125: if ( !e )
126: *c = (Num)ONE;
127: else if ( !a )
128: *c = 0;
129: else if ( !INT(e) ) {
130: #if PARI
131: GEN pa,pe,z;
132: int ltop,lbot;
133:
134: ltop = avma; ritopa((Obj)a,&pa); ritopa((Obj)e,&pe); lbot = avma;
135: z = gerepile(ltop,lbot,gpui(pa,pe,prec));
136: patori(z,(Obj *)c); cgiv(z);
137: #else
138: error("pwrcplx : can't calculate a fractional power");
139: #endif
140: } else {
141: ei = SGN((Q)e)*QTOS((Q)e);
142: pwrcplx0(a,ei,&t);
143: if ( SGN((Q)e) < 0 )
144: divnum(0,(Num)ONE,t,c);
145: else
146: *c = t;
147: }
148: }
149:
150: void pwrcplx0(a,e,c)
151: Num a;
152: int e;
153: Num *c;
154: {
155: Num t,s;
156:
157: if ( e == 1 )
158: *c = a;
159: else {
160: pwrcplx0(a,e/2,&t);
161: if ( !(e%2) )
162: mulnum(0,t,t,c);
163: else {
164: mulnum(0,t,t,&s); mulnum(0,s,a,c);
165: }
166: }
167: }
168:
169: void chsgncplx(a,c)
170: Num a,*c;
171: {
172: Num r,i;
173:
174: if ( !a )
175: *c = 0;
176: else if ( NID(a) <= N_B )
177: chsgnnum(a,c);
178: else {
179: chsgnnum(((C)a)->r,&r); chsgnnum(((C)a)->i,&i);
180: reimtocplx(r,i,c);
181: }
182: }
183:
184: int cmpcplx(a,b)
185: Num a,b;
186: {
187: Num ar,ai,br,bi;
188: int s;
189:
190: if ( !a ) {
191: if ( !b || (NID(b)<=N_B) )
192: return compnum(0,a,b);
193: else
194: return -1;
195: } else if ( !b ) {
196: if ( !a || (NID(a)<=N_B) )
197: return compnum(0,a,b);
198: else
199: return 1;
200: } else {
201: toreim(a,&ar,&ai); toreim(b,&br,&bi);
202: s = compnum(0,ar,br);
203: if ( s )
204: return s;
205: else
206: return compnum(0,ai,bi);
207: }
208: }