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