home *** CD-ROM | disk | FTP | other *** search
/ Il CD di internet / CD.iso / SOURCE / D / CLISP / CLISPSRC.TAR / clisp-1995-01-01 / src / realtran.d < prev    next >
Encoding:
Text File  |  1993-10-24  |  52.4 KB  |  1,140 lines
  1. # Transzendente Funktionen fⁿr reelle Zahlen
  2.  
  3. # pi_F_float_F(f) liefert die Zahl pi im selben Float-Format wie f.
  4. # kann GC ausl÷sen
  5.   local object pi_F_float_F (object f);
  6.   local object pi_F_float_F(f)
  7.     var reg5 object f;
  8.     { floatcase(f,
  9.                 { return O(SF_pi); },
  10.                 { return O(FF_pi); },
  11.                 { return O(DF_pi); },
  12.                 ;
  13.                );
  14.      {var reg6 object pi = O(LF_pi);
  15.       var reg4 uintC f_len = TheLfloat(f)->len; # gewⁿnschte LΣnge von Pi
  16.       {var reg1 uintC len = TheLfloat(pi)->len; # vorhandene LΣnge
  17.        if (f_len < len) { return LF_shorten_LF(pi,f_len); }
  18.        if (f_len == len) { return pi; }
  19.       }
  20.       # gewⁿnschte > vorhandene LΣnge -> mu▀ nachberechnen:
  21.       # Methode:
  22.       # [Richard P. Brent: Fast multiple-precision evaluation of elementary
  23.       #  functions. J. ACM 23(1976), 242-251.]
  24.       # d=f_len, n:=16*d. Verwende Long-Floats mit 16*(d+1) Mantissenbits.
  25.       # (let* ((a (coerce 1 'long-float)) ; 1
  26.       #        (b (sqrt (scale-float a -1))) ; 2^-(1/2)
  27.       #        (eps (scale-float a (- n))) ; 2^-n
  28.       #        (t (scale-float a -2)) ; 1/4
  29.       #        (x 0)
  30.       #       )
  31.       #   (loop
  32.       #     (when (< (- a b) eps) (return))
  33.       #     (let ((y a))
  34.       #       (setq a (scale-float (+ a b) -1))
  35.       #       (setq b (sqrt (* b y)))
  36.       #       (setq t (- t (scale-float (expt (- a y) 2) x)))
  37.       #     )
  38.       #     (incf x)
  39.       #   )
  40.       #   (/ (expt a 2) t)
  41.       # )
  42.       {var reg3 uintC len = f_len + 1; # Arbeite mit Long-Floats mit len Digits
  43.        if (uintCoverflow(len)) { fehler_LF_toolong(); }
  44.        {var reg2 uintL uexp_limit = LF_exp_mid - intDsize*(uintL)f_len; # LF_exp_mid - n
  45.         # Ein Long-Float ist genau dann betragsmΣ▀ig <2^-n, wenn
  46.         # sein Exponent < LF_exp_mid-n = uexp_limit ist.
  47.         {var reg1 object temp = I_to_LF(Fixnum_1,len); # 1 als Long-Float
  48.          pushSTACK(temp); # =: a
  49.          temp = LF_I_scale_float_LF(temp,Fixnum_minus1); # (scale-float a -1)
  50.          pushSTACK(LF_sqrt_LF(temp)); # daraus die Wurzel, =: b
  51.          pushSTACK(Fixnum_0); # x:=0
  52.          temp = LF_I_scale_float_LF(STACK_2,sfixnum(-2)); # (scale-float a -2)
  53.          pushSTACK(temp); # =: t
  54.         }# Stackaufbau: a, b, x, t.
  55.         loop
  56.           {{var reg1 object temp;
  57.             temp = LF_LF_minus_LF(STACK_3,STACK_2); # (- a b)
  58.             if (TheLfloat(temp)->expo < uexp_limit) # Exponent < uexp_limit
  59.               break; } # ja -> |a-b| < 2^-n -> fertig
  60.            {var reg1 object temp;
  61.             temp = LF_LF_plus_LF(STACK_3,STACK_2); # a+b
  62.             temp = LF_I_scale_float_LF(temp,Fixnum_minus1); # (a+b)/2
  63.             pushSTACK(temp); } # neues a
  64.             STACK_(2+1) = LF_sqrt_LF(LF_LF_mal_LF(STACK_(3+1),STACK_(2+1))); # b := sqrt(a*b)
  65.            {var reg1 object temp;
  66.             temp = STACK_(3+1); # altes a
  67.             temp = LF_LF_minus_LF(STACK_(3+1) = STACK_0, temp); # neues a - altes a
  68.             temp = LF_LF_mal_LF(temp,temp); # quadieren
  69.             temp = LF_I_scale_float_LF(temp,STACK_(1+1)); # mal 2^x
  70.             skipSTACK(1);
  71.             STACK_0 = LF_LF_minus_LF(STACK_0,temp); # von t subtrahieren
  72.             STACK_1 = fixnum_inc(STACK_1,1); # x:=x+1
  73.           }}
  74.         {var reg1 object temp;
  75.          temp = STACK_3; temp = LF_LF_mal_LF(temp,temp); # a quadieren
  76.          temp = LF_LF_durch_LF(temp,STACK_0); # durch t dividieren
  77.          skipSTACK(4);
  78.          # temp = Pi ist fertig.
  79.          return O(LF_pi) = LF_shorten_LF(temp,f_len); # wieder verkⁿrzen, als LF_pi abspeichern
  80.     }}}}}
  81.  
  82. # pi() liefert die Zahl pi im Default-Float-Format.
  83. # kann GC ausl÷sen
  84.   local object pi (void);
  85.   local object pi()
  86.     { defaultfloatcase(S(default_float_format),
  87.                        return O(SF_pi); , # pi als SF
  88.                        return O(FF_pi); , # pi als FF
  89.                        return O(DF_pi); , # pi als DF
  90.                        return O(pi); , # pi als LF der DefaultlΣnge
  91.                        ,); # nichts zu retten
  92.     }
  93.  
  94. # F_atanhx_F(x) liefert zu einem Float x (betragsmΣ▀ig <1/2) atanh(x) als Float.
  95. # kann GC ausl÷sen
  96.   local object F_atanhx_F (object x);
  97. # Methode:
  98. # e := Exponent aus (decode-float x), d := (float-digits x)
  99. # Bei x=0.0 oder e<=-d/2 liefere x
  100. #   (denn bei e<=-d/2 ist x^2 < 2^(-d), also
  101. #   1 <= atanh(x)/x = 1+x^2/3+x^4/5+... < 1+x^2/2 < 1+2^(-d-1) < 1+2^(-d),
  102. #   also ist atanh(x)/x, auf d Bits gerundet, gleich 1.0).
  103. # Bei e<=-sqrt(d) verwende die Potenzreihe
  104. #   atanh(x)/x = sum(j=0..inf,(x^2)^j/(2j+1)):
  105. #   a:=x^2, b:=1, i:=1, sum:=0,
  106. #   while (/= sum (setq sum (+ sum (/ b i)))) do i:=i+2, b:=b*a.
  107. #   Ergebnis x*sum.
  108. # Sonst setze y := x/(1+sqrt(1-x^2)), berechne rekursiv z:=atanh(y)
  109. #   und liefere 2*z = (scale-float z 1).
  110. # Diese Rekursion wird entrekursiviert. Statt k mal hintereinander
  111. #   x := x/(1+sqrt(1-x^2)) zu bilden, arbeitet man lieber mit den Kehrwerten,
  112. #   setzt also x := 1/|x|, dann k mal x := x+sqrt(x^2-1), dann x := +- 1/x.
  113. # Aufwand: asymptotisch d^2.5 .
  114.  
  115. # F_atanx_F(x) liefert zu einem Float x (betragsmΣ▀ig <=1) atan(x) als Float.
  116. # kann GC ausl÷sen
  117.   local object F_atanx_F (object x);
  118. # Methode:
  119. # e := Exponent aus (decode-float x), d := (float-digits x)
  120. # Bei x=0.0 oder e<=-d/2 liefere x
  121. #   (denn bei e<=-d/2 ist x^2/3 < x^2/2 < 2^(-d)/2 = 2^(-d-1), also
  122. #   1 >= atan(x)/x > 1-x^2/3 > 1-2^(-d-1),
  123. #   also ist atan(x)/x, auf d Bits gerundet, gleich 1.0).
  124. # Bei e<=-sqrt(d) verwende die Potenzreihe
  125. #   atan(x)/x = sum(j=0..inf,(-x^2)^j/(2j+1)):
  126. #   a:=-x^2, b:=1, i:=1, sum:=0,
  127. #   while (/= sum (setq sum (+ sum (/ b i)))) do i:=i+2, b:=b*a.
  128. #   Ergebnis x*sum.
  129. # Sonst setze y := x/(1+sqrt(1+x^2)), berechne rekursiv z:=atan(y)
  130. #   und liefere 2*z = (scale-float z 1).
  131. # Diese Rekursion wird entrekursiviert. Statt k mal hintereinander
  132. #   x := x/(1+sqrt(1+x^2)) zu bilden, arbeitet man lieber mit den Kehrwerten,
  133. #   setzt also x := 1/|x|, dann k mal x := x+sqrt(x^2+1), dann x := +- 1/x.
  134. # Aufwand: asymptotisch d^2.5 .
  135.  
  136. # Generiert eine Funktion wie F_atanx_F
  137.   #define GEN_F_atanx(name,Fixnum_plusminus1,F_plusminus_F)                            \
  138.     local object CONCAT3(F_,name,_F) PARM1(x,                                          \
  139.       var reg3 object x)                                                               \
  140.       { if (R_zerop(x)) { return x; }                                                  \
  141.        {var reg8 uintL d = F_float_digits(x);                                          \
  142.         var reg7 sintL e = F_exponent_L(x);                                            \
  143.         if (e <= (sintL)(-d)>>1) # e <= -d/2 <==> e <= -ceiling(d/2)                   \
  144.           { return x; } # ja -> x als Ergebnis                                         \
  145.         pushSTACK(x);                                                                  \
  146.         # Stackaufbau: x.                                                              \
  147.         {var reg4 object k = Fixnum_0; # RekursionszΣhler k:=0                         \
  148.          var reg6 uintL sqrt_d = UL_sqrt_UW(d); # floor(sqrt(d))                       \
  149.          # Bei e <= -1-floor(sqrt(d)) kann die Potenzreihe angewandt werden.           \
  150.          if (e >= (sintL)(-sqrt_d))                                                             \
  151.            # e > -1-floor(sqrt(d)) -> mu▀ |x| verkleinern.                             \
  152.            { var reg5 sintL e_limit = 1+sqrt_d; # 1+floor(sqrt(d))                     \
  153.              pushSTACK(x = F_durch_F(F_abs_F(x))); # 1/|x|                             \
  154.              # Stackaufbau: originales x, neues x.                                     \
  155.              loop                                                                      \
  156.                { # nΣchstes x nach der Formel x := x+sqrt(x^2 +- 1) berechnen:         \
  157.                  x = F_sqrt_F(R_R_plus_R(F_F_mal_F(x,x),Fixnum_plusminus1));           \
  158.                  STACK_0 = x = F_F_plus_F(STACK_0,x);                                  \
  159.                  k = fixnum_inc(k,1); # k:=k+1                                         \
  160.                  if (F_exponent_L(x) > e_limit) break;                                 \
  161.                }                                                                       \
  162.              # Schleifenende mit Exponent(x) > 1+floor(sqrt(d)), also                  \
  163.              # x >= 2^(1+floor(sqrt(d))), also 1/x <= 2^(-1-floor(sqrt(d))).           \
  164.              # Nun kann die Potenzreihe auf 1/x angewandt werden.                      \
  165.             {var reg1 object x = F_durch_F(popSTACK());                                \
  166.              if (R_mminusp(STACK_0)) { x = F_minus_F(x); } # Vorzeichen wieder rein    \
  167.              STACK_0 = x; # neues x ersetzt altes x                                    \
  168.            }}                                                                          \
  169.          # Stackaufbau: neues x.                                                       \
  170.          # Potenzreihe anwenden:                                                       \
  171.          {var reg2 object i = Fixnum_1;                                                \
  172.           pushSTACK(F_plusminus_F(F_F_mal_F(STACK_0,STACK_0))); # a := -x^2 bzw. x^2   \
  173.           pushSTACK(I_F_float_F(Fixnum_1,STACK_1)); # b := (float 1 x)                 \
  174.           pushSTACK(I_F_float_F(Fixnum_0,STACK_2)); # sum := (float 0 x)               \
  175.           # Stackaufbau: x, a, b, sum.                                                 \
  176.           loop                                                                         \
  177.             { var reg1 object temp;                                                    \
  178.               temp = R_R_durch_R(STACK_1,i); # (/ b i)                                 \
  179.               temp = F_F_plus_F(STACK_0,temp); # (+ sum (/ b i))                       \
  180.               if (eql(STACK_0,temp)) # = sum ?                                         \
  181.                 break; # ja -> Potenzreihe abbrechen                                   \
  182.               STACK_0 = temp;                                                          \
  183.               STACK_1 = F_F_mal_F(STACK_1,STACK_2); # b := b*a                         \
  184.               i = fixnum_inc(i,2); # i := i+2                                          \
  185.          }  }                                                                          \
  186.          {var reg1 object erg = F_F_mal_F(STACK_0,STACK_3); # sum*x als Ergebnis       \
  187.           skipSTACK(4);                                                                \
  188.           return F_I_scale_float_F(erg,k); # wegen Rekursion noch mal 2^k              \
  189.       }}}}
  190. # F_atanx_F : mit x -> x+sqrt(x^2-1), a = -x^2
  191.   GEN_F_atanx(atanx,Fixnum_1,F_minus_F)
  192. # F_atanhx_F : mit x -> x+sqrt(x^2+1), a = x^2
  193.   GEN_F_atanx(atanhx,Fixnum_minus1,)
  194.  
  195. # R_R_atan_R(x,y) liefert zu zwei reellen Zahlen x, y den Winkel von (x,y)
  196. # in Polarkoordinaten. Ergebnis rational nur, wenn x>0 und y=0.
  197. # kann GC ausl÷sen
  198.   local object R_R_atan_R (object x, object y);
  199. # Methode:
  200. # y=0 -> bei x>0: 0 als Ergebnis,
  201. #        bei x<0: pi als Ergebnis.
  202. #        bei x=0: Error.
  203. # x=0 -> bei y>0: pi/2 als Ergebnis.
  204. #        bei y<0: -pi/2 als Ergebnis.
  205. #        bei y=0: Error.
  206. # Falls x und y beide rational: beide in Floats umwandeln.
  207. # 0 <= |y| <= x  ->  atan(y/x)
  208. # 0 <= |x| <= y  ->  pi/2 - atan(x/y)
  209. # 0 <= |x| <= -y  ->  -pi/2 - atan(x/y)
  210. # 0 <= |y| <= -x  ->  fⁿr y>=0: pi + atan(y/x), fⁿr y<0: -pi + atan(y/x)
  211.   local object R_R_atan_R(x,y)
  212.     var reg2 object x;
  213.     var reg3 object y;
  214.     { if (eq(y,Fixnum_0))
  215.         # y=0 (exakt)
  216.         { if (R_zerop(x)) { divide_0(); } # x=0 -> Error
  217.           if (R_minusp(x)) { return pi(); } # x<0 -> pi in Default-Float-Genauigkeit
  218.           else { return Fixnum_0; } # x>0 -> 0
  219.         }
  220.       elif (eq(x,Fixnum_0))
  221.         # x=0 (exakt)
  222.         { if (R_zerop(y)) { divide_0(); } # y=0 -> Error
  223.           if (R_minusp(x)) # x<0 -> -pi/2
  224.             { return F_minus_F(F_I_scale_float_F(pi(),Fixnum_minus1)); }
  225.           else # x>0 -> pi/2
  226.             { return F_I_scale_float_F(pi(),Fixnum_minus1); }
  227.         }
  228.       pushSTACK(x); pushSTACK(y);
  229.       # Stackaufbau: x, y.
  230.       if (R_rationalp(x) && R_rationalp(y))
  231.         # x,y in Floats umwandeln:
  232.         { STACK_1 = RA_float_F(x); STACK_0 = RA_float_F(STACK_0); }
  233.       # x,y nicht exakt =0, x/y und y/x werden Floats sein.
  234.       pushSTACK(R_abs_R(STACK_1)); y = R_abs_R(STACK_(0+1)); x = popSTACK();
  235.       if (R_R_comp(x,y) >= 0) # (abs x) und (abs y) vergleichen
  236.         # |x| >= |y|
  237.         { var reg1 object z = F_atanx_F(R_R_durch_R(STACK_0,STACK_1)); # atan(y/x)
  238.           # Division war erfolgreich, also x/=0.
  239.           if (R_mminusp(STACK_1))
  240.             # x<0 -> pi bzw. -pi addieren:
  241.             { STACK_1 = z; # atan(y/x) retten
  242.               z = pi_F_float_F(z); # pi im selben Float-Format
  243.               if (!R_mminusp(STACK_0))
  244.                 { z = F_F_plus_F(STACK_1,z); } # y>=0 -> atan(y/x) + pi
  245.                 else
  246.                 { z = F_F_minus_F(STACK_1,z); } # y<0 -> atan(y/x) - pi
  247.             }
  248.           skipSTACK(2);
  249.           return z;
  250.         }
  251.         else
  252.         # |x| < |y|
  253.         { var reg1 object z = F_atanx_F(R_R_durch_R(STACK_1,STACK_0)); # atan(x/y)
  254.           # von pi/2 bzw. -pi/2 subtrahieren:
  255.           STACK_1 = z; # atan(x/y) retten
  256.           z = pi_F_float_F(z); # pi im selben Float-Format
  257.           z = F_I_scale_float_F(z,Fixnum_minus1); # pi/2
  258.           if (R_mminusp(STACK_0)) { z = F_minus_F(z); } # y<0 -> -pi/2 statt pi/2
  259.           z = F_F_minus_F(z,STACK_1); # +-pi/2 - atan(x/y)
  260.           skipSTACK(2);
  261.           return z;
  262.         }
  263.     }
  264.  
  265. # R_atan_R(x) liefert den Arctan einer reellen Zahl x.
  266. # Ergebnis rational nur, wenn x=0.
  267. # kann GC ausl÷sen
  268.   local object R_atan_R (object x);
  269. # Methode:
  270. # arctan(x) = arctan(X=1,Y=x).
  271. #if 0
  272.   local object R_atan_R(x)
  273.     var reg1 object x;
  274.     { return R_R_atan_R(Fixnum_1,x); }
  275. #else # Macro spart Code
  276.   #define R_atan_R(x)  R_R_atan_R(Fixnum_1,x)
  277. #endif
  278.  
  279. # F_sinx_F(x) liefert zu einem Float x (betragsmΣ▀ig <2) (sin(x)/x)^2 als Float.
  280. # kann GC ausl÷sen
  281.   local object F_sinx_F (object x);
  282. # Methode:
  283. # e := Exponent aus (decode-float x), d := (float-digits x)
  284. # Bei x=0.0 oder e<=-d/2 liefere 1.0
  285. #   (denn bei e<=-d/2 ist x^2/6 < x^2/4 < 2^(-d)/4 = 2^(-d-2), also
  286. #   1 >= sin(x)/x > 1-x^2/6 > 1-2^(-d-2), also 1 >= (sin(x)/x)^2 > 1-2^(-d-1),
  287. #   also ist (sin(x)/x)^2, auf d Bits gerundet, gleich 1.0).
  288. # Bei e<=-sqrt(d) verwende die Potenzreihe
  289. #   sin(x)/x = sum(j=0..inf,(-x^2)^j/(2j+1)!):
  290. #   a:=-x^2, b:=1, i:=1, sum:=0,
  291. #   while (/= sum (setq sum (+ sum b))) do b:=b*a/((i+1)*(i+2)), i:=i+2.
  292. #   Ergebnis sum^2.
  293. # Sonst setze y := x/2 = (scale-float x -1),
  294. #   berechne rekursiv z:=(sin(y)/y)^2 und liefere z*(1-y^2*z).
  295. # [Die Grenze sqrt(d) ergibt sich so:
  296. #  Man braucht bei der Potenzreihe mit x=2^-k etwa j Glieder, mit
  297. #  k*j*ln 2 + j*(ln j - 1) = d, und der Aufwand betrΣgt etwa 2.8*(j/2)
  298. #  Multiplikationen von d-Bit-Zahlen. Bei Halbierungen bis x=2^-k ist der
  299. #  Gesamtaufwand etwa 2*(k+e)+1.4*j(k). Dieses minimieren nach k: Soll sein
  300. #  -1.4 = d/dk j(k) = (d/dj k(j))^-1 = - j^2/(d+j)*ln 2, also j^2=2(d+j),
  301. #  grob j=sqrt(2d) und damit k=sqrt(d).]
  302. # Aufwand: asymptotisch d^2.5 .
  303.  
  304. # F_sinhx_F(x) liefert zu einem Float x (betragsmΣ▀ig <2) (sinh(x)/x)^2 als Float.
  305. # kann GC ausl÷sen
  306.   local object F_sinhx_F (object x);
  307. # Methode:
  308. # e := Exponent aus (decode-float x), d := (float-digits x)
  309. # Bei x=0.0 oder e<=(1-d)/2 liefere 1.0
  310. #   (denn bei e<=(1-d)/2 ist x^2/6 < x^2/4 < 2^(1-d)/4 = 2^(-d-1), also
  311. #   1 <= sinh(x)/x = 1+x^2/6+... < 1+2^(-d-1), also 1 <= (sinh(x)/x)^2 < 1+2^(-d),
  312. #   also ist (sinh(x)/x)^2, auf d Bits gerundet, gleich 1.0).
  313. # Bei e<=-sqrt(d) verwende die Potenzreihe
  314. #   sinh(x)/x = sum(j=0..inf,(x^2)^j/(2j+1)!):
  315. #   a:=x^2, b:=1, i:=1, sum:=0,
  316. #   while (/= sum (setq sum (+ sum b))) do b:=b*a/((i+1)*(i+2)), i:=i+2.
  317. #   Ergebnis sum^2.
  318. # Sonst setze y := x/2 = (scale-float x -1),
  319. #   berechne rekursiv z:=(sinh(y)/y)^2 und liefere z*(1+y^2*z).
  320. # [Die Grenze sqrt(d) ergibt sich so:
  321. #  Man braucht bei der Potenzreihe mit x=2^-k etwa j Glieder, mit
  322. #  k*j*ln 2 + j*(ln j - 1) = d, und der Aufwand betrΣgt etwa 2.8*(j/2)
  323. #  Multiplikationen von d-Bit-Zahlen. Bei Halbierungen bis x=2^-k ist der
  324. #  Gesamtaufwand etwa 2*(k+e)+1.4*j(k). Dieses minimieren nach k: Soll sein
  325. #  -1.4 = d/dk j(k) = (d/dj k(j))^-1 = - j^2/(d+j)*ln 2, also j^2=2(d+j),
  326. #  grob j=sqrt(2d) und damit k=sqrt(d).]
  327. # Aufwand: asymptotisch d^2.5 .
  328.  
  329. # Generiert eine Funktion wie F_sinx_F
  330.   #define GEN_F_sinx(name,f,flag,R_R_plusminus_R)                                               \
  331.     local object CONCAT3(F_,name,_F) PARM1(x,                                                   \
  332.       var reg4 object x)                                                                        \
  333.       { if (R_zerop(x)) { return I_F_float_F(Fixnum_1,x); }                                     \
  334.        {var reg6 uintL d = F_float_digits(x);                                                   \
  335.         var reg5 sintL e = F_exponent_L(x);                                                     \
  336.         if (e <= (sintL)(f-d)>>1) # e <= (f-d)/2 <==> e <= -ceiling((d-f)/2) ?                  \
  337.           { return I_F_float_F(Fixnum_1,x); } # ja -> 1.0 als Ergebnis                          \
  338.         pushSTACK(x); pushSTACK(NIL); pushSTACK(F_F_mal_F(x,x));                                \
  339.         # Stackaufbau: x, y^2-Liste, x^2.                                                       \
  340.         {# Bei e <= -1-floor(sqrt(d)) kann die Potenzreihe angewandt werden.                    \
  341.          var reg3 sintL e_limit = -1-UL_sqrt_UW(d); # -1-floor(sqrt(d))                         \
  342.          while (e > e_limit)                                                                    \
  343.            # e > -1-floor(sqrt(d)) -> mu▀ |x| verkleinern.                                      \
  344.            { STACK_2 = F_I_scale_float_F(STACK_2,Fixnum_minus1); # x := x/2                     \
  345.              STACK_0 = F_I_scale_float_F(STACK_0,sfixnum(-2)); # x^2 := x^2/4                   \
  346.              e = e-1; # e := e-1                                                                \
  347.              # (push x^2 y^2-Liste) :                                                           \
  348.             {var reg1 object new_cons = allocate_cons();                                        \
  349.              Car(new_cons) = STACK_0; Cdr(new_cons) = STACK_1;                                  \
  350.              STACK_1 = new_cons;                                                                \
  351.         }  }}                                                                                   \
  352.         # Potenzreihe anwenden:                                                                 \
  353.         if (flag) { STACK_0 = F_minus_F(STACK_0); } # a := -x^2 bzw. x^2                        \
  354.         {var reg2 object i = Fixnum_1;                                                          \
  355.          pushSTACK(I_F_float_F(Fixnum_1,STACK_2)); # b := (float 1 x)                           \
  356.          pushSTACK(I_F_float_F(Fixnum_0,STACK_3)); # sum := (float 0 x)                         \
  357.          # Stackaufbau: x, y^2-Liste, a, b, sum.                                                \
  358.          loop                                                                                   \
  359.            { var reg1 object temp;                                                              \
  360.              temp = F_F_plus_F(STACK_0,STACK_1); # (+ sum b)                                    \
  361.              if (eql(STACK_0,temp)) # = sum ?                                                   \
  362.                break; # ja -> Potenzreihe abbrechen                                             \
  363.              STACK_0 = temp;                                                                    \
  364.              STACK_1 = F_F_mal_F(STACK_1,STACK_2); # b := b*a                                   \
  365.              temp = I_I_mal_I(fixnum_inc(i,1),fixnum_inc(i,2)); # (i+1)*(i+2)                   \
  366.              i = fixnum_inc(i,2); # i := i+2                                                    \
  367.              STACK_1 = R_R_durch_R(STACK_1,temp); # b := b/((i+1)*(i+2))                        \
  368.         }  }                                                                                    \
  369.         {var reg1 object z = F_F_mal_F(STACK_0,STACK_0); # sum^2 als Ergebnis                   \
  370.          # Stackaufbau: x, y^2-Liste, -, -, -.                                                  \
  371.          # Wegen Rekursion noch die y^2-Liste abarbeiten:                                       \
  372.          while (mconsp(STACK_3))                                                                \
  373.            { var reg2 object temp;                                                              \
  374.              temp = STACK_3; STACK_3 = Cdr(temp); temp = Car(temp); # nΣchstes y^2              \
  375.              STACK_4 = z; # z retten                                                            \
  376.              z = R_R_plusminus_R(Fixnum_1,F_F_mal_F(temp,z)); # 1 +- y^2*z                      \
  377.              z = F_F_mal_F(STACK_4,z); # mit z multiplizieren                                   \
  378.            }                                                                                    \
  379.          skipSTACK(5);                                                                          \
  380.          return z;                                                                              \
  381.       }}}
  382. # F_sinx_F : mit z -> z*(1-y^2*z), a = -x^2, -d/2
  383.   GEN_F_sinx(sinx,0,TRUE,R_R_minus_R)
  384. # F_sinhx_F : mit z -> z*(1+y^2*z), a = x^2, (1-d)/2
  385.   GEN_F_sinx(sinhx,1,FALSE,R_R_plus_R)
  386.  
  387. # F_pi_round_I_F(x) dividiert ein Float x mit Rest durch pi.
  388. # Beide Werte von (round x (float pi x)) auf den Stack.
  389. # kann GC ausl÷sen
  390.   local void F_pi_round_I_F (object x);
  391.   local void F_pi_round_I_F(x)
  392.     var reg1 object x;
  393.     { if (F_exponent_L(x) <= 0)
  394.         # Exponent <=0 -> |x|<1 -> |x/pi| < 1/2, also Division unn÷tig
  395.         { pushSTACK(Fixnum_0); pushSTACK(x); } # Quotient 0, Rest x
  396.         else
  397.         { pushSTACK(x); # x retten
  398.          {var reg2 object pi = pi_F_float_F(x); # pi mit hinreichender Genauigkeit
  399.           R_R_round_I_R(popSTACK(),pi); # x durch pi dividieren
  400.         }}
  401.     }
  402.  
  403. # F_pi2_round_I_F(x) dividiert ein Float x mit Rest durch pi/2.
  404. # Beide Werte von (round x (float pi/2 x)) auf den Stack.
  405. # kann GC ausl÷sen
  406.   local void F_pi2_round_I_F (object x);
  407.   local void F_pi2_round_I_F(x)
  408.     var reg1 object x;
  409.     { if (F_exponent_L(x) < 0)
  410.         # Exponent <0 -> |x|<1/2 -> |x/(pi/2)| < 1/2, also Division unn÷tig
  411.         { pushSTACK(Fixnum_0); pushSTACK(x); } # Quotient 0, Rest x
  412.         else
  413.         { pushSTACK(x); # x retten
  414.          {var reg2 object pi = pi_F_float_F(x); # pi mit hinreichender Genauigkeit
  415.           pi = F_I_scale_float_F(pi,Fixnum_minus1); # pi/2 mit hinreichender Genauigkeit
  416.           R_R_round_I_R(popSTACK(),pi); # x durch pi/2 dividieren
  417.         }}
  418.     }
  419.  
  420. # R_sin_R(x) liefert den Sinus (sin x) einer reellen Zahl x.
  421. # kann GC ausl÷sen
  422.   local object R_sin_R (object x);
  423. # Methode:
  424. # x rational -> bei x=0 0 als Ergebnis, sonst x in Float umwandeln.
  425. # x Float -> Genauigkeit erh÷hen,
  426. #   (q,r) := (round x (float pi x)), so da▀ |r|<=pi/2.
  427. #   (sin(r)/r)^2 errechnen, Wurzel ziehen, mit r multiplizieren
  428. #   und - falls q ungerade - Vorzeichenwechsel.
  429.   local object R_sin_R(x)
  430.     var reg1 object x;
  431.     { if (R_rationalp(x))
  432.         { if (eq(x,Fixnum_0)) { return x; } # x=0 -> 0 als Ergebnis
  433.           x = RA_float_F(x); # sonst in Float umwandeln
  434.         }
  435.       # x Float
  436.       pushSTACK(x); # x retten
  437.       x = F_extend_F(x); # Rechengenauigkeit erh÷hen
  438.       F_pi_round_I_F(x); # durch pi dividieren
  439.       # Stackaufbau: Argument, q, r.
  440.      {var reg2 object x;
  441.       x = F_sqrt_F(F_sinx_F(STACK_0)); # Wurzel aus (sin(r)/r)^2
  442.       x = F_F_mal_F(x,STACK_0); # mit r multiplizieren
  443.       x = F_F_float_F(x,STACK_2); # und wieder runden
  444.       if (I_oddp(STACK_1)) { x = F_minus_F(x); } # q ungerade -> mal -1
  445.       skipSTACK(3);
  446.       return x;
  447.     }}
  448.  
  449. # R_cos_R(x) liefert den Cosinus (cos x) einer reellen Zahl x.
  450. # kann GC ausl÷sen
  451.   local object R_cos_R (object x);
  452. # Methode:
  453. # x rational -> bei x=0 1 als Ergebnis, sonst x in Float umwandeln.
  454. # x Float -> Genauigkeit erh÷hen,
  455. #   (q,r) := (round x (float pi x)), so da▀ |r|<=pi/2.
  456. #   e := Exponent aus (decode-float r), d := (float-digits r)
  457. #   Bei r=0.0 oder e<=-d/2 liefere 1.0
  458. #     (denn bei e<=-d/2 ist r^2/2 < 2^(-d)/2 = 2^(-d-1), also
  459. #     1 >= cos(r) > 1-r^2/2 > 1-2^(-d-1),
  460. #     also ist cos(r), auf d Bits gerundet, gleich 1.0).
  461. #   Sonst s := r/2 = (scale-float r -1),
  462. #     (sin(s)/s)^2 errechnen, cos(r) = 1-r*s*(sin(s)/s)^2 errechnen.
  463. #   Falls q ungerade: Vorzeichenwechsel.
  464.   local object R_cos_R(x)
  465.     var reg1 object x;
  466.     { if (R_rationalp(x))
  467.         { if (eq(x,Fixnum_0)) { return Fixnum_1; } # x=0 -> 1 als Ergebnis
  468.           x = RA_float_F(x); # sonst in Float umwandeln
  469.         }
  470.       # x Float
  471.       pushSTACK(x); # x retten
  472.       x = F_extend_F(x); # Rechengenauigkeit erh÷hen
  473.       F_pi_round_I_F(x); # durch pi dividieren
  474.       # Stackaufbau: Argument, q, r.
  475.      {var reg2 object x;
  476.       x = STACK_0;
  477.       if (R_zerop(x) # r=0.0 -> 1.0
  478.           || (F_exponent_L(x) <= (sintL)(-F_float_digits(x))>>1) # e <= -d/2 <==> e <= -ceiling(d/2) ?
  479.          )
  480.         { x = I_F_float_F(Fixnum_1,STACK_2); } # (cos r) = 1.0
  481.         else
  482.         { var reg1 object s = F_I_scale_float_F(STACK_0,Fixnum_minus1); # s := r/2
  483.           pushSTACK(s);
  484.           s = F_sinx_F(s); # (sin(s)/s)^2
  485.           s = F_F_mal_F(popSTACK(),s); # mit s multiplizieren
  486.           s = F_F_mal_F(STACK_0,s); # mit r multiplizieren
  487.           s = R_R_minus_R(Fixnum_1,s); # von 1 subtrahieren
  488.           x = F_F_float_F(s,STACK_2); # und wieder runden
  489.         }
  490.       if (I_oddp(STACK_1)) { x = F_minus_F(x); } # q ungerade -> mal -1
  491.       skipSTACK(3);
  492.       return x;
  493.     }}
  494.  
  495. # R_cos_sin_R_R(x) liefert ((cos x),(sin x)), beide Werte auf dem Stack.
  496. # kann GC ausl÷sen
  497.   local void R_cos_sin_R_R (object x);
  498. # Methode:
  499. # x rational -> bei x=0 (1,0) als Ergebnis, sonst x in Float umwandeln.
  500. # x Float -> Genauigkeit erh÷hen,
  501. #   (q,r) := (round x (float pi/2 x)), so da▀ |r|<=pi/4.
  502. #   y:=(sin(r)/r)^2 errechnen.
  503. #   (cos r) berechnen:
  504. #     e := Exponent aus (decode-float r), d := (float-digits r)
  505. #     Bei r=0.0 oder e<=-d/2 liefere 1.0
  506. #       (denn bei e<=-d/2 ist r^2/2 < 2^(-d)/2 = 2^(-d-1), also
  507. #       1 >= cos(r) > 1-r^2/2 > 1-2^(-d-1),
  508. #       also ist cos(r), auf d Bits gerundet, gleich 1.0).
  509. #     Sonst sqrt(1-r^2*y).
  510. #   (sin r) berechnen: r*sqrt(y).
  511. #   Genauigkeit wieder verringern.
  512. #   Falls q = 0 mod 4: ((cos r), (sin r))
  513. #   Falls q = 1 mod 4: ((- (sin r)), (cos r))
  514. #   Falls q = 2 mod 4: ((- (cos r)), (- (sin r)))
  515. #   Falls q = 3 mod 4: ((sin r), (- (cos r)))
  516.   local void R_cos_sin_R_R(x)
  517.     var reg3 object x;
  518.     { if (R_rationalp(x))
  519.         { if (eq(x,Fixnum_0)) # x=0 -> (1,0) als Ergebnis
  520.             { pushSTACK(Fixnum_1); pushSTACK(Fixnum_0); return; }
  521.           x = RA_float_F(x); # sonst in Float umwandeln
  522.         }
  523.       # x Float
  524.       pushSTACK(x); # x retten
  525.       x = F_extend_F(x); # Rechengenauigkeit erh÷hen
  526.       F_pi2_round_I_F(x); # durch pi/2 dividieren
  527.       # Stackaufbau: Argument, q, r.
  528.       pushSTACK(F_sinx_F(STACK_0)); # y := (sin(r)/r)^2
  529.       # Stackaufbau: Argument, q, r, y.
  530.       # erste Komponente cos(r) berechnen:
  531.       {var reg2 object x;
  532.        x = STACK_0;
  533.        if (R_zerop(x) # r=0.0 -> 1.0
  534.            || (F_exponent_L(x) <= (sintL)(-F_float_digits(x))>>1) # e <= -d/2 <==> e <= -ceiling(d/2) ?
  535.           )
  536.          { x = I_F_float_F(Fixnum_1,STACK_3); } # (cos r) = 1.0
  537.          else
  538.          { var reg1 object r = STACK_1;
  539.            r = F_F_mal_F(r,r); # r^2
  540.            r = F_F_mal_F(r,STACK_0); # r^2*y
  541.            r = R_R_minus_R(Fixnum_1,r); # 1-r^2*y
  542.            r = F_sqrt_F(r); # sqrt(1-r^2*y)
  543.            x = F_F_float_F(r,STACK_3); # und wieder runden
  544.          }
  545.        pushSTACK(x); # cos(r) retten
  546.       }
  547.       # zweite Komponente sin(r) berechnen:
  548.       {var reg1 object x = F_sqrt_F(STACK_(0+1)); # Wurzel aus y
  549.        x = F_F_mal_F(STACK_(1+1),x); # mit r multiplizieren
  550.        x = F_F_float_F(x,STACK_(3+1)); # und wieder runden
  551.       # evtl. Vorzeichenwechsel oder Vertauschen:
  552.        {var reg2 object q = STACK_(2+1);
  553.         switch (I_fixnump(q) # q mod 4, je nachdem ob q Fixnum oder Bignum
  554.                 ? (as_oint(q) >> oint_addr_shift) & (bit(2)-1)
  555.                 : TheBignum(q)->data[(uintP)(TheBignum(q)->length)-1] & (bit(2)-1)
  556.                )
  557.           { case 0:
  558.               STACK_(2+1) = x; STACK_(3+1) = STACK_0; break;
  559.             case 1:
  560.               STACK_(3+1) = F_minus_F(x); STACK_(2+1) = STACK_0; break;
  561.             case 2:
  562.               STACK_(2+1) = F_minus_F(x); STACK_(3+1) = F_minus_F(STACK_0); break;
  563.             case 3:
  564.               STACK_(3+1) = x; STACK_(2+1) = F_minus_F(STACK_0); break;
  565.       }}  }
  566.       skipSTACK(2+1);
  567.     }
  568.  
  569. # F_lnx_F(x) liefert zu einem Float x (>=1/2, <=2) ln(x) als Float.
  570. # kann GC ausl÷sen
  571.   local object F_lnx_F (object x);
  572. # Methode:
  573. # y:=x-1, e := Exponent aus (decode-float y), d := (float-digits y)
  574. # Bei y=0.0 oder e<=-d liefere y
  575. #   (denn bei e<=-d ist y/2 < 2^(-d)/2 = 2^(-d-1), also
  576. #   0 <= y - ln(x) < y^2/2 < 2^(-d-1)*y
  577. #   also ist ln(x)/y, auf d Bits gerundet, gleich y).
  578. # Bei e<=-sqrt(d) verwende die Potenzreihe
  579. #   ln(x) = sum(j=0..inf,(-1)^j*y^(j+1)/(j+1)):
  580. #   a:=-y, b:=y, i:=1, sum:=0,
  581. #   while (/= sum (setq sum (+ sum (/ b i)))) do i:=i+1, b:=b*a.
  582. #   Ergebnis sum.
  583. # Sonst setze y := sqrt(x), berechne rekursiv z:=ln(y)
  584. #   und liefere 2*z = (scale-float z 1).
  585. # Aufwand: asymptotisch d^2.5 .
  586.   local object F_lnx_F(x)
  587.     var reg4 object x;
  588.     { pushSTACK(x);
  589.       x = R_R_minus_R(x,Fixnum_1); # y := (- x 1)
  590.       if (R_zerop(x)) { skipSTACK(1); return x; } # y=0.0 -> y als Ergebnis
  591.       pushSTACK(x);
  592.       # Stackaufbau: x, y.
  593.      {var reg6 uintL d = F_float_digits(x);
  594.       var reg5 sintL e = F_exponent_L(x);
  595.        if (e <= (sintL)(-d)) # e <= -d ?
  596.          { x = STACK_0; skipSTACK(2); return x; } # ja -> y als Ergebnis
  597.       { var reg2 object k = Fixnum_0; # RekursionszΣhler k:=0
  598.        {# Bei e <= -1-floor(sqrt(d)) kann die Potenzreihe angewandt werden.
  599.         var reg3 sintL e_limit = -1-UL_sqrt_UW(d); # -1-floor(sqrt(d))
  600.         while (e > e_limit)
  601.           # e > -1-floor(sqrt(d)) -> mu▀ |y| verkleinern.
  602.           { var reg1 object x = F_sqrt_F(STACK_1); STACK_1 = x; # x := (sqrt x)
  603.             x = R_R_minus_R(x,Fixnum_1); STACK_0 = x; # y := (- x 1) und
  604.             e = F_exponent_L(x); # e neu berechnen
  605.             k = fixnum_inc(k,1); # k:=k+1
  606.        }  }
  607.        # Stackaufbau: x, y.
  608.        # Potenzreihe anwenden:
  609.        {var reg2 object i = Fixnum_1;
  610.         pushSTACK(I_F_float_F(Fixnum_0,STACK_1)); # sum := (float 0 x)
  611.         STACK_2 = F_minus_F(STACK_1); # a := -y, b := y
  612.         # Stackaufbau: a, b, sum.
  613.         loop
  614.           { var reg1 object temp;
  615.             temp = R_R_durch_R(STACK_1,i); # (/ b i)
  616.             temp = F_F_plus_F(STACK_0,temp); # (+ sum (/ b i))
  617.             if (eql(STACK_0,temp)) # = sum ?
  618.               break; # ja -> Potenzreihe abbrechen
  619.             STACK_0 = temp;
  620.             STACK_1 = F_F_mal_F(STACK_1,STACK_2); # b := b*a
  621.             i = fixnum_inc(i,1); # i := i+1
  622.        }  }
  623.        {var reg1 object erg = STACK_0; # sum als Ergebnis
  624.         skipSTACK(3);
  625.         return F_I_scale_float_F(erg,k); # wegen Rekursion noch mal 2^k
  626.     }}}}
  627.  
  628. # ln2_F_float_F(f) liefert die Zahl ln(2) im selben Float-Format wie f.
  629. # kann GC ausl÷sen
  630.   local object ln2_F_float_F (object f);
  631.   local object ln2_F_float_F(f)
  632.     var reg3 object f;
  633.     { var reg4 object ln2 = O(LF_ln2);
  634.       floatcase(f,
  635.                 { return LF_to_SF(ln2); },
  636.                 { return LF_to_FF(ln2); },
  637.                 { return LF_to_DF(ln2); },
  638.                 ;
  639.                );
  640.      {var reg2 uintC f_len = TheLfloat(f)->len; # gewⁿnschte LΣnge von ln(2)
  641.       {var reg1 uintC len = TheLfloat(ln2)->len; # vorhandene LΣnge
  642.        if (f_len < len) { return LF_shorten_LF(ln2,f_len); }
  643.        if (f_len == len) { return ln2; }
  644.       }
  645.       # gewⁿnschte > vorhandene LΣnge -> mu▀ nachberechnen:
  646.       {var reg4 uintC len = lf_len_extend(f_len); # einige Digits mehr verlangen
  647.        var reg1 object temp = F_lnx_F(I_to_LF(fixnum(2),len)); # (ln 2.0)
  648.        # temp = ln(2) ist fertig.
  649.        return O(LF_ln2) = LF_shorten_LF(temp,f_len); # wieder verkⁿrzen, als LF_ln2 abspeichern
  650.     }}}
  651.  
  652. # Vergr÷▀ert eine Long-Float-LΣnge n, so da▀ aus d = intDsize*n
  653. # mindestens d+sqrt(d)+2+(LF_exp_len-1) wird.
  654. # Allgemein: intDsize*n + sqrt(intDsize*n) + 2 + 31 < intDsize*(n+inc)
  655. # <==>       sqrt(intDsize*n) + 33 < intDsize*inc
  656. # <==>       sqrt(intDsize*n) < intDsize*inc - 33
  657. # <==>       intDsize*n < intDsize^2*inc^2 - 66*intDsize*inc + 1089
  658. # <==>       n <= intDsize*inc^2 - 66*inc + floor(1089/intDsize)
  659.   local uintC lf_len_extend2 (uintC n);
  660.   local uintC lf_len_extend2(n)
  661.     var reg1 uintC n;
  662.     { var reg2 uintC inc =
  663.         #define FITS(n,k)  ((n) <= (uintL)((intDsize*(k)-66)*(k)+floor(1089,intDsize)))
  664.         #define n_max  (uintL)(bitm(intCsize)-1)
  665.         #define TEST(i)  FITS(n_max,1UL<<i) || FITS(n,1UL<<i) ? 1UL<<i :
  666.         TEST(0) TEST(1) TEST(2) TEST(3) TEST(4) TEST(5) TEST(6) TEST(7)
  667.         TEST(8) TEST(9) TEST(10) TEST(11) TEST(12) TEST(13)
  668.         (fehler_LF_toolong(),0);
  669.         #undef TEST
  670.         #undef n_max
  671.         #undef FITS
  672.       if ((n = n+inc) < inc) { fehler_LF_toolong(); }
  673.       return n;
  674.     }
  675.  
  676. # F_extend2_F(x) erweitert die Genauigkeit eines Floats x um eine Stufe
  677. # SF -> FF -> DF -> LF(4) -> LF(5) -> LF(6) -> ...
  678. # Ein Float mit d Mantissenbits und l Exponentenbits wird so zu einem Float
  679. # mit mindestens d+sqrt(d)+2+(l-1) Mantissenbits.
  680. # SF -> DF wegen 17+sqrt(17)+2+7 = 30.2 < 53
  681. # FF -> DF wegen 24+sqrt(24)+2+7 = 37.9 < 53
  682. # DF -> LF(5) wegen 53+sqrt(53)+2+10 = 72.3 < 80
  683. # kann GC ausl÷sen
  684.   local object F_extend2_F (object x);
  685.   local object F_extend2_F(x)
  686.     var reg1 object x;
  687.     { floatcase(x,
  688.                 { return SF_to_DF(x); }, # 17+sqrt(17)+2+7 = 30.2 < 53
  689.                 { return FF_to_DF(x); }, # 24+sqrt(24)+2+7 = 37.9 < 53
  690.                 { return DF_to_LF(x,ceiling(73,intDsize)); }, # 53+sqrt(53)+2+10 = 72.3 < 73
  691.                 { return LF_extend_LF(x,lf_len_extend2(TheLfloat(x)->len)); }
  692.                );
  693.     }
  694.  
  695. # R_ln_R(x) liefert zu einer reellen Zahl x>0 die Zahl ln(x).
  696. # kann GC ausl÷sen
  697.   local object R_ln_R (object x);
  698. # Methode:
  699. # x rational -> bei x=1 0 als Ergebnis, sonst x in Float umwandeln.
  700. # x Float ->
  701. #   d := (float-digits x),
  702. #   Genauigkeit um sqrt(d)+max(integer-length(e)) Bits erh÷hen,
  703. #   (m,e) := (decode-float x), so da▀ 1/2 <= m < 1.
  704. #   m<2/3 -> m:=2m, e:=e-1, so da▀ 2/3 <= m <= 4/3.
  705. #   ln(m) errechnen, ln(x)=ln(m)+e*ln(2) als Ergebnis.
  706.   local object R_ln_R(x)
  707.     var reg2 object x;
  708.     { if (R_rationalp(x))
  709.         { if (eq(x,Fixnum_1)) { return Fixnum_0; } # x=1 -> 0 als Ergebnis
  710.           x = RA_float_F(x); # sonst in Float umwandeln
  711.         }
  712.       # x Float
  713.       pushSTACK(x); # x retten
  714.       x = F_extend2_F(x); # Rechengenauigkeit erh÷hen
  715.       F_decode_float_F_I_F(x); # m,e,s bestimmen
  716.       # Stackaufbau: x, m, e, s.
  717.       if (F_F_comp(STACK_2,
  718.                    make_SF(0,0+SF_exp_mid,floor(bit(SF_mant_len+2),3)) # Short-Float 2/3
  719.                   )
  720.           < 0
  721.          ) # m < 2/3 ->
  722.         { STACK_2 = F_I_scale_float_F(STACK_2,Fixnum_1); # m verdoppeln
  723.           STACK_1 = I_minus1_plus_I(STACK_1); # e decrementieren
  724.         }
  725.       STACK_2 = F_lnx_F(STACK_2); # ln(m) im genaueren Float-Format errechnen
  726.      {var reg1 object temp;
  727.       temp = ln2_F_float_F(STACK_0); # ln(2) im genaueren Float-Format
  728.       temp = R_R_mal_R(STACK_1,temp); # e*ln(2)
  729.       temp = R_R_plus_R(STACK_2,temp); # ln(m)+e*ln(2)
  730.       temp = F_F_float_F(temp,STACK_3); # (float ... x)
  731.       skipSTACK(4);
  732.       return temp;
  733.     }}
  734.   #define F_ln_F  R_ln_R
  735.  
  736. # I_I_log_RA(a,b) liefert zu Integers a>0, b>1 den Logarithmus log(a,b)
  737. # als exakte rationale Zahl, oder nullobj wenn er irrational ist.
  738. # kann GC ausl÷sen
  739.   local object I_I_log_RA (object a, object b);
  740. # Methode:
  741. # log(a,b) soll Bruch c/d mit teilerfremdem c>=0,d>0 ergeben.
  742. # a=1 -> c=0, d=1.
  743. # a>=b -> Dividiere a durch b. Rest da -> geht nicht.
  744. #         Sonst log(a,b) = 1+log(a/b,b).
  745. #         Berechne also c/d := log(a/b,b) und setze c:=c+d.
  746. # 1<a<b -> log(a,b) = 1/log(b,a).
  747. #         Berechne also c/d := log(b,a) und vertausche c und d.
  748. # Man konstruiert hierbei eigentlich die Kettenbruchentwicklung von c/d.
  749. # Wegen a>=2^c, b>=2^d sind c,d < (integer-length a,b) < intDsize*2^intCsize.
  750. # Entrekursiviert:
  751. # Wir werden (a,b) und damit auch c/d = log(a/b) verΣndern.
  752. # Invariante: Statt (c,d) wollen wir (uc*c+ud*d,vc*c+vd*d) zurⁿckliefern.
  753. # uc:=1, ud:=0, vc:=0, vd:=1.
  754. # Solange a>1,
  755. #   a>=b -> Dividiere a durch b. Rest da -> geht nicht.
  756. #           Sonst a:=a/b, und (fⁿr spΣter c:=c+d) uc:=uc+vc, ud:=ud+vd.
  757. #   1<a<b -> vertausche a und b, uc und vc, ud und vd.
  758. # Liefere (ud,vd), der Bruch ud/vd ist gekⁿrzt.
  759. # Offenbar braucht man uc und vc nicht.
  760.   local object I_I_log_RA(a,b)
  761.     var reg1 object a;
  762.     var reg2 object b;
  763.     { var reg3 uintL u = 0; # ud
  764.       var reg4 uintL v = 1; # vd
  765.       loop
  766.         { if (eq(a,Fixnum_1)) break; # a=1 -> Rekursion zu Ende
  767.           if (I_I_comp(a,b) >=0)
  768.             # a>=b
  769.             { pushSTACK(b);
  770.               I_I_divide_I_I(a,b); # a durch b dividieren
  771.               if (!eq(STACK_0,Fixnum_0)) # Rest /=0 ?
  772.                 { skipSTACK(3); return nullobj; } # -> fertig
  773.               a = STACK_1; b = STACK_2; skipSTACK(3); # a:=a/b
  774.               u = u + v; # ud:=ud+vd
  775.             }
  776.             else
  777.             # 1<a<b -> a und b vertauschen
  778.             {{var reg1 object temp = a; a = b; b = temp; }
  779.              {var reg1 uintL temp = u; u = v; v = temp; } # ud und vd vertauschen
  780.             }
  781.         }
  782.       # a=1 -> c=0,d=1 -> Ergebnis ud/vd
  783.       pushSTACK(UL_to_I(u)); # u als Integer
  784.      {var reg5 object y = UL_to_I(v); # v als Integer
  785.       var reg6 object x = popSTACK();
  786.       return I_I_to_RA(x,y);
  787.     }}
  788.  
  789. # R_R_log_R(a,b) liefert zu reellen Zahlen a>0, b>0 die Zahl
  790. # log(a,b)=ln(a)/ln(b).
  791. # Ergebnis rational nur, wenn a=1 oder a und b rational.
  792. # kann GC ausl÷sen
  793.   local object R_R_log_R (object a, object b);
  794. # Methode:
  795. # a und b rational:
  796. #   b=1 -> Error
  797. #   a=1 -> Ergebnis 0
  798. #   b Integer:
  799. #     a Integer: log(a,b) rational errechenbar -> liefern
  800. #     a Ratio: a=a1/a2 mit a1>0, a2>1.
  801. #              a1=1 und log(a2,b) rational errechenbar -> -log(a2,b) liefern
  802. #   b Ratio: a=a1/a2, b=b1/b2 mit a1>0, a2>0, b1>0, b2>1.
  803. #            log(a2,b2) rational errechenbar ->
  804. #               b1=1 -> bei a1=1 liefern, sonst nicht.
  805. #               b1>1 -> log(a1,b1) rational errechenbar und
  806. #                       log(a1,b1)=log(a2,b2) -> liefern, sonst nicht.
  807. #            sonst a1,a2 vertauschen:
  808. #              log(a2/a1,b1/b2) versuchen (wie oben) ->
  809. #                -log(a2/a1,b1/b2) liefern
  810. #   Sonst a und b in Floats umwandeln.
  811. # a Float, b rational -> bei b=1 Error, sonst b := (float b a)
  812. # a rational, b Float -> bei a=1 Ergebnis 0, sonst a := (float a b)
  813. # a,b Floats -> log(a,b) = ln(a)/ln(b)
  814.   local object R_R_log_R(a,b)
  815.     var reg2 object a;
  816.     var reg1 object b;
  817.     { pushSTACK(a); pushSTACK(b);
  818.       # Stackaufbau: a, b.
  819.       if (R_rationalp(b))
  820.         # b rational
  821.         { if (eq(b,Fixnum_1)) { divide_0(); } # b=1 -> Error
  822.           if (R_rationalp(a))
  823.             # a,b beide rational
  824.             { if (eq(a,Fixnum_1)) { skipSTACK(2); return Fixnum_0; } # a=1 -> Ergebnis 0
  825.               if (RA_integerp(b))
  826.                 # b Integer
  827.                 { if (RA_integerp(a))
  828.                     # a,b beide Integers
  829.                     { var reg3 object x = I_I_log_RA(a,b); # rationalen log(a,b) versuchen
  830.                       if (!eq(x,nullobj)) { skipSTACK(2); return x; }
  831.                     }
  832.                     else
  833.                     # a Ratio, b Integer
  834.                     { if (eq(TheRatio(a)->rt_num,Fixnum_1)) # a1=1
  835.                         { var reg3 object x = I_I_log_RA(TheRatio(a)->rt_den,b); # rationalen log(a2,b) versuchen
  836.                           if (!eq(x,nullobj)) { skipSTACK(2); return RA_minus_RA(x); }
  837.                 }   }   }
  838.                 else
  839.                 # a rational, b Ratio
  840.                 { if (RA_integerp(a))
  841.                     { pushSTACK(a); pushSTACK(a=Fixnum_1); }
  842.                     else
  843.                     { pushSTACK(TheRatio(a)->rt_num); pushSTACK(a=TheRatio(a)->rt_den); }
  844.                   pushSTACK(TheRatio(b)->rt_num); pushSTACK(b=TheRatio(b)->rt_den);
  845.                   # Stackaufbau: a, b, a1>0, a2>0, b1>0, b2>1.
  846.                   {var reg3 object x = I_I_log_RA(a,b); # rationalen log(a2,b2) versuchen
  847.                    if (!eq(x,nullobj))
  848.                      { if (eq(STACK_1,Fixnum_1)) # b1=1 ?
  849.                          { if (eq(STACK_3,Fixnum_1)) # a1=1 ?
  850.                              { skipSTACK(6); return x; } # ja -> x liefern
  851.                          }
  852.                          else
  853.                          { pushSTACK(x);
  854.                           {var reg4 object y = I_I_log_RA(STACK_(3+1),STACK_(1+1)); # rationalen log(a1,b1) versuchen
  855.                            if ((!eq(y,nullobj)) && eql(STACK_0,y)) # x=y ?
  856.                              { x = STACK_0; skipSTACK(6+1); return x; }
  857.                            skipSTACK(1);
  858.                   }  }   }}
  859.                   {var reg3 object x = I_I_log_RA(STACK_3,STACK_0); # rationalen log(a1,b2) versuchen
  860.                    if (!eq(x,nullobj))
  861.                      { if (eq(STACK_1,Fixnum_1)) # b1=1 ?
  862.                          { if (eq(STACK_2,Fixnum_1)) # a2=1 ?
  863.                              { skipSTACK(6); return RA_minus_RA(x); } # ja -> -x liefern
  864.                          }
  865.                          else
  866.                          { pushSTACK(x);
  867.                           {var reg4 object y = I_I_log_RA(STACK_(2+1),STACK_(1+1)); # rationalen log(a2,b1) versuchen
  868.                            if ((!eq(y,nullobj)) && eql(STACK_0,y)) # x=y ?
  869.                              { x = STACK_0; skipSTACK(6+1); return RA_minus_RA(x); }
  870.                            skipSTACK(1);
  871.                   }  }   }}
  872.                   skipSTACK(4);
  873.                 }
  874.               # a,b beide in Floats umwandeln:
  875.               STACK_1 = RA_float_F(STACK_1); STACK_0 = RA_float_F(STACK_0);
  876.             }
  877.             else
  878.             # a Float
  879.             { STACK_0 = RA_F_float_F(b,a); } # b := (float b a)
  880.         }
  881.         else
  882.         # b Float
  883.         { if (R_rationalp(a))
  884.             # a rational
  885.             { if (eq(a,Fixnum_1)) { skipSTACK(2); return Fixnum_0; } # a=1 -> Ergebnis 0
  886.               STACK_1 = RA_F_float_F(a,b); # a := (float a b)
  887.         }   }
  888.       # Nun a,b beide Floats.
  889.       STACK_1 = R_ln_R(STACK_1); # (ln a) errechnen
  890.      {var reg3 object lnb = R_ln_R(popSTACK()); # (ln b) errechnen
  891.       return F_F_durch_F(popSTACK(),lnb); # (/ (ln a) (ln b)) als Ergebnis
  892.     }}
  893.  
  894. # F_expx_F(x) liefert zu einem Float x (betragsmΣ▀ig <1) exp(x) als Float.
  895. # kann GC ausl÷sen
  896.   local object F_expx_F (object x);
  897. # Methode:
  898. # e := Exponent aus (decode-float x), d := (float-digits x)
  899. # Bei x=0.0 oder e<-d liefere 1.0
  900. #   (denn bei e<=-d-1 ist abs(exp(x)-1) = abs(x)+O(x^2) < 2^(-d-1),
  901. #    also ist exp(x), auf d Bits gerundet, gleich 1.0).
  902. # Bei e<=-sqrt(d) verwende die Potenzreihe
  903. #   exp(x) = sum(j=0..inf,x^j/j!):
  904. #   b:=1, i:=0, sum:=0,
  905. #   while (/= sum (setq sum (+ sum b))) do b:=b*x/(i+1), i:=i+1.
  906. #   Ergebnis sum.
  907. # Sonst setze y := x/2 = (scale-float x -1),
  908. #   berechne rekursiv z:=exp(y) und liefere z^2.
  909. # Aufwand: asymptotisch d^2.5 .
  910.   local object F_expx_F (x)
  911.     var reg4 object x;
  912.     { if (R_zerop(x)) { return I_F_float_F(Fixnum_1,x); }
  913.      {var reg6 uintL d = F_float_digits(x);
  914.       var reg5 sintL e = F_exponent_L(x);
  915.       if (e < (sintL)(-d)) # e < -d ?
  916.         { return I_F_float_F(Fixnum_1,x); } # ja -> 1.0 als Ergebnis
  917.       pushSTACK(x);
  918.       # Stackaufbau: x.
  919.       {var reg3 uintL k = 0; # RekursionszΣhler k:=0
  920.        {# Bei e <= -1-floor(sqrt(d)) kann die Potenzreihe angewandt werden.
  921.         var reg1 sintL e_limit = -1-UL_sqrt_UW(d); # -1-floor(sqrt(d))
  922.         if (e > e_limit)
  923.           # e > -1-floor(sqrt(d)) -> mu▀ |x| verkleinern.
  924.           { k = e - e_limit;
  925.            {var reg1 object temp = L_to_I((sintL)(-k));
  926.             STACK_0 = F_I_scale_float_F(STACK_0,temp); # x := x/2^k
  927.             # Neuer Exponent = e-k = e_limit.
  928.        }  }}
  929.        # Potenzreihe anwenden:
  930.        {var reg2 object i = Fixnum_0;
  931.         pushSTACK(I_F_float_F(Fixnum_1,STACK_0)); # b := (float 1 x)
  932.         pushSTACK(I_F_float_F(Fixnum_0,STACK_1)); # sum := (float 0 x)
  933.         # Stackaufbau: x, b, sum.
  934.         loop
  935.           { var reg1 object temp;
  936.             temp = F_F_plus_F(STACK_0,STACK_1); # (+ sum b)
  937.             if (eql(STACK_0,temp)) # = sum ?
  938.               break; # ja -> Potenzreihe abbrechen
  939.             STACK_0 = temp;
  940.             temp = F_F_mal_F(STACK_1,STACK_2); # b*x
  941.             i = fixnum_inc(i,1); # i := i+1
  942.             STACK_1 = R_R_durch_R(temp,i); # b := b*x/i
  943.        }  }
  944.        {var reg1 object z = STACK_0; # sum als Ergebnis
  945.         skipSTACK(3);
  946.         # Wegen Rekursion noch k mal quadrieren:
  947.         dotimesL(k,k, { z = F_F_mal_F(z,z); } );
  948.         return z;
  949.     }}}}
  950.  
  951. # R_exp_R(x) liefert zu einer reellen Zahl x die Zahl exp(x).
  952. # kann GC ausl÷sen
  953.   local object R_exp_R (object x);
  954. # Methode:
  955. # x rational -> bei x=0 1 als Ergebnis, sonst x in Float umwandeln.
  956. # x Float ->
  957. #   d := (float-digits x),
  958. #   Genauigkeit um sqrt(d)+max(integer-length(e)) Bits erh÷hen,
  959. #   (q,r) := (floor x ln(2))
  960. #   Ergebnis ist exp(q*ln(2)+r) = (scale-float exp(r) q).
  961.   local object R_exp_R(x)
  962.     var reg1 object x;
  963.     { if (R_rationalp(x))
  964.         # x rational
  965.         { if (eq(x,Fixnum_0)) { return Fixnum_1; } # x=0 -> 1 als Ergebnis
  966.           x = RA_float_F(x); # sonst in Float umwandeln
  967.         }
  968.       # x Float
  969.       pushSTACK(x); # x retten
  970.       x = F_extend2_F(x); # Genauigkeit vergr÷▀ern
  971.       # durch ln(2) dividieren (bei 0<=x<1/2 kann man sofort q:=0 setzen)
  972.       if ((!R_minusp(x)) && (F_exponent_L(x)<0))
  973.         { pushSTACK(Fixnum_0); pushSTACK(x); } # x>=0, Exponent <0 -> 0<=x<1/2 -> Division unn÷tig
  974.         else
  975.         { pushSTACK(x);
  976.          {var reg1 object ln2 = ln2_F_float_F(x); # ln(2) mit hinreichender Genauigkeit
  977.           x = popSTACK();
  978.           F_F_floor_I_F(x,ln2); # x durch ln(2) dividieren
  979.         }}
  980.       # Stackaufbau: originales x, q, r.
  981.      {var reg2 object temp = F_expx_F(STACK_0); # exp(r)
  982.       temp = F_I_scale_float_F(temp,STACK_1); # mal 2^q
  983.       temp = F_F_float_F(temp,STACK_2); # (float ... x) als Ergebnis
  984.       skipSTACK(3);
  985.       return temp;
  986.     }}
  987.  
  988. # R_sinh_R(x) liefert zu einer reellen Zahl x die Zahl sinh(x).
  989. # kann GC ausl÷sen
  990.   local object R_sinh_R (object x);
  991. # Methode:
  992. # x rational -> bei x=0 0 als Ergebnis, sonst x in Float umwandeln.
  993. # x Float -> Genauigkeit erh÷hen,
  994. #   e := Exponent aus (decode-float x)
  995. #   falls e<=0: (sinh(x)/x)^2 errechnen, Wurzel ziehen, mit x multiplizieren.
  996. #   falls e>0: y:=exp(x) errechnen, (scale-float (- y (/ y)) -1) bilden.
  997.   local object R_sinh_R(x)
  998.     var reg2 object x;
  999.     { if (R_rationalp(x))
  1000.         # x rational
  1001.         { if (eq(x,Fixnum_0)) { return x; } # x=0 -> 0 als Ergebnis
  1002.           x = RA_float_F(x); # sonst in Float umwandeln
  1003.         }
  1004.       # x Float
  1005.       if (F_exponent_L(x)<=0) # Exponent e abtesten
  1006.         # e<=0
  1007.         { var reg1 object temp;
  1008.           pushSTACK(x);
  1009.           pushSTACK(temp = F_extend_F(x)); # Rechengenauigkeit erh÷hen
  1010.           temp = F_sqrt_F(F_sinhx_F(x)); # Wurzel aus (sinh(x)/x)^2
  1011.           temp = F_F_mal_F(temp,STACK_0); # mit genauerem x multiplizieren
  1012.           temp = F_F_float_F(temp,STACK_1); # und wieder runden
  1013.           skipSTACK(2);
  1014.           return temp;
  1015.         }
  1016.         else
  1017.         # e>0 -> verwende exp(x)
  1018.         { var reg1 object temp;
  1019.           pushSTACK(temp = R_exp_R(x)); # y:=exp(x)
  1020.           temp = F_durch_F(temp); # (/ y)
  1021.           temp = F_F_minus_F(popSTACK(),temp); # von y subtrahieren
  1022.           return F_I_scale_float_F(temp,Fixnum_minus1); # (scale-float ... -1)
  1023.     }   }
  1024.  
  1025. # R_cosh_R(x) liefert zu einer reellen Zahl x die Zahl cosh(x).
  1026. # kann GC ausl÷sen
  1027.   local object R_cosh_R (object x);
  1028. # Methode:
  1029. # x rational -> bei x=0 1 als Ergebnis, sonst x in Float umwandeln.
  1030. # x Float -> Genauigkeit erh÷hen,
  1031. #   e := Exponent aus (decode-float x), d := (float-digits x)
  1032. #   falls x=0.0 oder e<=(1-d)/2 liefere 1.0
  1033. #     (denn bei e<=(1-d)/2 ist 1 <= cosh(x) = 1+x^2/2+... < 1+2^(-d),
  1034. #      also ist cosh(x), auf d Bits gerundet, gleich 1.0).
  1035. #   falls e<=0:
  1036. #     y := x/2 = (scale-float x -1), (sinh(y)/y)^2 errechnen,
  1037. #     cosh(x) = 1+x*y*(sinh(y)/y)^2 errechnen.
  1038. #   falls e>0: y:=exp(x) errechnen, (scale-float (+ y (/ y)) -1) bilden.
  1039.   local object R_cosh_R(x)
  1040.     var reg2 object x;
  1041.     { if (R_rationalp(x))
  1042.         # x rational
  1043.         { if (eq(x,Fixnum_0)) { return Fixnum_1; } # x=0 -> 1 als Ergebnis
  1044.           x = RA_float_F(x); # sonst in Float umwandeln
  1045.         }
  1046.       # x Float
  1047.       {var reg3 sintL e = F_exponent_L(x);
  1048.        if (e > 0)
  1049.          # e>0 -> verwende exp(x)
  1050.          { var reg1 object temp;
  1051.            pushSTACK(temp = R_exp_R(x)); # y:=exp(x)
  1052.            temp = F_durch_F(temp); # (/ y)
  1053.            temp = F_F_plus_F(popSTACK(),temp); # zu y addieren
  1054.            return F_I_scale_float_F(temp,Fixnum_minus1); # (scale-float ... -1)
  1055.          }
  1056.          else
  1057.          # e<=0
  1058.          { if (R_zerop(x)) { return I_F_float_F(Fixnum_1,x); }
  1059.           {var reg1 uintL d = F_float_digits(x);
  1060.            if (e <= (sintL)(1-d)>>1) # e <= (1-d)/2 <==> e <= -ceiling((d-1)/2) ?
  1061.              { return I_F_float_F(Fixnum_1,x); } # ja -> 1.0 als Ergebnis
  1062.           }
  1063.           {var reg1 object temp;
  1064.            pushSTACK(x);
  1065.            pushSTACK(temp = F_extend_F(x)); # Rechengenauigkeit erh÷hen
  1066.            pushSTACK(temp = F_I_scale_float_F(temp,Fixnum_minus1)); # y=(scale-float x -1)
  1067.            temp = F_sinhx_F(temp); # (sinh(y)/y)^2
  1068.            temp = F_F_mal_F(STACK_0,temp); # mit y multiplizieren
  1069.            temp = F_F_mal_F(STACK_1,temp); # mit x multiplizieren
  1070.            temp = R_R_plus_R(Fixnum_1,temp); # 1 addieren
  1071.            temp = F_F_float_F(temp,STACK_2); # und wieder runden
  1072.            skipSTACK(3);
  1073.            return temp;
  1074.          }}
  1075.     } }
  1076.  
  1077. # R_cosh_sinh_R_R(x) liefert ((cosh x),(sinh x)), beide Werte auf dem Stack.
  1078. # kann GC ausl÷sen
  1079.   local void R_cosh_sinh_R_R (object x);
  1080. # Methode:
  1081. # x rational -> bei x=0 (1,0) als Ergebnis, sonst x in Float umwandeln.
  1082. # x Float -> Genauigkeit erh÷hen,
  1083. #   e := Exponent aus (decode-float x), d := (float-digits x)
  1084. #   falls x=0.0 oder e<=(1-d)/2 liefere (1.0,x)
  1085. #     (denn bei e<=(1-d)/2 ist
  1086. #      1 <= sinh(x)/x < cosh(x) = 1+x^2/2+... < 1+2^(-d),
  1087. #      also ist cosh(x), auf d Bits gerundet, gleich 1.0
  1088. #      und sinh(x), auf d Bits gerundet, gleich x).
  1089. #   falls e<=0:
  1090. #     y:=(sinh(x)/x)^2 errechnen,
  1091. #     cosh(x) = sqrt(1+x^2*y) und sinh(x) = x*sqrt(y) errechnen.
  1092. #   falls e>0: y:=exp(x) errechnen,
  1093. #     (scale-float (+ y (/ y)) -1) und (scale-float (- y (/ y)) -1) bilden.
  1094. #   Genauigkeit wieder verringern.
  1095.   local void R_cosh_sinh_R_R(x)
  1096.     var reg3 object x;
  1097.     { if (R_rationalp(x))
  1098.         # x rational
  1099.         { if (eq(x,Fixnum_0)) { pushSTACK(Fixnum_1); pushSTACK(Fixnum_0); return; } # x=0 -> (1,0) als Ergebnis
  1100.           x = RA_float_F(x); # sonst in Float umwandeln
  1101.         }
  1102.       # x Float
  1103.       {var reg4 sintL e = F_exponent_L(x);
  1104.        if (e > 0)
  1105.          # e>0 -> verwende exp(x)
  1106.          { var reg1 object temp;
  1107.            pushSTACK(temp = R_exp_R(x)); # y:=exp(x)
  1108.            pushSTACK(temp = F_durch_F(temp)); # (/ y)
  1109.            # Stackaufbau: exp(x), exp(-x).
  1110.            temp = F_F_minus_F(STACK_1,temp); # von y subtrahieren
  1111.            temp = F_I_scale_float_F(temp,Fixnum_minus1); # (scale-float ... -1)
  1112.           {var reg2 object temp2 = STACK_0;
  1113.            STACK_0 = temp;
  1114.            temp = F_F_plus_F(STACK_1,temp2); # zu y addieren
  1115.            STACK_1 = F_I_scale_float_F(temp,Fixnum_minus1); # (scale-float ... -1)
  1116.            return;
  1117.          }}
  1118.          else
  1119.          # e<=0
  1120.          { if (R_zerop(x)
  1121.                || (e <= (sintL)(1-F_float_digits(x))>>1) # e <= (1-d)/2 <==> e <= -ceiling((d-1)/2) ?
  1122.               )
  1123.              { pushSTACK(x); pushSTACK(x); STACK_1 = I_F_float_F(Fixnum_1,x); return; }
  1124.           {var reg1 object temp;
  1125.            pushSTACK(x);
  1126.            pushSTACK(temp = F_extend_F(x)); # Rechengenauigkeit erh÷hen
  1127.            pushSTACK(F_F_mal_F(temp,temp)); # x*x
  1128.            pushSTACK(temp = F_sinhx_F(STACK_1)); # y:=(sinh(x)/x)^2
  1129.            # Stackaufbau: originales x, x, x^2, y.
  1130.            temp = F_sqrt_F(temp); # sqrt(y) = sinh(x)/x
  1131.            temp = F_F_mal_F(STACK_2,temp); # x*sqrt(y) = sinh(x)
  1132.            STACK_2 = F_F_float_F(temp,STACK_3); # und wieder runden
  1133.            temp = F_F_mal_F(STACK_1,STACK_0); # x^2*y
  1134.            temp = F_sqrt_F(R_R_plus_R(Fixnum_1,temp)); # sqrt(1+x^2*y)
  1135.            STACK_3 = F_F_float_F(temp,STACK_3); # und wieder runden
  1136.            skipSTACK(2); return;
  1137.          }}
  1138.     } }
  1139.  
  1140.