src/complex/elem/division/cl_C_recip.cc

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