home *** CD-ROM | disk | FTP | other *** search
/ Amiga MA Magazine 1998 #6 / amigamamagazinepolishissue1998.iso / coders / jËzyki_programowania / clisp / src / archive / clisp.src.lha / src / comptran.d < prev    next >
Text File  |  1996-04-15  |  43KB  |  1,008 lines
  1. # Transzendente Funktionen für komplexe Zahlen
  2.  
  3. # N_phase_R(x) liefert (phase x), wo x eine Zahl ist.
  4. # Ergebnis rational nur wenn (= x 0) oder wenn x reell und >0.
  5. # kann GC auslösen
  6.   local object N_phase_R (object x);
  7. # Methode:
  8. # (= x 0) -> willkürliches Ergebnis 0
  9. # x reell -> Winkel von (x,0) in Polarkoordinaten
  10. # x komplex -> Winkel von ((realpart x),(imagpart x)) in Polarkoordinaten
  11.   local object N_phase_R(x)
  12.     var reg1 object x;
  13.     { if (N_zerop(x)) { return Fixnum_0; }
  14.       elif (N_realp(x)) { return R_R_atan_R(x,Fixnum_0); }
  15.       else { return R_R_atan_R(TheComplex(x)->c_real,TheComplex(x)->c_imag); }
  16.     }
  17.  
  18. # N_exp_N(x) liefert (exp x), wo x eine Zahl ist.
  19. # kann GC auslösen
  20.   local object N_exp_N (object x);
  21. # Methode:
  22. # x reell -> klar.
  23. # x = a+bi -> (exp a) mit (cos b) + i (sin b) multiplizieren:
  24. #             (complex (* (exp a) (cos b)) (* (exp a) (sin b)))
  25.   local object N_exp_N(x)
  26.     var reg1 object x;
  27.     { if (N_realp(x))
  28.         { return R_exp_R(x); }
  29.         else
  30.         # x=a+bi
  31.         { pushSTACK(TheComplex(x)->c_real); # a retten
  32.           R_cos_sin_R_R(TheComplex(x)->c_imag); # (cos b) und (sin b) errechnen
  33.           # Stackaufbau: a, cos(b), sin(b).
  34.           # Da b nicht = Fixnum 0 ist, ist auch sin(b) nicht = Fixnum 0.
  35.           STACK_2 = x = R_exp_R(STACK_2); # (exp a)
  36.           # Stackaufbau: exp(a), cos(b), sin(b).
  37.           STACK_0 = R_R_mal_R(x,STACK_0); # (* (exp a) (sin b)), nicht = Fixnum 0
  38.           x = R_R_mal_R(STACK_2,STACK_1); # (* (exp a) (cos b))
  39.           x = R_R_complex_C(x,STACK_0); # (complex ... ...)
  40.           skipSTACK(3); return x;
  41.         }
  42.     }
  43.  
  44. # N_log_N(x) liefert (log x), wo x eine Zahl ist.
  45. # kann GC auslösen
  46.   local object N_log_N (object x);
  47. # Methode:
  48. # (complex (log (abs x)) (phase x))
  49.   local object N_log_N(x)
  50.     var reg1 object x;
  51.     { pushSTACK(x); # x retten
  52.      {var reg2 object r = N_abs_R(x); # (abs x)
  53.       if (R_zerop(r)) { divide_0(); } # (abs x) = 0 -> Error
  54.       r = R_ln_R(r); # (log (abs x))
  55.       x = STACK_0; STACK_0 = r;
  56.       x = N_phase_R(x); # (phase x)
  57.       return R_R_complex_N(popSTACK(),x); # (complex r (phase x))
  58.     }}
  59.  
  60. # N_N_log_N(a,b) liefert (log a b), wo a und b Zahlen sind.
  61. # kann GC auslösen
  62.   local object N_N_log_N (object a, object b);
  63. # Methode:
  64. # (log a b) =
  65. #   falls b reell, >0:
  66. #     (complex (/ (log (abs a)) (log b)) (/ (phase a) (log b))), genauer:
  67. #     falls a reell, >0: bekannt
  68. #     falls (= a 0): Error
  69. #     sonst: (phase a) errechnen, ein Float.
  70. #            b (falls rational) ins selbe Float-Format umwandeln,
  71. #            Imaginärteil := (/ (phase a) (log dieses_b)).
  72. #            Falls a rational: (log (abs a) b).
  73. #            Falls a komplex mit rationalem Real- und Imaginärteil,
  74. #              Betragsquadrat  (expt (abs a) 2)  exakt ausrechnen als
  75. #              (+ (expt (realpart a) 2) (expt (imagpart a) 2)).
  76. #              Setze  Realteil := (/ (log Betragsquadrat b) 2).
  77. #              [Eventuell wird hierbei (log b) ein zweites Mal ausgerechnet,
  78. #               aber dies sowieso nur in Single-Precision.]
  79. #            Sonst bilde (abs a), ein Float, und (log (abs a)), ein Float,
  80. #              wandle b (falls rational) ins selbe Float-Format um,
  81. #              setze  Realteil := (/ (log (abs a)) (log dieses_b)).
  82. #   sonst: (/ (log a) (log b))
  83.   local object N_N_log_N(a,b)
  84.     var reg1 object a;
  85.     var reg2 object b;
  86.     { if (N_realp(b) && R_plusp(b))
  87.         # b ist reell und >0
  88.         { if (N_realp(a) && R_plusp(a))
  89.             # a und b sind beide reell und >0
  90.             { return R_R_log_R(a,b); }
  91.             else
  92.             # b ist reell und >0, a aber nicht
  93.             { pushSTACK(a); pushSTACK(b); # a,b retten
  94.               # Imaginärteil (/ (phase a) (log b)) errechnen:
  95.               {var reg3 object angle = N_phase_R(a); # (phase a)
  96.                if (eq(angle,Fixnum_0)) # = Fixnum 0 <==> (= a 0) -> Error
  97.                  { divide_0(); }
  98.                # durch (log b) dividieren, liefert den Imaginärteil:
  99.                pushSTACK(angle);
  100.                b = STACK_1; if (R_rationalp(b)) { b = RA_F_float_F(b,angle); }
  101.                b = F_ln_F(b); STACK_0 = F_F_durch_F(STACK_0,b);
  102.               }
  103.               # Stackaufbau: a, b, Imaginärteil.
  104.               # Realteil (/ (log (abs a)) (log b)) errechnen:
  105.               a = STACK_2;
  106.               if (N_realp(a))
  107.                 { if (R_rationalp(a))
  108.                     # a rational -> (log (abs a) b) errechnen:
  109.                     { a = R_abs_R(a); # Betrag (>0)
  110.                       pushSTACK(R_R_log_R(a,STACK_1));
  111.                       goto real_ok;
  112.                 }   }
  113.                 else
  114.                 { if (R_mrationalp(TheComplex(a)->c_real)
  115.                       && R_mrationalp(TheComplex(a)->c_imag)
  116.                      )
  117.                     # a komplex mit rationalem Real- und Imaginärteil a1,a2
  118.                     { # Betragsquadrat a1^2+a2^2 errechnen:
  119.                       pushSTACK(TheComplex(a)->c_imag);
  120.                       {var reg3 object a1 = TheComplex(a)->c_real;
  121.                        a1 = RA_RA_mal_RA(a1,a1); # a1*a1
  122.                        {var reg4 object a2 = STACK_0; STACK_0 = a1;
  123.                         a1 = RA_RA_mal_RA(a2,a2); # a2*a2
  124.                         a = RA_RA_plus_RA(STACK_0,a1);
  125.                       }}
  126.                       # davon der Logarithmus zur Basis b, durch 2:
  127.                       STACK_0 = R_R_durch_R(R_R_log_R(a,STACK_2),fixnum(2));
  128.                       goto real_ok;
  129.                 }   }
  130.               # Keine Chance für rationalen Realteil
  131.               pushSTACK(F_ln_F(N_abs_R(a))); # (log (abs a)), ein Float
  132.               # durch (log b) dividieren, liefert den Realteil:
  133.               b = STACK_2; if (R_rationalp(b)) { b = RA_F_float_F(b,STACK_0); }
  134.               b = F_ln_F(b); STACK_0 = F_F_durch_F(STACK_0,b);
  135.               real_ok:
  136.               # Stackaufbau: a, b, Imaginärteil, Realteil.
  137.              {var reg3 object erg = R_R_complex_C(STACK_0,STACK_1);
  138.               skipSTACK(4); return erg;
  139.             }}
  140.         }
  141.         else
  142.         # normaler komplexer Fall
  143.         { pushSTACK(b); a = N_log_N(a); b = STACK_0; # (log a)
  144.           STACK_0 = a; b = N_log_N(b); a = popSTACK(); # (log b)
  145.           return N_N_durch_N(a,b); # dividieren
  146.         }
  147.     }
  148.  
  149. # N_I_expt_N(x,y) = (expt x y), wo x eine Zahl und y ein Integer ist.
  150. # kann GC auslösen
  151.   local object N_I_expt_N (object x, object y);
  152.   # Methode:
  153.   # Für y>0:
  154.   #   a:=x, b:=y.
  155.   #   Solange b gerade, setze a:=a*a, b:=b/2. [a^b bleibt invariant, = x^y.]
  156.   #   c:=a.
  157.   #   Solange b:=floor(b/2) >0 ist,
  158.   #     setze a:=a*a, und falls b ungerade, setze c:=a*c.
  159.   #   Ergebnis c.
  160.   # Für y=0: Ergebnis 1.
  161.   # Für y<0: (/ (expt x (- y))).
  162.   local object N_I_expt_N(x,y)
  163.     var reg3 object x;
  164.     var reg1 object y;
  165.     { if (N_realp(x)) { return R_I_expt_R(x,y); } # x reell -> schnellere Routine
  166.       if (eq(y,Fixnum_0)) { return Fixnum_1; } # y=0 -> Ergebnis 1
  167.       pushSTACK(x);
  168.      {var reg4 boolean y_negative = FALSE;
  169.       if (R_minusp(y)) { y = I_minus_I(y); y_negative = TRUE; } # Betrag von y nehmen
  170.       # Nun ist y>0.
  171.       pushSTACK(y);
  172.       # Stackaufbau: a, b.
  173.       while (!I_oddp(y))
  174.         { var reg2 object a = STACK_1; STACK_1 = N_N_mal_N(a,a); # a:=a*a
  175.           STACK_0 = y = I_I_ash_I(STACK_0,Fixnum_minus1); # b := (ash b -1)
  176.         }
  177.       pushSTACK(STACK_1); # c:=a
  178.       # Stackaufbau: a, b, c.
  179.       until (eq(y=STACK_1,Fixnum_1)) # Solange b/=1
  180.         { STACK_1 = I_I_ash_I(y,Fixnum_minus1); # b := (ash b -1)
  181.          {var reg2 object a = STACK_2; STACK_2 = a = N_N_mal_N(a,a); # a:=a*a
  182.           if (I_oddp(STACK_1)) { STACK_0 = N_N_mal_N(a,STACK_0); } # evtl. c:=a*c
  183.         }}
  184.       x = STACK_0; skipSTACK(3);
  185.       # (expt x (abs y)) ist jetzt in x.
  186.       return (y_negative ? N_durch_N(x) : x); # evtl. noch Kehrwert nehmen
  187.     }}
  188.  
  189. # N_N_expt_N(x,y) = (expt x y), wo x und y Zahlen sind.
  190. # kann GC auslösen
  191.   local object N_N_expt_N (object x, object y);
  192.   # Methode:
  193.   # Falls y rational:
  194.   #   Falls y Integer:
  195.   #     Falls y=0: Ergebnis 1,
  196.   #       [Nach CLTL folgendermaßen:
  197.   #         x reell:
  198.   #           x rational -> Fixnum 1
  199.   #           x Float -> (float 1 x)
  200.   #         x komplex:
  201.   #           x komplex rational -> Fixnum 1
  202.   #           sonst: #C(1.0 0.0) im Float-Format des Real- bzw. Imaginärteils von x
  203.   #       ]
  204.   #     Falls x rational oder komplex rational oder |y| klein:
  205.   #       x^|y| durch wiederholtes Quadrieren und Multiplizieren und evtl.
  206.   #       Kehrwert-Bilden ermitteln.
  207.   #     Sonst wie bei 'y Float'.
  208.   #   Falls y Ratio m/n:
  209.   #     Es gilt (expt x m/n) = (expt (expt x 1/n) m).
  210.   #     Falls x in Q(i) liegt (also rational oder komplex rational ist):
  211.   #       Sollte x^(m/n) in Q(i) liegen, so auch eine n-te Wurzel x^(1/n)
  212.   #       (und bei n=2 oder n=4 damit auch alle n-ten Wurzeln x^(1/n) ).
  213.   #       Falls x rational >=0: n-te Wurzel aus x nehmen. Ist sie rational,
  214.   #         deren m-te Potenz als Ergebnis.
  215.   #       Falls x rational <=0 oder komplex rational und n Zweierpotenz:
  216.   #         n-te Wurzel aus x nehmen (mehrfaches sqrt). Ist sie rational oder
  217.   #         komplex rational, deren m-te Potenz als Ergebnis.
  218.   #         [Beliebige n betrachten!??]
  219.   #     Falls n Zweierpotenz und |m|,n klein: n-te Wurzel aus x nehmen
  220.   #       (mehrfaches sqrt), davon die m-te Potenz durch wiederholtes
  221.   #       Quadrieren und Multiplizieren und evtl. Kehrwert-Bilden.
  222.   #     Sonst wie bei 'y Float'.
  223.   # Falls y Float oder komplex:
  224.   #   Falls (zerop x):
  225.   #     Falls Realteil von y >0 :
  226.   #       liefere 0.0 falls x und y reell, #C(0.0 0.0) sonst.
  227.   #     Sonst Error.
  228.   #   Falls y=0.0:
  229.   #     liefere 1.0 falls x und y reell, #C(1.0 0.0) sonst.
  230.   #   Sonst: (exp (* (log x) y))
  231.   # Das Ergebnis liegt in Q(i), falls x in Q(i) liegt und 4y ein Integer ist.??
  232.   # Genauigkeit erhöhen, log2(|y|) Bits mehr??
  233.   # Bei x oder y rational und der andere Long-Float: bitte kein Single-Float!??
  234.   local object N_N_expt_N(x,y)
  235.     var reg2 object x;
  236.     var reg1 object y;
  237.     { if (N_realp(y) && R_rationalp(y))
  238.         # y rational
  239.         { if (RA_integerp(y))
  240.             # y Integer
  241.             { if (eq(y,Fixnum_0))
  242.                 #ifdef STRICT_CLTL
  243.                 # y=0 -> 1 im Format von x.
  244.                 { if (N_realp(x))
  245.                     { if (R_rationalp(x))
  246.                         { return Fixnum_1; }
  247.                         else
  248.                         { return I_F_float_F(Fixnum_1,x); }
  249.                     }
  250.                     else
  251.                     { x = R_R_contagion_R(TheComplex(x)->c_real,TheComplex(x)->c_imag);
  252.                       if (R_rationalp(x))
  253.                         { return Fixnum_1; }
  254.                         else
  255.                         { x = I_F_float_F(Fixnum_0,x); # 0.0
  256.                           pushSTACK(x);
  257.                           x = I_F_float_F(Fixnum_1,x); # 1.0
  258.                           return R_R_complex_C(x,popSTACK()); # #C(1.0 0.0)
  259.                     }   }
  260.                 }
  261.                 #else
  262.                 # Exponent exakt 0 -> Ergebnis exakt 1
  263.                 { return Fixnum_1; }
  264.                 #endif
  265.               if (I_fixnump(y)) { return N_I_expt_N(x,y); } # |y| klein ?
  266.               if (N_realp(x))
  267.                 { if (R_rationalp(x)) { return R_I_expt_R(x,y); } }
  268.                 else
  269.                 { if (R_rationalp(TheComplex(x)->c_real) && R_rationalp(TheComplex(x)->c_imag))
  270.                     { return N_I_expt_N(x,y); }
  271.                 }
  272.             }
  273.             else
  274.             # y Ratio
  275.             { if (N_realp(x))
  276.                 { if (R_rationalp(x))
  277.                     { if (R_minusp(x)) goto complex_rational;
  278.                       # x rational >=0
  279.                       pushSTACK(x); pushSTACK(y);
  280.                       {var reg1 object temp = RA_rootp(x,TheRatio(y)->rt_den); # n-te Wurzel versuchen
  281.                        if (!eq(temp,nullobj)) # Wurzel rational?
  282.                          {var reg1 object m = TheRatio(STACK_0)->rt_num;
  283.                           skipSTACK(2);
  284.                           return R_I_expt_R(temp,m); # (x^(1/n))^m
  285.                       }  }
  286.                       y = popSTACK(); x = popSTACK();
  287.                 }   }
  288.                 else
  289.                 { if (R_rationalp(TheComplex(x)->c_real) && R_rationalp(TheComplex(x)->c_imag))
  290.                     complex_rational: # x in Q(i)
  291.                     {var reg1 uintL k = I_power2p(TheRatio(y)->rt_den);
  292.                      if (!(k==0))
  293.                        # n Zweierpotenz = 2^(k-1). n>1, also k>1
  294.                        { pushSTACK(TheRatio(y)->rt_num); # m retten
  295.                          dotimespL(k,k-1, { x = N_sqrt_N(x); } ); # k-1 mal Quadratwurzel
  296.                          return N_I_expt_N(x,popSTACK()); # dann hoch m
  297.                 }   }  }
  298.               if (I_fixnump(TheRatio(y)->rt_num) # |m| klein
  299.                   && I_fixnump(TheRatio(y)->rt_den) # n klein
  300.                  )
  301.                 {var reg1 uintL n = posfixnum_to_L(TheRatio(y)->rt_den);
  302.                  if ((n & (n-1)) == 0) # n Zweierpotenz?
  303.                    { pushSTACK(TheRatio(y)->rt_num); # m retten
  304.                      until ((n = n>>1) ==0) { x = N_sqrt_N(x); } # n-te Wurzel ziehen
  305.                      return N_I_expt_N(x,popSTACK()); # dann hoch m
  306.                 }  }
  307.             }
  308.         }
  309.       # allgemeiner Fall (z.B. y Float oder komplex):
  310.       if (N_zerop(x)) # x=0.0 ?
  311.         { if (!R_plusp(N_realpart_R(y))) # Realteil von y <=0 ?
  312.             { divide_0(); } # ja -> Error
  313.           if (N_realp(x) && N_realp(y))
  314.             { x = R_R_contagion_R(x,y); # ein Float, da sonst x = Fixnum 0 gewesen wäre
  315.               return I_F_float_F(Fixnum_0,x); # 0.0
  316.             }
  317.             else
  318.             { if (!N_realp(x)) { x = R_R_contagion_R(TheComplex(x)->c_real,TheComplex(x)->c_imag); }
  319.               if (!N_realp(y)) { y = R_R_contagion_R(TheComplex(y)->c_real,TheComplex(y)->c_imag); }
  320.               x = R_R_contagion_R(x,y); # ein Float, da sonst x = Fixnum 0 gewesen wäre
  321.               x = I_F_float_F(Fixnum_0,x); # 0.0
  322.               return R_R_complex_C(x,x); # #C(0.0 0.0)
  323.             }
  324.         }
  325.       if (N_zerop(y)) # y=0.0 ?
  326.         { if (N_realp(x) && N_realp(y))
  327.             { x = R_R_contagion_R(x,y); # ein Float, da sonst y = Fixnum 0 gewesen wäre
  328.               return I_F_float_F(Fixnum_1,x); # 1.0
  329.             }
  330.             else
  331.             { if (!N_realp(x)) { x = R_R_contagion_R(TheComplex(x)->c_real,TheComplex(x)->c_imag); }
  332.               if (!N_realp(y)) { y = R_R_contagion_R(TheComplex(y)->c_real,TheComplex(y)->c_imag); }
  333.               x = R_R_contagion_R(x,y); # ein Float, da sonst y = Fixnum 0 gewesen wäre
  334.               x = I_F_float_F(Fixnum_0,x); # 0.0
  335.               pushSTACK(x);
  336.               x = I_F_float_F(Fixnum_1,x); # 1.0
  337.               return R_R_complex_C(x,popSTACK()); # #C(1.0 0.0)
  338.             }
  339.         }
  340.       pushSTACK(y);
  341.      {var reg1 object temp = N_log_N(x); # (log x)
  342.       temp = N_N_mal_N(temp,popSTACK()); # mal y
  343.       return N_exp_N(temp); # exp davon
  344.     }}
  345.  
  346. # N_sin_N(x) liefert (sin x), wo x eine Zahl ist.
  347. # kann GC auslösen
  348.   local object N_sin_N (object x);
  349. # Methode:
  350. # x reell -> klar
  351. # x = a+bi -> (complex (* (sin a) (cosh b)) (* (cos a) (sinh b)))
  352.   local object N_sin_N(x)
  353.     var reg1 object x;
  354.     { if (N_realp(x))
  355.         { return R_sin_R(x); }
  356.         else
  357.         # x=a+bi
  358.         { pushSTACK(TheComplex(x)->c_real); # a retten
  359.           R_cosh_sinh_R_R(TheComplex(x)->c_imag); # cosh(b), sinh(b) errechnen
  360.           # Stackaufbau: a, cosh(b), sinh(b).
  361.           # Da b nicht = Fixnum 0 ist, ist auch sinh(b) nicht = Fixnum 0.
  362.           R_cos_sin_R_R(STACK_2); # cos(a), sin(a) errechnen, cos(a) /= Fixnum 0
  363.           # Stackaufbau: a, cosh(b), sinh(b), cos(a), sin(a).
  364.           STACK_0 = R_R_mal_R(STACK_0,STACK_3); # sin(a)*cosh(b)
  365.          {var reg2 object erg = R_R_mal_R(STACK_1,STACK_2); # cos(a)*sinh(b), nicht Fixnum 0
  366.           erg = R_R_complex_C(STACK_0,erg); skipSTACK(5); return erg;
  367.         }}
  368.     }
  369.  
  370. # N_cos_N(x) liefert (cos x), wo x eine Zahl ist.
  371. # kann GC auslösen
  372.   local object N_cos_N (object x);
  373. # Methode:
  374. # x reell -> klar
  375. # x = a+bi -> (complex (* (cos a) (cosh b)) (- (* (sin a) (sinh b))))
  376.   local object N_cos_N(x)
  377.     var reg1 object x;
  378.     { if (N_realp(x))
  379.         { return R_cos_R(x); }
  380.         else
  381.         # x=a+bi
  382.         { pushSTACK(TheComplex(x)->c_real); # a retten
  383.           R_cosh_sinh_R_R(TheComplex(x)->c_imag); # cosh(b), sinh(b) errechnen
  384.           # Stackaufbau: a, cosh(b), sinh(b).
  385.           R_cos_sin_R_R(STACK_2); # cos(a), sin(a) errechnen
  386.           # Stackaufbau: a, cosh(b), sinh(b), cos(a), sin(a).
  387.           STACK_0 = R_minus_R(R_R_mal_R(STACK_0,STACK_2)); # -sin(a)*sinh(b)
  388.          {var reg2 object erg = R_R_mal_R(STACK_1,STACK_3); # cos(a)*cosh(b)
  389.           erg = R_R_complex_N(erg,STACK_0); skipSTACK(5); return erg;
  390.         }}
  391.     }
  392.  
  393. # N_tan_N(x) liefert (tan x), wo x eine Zahl ist.
  394. # kann GC auslösen
  395.   local object N_tan_N (object x);
  396. # Methode:
  397. # x reell -> (/ (sin x) (cos x))
  398. # x = a+bi -> (/ (complex (* (sin a) (cosh b)) (* (cos a) (sinh b)))
  399. #                (complex (* (cos a) (cosh b)) (- (* (sin a) (sinh b)))) )
  400.   local object N_tan_N(x)
  401.     var reg2 object x;
  402.     { if (N_realp(x))
  403.         { R_cos_sin_R_R(x);
  404.           # Stackaufbau: cos(x), sin(x).
  405.          {var reg1 object erg = R_R_durch_R(STACK_0,STACK_1);
  406.           skipSTACK(2); return erg;
  407.         }}
  408.         else
  409.         # x=a+bi
  410.         { pushSTACK(TheComplex(x)->c_real); # a retten
  411.           R_cosh_sinh_R_R(TheComplex(x)->c_imag); # cosh(b), sinh(b) errechnen
  412.           # Stackaufbau: a, cosh(b), sinh(b).
  413.           R_cos_sin_R_R(STACK_2); # cos(a), sin(a) errechnen
  414.           # Stackaufbau: a, cosh(b), sinh(b), cos(a), sin(a).
  415.          {var reg1 object temp;
  416.           STACK_4 = R_R_mal_R(STACK_0,STACK_3); # sin(a)*cosh(b)
  417.           temp = R_R_mal_R(STACK_1,STACK_2); # cos(a)*sinh(b) /= Fixnum 0
  418.           STACK_4 = R_R_complex_C(STACK_4,temp); # Zähler
  419.           # Stackaufbau: Zähler, cosh(b), sinh(b), cos(a), sin(a).
  420.           STACK_3 = R_R_mal_R(STACK_1,STACK_3); # cos(a)*cosh(b)
  421.           temp = R_minus_R(R_R_mal_R(STACK_0,STACK_2)); # -sin(a)*sinh(b)
  422.           temp = R_R_complex_N(STACK_3,temp); # Nenner
  423.           temp = N_N_durch_N(STACK_4,temp); # Zähler/Nenner
  424.           skipSTACK(5); return temp;
  425.         }}
  426.     }
  427.  
  428. # N_cis_N(x) liefert (cis x), wo x eine Zahl ist.
  429. # kann GC auslösen
  430.   local object N_cis_N (object x);
  431. # Methode:
  432. # x reell -> (complex (cos x) (sin x))
  433. # x = a+bi -> (complex (* (exp (- b)) (cos a)) (* (exp (- b)) (sin a)))
  434.   local object N_cis_N(x)
  435.     var reg1 object x;
  436.     { if (N_realp(x))
  437.         { R_cos_sin_R_R(x);
  438.           # Stackaufbau: cos(x), sin(x).
  439.          {var reg1 object erg = R_R_complex_N(STACK_1,STACK_0);
  440.           skipSTACK(2); return erg;
  441.         }}
  442.         else
  443.         # x=a+bi
  444.         { pushSTACK(TheComplex(x)->c_imag); # b retten
  445.           R_cos_sin_R_R(TheComplex(x)->c_real); # (cos a) und (sin a) errechnen
  446.           # Stackaufbau: b, cos(a), sin(a).
  447.           STACK_2 = x = R_exp_R(R_minus_R(STACK_2)); # (exp (- b))
  448.           # Stackaufbau: exp(-b), cos(a), sin(a).
  449.           STACK_0 = R_R_mal_R(x,STACK_0); # (* (exp (- b)) (sin a))
  450.           x = R_R_mal_R(STACK_2,STACK_1); # (* (exp (- b)) (cos a))
  451.           x = R_R_complex_N(x,STACK_0); # (complex ... ...)
  452.           skipSTACK(3); return x;
  453.         }
  454.     }
  455.  
  456. # N_sinh_N(x) liefert (sinh x), wo x eine Zahl ist.
  457. # kann GC auslösen
  458.   local object N_sinh_N (object x);
  459. # Methode:
  460. # x reell -> klar
  461. # x = a+bi -> (complex (* (sinh a) (cos b)) (* (cosh a) (sin b)))
  462.   local object N_sinh_N(x)
  463.     var reg1 object x;
  464.     { if (N_realp(x))
  465.         { return R_sinh_R(x); }
  466.         else
  467.         # x=a+bi
  468.         { pushSTACK(TheComplex(x)->c_real); # a retten
  469.           R_cos_sin_R_R(TheComplex(x)->c_imag); # cos(b), sin(b) errechnen
  470.           # Stackaufbau: a, cos(b), sin(b).
  471.           # Da b nicht = Fixnum 0 ist, ist auch sin(b) nicht = Fixnum 0.
  472.           R_cosh_sinh_R_R(STACK_2); # cosh(a), sinh(a) errechnen, cosh(a) /= Fixnum 0
  473.           # Stackaufbau: a, cos(b), sin(b), cosh(a), sinh(a).
  474.           STACK_0 = R_R_mal_R(STACK_0,STACK_3); # sinh(a)*cos(b)
  475.          {var reg2 object erg = R_R_mal_R(STACK_1,STACK_2); # cosh(a)*sin(b), nicht Fixnum 0
  476.           erg = R_R_complex_C(STACK_0,erg); skipSTACK(5); return erg;
  477.         }}
  478.     }
  479.  
  480. # N_cosh_N(x) liefert (cosh x), wo x eine Zahl ist.
  481. # kann GC auslösen
  482.   local object N_cosh_N (object x);
  483. # Methode:
  484. # x reell -> klar
  485. # x = a+bi -> (complex (* (cosh a) (cos b)) (* (sinh a) (sin b)))
  486.   local object N_cosh_N(x)
  487.     var reg1 object x;
  488.     { if (N_realp(x))
  489.         { return R_cosh_R(x); }
  490.         else
  491.         # x=a+bi
  492.         { pushSTACK(TheComplex(x)->c_real); # a retten
  493.           R_cos_sin_R_R(TheComplex(x)->c_imag); # cos(b), sin(b) errechnen
  494.           # Stackaufbau: a, cos(b), sin(b).
  495.           R_cosh_sinh_R_R(STACK_2); # cosh(a), sinh(a) errechnen
  496.           # Stackaufbau: a, cos(b), sin(b), cosh(a), sinh(a).
  497.           STACK_0 = R_R_mal_R(STACK_0,STACK_2); # sinh(a)*sin(b)
  498.          {var reg2 object erg = R_R_mal_R(STACK_1,STACK_3); # cosh(a)*cos(b)
  499.           erg = R_R_complex_N(erg,STACK_0); skipSTACK(5); return erg;
  500.         }}
  501.     }
  502.  
  503. # N_tanh_N(x) liefert (tanh x), wo x eine Zahl ist.
  504. # kann GC auslösen
  505.   local object N_tanh_N (object x);
  506. # Methode:
  507. # x reell -> (/ (sinh x) (cosh x))
  508. # x = a+bi -> (/ (complex (* (sinh a) (cos b)) (* (cosh a) (sin b)))
  509. #                (complex (* (cosh a) (cos b)) (* (sinh a) (sin b))) )
  510.   local object N_tanh_N(x)
  511.     var reg2 object x;
  512.     { if (N_realp(x))
  513.         { R_cosh_sinh_R_R(x);
  514.           # Stackaufbau: cosh(x), sinh(x).
  515.          {var reg1 object erg = R_R_durch_R(STACK_0,STACK_1);
  516.           skipSTACK(2); return erg;
  517.         }}
  518.         else
  519.         # x=a+bi
  520.         { pushSTACK(TheComplex(x)->c_real); # a retten
  521.           R_cos_sin_R_R(TheComplex(x)->c_imag); # cos(b), sin(b) errechnen
  522.           # Stackaufbau: a, cos(b), sin(b).
  523.           R_cosh_sinh_R_R(STACK_2); # cosh(a), sinh(a) errechnen
  524.           # Stackaufbau: a, cos(b), sin(b), cosh(a), sinh(a).
  525.          {var reg1 object temp;
  526.           STACK_4 = R_R_mal_R(STACK_0,STACK_3); # sinh(a)*cos(b)
  527.           temp = R_R_mal_R(STACK_1,STACK_2); # cosh(a)*sin(b) /= Fixnum 0
  528.           STACK_4 = R_R_complex_C(STACK_4,temp); # Zähler
  529.           # Stackaufbau: Zähler, cos(b), sin(b), cosh(a), sinh(a).
  530.           STACK_3 = R_R_mal_R(STACK_1,STACK_3); # cosh(a)*cos(b)
  531.           temp = R_R_mal_R(STACK_0,STACK_2); # sinh(a)*sin(b)
  532.           temp = R_R_complex_N(STACK_3,temp); # Nenner
  533.           temp = N_N_durch_N(STACK_4,temp); # Zähler/Nenner
  534.           skipSTACK(5); return temp;
  535.         }}
  536.     }
  537.  
  538. # N_atanh_N(z) liefert den Artanh einer Zahl z.
  539. # kann GC auslösen
  540.   local object N_atanh_N (object z);
  541. # Methode:
  542. # Wert und Branch Cuts nach der Formel CLTL2, S. 315:
  543. #   artanh(z) = (log(1+z)-log(1-z)) / 2
  544. # Sei z=x+iy, Ergebnis u+iv.
  545. # Falls x=0 und y=0: u=0, v=0.
  546. # Falls x=0: u = 0, v = atan(X=1,Y=y).
  547. # Falls y=0:
  548. #   x rational -> x in Float umwandeln.
  549. #   |x|<1/2: u = atanh(x), v = 0.
  550. #   |x|>=1/2: (1+x)/(1-x) errechnen,
  551. #             =0 -> Error,
  552. #             >0 (also |x|<1) -> u = 1/2 log((1+x)/(1-x)), v = 0.
  553. #             <0 (also |x|>1) -> u = 1/2 log(-(1+x)/(1-x)),
  554. #                                v = (-pi/2 für x>1, pi/2 für x<-1).
  555. # Sonst:
  556. #   1+x und 1-x errechnen.
  557. #   x und y in Floats umwandeln.
  558. #   |4x| und 1+x^2+y^2 errechnen,
  559. #   |4x| < 1+x^2+y^2 -> u = 1/2 atanh(2x/(1+x^2+y^2)),
  560. #   |4x| >= 1+x^2+y^2 -> u = 1/4 ln ((1+x^2+y^2)+2x)/((1+x^2+y^2)-2x)
  561. #                        oder besser (an der Singularität: |x|-1,|y| klein):
  562. #                        u = 1/4 ln ((1+x)^2+y^2)/((1-x)^2+y^2).
  563. #   v = 1/2 atan(X=(1-x)(1+x)-y^2,Y=2y) * (-1 falls Y=0.0 und X<0.0 und x>=0.0,
  564. #                                          1 sonst)
  565. # Ergebnis ist reell nur, wenn z reell.
  566. # Real- und Imaginärteil des Ergebnisses sind Floats, außer wenn z reell oder
  567. # rein imaginär ist.
  568.  
  569. # N_atan_N(z) liefert den Arctan einer Zahl z.
  570. # kann GC auslösen
  571.   local object N_atan_N (object z);
  572. # Methode:
  573. # Wert und Branch Cuts nach der Formel CLTL2, S. 307/312/313:
  574. #   arctan(z) = (log(1+iz)-log(1-iz)) / 2i
  575. # Sei z=x+iy, errechne u+iv = artanh(-y+ix) wie oben, Ergebnis v-iu.
  576. # Real- und Imaginärteil des Ergebnisses sind Floats, außer wenn z reell oder
  577. # rein imaginär ist.
  578.  
  579. # Hilfsfunktion für beide: u+iv := artanh(x+iy), u,v beide auf den Stack.
  580.   local void R_R_atanh_R_R (object x, object y);
  581.   local void R_R_atanh_R_R(x,y)
  582.     var reg2 object x;
  583.     var reg3 object y;
  584.     { if (eq(x,Fixnum_0)) # x=0 -> u=0, v=atan(X=1,Y=y) (Fall y=0 ist inbegriffen)
  585.         { pushSTACK(x); pushSTACK(R_R_atan_R(Fixnum_1,y)); return; }
  586.       if (eq(y,Fixnum_0))
  587.         { if (R_rationalp(x)) { x = RA_float_F(x); } # x in Float umwandeln
  588.           # x Float
  589.           if (R_zerop(x)) # x=0.0 -> x als Ergebnis
  590.             { pushSTACK(x); pushSTACK(Fixnum_0); return; }
  591.           if (F_exponent_L(x) < 0)
  592.             # Exponent e<0, also |x|<1/2
  593.             { pushSTACK(F_atanhx_F(x)); pushSTACK(Fixnum_0); return; }
  594.           # e>=0, also |x|>=1/2
  595.           pushSTACK(x);
  596.           pushSTACK(R_R_minus_R(Fixnum_1,x)); # 1-x
  597.           # Stackaufbau: x, 1-x.
  598.           {var reg1 object temp;
  599.            temp = R_R_plus_R(Fixnum_1,STACK_1); # 1+x
  600.            temp = F_F_durch_F(temp,STACK_0); # (1+x)/(1-x)
  601.            if (!R_minusp(temp))
  602.              { STACK_1 = temp; STACK_0 = Fixnum_0; # Imaginärteil:=0
  603.                if (R_zerop(temp)) { divide_0(); } # x = -1 -> Error
  604.              }
  605.              else
  606.              # (1+x)/(1-x) < 0 -> Betrag nehmen, Imaginärteil berechnen:
  607.              { STACK_1 = F_minus_F(temp);
  608.                temp = F_I_scale_float_F(pi(),Fixnum_minus1); # (scale-float pi -1) = pi/2
  609.                if (R_mminusp(STACK_0)) { temp = F_minus_F(temp); } # 1-x<0 -> dann -pi/2
  610.                STACK_0 = temp;
  611.              }
  612.            # Stackaufbau: |(1+x)/(1-x)| (>0), Imaginärteil.
  613.            STACK_1 = F_I_scale_float_F(R_ln_R(STACK_1),Fixnum_minus1); # ln bilden, durch 2
  614.            return;
  615.         } }
  616.       pushSTACK(x); pushSTACK(y);
  617.       pushSTACK(R_R_plus_R(Fixnum_1,x)); # 1+x
  618.       pushSTACK(R_R_minus_R(Fixnum_1,STACK_2)); # 1-x
  619.       # Stackaufbau: x, y, 1+x, 1-x.
  620.       # x und y in Floats umwandeln:
  621.       if (R_mrationalp(STACK_3))
  622.         { if (R_mrationalp(STACK_2)) { STACK_2 = RA_float_F(STACK_2); }
  623.           STACK_3 = RA_F_float_F(STACK_3,STACK_2);
  624.         }
  625.         else
  626.         { if (R_mrationalp(STACK_2)) { STACK_2 = RA_F_float_F(STACK_2,STACK_3); } }
  627.       {var reg1 object temp = STACK_2; pushSTACK(R_R_mal_R(temp,temp)); } # y^2
  628.       {var reg1 object temp = STACK_4; temp = R_R_mal_R(temp,temp); # x^2
  629.        pushSTACK(R_R_plus_R(Fixnum_1,R_R_plus_R(temp,STACK_0))); # 1+x^2+y^2
  630.       }
  631.       # Stackaufbau: x, y, 1+x, 1-x, y^2, 1+x^2+y^2.
  632.       {var reg1 object temp = # |4x|
  633.          F_abs_F(F_I_scale_float_F(STACK_5,fixnum(2)));
  634.        if (F_F_comp(temp,STACK_0) < 0) # |4x| < 1+x^2+y^2 ?
  635.          { temp = F_I_scale_float_F(STACK_5,Fixnum_1); # 2x
  636.            temp = F_F_durch_F(temp,STACK_0); # 2x/(1+x^2+y^2)
  637.            temp = F_atanhx_F(temp); # atanh davon
  638.            temp = F_I_scale_float_F(temp,Fixnum_minus1); # .../2 =: u
  639.          }
  640.          else
  641.          { temp = STACK_3; temp = R_R_mal_R(temp,temp); # (1+x)^2
  642.            STACK_0 = R_R_plus_R(temp,STACK_1); # (1+x)^2+y^2, ein Float >=0
  643.            temp = STACK_2; temp = R_R_mal_R(temp,temp); # (1-x)^2
  644.            temp = R_R_plus_R(temp,STACK_1); # (1-x)^2+y^2, ein Float >=0
  645.            temp = F_F_durch_F(STACK_0,temp); # ((1+x)^2+y^2)/((1-x)^2+y^2), ein Float >=0
  646.            if (R_zerop(temp)) { divide_0(); } # sollte >0 sein
  647.            temp = R_ln_R(temp); # ln davon, ein Float
  648.            temp = F_I_scale_float_F(temp,sfixnum(-2)); # .../4 =: u
  649.          }
  650.        x = STACK_5; # x retten (Vorzeichen ist GC-invariant!)
  651.        STACK_5 = temp;
  652.        # Stackaufbau: u, y, 1+x, 1-x, y^2, -.
  653.        temp = R_R_mal_R(STACK_3,STACK_2); # (1+x)(1-x)
  654.        STACK_0 = R_R_minus_R(temp,STACK_1); # (1+x)(1-x)-y^2, ein Float
  655.        temp = F_I_scale_float_F(STACK_4,Fixnum_1); # 2y, ein Float
  656.        temp = R_R_atan_R(STACK_0,temp); # atan(X=(1-x)(1+x)-y^2,Y=2y), ein Float
  657.        if (R_mminusp(STACK_0) && !R_minusp(x) && R_zerop(STACK_4)) # X<0.0 und x>=0.0 und Y=0.0 ?
  658.          { temp = F_minus_F(temp); } # ja -> Vorzeichenwechsel
  659.        STACK_4 = F_I_scale_float_F(temp,Fixnum_minus1); # .../2 =: v
  660.       }
  661.       # Stackaufbau: u, v, 1+x, 1-x, y^2, -.
  662.       skipSTACK(4); return;
  663.     }
  664.   #
  665.   local object N_atanh_N(z)
  666.     var reg1 object z;
  667.     { if (N_realp(z))
  668.         { R_R_atanh_R_R(z,Fixnum_0); }
  669.         else
  670.         { R_R_atanh_R_R(TheComplex(z)->c_real,TheComplex(z)->c_imag); }
  671.       # Stackaufbau: u, v.
  672.       z = R_R_complex_N(STACK_1,STACK_0); skipSTACK(2); return z;
  673.     }
  674.   #
  675.   local object N_atan_N(z)
  676.     var reg1 object z;
  677.     { # atanh(iz) errechnen:
  678.       if (N_realp(z))
  679.         { R_R_atanh_R_R(Fixnum_0,z); }
  680.         else
  681.         { pushSTACK(TheComplex(z)->c_real);
  682.           z = R_minus_R(TheComplex(z)->c_imag);
  683.           R_R_atanh_R_R(z,popSTACK());
  684.         }
  685.       # Stackaufbau: u, v.
  686.       z = R_minus_R(STACK_1); z = R_R_complex_N(STACK_0,z); # z := v-iu
  687.       skipSTACK(2); return z;
  688.     }
  689.  
  690. # Um für zwei Zahlen u,v mit u^2-v^2=1 und u,v beide in Bild(sqrt)
  691. # (d.h. Realteil>0.0 oder Realteil=0.0 und Imaginärteil>=0.0)
  692. # log(u+v) zu berechnen:
  693. #               log(u+v) = 2 artanh(v/(u+1))                            (!)
  694. # (Beweis: 2 artanh(v/(u+1)) = log(1+(v/(u+1))) - log(1-(v/(u+1)))
  695. #  = log((1+u+v)/(u+1)) - log((1+u-v)/(u+1)) == log((1+u+v)/(1+u-v))
  696. #  = log(u+v) mod 2 pi i, und beider Imaginärteil ist > -pi und <= pi.)
  697.  
  698. # N_asinh_N(z) liefert den Arsinh einer Zahl z.
  699. # kann GC auslösen
  700.   local object N_asinh_N (object z);
  701. # Methode:
  702. # Wert und Branch Cuts nach der Formel CLTL2, S. 313:
  703. #   arsinh(z) = log(z+sqrt(1+z^2))
  704. # z=x+iy, Ergebnis u+iv.
  705. # Falls x=0 und y=0: u=0, v=0.
  706. # Falls x=0: arsinh(iy) = i arcsin(y).
  707. #   y rational ->
  708. #     Bei y=1: u = 0, v = pi/2.
  709. #     Bei y=1/2: u = 0, v = pi/6.
  710. #     Bei y=0: u = 0, v = 0.
  711. #     Bei y=-1/2: u = 0, v = -pi/6.
  712. #     Bei y=-1: u = 0, v = -pi/2.
  713. #     Sonst y in Float umwandeln.
  714. #   e := Exponent aus (decode-float y), d := (float-digits y)
  715. #   Bei y=0.0 oder e<=-d/2 liefere u = 0, v = y
  716. #     (denn bei e<=-d/2 ist y^2/3 < y^2/2 < 2^(-d)/2 = 2^(-d-1), also
  717. #     1 <= asin(y)/y < 1+y^2/3 < 1+2^(-d-1) < 1+2^(-d),
  718. #     also ist asin(y)/y, auf d Bits gerundet, gleich 1.0).
  719. #   Berechne 1-y^2.
  720. #   Bei y>1 liefere  u = ln(y+sqrt(y^2-1)), v = pi/2.
  721. #   Bei y<-1 liefere  u = -ln(|y|+sqrt(|y|^2-1)), v = -pi/2.
  722. #   Bei |y|<=1 liefere  u = 0, v = atan(X=sqrt(1-y^2),Y=y).
  723. # Falls y=0:
  724. #   x rational -> x in Float umwandeln.
  725. #   |x|<1/2: u = atanh(x/sqrt(1+x^2)),
  726. #   x>=1/2: u = ln(x+sqrt(1+x^2)),
  727. #   x<=-1/2: u = -ln(-x+sqrt(1+x^2)).
  728. #   v = 0.
  729. # Sonst:
  730. #   z in Bild(sqrt) -> log(sqrt(1+z^2)+z) = (!) = 2 artanh(z/(1+sqrt(1+z^2))).
  731. #   z nicht in Bild(sqrt) ->
  732. #     arsinh(z) = -arsinh(-z).
  733. #     (Denn arsinh(z)+arsinh(-z) == log((z+sqrt(1+z^2))(-z+sqrt(1+z^2)))
  734. #           = log((1+z^2)-z^2) = log(1) = 0 mod 2 pi i, und links ist
  735. #      der Imaginärteil betragsmäßig <=pi.)
  736. #     Also arsinh(z) = -arsinh(-z) = - 2 artanh(-z/(1+sqrt(1+z^2)))
  737. #          = (wegen -artanh(-w) = artanh(w)) = 2 artanh(z/(1+sqrt(1+z^2))).
  738. # Real- und Imaginärteil des Ergebnisses sind Floats, außer wenn z reell oder
  739. # rein imaginär ist.
  740.  
  741. # N_asin_N(z) liefert den Arcsin einer Zahl z.
  742. # kann GC auslösen
  743.   local object N_asin_N (object z);
  744. # Methode:
  745. # Wert und Branch Cuts nach der Formel CLTL2, S. 311:
  746. #   arcsin(z) = log(iz+sqrt(1-z^2))/i
  747. # Sei z=x+iy, errechne u+iv = arsinh(-y+ix) wie oben, Ergebnis v-iu.
  748. # Real- und Imaginärteil des Ergebnisses sind Floats, außer wenn z reell oder
  749. # rein imaginär ist.
  750.  
  751. # Hilfsfunktion für beide: u+iv := arsinh(x+iy), u,v beide auf den Stack.
  752.   local void R_R_asinh_R_R (object x, object y);
  753.   local void R_R_asinh_R_R(x,y)
  754.     var reg2 object x;
  755.     var reg3 object y;
  756.     { if (eq(x,Fixnum_0)) # x=0 ?
  757.         { pushSTACK(x); pushSTACK(y);
  758.           if (R_rationalp(y))
  759.             # y rational
  760.             { if (eq(y,Fixnum_0)) { return; } # x=0, y=0 -> u=0, v=0 bereits im Stack
  761.               if (RA_integerp(y))
  762.                 # y Integer
  763.                 { if (eq(y,Fixnum_1)) # x=0, y=1 -> v = pi/2
  764.                     { STACK_0 = F_I_scale_float_F(pi(),Fixnum_minus1); return; }
  765.                   if (eq(y,Fixnum_minus1)) # x=0, y=-1 -> v = -pi/2
  766.                     { STACK_0 = F_minus_F(F_I_scale_float_F(pi(),Fixnum_minus1)); return; }
  767.                   STACK_0 = y = I_float_F(y); # y in Float umwandeln
  768.                 }
  769.                 else
  770.                 # y Ratio
  771.                 { if (eq(TheRatio(y)->rt_den,fixnum(2))) # Nenner = 2 ?
  772.                     { var reg1 object temp = TheRatio(y)->rt_num; # Zähler
  773.                       if (eq(temp,Fixnum_1)) # x=0, y=1/2 -> v = pi/6
  774.                         { STACK_0 = R_R_durch_R(pi(),fixnum(6)); return; }
  775.                       if (eq(temp,Fixnum_minus1)) # x=0, y=-1/2 -> v = -pi/6
  776.                         { STACK_0 = F_minus_F(R_R_durch_R(pi(),fixnum(6))); return; }
  777.                     }
  778.                   STACK_0 = y = RA_float_F(y); # y in Float umwandeln
  779.             }   }
  780.           # y Float
  781.           if (R_zerop(y) # y=0.0 -> arcsin(y) = y als Ergebnis
  782.               || (F_exponent_L(y) <= (sintL)(-F_float_digits(y))>>1) # e <= -d/2 <==> e <= -ceiling(d/2) ?
  783.              )
  784.             { return; } # u=0, v=y bereits im Stack
  785.           # Stackaufbau: 0, y.
  786.           {var reg1 object temp = R_R_minus_R(Fixnum_1,F_F_mal_F(y,y)); # 1-y*y
  787.            if (!R_minusp(temp))
  788.              # 1-y*y>=0, also |y|<=1
  789.              { temp = F_sqrt_F(temp); # sqrt(1-y*y)
  790.                STACK_0 = R_R_atan_R(temp,STACK_0); # v = atan(X=sqrt(1-y*y),Y=y)
  791.              }
  792.              else
  793.              # 1-y*y<0, also |y|>1
  794.              { temp = F_sqrt_F(F_minus_F(temp)); # sqrt(y*y-1)
  795.                y = STACK_0; # |y| zu temp addieren:
  796.                if (R_minusp(y)) { temp = F_F_minus_F(temp,y); } else { temp = F_F_plus_F(temp,y); }
  797.                # temp = sqrt(y^2-1)+|y|, ein Float >1
  798.                STACK_1 = R_ln_R(temp); # ln(|y|+sqrt(y^2-1)), ein Float >0
  799.                temp = F_I_scale_float_F(pi(),Fixnum_minus1); # (scale-float pi -1) = pi/2
  800.                if (!R_mminusp(STACK_0)) # Vorzeichen von y
  801.                  # y>1 -> v = pi/2
  802.                  { STACK_0 = temp; }
  803.                  else
  804.                  # y<-1 -> v = -pi/2, u = -ln(...)
  805.                  { STACK_0 = F_minus_F(temp); STACK_1 = F_minus_F(STACK_1); }
  806.              }
  807.            return;
  808.         } }
  809.       if (eq(y,Fixnum_0)) # y=0 ?
  810.         { if (R_rationalp(x)) { x = RA_float_F(x); } # x in Float umwandeln
  811.           # x Float
  812.           pushSTACK(x); pushSTACK(Fixnum_0); # x retten, v = 0
  813.           if (R_zerop(x)) { return; } # x=0.0 -> u=x, v=0.
  814.           {var reg1 object temp = # sqrt(1+x^2)
  815.              F_sqrt_F(R_R_plus_R(Fixnum_1,F_F_mal_F(x,x)));
  816.            x = STACK_1;
  817.            if (F_exponent_L(x) < 0) # Exponent e (von x/=0) <0 ?
  818.              # |x|<1/2
  819.              { STACK_1 = F_atanhx_F(F_F_durch_F(x,temp)); } # u = atanh(x/sqrt(1+x^2))
  820.              else
  821.              # |x|>=1/2
  822.              { if (!R_minusp(x))
  823.                  # x>=1
  824.                  { STACK_1 = R_ln_R(F_F_plus_F(temp,x)); } # u = ln(x+sqrt(1+x^2))
  825.                  else
  826.                  # x<=-1
  827.                  { STACK_1 = F_minus_F(R_ln_R(F_F_minus_F(temp,x))); } # u = -ln(-x+sqrt(1+x^2))
  828.              }
  829.            return;
  830.         } }
  831.      {var reg1 object z = R_R_complex_C(x,y); # z=x+iy
  832.       pushSTACK(z);
  833.       z = N_1_plus_N(N_sqrt_N(N_1_plus_N(N_N_mal_N(z,z)))); # 1+sqrt(1+z^2)
  834.       z = N_N_durch_N(popSTACK(),z); # z/(1+sqrt(1+z^2))
  835.       # Da z=x+iy weder reell noch rein imaginär ist, ist auch
  836.       # w := z/(1+sqrt(1+z^2)) weder reell noch rein imaginär.
  837.       # (Beweis: Sollte sqrt(1+z^2) rationalen Real- und Imaginärteil haben,
  838.       # so auch z, also auch w, und die Formel z = 2w/(1-w^2) zeigt, daß dann
  839.       # z reell oder rein imaginär sein müßte. Also hat sqrt(1+z^2) ein
  840.       # Float als Real- oder Imaginärteil, das Betragsquadrat des Nenners
  841.       # ist also ein Float, und da Real- und Imaginärteil von z /=0 sind,
  842.       # sind Real- und Imaginärteil von w Floats.)
  843.       # Daher hat dann atanh(...) Floats als Realteil u und Imaginärteil v.
  844.       R_R_atanh_R_R(TheComplex(z)->c_real,TheComplex(z)->c_imag); # atanh nehmen
  845.       # u und v mit 2 multiplizieren:
  846.       STACK_1 = F_I_scale_float_F(STACK_1,Fixnum_1); # u:=2*u
  847.       STACK_0 = F_I_scale_float_F(STACK_0,Fixnum_1); # v:=2*v
  848.       return;
  849.     }}
  850.   #
  851.   local object N_asinh_N(z)
  852.     var reg1 object z;
  853.     { if (N_realp(z))
  854.         { R_R_asinh_R_R(z,Fixnum_0); }
  855.         else
  856.         { R_R_asinh_R_R(TheComplex(z)->c_real,TheComplex(z)->c_imag); }
  857.       # Stackaufbau: u, v.
  858.       z = R_R_complex_N(STACK_1,STACK_0); skipSTACK(2); return z;
  859.     }
  860.   #
  861.   local object N_asin_N(z)
  862.     var reg1 object z;
  863.     { # asinh(iz) errechnen:
  864.       if (N_realp(z))
  865.         { R_R_asinh_R_R(Fixnum_0,z); }
  866.         else
  867.         { pushSTACK(TheComplex(z)->c_real);
  868.           z = R_minus_R(TheComplex(z)->c_imag);
  869.           R_R_asinh_R_R(z,popSTACK());
  870.         }
  871.       # Stackaufbau: u, v.
  872.       z = R_minus_R(STACK_1); z = R_R_complex_N(STACK_0,z); # z := v-iu
  873.       skipSTACK(2); return z;
  874.     }
  875.  
  876. # N_acos_N(z) liefert den Arccos einer Zahl z.
  877. # kann GC auslösen
  878.   local object N_acos_N (object z);
  879. # Methode:
  880. # Wert und Branch Cuts nach der Formel CLTL2, S. 312:
  881. #   arccos(z) = log(z+i*sqrt(1-z^2))/i = pi/2 - arcsin(z)
  882. # Sei z=x+iy.
  883. # Falls y=0:
  884. #   Falls x rational:
  885. #     Bei x=1: Ergebnis 0.
  886. #     Bei x=1/2: Ergebnis pi/3.
  887. #     Bei x=0: Ergebnis pi/2.
  888. #     Bei x=-1/2: Ergebnis 2pi/3.
  889. #     Bei x=-1: Ergebnis pi.
  890. #     Sonst x in Float umwandeln.
  891. #   Falls x>1: Ergebnis i ln(x+sqrt(x^2-1)).
  892. # Sonst errechne u+iv = arsinh(-y+ix) wie oben, Ergebnis (pi/2-v)+iu.
  893.   local object N_acos_N(z)
  894.     var reg1 object z;
  895.     { if (N_realp(z)) # y=0 ?
  896.         { if (R_rationalp(z))
  897.             # z rational
  898.             { if (RA_integerp(z))
  899.                 # z Integer
  900.                 { if (eq(z,Fixnum_0)) # x=0 -> Ergebnis pi/2
  901.                     { return F_I_scale_float_F(pi(),Fixnum_minus1); }
  902.                   if (eq(z,Fixnum_1)) # x=1 -> Ergebnis 0
  903.                     { return Fixnum_0; }
  904.                   if (eq(z,Fixnum_minus1)) # x=-1 -> Ergebnis pi
  905.                     { return pi(); }
  906.                   z = I_float_F(z); # z in Float umwandeln
  907.                 }
  908.                 else
  909.                 # z Ratio
  910.                 { if (eq(TheRatio(z)->rt_den,fixnum(2))) # Nenner = 2 ?
  911.                     { var reg2 object temp = TheRatio(z)->rt_num; # Zähler
  912.                       if (eq(temp,Fixnum_1)) # x=1/2 -> Ergebnis pi/3
  913.                         { return R_R_durch_R(pi(),fixnum(3)); }
  914.                       if (eq(temp,Fixnum_minus1)) # x=-1/2 -> Ergebnis 2pi/3
  915.                         { return R_R_durch_R(F_I_scale_float_F(pi(),Fixnum_1),fixnum(3)); }
  916.                     }
  917.                   z = RA_float_F(z); # z in Float umwandeln
  918.             }   }
  919.           # z Float
  920.           pushSTACK(z);
  921.           if (R_R_comp(Fixnum_1,z)<0) # 1<z ?
  922.             { var reg1 object temp = STACK_0; # z
  923.               temp = R_R_minus_R(F_F_mal_F(temp,temp),Fixnum_1); # z^2-1, ein Float >=0
  924.               temp = F_sqrt_F(temp); # sqrt(z^2-1), ein Float >=0
  925.               temp = F_F_plus_F(popSTACK(),temp); # z+sqrt(z^2-1), ein Float >1
  926.               temp = R_ln_R(temp); # ln(z+sqrt(z^2-1)), ein Float >=0
  927.               return R_R_complex_C(Fixnum_0,temp);
  928.             }
  929.           R_R_asinh_R_R(Fixnum_0,popSTACK());
  930.         }
  931.         else
  932.         { pushSTACK(TheComplex(z)->c_real);
  933.           z = R_minus_R(TheComplex(z)->c_imag);
  934.           R_R_asinh_R_R(z,popSTACK());
  935.         }
  936.       # Stackaufbau: u, v.
  937.       # Bilde pi/2-v :
  938.       z = STACK_0;
  939.       z = (R_rationalp(z) ? pi() : pi_F_float_F(z)); # pi im Float-Format von v
  940.       z = F_I_scale_float_F(z,Fixnum_minus1); # pi/2
  941.       z = R_R_minus_R(z,STACK_0); # pi/2-v
  942.       z = R_R_complex_N(z,STACK_1); # (pi/2-v)+iu
  943.       skipSTACK(2); return z;
  944.     }
  945.  
  946. # N_acosh_N(z) liefert den Arcosh einer Zahl z.
  947. # kann GC auslösen
  948.   local object N_acosh_N (object z);
  949. # Methode:
  950. # Wert und Branch Cuts nach der Formel CLTL2, S. 314:
  951. #   arcosh(z) = 2 log(sqrt((z+1)/2)+sqrt((z-1)/2))
  952. # Sei z=x+iy.
  953. # Falls y=0:
  954. #   Falls x rational:
  955. #     Bei x=1: Ergebnis 0.
  956. #     Bei x=1/2: Ergebnis pi/3 i.
  957. #     Bei x=0: Ergebnis pi/2 i.
  958. #     Bei x=-1/2: Ergebnis 2pi/3 i.
  959. #     Bei x=-1: Ergebnis pi i.
  960. #   Falls x<-1:
  961. #     x in Float umwandeln, Ergebnis log(sqrt(x^2-1)-x) + i pi.
  962. # Sonst nach (!) mit u = sqrt((z+1)/2) und v = sqrt((z-1)/2) :
  963. # arcosh(z) = 4 artanh(v/(u+1)) = 4 artanh(sqrt((z-1)/2)/(1+sqrt((z+1)/2)))
  964.   local object N_acosh_N(z)
  965.     var reg1 object z;
  966.     { if (N_realp(z)) # y=0 ?
  967.         { if (R_rationalp(z))
  968.             # z rational
  969.             { if (RA_integerp(z))
  970.                 # z Integer
  971.                 { if (eq(z,Fixnum_0)) # x=0 -> Ergebnis pi/2 i
  972.                     { return R_R_complex_C(Fixnum_0,F_I_scale_float_F(pi(),Fixnum_minus1)); }
  973.                   if (eq(z,Fixnum_1)) # x=1 -> Ergebnis 0
  974.                     { return Fixnum_0; }
  975.                   if (eq(z,Fixnum_minus1)) # x=-1 -> Ergebnis pi i
  976.                     { return R_R_complex_C(Fixnum_0,pi()); }
  977.                 }
  978.                 else
  979.                 # z Ratio
  980.                 { if (eq(TheRatio(z)->rt_den,fixnum(2))) # Nenner = 2 ?
  981.                     { var reg2 object temp = TheRatio(z)->rt_num; # Zähler
  982.                       if (eq(temp,Fixnum_1)) # x=1/2 -> Ergebnis pi/3 i
  983.                         { return R_R_complex_C(Fixnum_0,R_R_durch_R(pi(),fixnum(3))); }
  984.                       if (eq(temp,Fixnum_minus1)) # x=-1/2 -> Ergebnis 2pi/3 i
  985.                         { return R_R_complex_C(Fixnum_0,R_R_durch_R(F_I_scale_float_F(pi(),Fixnum_1),fixnum(3))); }
  986.                     }
  987.             }   }
  988.           pushSTACK(z);
  989.           if (R_R_comp(z,Fixnum_minus1)<0) # z<-1 ?
  990.             { z = STACK_0;
  991.               if (R_rationalp(z)) { STACK_0 = z = RA_float_F(z); }
  992.               # z Float <= -1
  993.               z = F_sqrt_F(R_R_minus_R(F_F_mal_F(z,z),Fixnum_1)); # sqrt(z^2-1), ein Float >=0
  994.               STACK_0 = R_ln_R(F_F_minus_F(z,STACK_0)); # log(sqrt(z^2-1)-z), ein Float >=0
  995.               z = pi(); # und Imaginärteil pi
  996.               return R_R_complex_C(popSTACK(),z);
  997.             }
  998.           z = popSTACK();
  999.         }
  1000.       pushSTACK(z);
  1001.      {var reg2 object temp;
  1002.       temp = N_sqrt_N(N_N_durch_N(N_minus1_plus_N(z),fixnum(2))); # Zähler
  1003.       z = STACK_0; STACK_0 = temp;
  1004.       temp = N_1_plus_N(N_sqrt_N(N_N_durch_N(N_1_plus_N(z),fixnum(2)))); # Nenner
  1005.       return N_N_mal_N(fixnum(4),N_atanh_N(N_N_durch_N(popSTACK(),temp)));
  1006.     }}
  1007.  
  1008.