![[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: }