home *** CD-ROM | disk | FTP | other *** search
/ Amiga ISO Collection / AmigaUtilCD2.iso / Programming / Misc / CLISP-1.LHA / CLISP960530-sr.lha / src / compelem.d < prev    next >
Encoding:
Text File  |  1996-04-15  |  33.9 KB  |  857 lines
  1. # Elementare Funktionen für komplexe Zahlen
  2.  
  3. # Liefert zu reellen Zahlen a und b /= Fixnum 0 die komplexe Zahl a+bi.
  4. # R_R_complex_C(a,b)
  5. # kann GC auslösen
  6.   #define R_R_complex_C  make_complex
  7.  
  8. # Liefert zu reellen Zahlen a und b die komplexe Zahl a+bi.
  9. # R_R_complex_N(a,b)
  10. # kann GC auslösen
  11.   local object R_R_complex_N (object a, object b);
  12. # Methode:
  13. # Falls b=0, nur a. sonst komplexe Zahl erzeugen.
  14.   local object R_R_complex_N(a,b)
  15.     var reg2 object a;
  16.     var reg1 object b;
  17.     { return (eq(b,Fixnum_0) ? a : R_R_complex_C(a,b)); }
  18.  
  19. # N_realpart_R(x) liefert den Realteil der Zahl x.
  20.   local object N_realpart_R (object x);
  21. #if 0
  22.   local object N_realpart_R(x)
  23.     var reg1 object x;
  24.     { return (N_realp(x) ? x : TheComplex(x)->c_real); }
  25. #else # Macro spart Code
  26.   #define N_realpart_R(x)  (N_realp(x) ? x : TheComplex(x)->c_real)
  27. #endif
  28.  
  29. # N_imagpart_R(x) liefert den Imaginärteil der Zahl x.
  30.   local object N_imagpart_R (object x);
  31. #if 0
  32.   local object N_imagpart_R(x)
  33.     var reg1 object x;
  34.     { return (N_realp(x) ? Fixnum_0 : TheComplex(x)->c_imag); }
  35. #else # Macro spart Code
  36.   #define N_imagpart_R(x)  (N_realp(x) ? Fixnum_0 : TheComplex(x)->c_imag)
  37. #endif
  38.  
  39. # N_conjugate_N(x) liefert die konjugiert komplexe Zahl zur Zahl x.
  40. # kann GC auslösen
  41.   local object N_conjugate_N (object x);
  42.   local object N_conjugate_N(x)
  43.     var reg1 object x;
  44.     { if (N_realp(x))
  45.         { return x; }
  46.         else
  47.         { pushSTACK(TheComplex(x)->c_real);
  48.          {var reg2 object b = TheComplex(x)->c_imag;
  49.           b = R_minus_R(b); # Vorzeichenwechsel beim Imaginärteil
  50.           return R_R_complex_C(popSTACK(),b);
  51.     }   }}
  52.  
  53. # N_minus_N(x) liefert (- x), wo x eine Zahl ist.
  54. # kann GC auslösen
  55.   local object N_minus_N (object x);
  56. # Methode:
  57. # x reell -> klar.
  58. # x=a+bi -> (-a) + (-b) i
  59.   local object N_minus_N(x)
  60.     var reg1 object x;
  61.     { if (N_realp(x))
  62.         { return R_minus_R(x); }
  63.         else
  64.         { pushSTACK(TheComplex(x)->c_real);
  65.          {var reg2 object b = R_minus_R(TheComplex(x)->c_imag);
  66.           var reg3 object a = STACK_0; STACK_0 = b;
  67.           a = R_minus_R(a);
  68.           return R_R_complex_C(a,popSTACK());
  69.     }   }}
  70.  
  71. # N_N_plus_N(x) liefert (+ x y), wo x und y Zahlen sind.
  72. # kann GC auslösen
  73.   local object N_N_plus_N (object x, object y);
  74. # Methode:
  75. # x,y beide reell -> klar.
  76. # x=a, y=b+ci -> (a+b)+ci
  77. # x=a+bi, y=c -> (a+c)+bi
  78. # x=a+bi, y=c+di -> (a+c)+(b+d)i
  79.   local object N_N_plus_N(x,y)
  80.     var reg2 object x;
  81.     var reg1 object y;
  82.     { if (N_realp(x))
  83.         { if (N_realp(y))
  84.             { return R_R_plus_R(x,y); }
  85.             else
  86.             # x=a, y=b+ci
  87.             { pushSTACK(TheComplex(y)->c_imag); # c
  88.              {var reg3 object re = R_R_plus_R(x,TheComplex(y)->c_real); # a+b
  89.               return R_R_complex_C(re,popSTACK());
  90.             }}
  91.         }
  92.         else
  93.         { if (N_realp(y))
  94.             # x=a+bi, y=c
  95.             { pushSTACK(TheComplex(x)->c_imag); # b
  96.              {var reg3 object re = R_R_plus_R(TheComplex(x)->c_real,y); # a+c
  97.               return R_R_complex_C(re,popSTACK());
  98.             }}
  99.             else
  100.             # x=a+bi, y=c+di
  101.             { pushSTACK(TheComplex(x)->c_real); # a
  102.               pushSTACK(TheComplex(y)->c_real); # c
  103.              {var reg3 object temp = # b+d
  104.                 R_R_plus_R(TheComplex(x)->c_imag,TheComplex(y)->c_imag);
  105.               y = STACK_0; STACK_0 = temp;
  106.               temp = R_R_plus_R(STACK_1,y); # a+c
  107.               temp = R_R_complex_N(temp,STACK_0);
  108.               skipSTACK(2); return temp;
  109.             }}
  110.     }   }
  111.  
  112. # N_N_minus_N(x) liefert (- x y), wo x und y Zahlen sind.
  113. # kann GC auslösen
  114.   local object N_N_minus_N (object x, object y);
  115. # Methode:
  116. # x,y beide reell -> klar.
  117. # x=a, y=b+ci -> (a-b)+(-c)i
  118. # x=a+bi, y=c -> (a-c)+bi
  119. # x=a+bi, y=c+di -> (a-c)+(b-d)i
  120.   local object N_N_minus_N(x,y)
  121.     var reg2 object x;
  122.     var reg1 object y;
  123.     { if (N_realp(x))
  124.         { if (N_realp(y))
  125.             { return R_R_minus_R(x,y); }
  126.             else
  127.             # x=a, y=b+ci
  128.             { pushSTACK(TheComplex(y)->c_imag); # c
  129.              {var reg3 object temp = R_R_minus_R(x,TheComplex(y)->c_real); # a-b
  130.               y = STACK_0; STACK_0 = temp;
  131.               y = R_minus_R(y); # -c
  132.               return R_R_complex_C(popSTACK(),y);
  133.             }}
  134.         }
  135.         else
  136.         { if (N_realp(y))
  137.             # x=a+bi, y=c
  138.             { pushSTACK(TheComplex(x)->c_imag); # b
  139.              {var reg3 object re = R_R_minus_R(TheComplex(x)->c_real,y); # a-c
  140.               return R_R_complex_C(re,popSTACK());
  141.             }}
  142.             else
  143.             # x=a+bi, y=c+di
  144.             { pushSTACK(TheComplex(x)->c_real); # a
  145.               pushSTACK(TheComplex(y)->c_real); # c
  146.              {var reg3 object temp = # b-d
  147.                 R_R_minus_R(TheComplex(x)->c_imag,TheComplex(y)->c_imag);
  148.               y = STACK_0; STACK_0 = temp;
  149.               temp = R_R_minus_R(STACK_1,y); # a-c
  150.               temp = R_R_complex_N(temp,STACK_0);
  151.               skipSTACK(2); return temp;
  152.             }}
  153.     }   }
  154.  
  155. # N_1_plus_N(x) liefert (1+ x), wo x eine Zahl ist.
  156. # kann GC auslösen
  157.   local object N_1_plus_N (object x);
  158. # Methode:
  159. # x reell -> klar.
  160. # x=a+bi -> (a+1)+bi
  161.   local object N_1_plus_N(x)
  162.     var reg1 object x;
  163.     { if (N_realp(x))
  164.         { return R_1_plus_R(x); }
  165.         else
  166.         { pushSTACK(TheComplex(x)->c_imag);
  167.          {var reg2 object re = R_1_plus_R(TheComplex(x)->c_real); # a+1
  168.           return R_R_complex_C(re,popSTACK());
  169.     }   }}
  170.  
  171. # N_minus1_plus_N(x) liefert (1- x), wo x eine Zahl ist.
  172. # kann GC auslösen
  173.   local object N_minus1_plus_N (object x);
  174. # Methode:
  175. # x reell -> klar.
  176. # x=a+bi -> (a-1)+bi
  177.   local object N_minus1_plus_N(x)
  178.     var reg1 object x;
  179.     { if (N_realp(x))
  180.         { return R_minus1_plus_R(x); }
  181.         else
  182.         { pushSTACK(TheComplex(x)->c_imag);
  183.          {var reg2 object re = R_minus1_plus_R(TheComplex(x)->c_real); # a-1
  184.           return R_R_complex_C(re,popSTACK());
  185.     }   }}
  186.  
  187. # N_N_mal_N(x) liefert (* x y), wo x und y Zahlen sind.
  188. # kann GC auslösen
  189.   local object N_N_mal_N (object x, object y);
  190. # Methode:
  191. # x,y beide reell -> klar.
  192. # x=a, y=b+ci -> (a*b)+(a*c)i
  193. # x=a+bi, y=c -> (a*c)+(b*c)i
  194. # x=a+bi, y=c+di -> (a*c-b*d)+(a*d+b*c)i
  195.   local object N_N_mal_N(x,y)
  196.     var reg2 object x;
  197.     var reg1 object y;
  198.     { if (N_realp(x))
  199.         { if (N_realp(y))
  200.             { return R_R_mal_R(x,y); }
  201.             else
  202.             # x=a, y=b+ci
  203.             { pushSTACK(x); # a
  204.               pushSTACK(TheComplex(y)->c_real); # b
  205.              {var reg3 object temp = R_R_mal_R(x,TheComplex(y)->c_imag); # a*c
  206.               y = popSTACK(); x = STACK_0; STACK_0 = temp;
  207.               y = R_R_mal_R(x,y); # a*b
  208.               return R_R_complex_N(y,popSTACK());
  209.             }}
  210.         }
  211.         else
  212.         { if (N_realp(y))
  213.             # x=a+bi, y=c
  214.             { pushSTACK(y); # c
  215.               pushSTACK(TheComplex(x)->c_real); # a
  216.              {var reg3 object temp = R_R_mal_R(TheComplex(x)->c_imag,y); # b*c
  217.               x = popSTACK(); y = STACK_0; STACK_0 = temp;
  218.               x = R_R_mal_R(x,y); # a*c
  219.               return R_R_complex_N(x,popSTACK());
  220.             }}
  221.             else
  222.             # x=a+bi, y=c+di -> (a*c-b*d)+(a*d+b*c)i
  223.             { pushSTACK(TheComplex(x)->c_real); # a
  224.               pushSTACK(TheComplex(y)->c_real); # c
  225.               pushSTACK(TheComplex(x)->c_imag); # b
  226.               pushSTACK(y = TheComplex(y)->c_imag); # d
  227.               # Stackaufbau: a, c, b, d.
  228.               y = R_R_mal_R(STACK_3,y); # a*d
  229.               x = STACK_3; STACK_3 = y;
  230.               # Stackaufbau: a*d, c, b, d.
  231.               x = R_R_mal_R(x,STACK_2); # a*c
  232.               y = STACK_2; STACK_2 = x;
  233.               # Stackaufbau: a*d, a*c, b, d.
  234.               y = R_R_mal_R(STACK_1,y); # b*c
  235.               STACK_3 = R_R_plus_R(STACK_3,y); # a*d+b*c
  236.               # Stackaufbau: a*d+b*c, a*c, b, d.
  237.               y = R_R_mal_R(STACK_1,STACK_0); # b*d
  238.               y = R_R_minus_R(STACK_2,y); # a*c-b*d
  239.               x = STACK_3; # a*d+b*c
  240.               skipSTACK(4);
  241.               return R_R_complex_N(y,x);
  242.             }
  243.     }   }
  244.  
  245. # N_zerop(x) stellt fest, ob (= x 0), wo x eine Zahl ist.
  246.   local boolean N_zerop (object x);
  247.   local boolean N_zerop(x)
  248.     var reg1 object x;
  249.     { if (N_realp(x))
  250.         { return R_zerop(x); }
  251.         else
  252.         # x komplex, teste ob Real- und Imaginärteil beide = 0 sind.
  253.         { return (R_zerop(TheComplex(x)->c_real) && R_zerop(TheComplex(x)->c_imag)); }
  254.     }
  255.  
  256. # N_durch_N(x) liefert (/ x), wo x eine Zahl ist.
  257. # kann GC auslösen
  258.   local object N_durch_N (object x);
  259. # Methode:
  260. # Falls x reell, klar.
  261. # Falls x=a+bi:
  262. #    Falls a=0: 0+(-1/b)i.
  263. #    Falls a und b beide rational sind:
  264. #      c:=a*a+b*b, c:=1/c, liefere a*c+(-b*c)i.
  265. #    Falls a oder b Floats sind:
  266. #      Falls einer von beiden rational ist, runde ihn zum selben Float-Typ
  267. #        wie der andere und führe das UP durch.
  268. #      Falls beide Floats sind, erweitere auf den genaueren, führe das UP
  269. #        durch und runde wieder auf den ungenaueren.
  270. #      Das Ergebnis ist eine komplexe Zahl, da beide Komponenten Floats sind.
  271. # UP: [a,b Floats vom selben Typ]
  272. #  a=0.0 -> liefere die Komponenten a=0.0 und -1/b.
  273. #  b=0.0 -> liefere die Komponenten 1/a und b=0.0.
  274. #  e:=max(exponent(a),exponent(b)).
  275. #  a':=a/2^e bzw. 0.0 bei Underflowmöglichkeit (beim Skalieren a':=a/2^e
  276. #      oder beim Quadrieren a'*a':  2*(e-exponent(a))>exp_mid-exp_low-1
  277. #      d.h. exponent(b)-exponent(a)>floor((exp_mid-exp_low-1)/2) ).
  278. #  b':=b/2^e bzw. 0.0 bei Underflowmöglichkeit (beim Skalieren b':=b/2^e
  279. #      oder beim Quadrieren b'*b':  2*(e-exponent(b))>exp_mid-exp_low-1
  280. #      d.h. exponent(a)-exponent(b)>floor((exp_mid-exp_low-1)/2) ).
  281. #  c':=a'*a'+b'*b',
  282. #  liefere die beiden Komponenten 2^(-e)*a'/c' und -2^(-e)*b'/c'.
  283.   local void SFC_durch_SFC (object a, object b);
  284.   local void FFC_durch_FFC (object a, object b);
  285.   local void DFC_durch_DFC (object a, object b);
  286.   local void LFC_durch_LFC (object a, object b);
  287.   local void SFC_durch_SFC(a,b)
  288.     var reg1 object a;
  289.     var reg2 object b;
  290.     { var reg5 sintWL a_exp;
  291.       var reg6 sintWL b_exp;
  292.       {# Exponenten von a holen:
  293.        var reg3 uintBWL uexp = SF_uexp(a);
  294.        if (uexp==0) # a=0.0 -> liefere (complex a (- (/ b))) :
  295.          { pushSTACK(a); pushSTACK(SF_minus_SF(SF_durch_SF(b))); return; }
  296.        a_exp = (sintWL)((uintWL)uexp - SF_exp_mid);
  297.       }
  298.       {# Exponenten von b holen:
  299.        var reg3 uintBWL uexp = SF_uexp(b);
  300.        if (uexp==0) # b=0.0 -> liefere (complex (/ a) b) :
  301.          { pushSTACK(SF_durch_SF(a)); pushSTACK(b); return; }
  302.        b_exp = (sintWL)((uintWL)uexp - SF_exp_mid);
  303.       }
  304.        # Nun a_exp = exponent(a), b_exp = exponent(b).
  305.       {var reg4 sintWL e = (a_exp > b_exp ? a_exp : b_exp); # Maximum der Exponenten
  306.        var reg3 object delta = L_to_FN(-(sintL)e); # -e als Fixnum
  307.        # a und b durch 2^e dividieren:
  308.        a = (b_exp-a_exp > floor(SF_exp_mid-SF_exp_low-1,2) ? SF_0 : SF_I_scale_float_SF(a,delta));
  309.        b = (a_exp-b_exp > floor(SF_exp_mid-SF_exp_low-1,2) ? SF_0 : SF_I_scale_float_SF(b,delta));
  310.        # c' := a'*a'+b'*b' berechnen:
  311.        {var reg4 object c = SF_SF_plus_SF(SF_SF_mal_SF(a,a),SF_SF_mal_SF(b,b));
  312.         a = SF_I_scale_float_SF(SF_SF_durch_SF(a,c),delta); # 2^(-e)*a'/c'
  313.         b = SF_I_scale_float_SF(SF_minus_SF(SF_SF_durch_SF(b,c)),delta); # -2^(-e)*b'/c'
  314.         pushSTACK(a); pushSTACK(b); return;
  315.     } }}
  316.   local void FFC_durch_FFC(a,b)
  317.     var reg1 object a;
  318.     var reg2 object b;
  319.     { var reg5 sintWL a_exp;
  320.       var reg6 sintWL b_exp;
  321.       {# Exponenten von a holen:
  322.        var reg3 uintBWL uexp = FF_uexp(ffloat_value(a));
  323.        if (uexp==0) # a=0.0 -> liefere (complex a (- (/ b))) :
  324.          { pushSTACK(a); pushSTACK(FF_minus_FF(FF_durch_FF(b))); return; }
  325.        a_exp = (sintWL)((uintWL)uexp - FF_exp_mid);
  326.       }
  327.       {# Exponenten von b holen:
  328.        var reg3 uintBWL uexp = FF_uexp(ffloat_value(b));
  329.        if (uexp==0) # b=0.0 -> liefere (complex (/ a) b) :
  330.          { pushSTACK(FF_durch_FF(a)); pushSTACK(FF_0); return; }
  331.        b_exp = (sintWL)((uintWL)uexp - FF_exp_mid);
  332.       }
  333.        pushSTACK(a); pushSTACK(b);
  334.        # Stackaufbau: a, b.
  335.        # Nun a_exp = exponent(a), b_exp = exponent(b).
  336.       {var reg4 sintWL e = (a_exp > b_exp ? a_exp : b_exp); # Maximum der Exponenten
  337.        var reg3 object delta = L_to_FN(-(sintL)e); # -e als Fixnum
  338.        # a und b durch 2^e dividieren:
  339.        STACK_1 = (b_exp-a_exp > floor(FF_exp_mid-FF_exp_low-1,2) ? FF_0 : FF_I_scale_float_FF(STACK_1,delta));
  340.        STACK_0 = (a_exp-b_exp > floor(FF_exp_mid-FF_exp_low-1,2) ? FF_0 : FF_I_scale_float_FF(STACK_0,delta));
  341.        # Stackaufbau: a', b'.
  342.        # c' := a'*a'+b'*b' berechnen:
  343.        {var reg4 object temp;
  344.         temp = STACK_1; pushSTACK(FF_FF_mal_FF(temp,temp)); # a'*a'
  345.         temp = STACK_1; temp = FF_FF_mal_FF(temp,temp); # b'*b'
  346.         STACK_0 = temp = FF_FF_plus_FF(STACK_0,temp); # c' = a'*a'+b'*b'
  347.         STACK_2 = FF_I_scale_float_FF(FF_FF_durch_FF(STACK_2,temp),delta); # 2^(-e)*a'/c'
  348.         STACK_1 = FF_I_scale_float_FF(FF_minus_FF(FF_FF_durch_FF(STACK_1,STACK_0)),delta); # -2^(-e)*b'/c'
  349.         skipSTACK(1); return;
  350.     } }}
  351.   local void DFC_durch_DFC(a,b)
  352.     var reg1 object a;
  353.     var reg2 object b;
  354.     { var reg5 sintWL a_exp;
  355.       var reg6 sintWL b_exp;
  356.       {# Exponenten von a holen:
  357.        var reg3 uintWL uexp = DF_uexp(TheDfloat(a)->float_value_semhi);
  358.        if (uexp==0) # a=0.0 -> liefere (complex a (- (/ b))) :
  359.          { pushSTACK(a); pushSTACK(DF_minus_DF(DF_durch_DF(b))); return; }
  360.        a_exp = (sintWL)(uexp - DF_exp_mid);
  361.       }
  362.       {# Exponenten von b holen:
  363.        var reg3 uintWL uexp = DF_uexp(TheDfloat(b)->float_value_semhi);
  364.        if (uexp==0) # b=0.0 -> liefere (complex (/ a) b) :
  365.          { pushSTACK(DF_durch_DF(a)); pushSTACK(DF_0); return; }
  366.        b_exp = (sintWL)(uexp - DF_exp_mid);
  367.       }
  368.        pushSTACK(a); pushSTACK(b);
  369.        # Stackaufbau: a, b.
  370.        # Nun a_exp = exponent(a), b_exp = exponent(b).
  371.       {var reg4 sintWL e = (a_exp > b_exp ? a_exp : b_exp); # Maximum der Exponenten
  372.        var reg3 object delta = L_to_FN(-(sintL)e); # -e als Fixnum
  373.        # a und b durch 2^e dividieren:
  374.        STACK_1 = (b_exp-a_exp > floor(DF_exp_mid-DF_exp_low-1,2) ? DF_0 : DF_I_scale_float_DF(STACK_1,delta));
  375.        STACK_0 = (a_exp-b_exp > floor(DF_exp_mid-DF_exp_low-1,2) ? DF_0 : DF_I_scale_float_DF(STACK_0,delta));
  376.        # Stackaufbau: a', b'.
  377.        # c' := a'*a'+b'*b' berechnen:
  378.        {var reg4 object temp;
  379.         temp = STACK_1; pushSTACK(DF_DF_mal_DF(temp,temp)); # a'*a'
  380.         temp = STACK_1; temp = DF_DF_mal_DF(temp,temp); # b'*b'
  381.         STACK_0 = temp = DF_DF_plus_DF(STACK_0,temp); # c' = a'*a'+b'*b'
  382.         STACK_2 = DF_I_scale_float_DF(DF_DF_durch_DF(STACK_2,temp),delta); # 2^(-e)*a'/c'
  383.         STACK_1 = DF_I_scale_float_DF(DF_minus_DF(DF_DF_durch_DF(STACK_1,STACK_0)),delta); # -2^(-e)*b'/c'
  384.         skipSTACK(1); return;
  385.     } }}
  386.   local void LFC_durch_LFC(a,b)
  387.     var reg1 object a;
  388.     var reg2 object b;
  389.     { var reg4 uintL a_exp;
  390.       var reg5 uintL b_exp;
  391.       {# Exponenten von a holen:
  392.        a_exp = TheLfloat(a)->expo;
  393.        if (a_exp==0) # a=0.0 -> liefere (complex a (- (/ b))) :
  394.          { pushSTACK(a); pushSTACK(LF_minus_LF(LF_durch_LF(b))); return; }
  395.       }
  396.       {# Exponenten von b holen:
  397.        b_exp = TheLfloat(b)->expo;
  398.        if (b_exp==0) # b=0.0 -> liefere (complex (/ a) b) :
  399.          { pushSTACK(b); pushSTACK(b); STACK_1 = LF_durch_LF(a); return; }
  400.       }
  401.       pushSTACK(a); pushSTACK(b);
  402.       # Stackaufbau: a, b.
  403.       # Nun a_exp = exponent(a)+LF_exp_mid, b_exp = exponent(b)+LF_exp_mid.
  404.       {var reg3 uintL e = (a_exp > b_exp ? a_exp : b_exp); # Maximum der Exponenten
  405.        pushSTACK(L_to_I(-(sintL)(e-LF_exp_mid))); # -e als Integer
  406.       }
  407.       # a und b durch 2^e dividieren:
  408.       if ((b_exp>a_exp) && (b_exp-a_exp > (uintL)floor((uintL)(LF_exp_mid-LF_exp_low-1),2)))
  409.         { encode_LF0(TheLfloat(STACK_2)->len, STACK_2=); }
  410.         else
  411.         { STACK_2 = LF_I_scale_float_LF(STACK_2,STACK_0); }
  412.       if ((a_exp>b_exp) && (a_exp-b_exp > (uintL)floor((uintL)(LF_exp_mid-LF_exp_low-1),2)))
  413.         { encode_LF0(TheLfloat(STACK_1)->len, STACK_1=); }
  414.         else
  415.         { STACK_1 = LF_I_scale_float_LF(STACK_1,STACK_0); }
  416.       # Stackaufbau: a', b', -e.
  417.       # c' := a'*a'+b'*b' berechnen:
  418.       {var reg3 object temp;
  419.        temp = STACK_2; pushSTACK(LF_LF_mal_LF(temp,temp)); # a'*a'
  420.        temp = STACK_2; temp = LF_LF_mal_LF(temp,temp); # b'*b'
  421.        STACK_0 = temp = LF_LF_plus_LF(STACK_0,temp); # c' = a'*a'+b'*b'
  422.        temp = LF_LF_durch_LF(STACK_3,temp); STACK_3 = LF_I_scale_float_LF(temp,STACK_1); # 2^(-e)*a'/c'
  423.        temp = LF_minus_LF(LF_LF_durch_LF(STACK_2,STACK_0)); STACK_2 = LF_I_scale_float_LF(temp,STACK_1); # -2^(-e)*b'/c'
  424.        skipSTACK(2); return;
  425.     } }
  426.   local object N_durch_N(x)
  427.     var reg2 object x;
  428.     { if (N_realp(x))
  429.         { return R_durch_R(x); }
  430.         else
  431.         { var reg3 object a = TheComplex(x)->c_real;
  432.           var reg4 object b = TheComplex(x)->c_imag;
  433.           # x = a+bi
  434.           if (R_rationalp(a))
  435.             { if (eq(a,Fixnum_0))
  436.                 { return R_R_complex_C(Fixnum_0,R_minus_R(R_durch_R(b))); } # (complex 0 (- (/ b)))
  437.               if (R_rationalp(b))
  438.                 # a,b beide rational
  439.                 { var reg1 object temp;
  440.                   pushSTACK(a); pushSTACK(b);
  441.                   temp = RA_RA_mal_RA(a,a); pushSTACK(temp); # a*a
  442.                   temp = STACK_1; temp = RA_RA_mal_RA(temp,temp); # b*b
  443.                   temp = RA_RA_plus_RA(STACK_0,temp); # a*a+b*b = c
  444.                   STACK_0 = temp = RA_durch_RA(temp); # c:=1/c
  445.                   STACK_2 = RA_RA_mal_RA(STACK_2,temp); # a*c
  446.                   STACK_1 = RA_minus_RA(RA_RA_mal_RA(STACK_1,STACK_0)); # -b*c
  447.                   temp = R_R_complex_C(STACK_2,STACK_1); # (a*c) + (-b*c) i
  448.                   skipSTACK(3);
  449.                   return temp;
  450.                 }
  451.                 else
  452.                 # a rational, b Float
  453.                 { pushSTACK(b);
  454.                   floatcase(b,
  455.                             { a = RA_to_SF(a); SFC_durch_SFC(a,popSTACK()); },
  456.                             { a = RA_to_FF(a); FFC_durch_FFC(a,popSTACK()); },
  457.                             { a = RA_to_DF(a); DFC_durch_DFC(a,popSTACK()); },
  458.                             { a = RA_to_LF(a,TheLfloat(b)->len); LFC_durch_LFC(a,popSTACK()); }
  459.                            );
  460.                 }
  461.             }
  462.             else
  463.             { if (R_rationalp(b))
  464.                 # a Float, b rational
  465.                 { pushSTACK(a);
  466.                   floatcase(a,
  467.                             { b = RA_to_SF(b); SFC_durch_SFC(popSTACK(),b); },
  468.                             { b = RA_to_FF(b); FFC_durch_FFC(popSTACK(),b); },
  469.                             { b = RA_to_DF(b); DFC_durch_DFC(popSTACK(),b); },
  470.                             { b = RA_to_LF(b,TheLfloat(a)->len); LFC_durch_LFC(popSTACK(),b); }
  471.                            );
  472.                 }
  473.                 else
  474.                 # a,b Floats
  475.                 GEN_F_op2(a,b,SFC_durch_SFC,FFC_durch_FFC,DFC_durch_DFC,LFC_durch_LFC,2,0,_EMA_)
  476.             }
  477.           # beide Komponenten zu einer komplexen Zahl zusammenfügen:
  478.          {var reg2 object a = STACK_1;
  479.           var reg1 object b = STACK_0;
  480.           skipSTACK(2); return R_R_complex_C(a,b);
  481.         }}
  482.     }
  483.   #define C_durch_C  N_durch_N
  484.  
  485. # N_N_durch_N(x,y) liefert (/ x y), wo x und y Zahlen sind.
  486. # kann GC auslösen
  487.   local object N_N_durch_N (object x, object y);
  488. # Methode:
  489. # x,y beide reell -> klar.
  490. # x=a+bi, y=c reell -> (a/c)+(b/c)i
  491. # y komplex -> (* x (/ y))
  492.   local object N_N_durch_N(x,y)
  493.     var reg2 object x;
  494.     var reg1 object y;
  495.     { if (N_realp(y))
  496.         # y reell
  497.         { if (N_realp(x))
  498.             # x,y beide reell
  499.             { return R_R_durch_R(x,y); }
  500.             else
  501.             # x komplex: x=a+bi, y=c
  502.             { pushSTACK(y); # c
  503.               pushSTACK(TheComplex(x)->c_real); # a
  504.              {var reg3 object temp = R_R_durch_R(TheComplex(x)->c_imag,y); # b/c
  505.               x = popSTACK(); y = STACK_0; STACK_0 = temp;
  506.               x = R_R_durch_R(x,y); # a/c
  507.               return R_R_complex_N(x,popSTACK());
  508.             }}
  509.         }
  510.         else
  511.         # y komplex
  512.         { pushSTACK(x); y = C_durch_C(y); # mit dem Kehrwert von y
  513.           return N_N_mal_N(popSTACK(),y); # multiplizieren
  514.         }
  515.     }
  516.  
  517. # R_R_hypot_R(a,b) liefert sqrt(a^2+b^2), wo a und b reelle Zahlen sind.
  518. # kann GC auslösen
  519.   local object R_R_hypot_R (object a, object b);
  520. # Methode:
  521. # Falls a=0: (abs b).
  522. # Falls b=0: (abs a).
  523. # Falls a und b beide rational sind:
  524. #   c:=a*a+b*b, liefere (sqrt c).
  525. # Falls a oder b Floats sind:
  526. #   Falls einer von beiden rational ist, runde ihn zum selben Float-Typ
  527. #     wie der andere und führe das UP durch.
  528. #   Falls beide Floats sind, erweitere auf den genaueren, führe das UP
  529. #     durch und runde wieder auf den ungenaueren.
  530. #   Das Ergebnis ist ein Float >=0.
  531. # UP: [a,b Floats vom selben Typ]
  532. #  a=0.0 -> liefere abs(b).
  533. #  b=0.0 -> liefere abs(a).
  534. #  e:=max(exponent(a),exponent(b)).
  535. #  a':=a/2^e bzw. 0.0 bei Underflowmöglichkeit (beim Skalieren a':=a/2^e
  536. #      oder beim Quadrieren a'*a':  2*(e-exponent(a))>exp_mid-exp_low-1
  537. #      d.h. exponent(b)-exponent(a)>floor((exp_mid-exp_low-1)/2) ).
  538. #  b':=b/2^e bzw. 0.0 bei Underflowmöglichkeit (beim Skalieren b':=b/2^e
  539. #      oder beim Quadrieren b'*b':  2*(e-exponent(b))>exp_mid-exp_low-1
  540. #      d.h. exponent(a)-exponent(b)>floor((exp_mid-exp_low-1)/2) ).
  541. #  c':=a'*a'+b'*b', c':=sqrt(c'), liefere 2^e*c'.
  542.   local object SF_SF_hypot_SF (object a, object b);
  543.   local object FF_FF_hypot_FF (object a, object b);
  544.   local object DF_DF_hypot_DF (object a, object b);
  545.   local object LF_LF_hypot_LF (object a, object b);
  546.   local object SF_SF_hypot_SF(a,b)
  547.     var reg1 object a;
  548.     var reg2 object b;
  549.     { var reg5 sintWL a_exp;
  550.       var reg6 sintWL b_exp;
  551.       {# Exponenten von a holen:
  552.        var reg3 uintBWL uexp = SF_uexp(a);
  553.        if (uexp==0) # a=0.0 -> liefere (abs b) :
  554.          { return (R_minusp(b) ? SF_minus_SF(b) : b); }
  555.        a_exp = (sintWL)((uintWL)uexp - SF_exp_mid);
  556.       }
  557.       {# Exponenten von b holen:
  558.        var reg3 uintBWL uexp = SF_uexp(b);
  559.        if (uexp==0) # b=0.0 -> liefere (abs a) :
  560.          { return (R_minusp(a) ? SF_minus_SF(a) : a); }
  561.        b_exp = (sintWL)((uintWL)uexp - SF_exp_mid);
  562.       }
  563.        # Nun a_exp = exponent(a), b_exp = exponent(b).
  564.       {var reg4 sintWL e = (a_exp > b_exp ? a_exp : b_exp); # Maximum der Exponenten
  565.        var reg3 object delta = L_to_FN(-(sintL)e); # -e als Fixnum
  566.        # a und b durch 2^e dividieren:
  567.        a = (b_exp-a_exp > floor(SF_exp_mid-SF_exp_low-1,2) ? SF_0 : SF_I_scale_float_SF(a,delta));
  568.        b = (a_exp-b_exp > floor(SF_exp_mid-SF_exp_low-1,2) ? SF_0 : SF_I_scale_float_SF(b,delta));
  569.        # c' := a'*a'+b'*b' berechnen:
  570.        {var reg4 object c = SF_SF_plus_SF(SF_SF_mal_SF(a,a),SF_SF_mal_SF(b,b));
  571.         c = SF_sqrt_SF(c); # c':=2^e*c'
  572.         return SF_I_scale_float_SF(c,L_to_FN((sintL)e)); # 2^e*c'
  573.     } }}
  574.   local object FF_FF_hypot_FF(a,b)
  575.     var reg1 object a;
  576.     var reg2 object b;
  577.     { var reg5 sintWL a_exp;
  578.       var reg6 sintWL b_exp;
  579.       {# Exponenten von a holen:
  580.        var reg3 uintBWL uexp = FF_uexp(ffloat_value(a));
  581.        if (uexp==0) # a=0.0 -> liefere (abs b) :
  582.          { return (R_minusp(b) ? FF_minus_FF(b) : b); }
  583.        a_exp = (sintWL)((uintWL)uexp - FF_exp_mid);
  584.       }
  585.       {# Exponenten von b holen:
  586.        var reg3 uintBWL uexp = FF_uexp(ffloat_value(b));
  587.        if (uexp==0) # b=0.0 -> liefere (abs a) :
  588.          { return (R_minusp(a) ? FF_minus_FF(a) : a); }
  589.        b_exp = (sintWL)((uintWL)uexp - FF_exp_mid);
  590.       }
  591.        pushSTACK(a); pushSTACK(b);
  592.        # Stackaufbau: a, b.
  593.        # Nun a_exp = exponent(a), b_exp = exponent(b).
  594.       {var reg4 sintWL e = (a_exp > b_exp ? a_exp : b_exp); # Maximum der Exponenten
  595.        var reg3 object delta = L_to_FN(-(sintL)e); # -e als Fixnum
  596.        # a und b durch 2^e dividieren:
  597.        STACK_1 = (b_exp-a_exp > floor(FF_exp_mid-FF_exp_low-1,2) ? FF_0 : FF_I_scale_float_FF(STACK_1,delta));
  598.        STACK_0 = (a_exp-b_exp > floor(FF_exp_mid-FF_exp_low-1,2) ? FF_0 : FF_I_scale_float_FF(STACK_0,delta));
  599.        # Stackaufbau: a', b'.
  600.        # c' := a'*a'+b'*b' berechnen:
  601.        {var reg4 object temp;
  602.         temp = STACK_1; pushSTACK(FF_FF_mal_FF(temp,temp)); # a'*a'
  603.         temp = STACK_1; temp = FF_FF_mal_FF(temp,temp); # b'*b'
  604.         temp = FF_FF_plus_FF(STACK_0,temp); # c' = a'*a'+b'*b'
  605.         skipSTACK(3);
  606.         temp = FF_sqrt_FF(temp); # c':=2^e*c'
  607.         return FF_I_scale_float_FF(temp,L_to_FN((sintL)e)); # 2^e*c'
  608.     } }}
  609.   local object DF_DF_hypot_DF(a,b)
  610.     var reg1 object a;
  611.     var reg2 object b;
  612.     { var reg5 sintWL a_exp;
  613.       var reg6 sintWL b_exp;
  614.       {# Exponenten von a holen:
  615.        var reg3 uintWL uexp = DF_uexp(TheDfloat(a)->float_value_semhi);
  616.        if (uexp==0) # a=0.0 -> liefere (abs b) :
  617.          { return (R_minusp(b) ? DF_minus_DF(b) : b); }
  618.        a_exp = (sintWL)(uexp - DF_exp_mid);
  619.       }
  620.       {# Exponenten von b holen:
  621.        var reg3 uintWL uexp = DF_uexp(TheDfloat(b)->float_value_semhi);
  622.        if (uexp==0) # b=0.0 -> liefere (abs a) :
  623.          { return (R_minusp(a) ? DF_minus_DF(a) : a); }
  624.        b_exp = (sintWL)(uexp - DF_exp_mid);
  625.       }
  626.        pushSTACK(a); pushSTACK(b);
  627.        # Stackaufbau: a, b.
  628.        # Nun a_exp = exponent(a), b_exp = exponent(b).
  629.       {var reg4 sintWL e = (a_exp > b_exp ? a_exp : b_exp); # Maximum der Exponenten
  630.        var reg3 object delta = L_to_FN(-(sintL)e); # -e als Fixnum
  631.        # a und b durch 2^e dividieren:
  632.        STACK_1 = (b_exp-a_exp > floor(DF_exp_mid-DF_exp_low-1,2) ? DF_0 : DF_I_scale_float_DF(STACK_1,delta));
  633.        STACK_0 = (a_exp-b_exp > floor(DF_exp_mid-DF_exp_low-1,2) ? DF_0 : DF_I_scale_float_DF(STACK_0,delta));
  634.        # Stackaufbau: a', b'.
  635.        # c' := a'*a'+b'*b' berechnen:
  636.        {var reg4 object temp;
  637.         temp = STACK_1; pushSTACK(DF_DF_mal_DF(temp,temp)); # a'*a'
  638.         temp = STACK_1; temp = DF_DF_mal_DF(temp,temp); # b'*b'
  639.         temp = DF_DF_plus_DF(STACK_0,temp); # c' = a'*a'+b'*b'
  640.         skipSTACK(3);
  641.         temp = DF_sqrt_DF(temp); # c':=2^e*c'
  642.         return DF_I_scale_float_DF(temp,L_to_FN((sintL)e)); # 2^e*c'
  643.     } }}
  644.   local object LF_LF_hypot_LF(a,b)
  645.     var reg1 object a;
  646.     var reg2 object b;
  647.     { var reg4 uintL a_exp;
  648.       var reg5 uintL b_exp;
  649.       {# Exponenten von a holen:
  650.        a_exp = TheLfloat(a)->expo;
  651.        if (a_exp==0) # a=0.0 -> liefere (abs b) :
  652.          { return (R_minusp(b) ? LF_minus_LF(b) : b); }
  653.       }
  654.       {# Exponenten von b holen:
  655.        b_exp = TheLfloat(b)->expo;
  656.        if (b_exp==0) # b=0.0 -> liefere (abs a) :
  657.          { return (R_minusp(a) ? LF_minus_LF(a) : a); }
  658.       }
  659.       pushSTACK(a); pushSTACK(b);
  660.       # Stackaufbau: a, b.
  661.       # Nun a_exp = exponent(a)+LF_exp_mid, b_exp = exponent(b)+LF_exp_mid.
  662.       {var reg6 uintL exp = (a_exp > b_exp ? a_exp : b_exp); # Maximum der Exponenten
  663.        var reg3 sintL e = (sintL)(exp-LF_exp_mid);
  664.        pushSTACK(L_to_I(e)); # e als Integer
  665.        pushSTACK(L_to_I(-e)); # -e als Integer
  666.       }
  667.       # a und b durch 2^e dividieren:
  668.       if ((b_exp>a_exp) && (b_exp-a_exp > (uintL)floor((uintL)(LF_exp_mid-LF_exp_low-1),2)))
  669.         { encode_LF0(TheLfloat(STACK_3)->len, STACK_3=); }
  670.         else
  671.         { STACK_3 = LF_I_scale_float_LF(STACK_3,STACK_0); }
  672.       if ((a_exp>b_exp) && (a_exp-b_exp > (uintL)floor((uintL)(LF_exp_mid-LF_exp_low-1),2)))
  673.         { encode_LF0(TheLfloat(STACK_2)->len, STACK_2=); }
  674.         else
  675.         { STACK_2 = LF_I_scale_float_LF(STACK_2,STACK_0); }
  676.       # Stackaufbau: a', b', e, -e.
  677.       # c' := a'*a'+b'*b' berechnen:
  678.       {var reg3 object temp;
  679.        temp = STACK_3; pushSTACK(LF_LF_mal_LF(temp,temp)); # a'*a'
  680.        temp = STACK_3; temp = LF_LF_mal_LF(temp,temp); # b'*b'
  681.        temp = LF_LF_plus_LF(STACK_0,temp); # c' = a'*a'+b'*b'
  682.        temp = LF_sqrt_LF(temp); # c':=2^e*c'
  683.        temp = LF_I_scale_float_LF(temp,STACK_2); # 2^e*c'
  684.        skipSTACK(5); return temp;
  685.     } }
  686.   local object R_R_hypot_R(a,b)
  687.     var reg3 object a;
  688.     var reg2 object b;
  689.     { if (R_rationalp(a))
  690.         { if (eq(a,Fixnum_0)) { return R_abs_R(b); } # a=0 -> (abs b)
  691.           if (R_rationalp(b))
  692.             # a,b beide rational
  693.             { var reg1 object temp;
  694.               if (eq(b,Fixnum_0)) { return R_abs_R(a); } # b=0 -> (abs a)
  695.               pushSTACK(a); pushSTACK(b);
  696.               temp = RA_RA_mal_RA(a,a); pushSTACK(temp); # a*a
  697.               temp = STACK_1; temp = RA_RA_mal_RA(temp,temp); # b*b
  698.               temp = RA_RA_plus_RA(STACK_0,temp); # a*a+b*b = c
  699.               skipSTACK(3);
  700.               return RA_sqrt_R(temp); # (sqrt c)
  701.             }
  702.             else
  703.             # a rational, b Float
  704.             { pushSTACK(b);
  705.               floatcase(b,
  706.                         { a = RA_to_SF(a); return SF_SF_hypot_SF(a,popSTACK()); },
  707.                         { a = RA_to_FF(a); return FF_FF_hypot_FF(a,popSTACK()); },
  708.                         { a = RA_to_DF(a); return DF_DF_hypot_DF(a,popSTACK()); },
  709.                         { a = RA_to_LF(a,TheLfloat(b)->len); return LF_LF_hypot_LF(a,popSTACK()); }
  710.                        );
  711.             }
  712.         }
  713.         else
  714.         { if (R_rationalp(b))
  715.             # a Float, b rational
  716.             { if (eq(b,Fixnum_0)) { return R_abs_R(a); } # b=0 -> (abs a)
  717.               pushSTACK(a);
  718.               floatcase(a,
  719.                         { b = RA_to_SF(b); return SF_SF_hypot_SF(popSTACK(),b); },
  720.                         { b = RA_to_FF(b); return FF_FF_hypot_FF(popSTACK(),b); },
  721.                         { b = RA_to_DF(b); return DF_DF_hypot_DF(popSTACK(),b); },
  722.                         { b = RA_to_LF(b,TheLfloat(a)->len); return LF_LF_hypot_LF(popSTACK(),b); }
  723.                        );
  724.             }
  725.             else
  726.             # a,b Floats
  727.             GEN_F_op2(a,b,SF_SF_hypot_SF,FF_FF_hypot_FF,DF_DF_hypot_DF,LF_LF_hypot_LF,1,0,return)
  728.         }
  729.     }
  730.  
  731. # N_abs_R(x) liefert (abs x), wo x eine Zahl ist.
  732. # kann GC auslösen
  733.   local object N_abs_R (object x);
  734. # Methode:
  735. # Falls x reell: klar
  736. # Falls x=a+bi: sqrt(a^2+b^2)
  737.   local object N_abs_R(x)
  738.     var reg1 object x;
  739.     { if (N_realp(x))
  740.         { return R_abs_R(x); }
  741.         else
  742.         { return R_R_hypot_R(TheComplex(x)->c_real,TheComplex(x)->c_imag); }
  743.     }
  744.   #define C_abs_R(x)  R_R_hypot_R(TheComplex(x)->c_real,TheComplex(x)->c_imag)
  745.  
  746. # N_signum_N(x) liefert (signum x), wo x eine Zahl ist.
  747. # kann GC auslösen
  748.   local object N_signum_N (object x);
  749. # Methode:
  750. # x reell -> klar.
  751. # x komplex -> falls (zerop x), x als Ergebnis, sonst (/ x (abs x)).
  752.   local object N_signum_N(x)
  753.     var reg1 object x;
  754.     { if (N_realp(x))
  755.         { return R_signum_R(x); }
  756.         else
  757.         { if (N_zerop(x)) return x;
  758.           pushSTACK(x); x = C_abs_R(x); # (abs x) errechnen
  759.           return N_N_durch_N(popSTACK(),x); # (/ x (abs x))
  760.     }   }
  761.  
  762. # N_sqrt_N(x) liefert (sqrt x), wo x eine beliebige Zahl ist.
  763. # kann GC auslösen
  764.   local object N_sqrt_N (object x);
  765. # Methode:
  766. # x reell -> Für x>=0 klar, für x<0: sqrt(-x)*i.
  767. # x=a+bi ->
  768. #   Bestimme r=abs(x)=sqrt(a*a+b*b).
  769. #   Falls a>=0: Setze c:=sqrt((r+a)/2), d:=(b/(2*c) falls c>0, c falls c=0).
  770. #   Falls a<0: Setze d:=sqrt((r-a)/2)*(1 falls b>=0, -1 falls b<0), c:=b/(2*d).
  771. #   Damit ist c>=0, 2*c*d=b, c*c=(r+a)/2, d*d=(r-a)/2, c*c-d*d=a, c*c+d*d=r,
  772. #   also c+di die gesuchte Wurzel.
  773.   local object N_sqrt_N(x)
  774.     var reg2 object x;
  775.     { if (N_realp(x))
  776.         # x reell
  777.         { if (!R_minusp(x))
  778.             { return R_sqrt_R(x); }
  779.             else
  780.             { return R_R_complex_C(Fixnum_0,R_sqrt_R(R_minus_R(x))); }
  781.         }
  782.         else
  783.         # x komplex
  784.         { var reg5 object a = TheComplex(x)->c_real;
  785.           var reg6 object b = TheComplex(x)->c_imag;
  786.           pushSTACK(b); pushSTACK(a);
  787.          {var reg4 object r = R_R_hypot_R(a,b); # r = (abs x)
  788.           var reg3 object a = STACK_0;
  789.           if (!R_minusp(a))
  790.             # a>=0
  791.             { var reg1 object c = # sqrt((r+a)/2)
  792.                 R_sqrt_R(R_R_durch_R(R_R_plus_R(r,a),fixnum(2)));
  793.               STACK_0 = c;
  794.               if (!R_zerop(c))
  795.                 { c = R_R_mal_R(fixnum(2),c); # 2*c
  796.                   c = R_R_durch_R(STACK_1,c); # d:=b/(2*c)
  797.                 }
  798.               c = R_R_complex_C(STACK_0,c); # c+di
  799.               skipSTACK(2); return c; # als Ergebnis
  800.             }
  801.             else
  802.             # a<0
  803.             { var reg1 object d = # sqrt((r-a)/2)
  804.                 R_sqrt_R(R_R_durch_R(R_R_minus_R(r,a),fixnum(2)));
  805.               if (R_mminusp(STACK_1)) { d = R_minus_R(d); } # bei b<0 Vorzeichenwechsel
  806.               STACK_0 = d;
  807.               d = R_R_mal_R(fixnum(2),d); # 2*d
  808.               d = R_R_durch_R(STACK_1,d); # c:=b/(2*d)
  809.               d = R_R_complex_C(d,STACK_0); # c+di
  810.               skipSTACK(2); return d; # als Ergebnis
  811.             }
  812.         }}
  813.     }
  814.  
  815. # N_N_gleich(x,y) vergleicht zwei reelle Zahlen x und y.
  816. # Ergebnis: TRUE falls x=y, FALSE sonst.
  817. # kann GC auslösen
  818.   global /* local */ boolean N_N_gleich (object x, object y);
  819. # Methode:
  820. # Falls beide reell, klar.
  821. # Falls x reell, y komplex: (= x (realpart y)) und (zerop (imagpart y)).
  822. # Falls x komplex, y reell: analog
  823. # Falls beide komplex: Realteile und Imaginärteile jeweils gleich?
  824.   global /* local */ boolean N_N_gleich(x,y)
  825.     var reg1 object x;
  826.     var reg2 object y;
  827.     { if (N_realp(x))
  828.         # x reell
  829.         { if (N_realp(y))
  830.             # x,y beide reell
  831.             { return (R_R_comp(x,y)==0 ? TRUE : FALSE); }
  832.             else
  833.             # x reell, y komplex
  834.             { if (!R_zerop(TheComplex(y)->c_imag)) return FALSE;
  835.               return (R_R_comp(x,TheComplex(y)->c_real)==0 ? TRUE : FALSE);
  836.             }
  837.         }
  838.         else
  839.         # x komplex
  840.         { if (N_realp(y))
  841.             # x komplex, y reell
  842.             { if (!R_zerop(TheComplex(x)->c_imag)) return FALSE;
  843.               return (R_R_comp(TheComplex(x)->c_real,y)==0 ? TRUE : FALSE);
  844.             }
  845.             else
  846.             # x,y beide komplex
  847.             { pushSTACK(TheComplex(x)->c_imag); pushSTACK(TheComplex(y)->c_imag);
  848.               x = TheComplex(x)->c_real; y = TheComplex(y)->c_real;
  849.               if (!(R_R_comp(x,y)==0)) { skipSTACK(2); return FALSE; } # Realteile vergleichen
  850.               y = popSTACK(); x = popSTACK();
  851.               if (!(R_R_comp(x,y)==0)) return FALSE; # Imaginärteile vergleichen
  852.               return TRUE;
  853.             }
  854.         }
  855.     }
  856.  
  857.