src/complex/elem/division/cl_C_recip.cc

   1 // recip().
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cl_complex.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cl_C.h"
  13 #include "cl_real.h"
  14 #include "cl_R.h"
  15 #include "cl_rational.h"
  16 #include "cl_RA.h"
  17 #include "cl_F.h"
  18 #include "cl_SF.h"
  19 #include "cl_FF.h"
  20 #include "cl_DF.h"
  21 #include "cl_LF.h"
  22
  23 ALL_cl_LF_OPERATIONS_SAME_PRECISION()
  24
  25 // for GEN_F_OP2:
  26 #define NOMAP2(F,EXPR)  \
  27   cl_C_##F _tmp = EXPR;                                                 \
  28   return complex_C(_tmp.realpart,_tmp.imagpart);
  29 #define MAP2(F,FN,EXPR)  \
  30   cl_C_##F _tmp = EXPR;                                                 \
  31   return complex_C(FN(_tmp.realpart),FN(_tmp.imagpart))
  32
  33 const cl_N recip (const cl_N& x)
  34 {
  35 // Methode:
  36 // Falls x reell, klar.
  37 // Falls x=a+bi:
  38 //    Falls a=0: 0+(-1/b)i.
  39 //    Falls a und b beide rational sind:
  40 //      c:=a*a+b*b, c:=1/c, liefere a*c+(-b*c)i.
  41 //    Falls a oder b Floats sind:
  42 //      Falls einer von beiden rational ist, runde ihn zum selben Float-Typ
  43 //        wie der andere und führe das UP durch.
  44 //      Falls beide Floats sind, erweitere auf den genaueren, führe das UP
  45 //        durch und runde wieder auf den ungenaueren.
  46 //      Das Ergebnis ist eine komplexe Zahl, da beide Komponenten Floats sind.
  47 // UP: [a,b Floats vom selben Typ]
  48 //  a=0.0 -> liefere die Komponenten a=0.0 und -1/b.
  49 //  b=0.0 -> liefere die Komponenten 1/a und b=0.0.
  50 //  e:=max(exponent(a),exponent(b)).
  51 //  a':=a/2^e bzw. 0.0 bei Underflowmöglichkeit (beim Skalieren a':=a/2^e
  52 //      oder beim Quadrieren a'*a':  2*(e-exponent(a))>exp_mid-exp_low-1
  53 //      d.h. exponent(b)-exponent(a)>floor((exp_mid-exp_low-1)/2) ).
  54 //  b':=b/2^e bzw. 0.0 bei Underflowmöglichkeit (beim Skalieren b':=b/2^e
  55 //      oder beim Quadrieren b'*b':  2*(e-exponent(b))>exp_mid-exp_low-1
  56 //      d.h. exponent(a)-exponent(b)>floor((exp_mid-exp_low-1)/2) ).
  57 //  c':=a'*a'+b'*b',
  58 //  liefere die beiden Komponenten 2^(-e)*a'/c' und -2^(-e)*b'/c'.
  59
  60         if (realp(x)) {
  61                 DeclareType(cl_R,x);
  62                 return recip(x);
  63         } else
  64     {
  65         DeclareType(cl_C,x);
  66         var const cl_R& a = realpart(x);
  67         var const cl_R& b = imagpart(x);
  68         // x = a+bi
  69         if (rationalp(a)) {
  70                 DeclareType(cl_RA,a);
  71                 if (eq(a,0))
  72                         // (complex 0 (- (/ b)))
  73                         return complex_C(0,-recip(b));
  74                 if (rationalp(b)) {
  75                         DeclareType(cl_RA,b);
  76                         // a,b beide rational
  77                         var cl_RA c = recip(square(a)+square(b));
  78                         return complex_C(a*c,-b*c);
  79                 } else {
  80                         DeclareType(cl_F,b);
  81                         // a rational, b Float
  82                         floatcase(b
  83                         ,       return complex_C(cl_C_recip(cl_RA_to_SF(a),b));
  84                         ,       return complex_C(cl_C_recip(cl_RA_to_FF(a),b));
  85                         ,       return complex_C(cl_C_recip(cl_RA_to_DF(a),b));
  86                         ,       return complex_C(cl_C_recip(cl_RA_to_LF(a,TheLfloat(b)->len),b));
  87                         );
  88                 }
  89         } else {
  90                 DeclareType(cl_F,a);
  91                 if (rationalp(b)) {
  92                         DeclareType(cl_RA,b);
  93                         // a Float, b rational
  94                         floatcase(a
  95                         ,       return complex_C(cl_C_recip(a,cl_RA_to_SF(b)));
  96                         ,       return complex_C(cl_C_recip(a,cl_RA_to_FF(b)));
  97                         ,       return complex_C(cl_C_recip(a,cl_RA_to_DF(b)));
  98                         ,       return complex_C(cl_C_recip(a,cl_RA_to_LF(b,TheLfloat(a)->len)));
  99                         );
 100                 } else {
 101                         DeclareType(cl_F,b);
 102                         // a,b Floats
 103                         #ifndef CL_LF_PEDANTIC
 104                         GEN_F_OP2(a,b,cl_C_recip,2,1,); // uses NOMAP2, MAP2.
 105                         #else
 106                         GEN_F_OP2(a,b,cl_C_recip,2,0,); // uses NOMAP2, MAP2.
 107                         #endif
 108                 }
 109         }
 110     }
 111 }