## clisp-devel

 Numeric inconsistency From: Raymond Toy - 2002-07-17 16:32:21 ```Consider this, using clisp 2.28 on a Solaris system: [1]> (tan (/ pi 2)) -3.986797629004264116L19 [2]> (cos (/ pi 2)) -2.5082788076048218878L-20 [3]> (/ (cos (/ pi 2))) -3.9867976277920596936L19 Note that tan(pi/2) differs quite a bit from 1/cos(pi/2). While I don't consider this a bug, I do think it's rather inconsistent since sin(pi/2) is essentially 1 so that tan(pi/2) should be essentially equal to 1/cos(pi/2). It almost seems as if some intermediate calculation of tan used single-precision since the results are the same to about 7-8 digits. Ray ```
 Re: Numeric inconsistency From: Sam Steingold - 2002-07-17 16:55:25 ```> * In message <4nfzyi2u45.fsf@...> > * On the subject of "Numeric inconsistency" > * Sent on 17 Jul 2002 12:32:10 -0400 > * Honorable Raymond Toy writes: > > Consider this, using clisp 2.28 on a Solaris system: > > [1]> (tan (/ pi 2)) > -3.986797629004264116L19 > [2]> (cos (/ pi 2)) > -2.5082788076048218878L-20 > [3]> (/ (cos (/ pi 2))) > -3.9867976277920596936L19 > > Note that tan(pi/2) differs quite a bit from 1/cos(pi/2). > While I don't consider this a bug, I do. > I do think it's rather inconsistent since sin(pi/2) is essentially 1 > so that tan(pi/2) should be essentially equal to 1/cos(pi/2). It > almost seems as if some intermediate calculation of tan used > single-precision since the results are the same to about 7-8 digits. I looked into this. the reason is that when CLISP computes TAN, it takes an SQRT and (SQRT 0) does not make sense (I wonder how many people would understand what I mean here :-) -- Sam Steingold (http://www.podval.org/~sds) running RedHat7.2 GNU/Linux ; ; ; ; ; (let((a'(list'let(list(list'a(list'quote a)))a)))`(let((a(quote ,a))),a)) ```
 Re: Numeric inconsistency From: Raymond Toy - 2002-07-17 17:08:58 ```>>>>> "Sam" == Sam Steingold writes: Sam> I looked into this. Sam> the reason is that when CLISP computes TAN, it takes an SQRT and Sam> (SQRT 0) does not make sense (I wonder how many people would understand Sam> what I mean here :-) Certainly not me. I might if I could actually read German and find out where tan was implemented and how. :-) Ray ```
 Re: Numeric inconsistency From: Sam Steingold - 2002-07-17 17:32:06 ```> * In message <4ny9ca1duj.fsf@...> > * On the subject of "Re: Numeric inconsistency" > * Sent on 17 Jul 2002 13:08:52 -0400 > * Honorable Raymond Toy writes: > > >>>>> "Sam" == Sam Steingold writes: > > Sam> the reason is that when CLISP computes TAN, it takes an SQRT > Sam> and (SQRT 0) does not make sense (I wonder how many people > Sam> would understand what I mean here :-) > > Certainly not me. I might if I could actually read German and find > out where tan was implemented and how. :-) sqrt has infinite derivative at 0, so approximate computations of sqrt of small number will have infinite relative errors. -- Sam Steingold (http://www.podval.org/~sds) running RedHat7.2 GNU/Linux ; ; ; ; ; Don't use force -- get a bigger hammer. ```
 Re: Numeric inconsistency From: Raymond Toy - 2002-07-17 18:58:27 ```>>>>> "Sam" == Sam Steingold writes: >> * In message <4ny9ca1duj.fsf@...> >> * On the subject of "Re: Numeric inconsistency" >> * Sent on 17 Jul 2002 13:08:52 -0400 >> * Honorable Raymond Toy writes: >> >> >>>>> "Sam" == Sam Steingold writes: >> Sam> the reason is that when CLISP computes TAN, it takes an SQRT Sam> and (SQRT 0) does not make sense (I wonder how many people Sam> would understand what I mean here :-) >> >> Certainly not me. I might if I could actually read German and find >> out where tan was implemented and how. :-) Sam> sqrt has infinite derivative at 0, so approximate computations of sqrt Sam> of small number will have infinite relative errors. Yes. If you could, can you tell me why the computation of tan requires the computation of sqrt? Don't bother if it's lengthy. Just curious. Ray ```
 Re: Numeric inconsistency From: Marco Baringer - 2002-07-17 17:50:47 ```Raymond Toy writes: > Note that tan(pi/2) differs quite a bit from 1/cos(pi/2). While I > don't consider this a bug, I do think it's rather inconsistent since > sin(pi/2) is essentially 1 so that tan(pi/2) should be essentially > equal to 1/cos(pi/2). It almost seems as if some intermediate you are asking for two infinte numbers, and CLISP is giving you two infinte numbers, or as close as CLISP can get to infinity. since it goes about calculating them in different ways it gives you different infinities, but this is what hapens when you ask for something which doesn't exist. the cool thing to do would be to implement "trigonometric" numbers ala complex and rationals...ie where pi was NOT 3.14159....but exactly 4arctan(1). -- -Marco Ring the bells that still can ring. Forget your perfect offering. There is a crack in everything. That's how the light gets in. -Leonard Cohen ```
 Re: Numeric inconsistency From: Raymond Toy - 2002-07-17 18:57:07 ```>>>>> "Marco" == Marco Baringer writes: Marco> Raymond Toy writes: >> Note that tan(pi/2) differs quite a bit from 1/cos(pi/2). While I >> don't consider this a bug, I do think it's rather inconsistent since >> sin(pi/2) is essentially 1 so that tan(pi/2) should be essentially >> equal to 1/cos(pi/2). It almost seems as if some intermediate Marco> you are asking for two infinte numbers, and CLISP is giving you two Marco> infinte numbers, or as close as CLISP can get to infinity. since it Marco> goes about calculating them in different ways it gives you different Marco> infinities, but this is what hapens when you ask for something which Marco> doesn't exist. No, I'm not asking for infinite numbers. I'm asking Clisp to compute tan (and cos) for the (mathematically) rational number (/ pi 2). And if Clisp really wanted to return the "closest" to infinity, most-positive-long-float would have been a better answer. Clisp returns it's best answer. That's fine. Note that I wasn't complaining that the numbers are wrong (CMUCL says the result is positive, not negative), just that tan and 1/cos gave significantly different numbers when I expected them to be much closer. They don't have to be exactly equal, but should be closer than the 7-8 digits it gets now. Note that CMUCL says (zerop (tan (/ pi 2)) (/ (cos (/ pi 2)))) is T, which actually exceeded my expectations. Marco> the cool thing to do would be to implement "trigonometric" numbers ala Marco> complex and rationals...ie where pi was NOT 3.14159....but exactly Marco> 4arctan(1). That's a different issue. :-) Maxima solves this pretty well. Ray ```
 Re[2]: Numeric inconsistency From: Arseny Slobodjuck - 2002-07-18 08:12:16 ```Hello Raymond, Thursday, July 18, 2002, 4:58:24 AM, you wrote: Raymond> Yes. If you could, can you tell me why the computation of tan Raymond> requires the computation of sqrt? Don't bother if it's lengthy. Just Raymond> curious. I have an supposition: it allows to compute only one of sin or cos, not both (supposing sqrt being calculated more quickly than trigonometric function). For example tg(a) = (+-)sqrt((1-cos(2a))/(1+cos(2a))). -- Best regards, Arseny ```
 Re: Re[2]: Numeric inconsistency From: Raymond Toy - 2002-07-18 17:11:41 ```>>>>> "Arseny" == Arseny Slobodjuck writes: Arseny> Hello Raymond, Arseny> Thursday, July 18, 2002, 4:58:24 AM, you wrote: Raymond> Yes. If you could, can you tell me why the computation of tan Raymond> requires the computation of sqrt? Don't bother if it's lengthy. Just Raymond> curious. Arseny> I have an supposition: it allows to compute only one of sin or cos, not both Arseny> (supposing sqrt being calculated more quickly than trigonometric function). Arseny> For example tg(a) = (+-)sqrt((1-cos(2a))/(1+cos(2a))). Sounds reasonable. Seems that this has lots of extra round-off problems than just computing cos(a), though. Ray ```
 Re[4]: Numeric inconsistency From: Arseny Slobodjuck - 2002-07-18 23:36:56 ```Hello, Raymond, Friday, July 19, 2002, 3:11:35 AM, you wrote: > Raymond> Yes. If you could, can you tell me why the computation of tan > Raymond> requires the computation of sqrt? Don't bother if it's lengthy. Just > Raymond> curious. > Arseny> I have an supposition: it allows to compute only one of sin or cos, not both > Arseny> (supposing sqrt being calculated more quickly than trigonometric function). > Arseny> For example tg(a) = (+-)sqrt((1-cos(2a))/(1+cos(2a))). > Sounds reasonable. Seems that this has lots of extra round-off > problems than just computing cos(a), though. This formula doesn't have sqrt(0) problem for tg(pi/2) (it haves sqrt(inf) one). If you mean sqrt(x - > 0) error, it can be fixed (if the code doesn't contain a workaround already). For example sqrt(x) = sqrt(n^2 * x) / n (There are possibly exist smarter workarounds, I don't know). -- Best regards, Arseny mailto:ampy@... ```
 Re: Re[4]: Numeric inconsistency From: Raymond Toy - 2002-07-19 13:14:51 ```>>>>> "Arseny" == Arseny Slobodjuck writes: Arseny> Hello, Raymond, Arseny> Friday, July 19, 2002, 3:11:35 AM, you wrote: Raymond> Yes. If you could, can you tell me why the computation of tan Raymond> requires the computation of sqrt? Don't bother if it's lengthy. Just Raymond> curious. Arseny> I have an supposition: it allows to compute only one of sin or cos, not both Arseny> (supposing sqrt being calculated more quickly than trigonometric function). Arseny> For example tg(a) = (+-)sqrt((1-cos(2a))/(1+cos(2a))). >> Sounds reasonable. Seems that this has lots of extra round-off >> problems than just computing cos(a), though. Arseny> This formula doesn't have sqrt(0) problem for tg(pi/2) (it haves Arseny> sqrt(inf) one). Yep. Arseny> If you mean sqrt(x - > 0) error, it can be fixed (if the code doesn't Arseny> contain a workaround already). For example sqrt(x) = sqrt(n^2 * x) / Arseny> n (There are possibly exist smarter workarounds, I don't know). Somehow these workarounds just seem to complicate things. Aren't they good algorithms that can comput sin and cos at the same time with only slightly more cost than just sin? I know the CORDIC algorithms do this. Ray ```
 Re: Numeric inconsistency From: Sam Steingold - 2002-09-09 21:45:12 ```the appended patch fixes the expt precision problem which you reported in another e-mail. I would greatly appreciate it if you could test it as much as possible. as far as TAN/COS are concerned, I am at a loss. I guess using R_cos_sin_R_R to compute just cos would be a solution... thanks. --=20 Sam Steingold (http://www.podval.org/~sds) running RedHat7.3 GNU/Linux ; ; ; ; ; As a computer, I find your faith in technology amusing. Index: src/realtran.d =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D RCS file: /cvsroot/clisp/clisp/src/realtran.d,v retrieving revision 1.8 diff -u -w -b -r1.8 realtran.d --- src/realtran.d 4 Sep 2002 19:48:24 -0000 1.8 +++ src/realtran.d 9 Sep 2002 21:41:43 -0000 @@ -91,7 +91,7 @@ return O(SF_pi), # pi as SF return O(FF_pi), # pi as FF return O(DF_pi), # pi as DF - return O(pi) , # pi as LF with default length + return pi_F_float_F(x), # pi as LF in the same format a= s x ,); # nothing to save } =20 @@ -494,6 +494,10 @@ }} =20 /* R_cos_sin_R_R(x) places ((cos x),(sin x)), on the Stack. + when start_p, this is a start of a computation, + so the accuracy will be raised + when end_p, this is also the end of the computation, + so the accuracy will be lowed back to this number can trigger GC Method: x rational -> if x=3D0 =3D=3D> (1,0), otherwise x =3D=3D> float. @@ -514,7 +518,7 @@ if q =3D 1 mod 4: ((- (sin r)), (cos r)) if q =3D 2 mod 4: ((- (cos r)), (- (sin r))) if q =3D 3 mod 4: ((sin r), (- (cos r))) */ -local void R_cos_sin_R_R (object x) { +local void R_cos_sin_R_R (object x,bool start_p,object *end_p) { if (R_rationalp(x)) { if (eq(x,Fixnum_0)) # x=3D0 -> return (1,0) { pushSTACK(Fixnum_1); pushSTACK(Fixnum_0); return; } @@ -522,14 +526,20 @@ } /* x -- float */ pushSTACK(x); /* save x */ - x =3D F_extend_F(x); /* increase computational accuracy */ + if (start_p) /* increase computational precision */ + x =3D F_extend_F(x); F_pi2_round_I_F(x); /* divide by pi/2 */ /* stack layout: Argument, q, r. */ x =3D STACK_0; if (R_zerop(x) /* r=3D0.0 -> cos=3D1.0+O(r^2), sin=3Dr+O(r^3) */ || (F_exponent_L(x) <=3D (sintL)(-F_float_digits(x))>>1)) { /* e <= =3D -d/2 <=3D=3D> e <=3D -ceiling(d/2) ? */ - pushSTACK(I_F_float_F(Fixnum_1,STACK_2)); /* cos=3D1 */ - pushSTACK(F_F_float_F(STACK_1,STACK_3)); /* sin=3Dr */ + if (end_p !=3D NULL) { + pushSTACK(RA_R_float_F(Fixnum_1,*end_p)); /* cos=3D1 */ + pushSTACK(F_R_float_F(STACK_1,*end_p)); /* sin=3Dr */ + } else { + pushSTACK(I_F_float_F(Fixnum_1,STACK_0)); /* cos=3D1 */ + pushSTACK(STACK_1); /* sin=3Dr */ + } } else { pushSTACK(F_I_scale_float_F(STACK_0,Fixnum_minus1)); /* s :=3D r/2 */ pushSTACK(F_sinx_F(STACK_0)); /* y :=3D (sin(s)/s)^2 */ @@ -537,7 +547,9 @@ x =3D F_F_mal_F(STACK_0,STACK_1); /* y*s */ x =3D F_F_mal_F(STACK_2,x); /* y*s*r */ x =3D R_R_minus_R(Fixnum_1,x); /* 1-y*s*r */ - pushSTACK(F_F_float_F(x,STACK_4)); /* scale and save cos(r) */ + if (end_p !=3D NULL) /* scale and save cos(r) */ + pushSTACK(F_R_float_F(x,*end_p)); + else pushSTACK(x); x =3D F_F_mal_F(STACK_1,STACK_2); /* y*s */ x =3D F_F_mal_F(x,STACK_2); /* y*s*s */ x =3D R_R_minus_R(Fixnum_1,x); /* 1-y*s*s =3D cos(s)^2 */ @@ -545,7 +557,8 @@ x =3D F_sqrt_F(x); /* cos(s)*sin(s)/s */ x =3D F_F_mal_F(x,STACK_2); /* cos(s)*sin(s) */ x =3D F_I_scale_float_F(x,Fixnum_1); /* 2*cos(s)*sin(s) =3D sin(r) */ - x =3D F_F_float_F(x,STACK_5); /* scale sin(r) */ + if (end_p !=3D NULL) /* scale sin(r) */ + x =3D F_R_float_F(x,*end_p); STACK_2 =3D STACK_0; STACK_1 =3D x; skipSTACK(1); @@ -699,7 +712,6 @@ =20 # R_ln_R(x) liefert zu einer reellen Zahl x>0 die Zahl ln(x). # can trigger GC - local object R_ln_R (object x); # Methode: # x rational -> bei x=3D1 0 als Ergebnis, sonst x in Float umwandeln. # x Float -> @@ -708,34 +720,35 @@ # (m,e) :=3D (decode-float x), so dass 1/2 <=3D m < 1. # m<2/3 -> m:=3D2m, e:=3De-1, so dass 2/3 <=3D m <=3D 4/3. # ln(m) errechnen, ln(x)=3Dln(m)+e*ln(2) als Ergebnis. - local object R_ln_R(x) - var object x; - { if (R_rationalp(x)) - { if (eq(x,Fixnum_1)) { return Fixnum_0; } # x=3D1 -> 0 als Ergebn= is - x =3D RA_float_F(x); # sonst in Float umwandeln +local object R_ln_R (object x, bool start_p, object* end_p) +{ + if (R_rationalp(x)) { + if (eq(x,Fixnum_1)) { return Fixnum_0; } /* x=3D1 -> return 0 */ + x =3D RA_float_F(x); /* convert to float */ } - # x Float - pushSTACK(x); # x retten - x =3D F_extend2_F(x); # Rechengenauigkeit erh=C3=B6hen - F_decode_float_F_I_F(x); # m,e,s bestimmen - # Stackaufbau: x, m, e, s. + /* x -- float */ + pushSTACK(x); /* save x */ + if (start_p) /* increase accuracy */ + x =3D F_extend2_F(x); + F_decode_float_F_I_F(x); /* compute m,e,s */ + /* Stack layout: x, m, e, s. */ if (F_F_comp(STACK_2, - make_SF(0,0+SF_exp_mid,floor(bit(SF_mant_len+2),3)) # S= hort-Float 2/3 - ) - < 0 - ) # m < 2/3 -> - { STACK_2 =3D F_I_scale_float_F(STACK_2,Fixnum_1); # m verdoppeln - STACK_1 =3D I_minus1_plus_I(STACK_1); # e decrementieren + make_SF(0,0+SF_exp_mid,floor(bit(SF_mant_len+2),3))) /* sho= rt-float 2/3 */ + < 0) { /* m < 2/3 -> */ + STACK_2 =3D F_I_scale_float_F(STACK_2,Fixnum_1); /* double m */ + STACK_1 =3D I_minus1_plus_I(STACK_1); /* decrement e */ } - STACK_2 =3D F_lnx_F(STACK_2); # ln(m) im genaueren Float-Format erre= chnen + STACK_2 =3D F_lnx_F(STACK_2); /* ln(m) in the more accurate float format= */ {var object temp; - temp =3D ln2_F_float_F(STACK_0); # ln(2) im genaueren Float-Format - temp =3D R_R_mal_R(STACK_1,temp); # e*ln(2) - temp =3D R_R_plus_R(STACK_2,temp); # ln(m)+e*ln(2) - temp =3D F_F_float_F(temp,STACK_3); # (float ... x) + temp =3D ln2_F_float_F(STACK_0); /* ln(2) in that float format */ + temp =3D R_R_mal_R(STACK_1,temp); /* e*ln(2) */ + temp =3D R_R_plus_R(STACK_2,temp); /* ln(m)+e*ln(2) */ + if (end_p !=3D NULL) /* (float ... x) */ + temp =3D F_R_float_F(temp,*end_p); skipSTACK(4); return temp; - }} + } +} #define F_ln_F R_ln_R =20 # I_I_log_RA(a,b) liefert zu Integers a>0, b>1 den Logarithmus log(a,b) @@ -902,9 +915,13 @@ STACK_1 =3D RA_F_float_F(a,b); # a :=3D (float a b) } } # Nun a,b beide Floats. - STACK_1 =3D R_ln_R(STACK_1); # (ln a) errechnen - {var object lnb =3D R_ln_R(popSTACK()); # (ln b) errechnen - return F_F_durch_F(popSTACK(),lnb); # (/ (ln a) (ln b)) als Ergebnis + pushSTACK(R_ln_R(STACK_1,true,NULL)); /* (ln a) */ + pushSTACK(R_ln_R(STACK_1,true,NULL)); /* (ln b) */ + STACK_0 =3D F_F_durch_F(STACK_1,STACK_0); /* (/ (ln a) (ln b)) */ + STACK_1 =3D R_R_contagion_R(STACK_2,STACK_3); + { var object ret =3D F_R_float_F(STACK_0,STACK_1); + skipSTACK(4); + return ret; }} =20 # F_expx_F(x) liefert zu einem Float x (betragsm=C3=A4=C3=C2=9Fig <1) exp(= x) als Float. @@ -966,7 +983,6 @@ =20 # R_exp_R(x) liefert zu einer reellen Zahl x die Zahl exp(x). # can trigger GC - local object R_exp_R (object x); # Methode: # x rational -> bei x=3D0 1 als Ergebnis, sonst x in Float umwandeln. # x Float -> @@ -974,32 +990,35 @@ # Genauigkeit um sqrt(d)+max(integer-length(e)) Bits erh=C3=B6hen, # (q,r) :=3D (floor x ln(2)) # Ergebnis ist exp(q*ln(2)+r) =3D (scale-float exp(r) q). - local object R_exp_R(x) - var object x; - { if (R_rationalp(x)) - # x rational - { if (eq(x,Fixnum_0)) { return Fixnum_1; } # x=3D0 -> 1 als Ergebn= is - x =3D RA_float_F(x); # sonst in Float umwandeln +local object R_exp_R (object x, bool start_p, object* end_p) +{ + if (R_rationalp(x)) { /* x rational */ + if (eq(x,Fixnum_0)) { return Fixnum_1; } /* x=3D0 -> return 1 */ + x =3D RA_float_F(x); /* =3D=3D> float */ } - # x Float - pushSTACK(x); # x retten - x =3D F_extend2_F(x); # Genauigkeit vergr=C3=B6=C3=C2=9Fern - # durch ln(2) dividieren (bei 0<=3Dx<1/2 kann man sofort q:=3D0 setz= en) - if ((!R_minusp(x)) && (F_exponent_L(x)<0)) - { pushSTACK(Fixnum_0); pushSTACK(x); } # x>=3D0, Exponent <0 -> 0<= =3Dx<1/2 -> Division unn=C3=B6tig - else - { pushSTACK(x); - {var object ln2 =3D ln2_F_float_F(x); # ln(2) mit hinreichender G= enauigkeit + /* x -- float */ + pushSTACK(x); /* save x */ + if (start_p) /* increase accuracy */ + x =3D F_extend2_F(x); + /* divide by ln(2) (if 0<=3Dx<1/2 can immediately set q:=3D0) */ + if ((!R_minusp(x)) && (F_exponent_L(x)<0)) { + /* x>=3D0, Exponent <0 -> 0<=3Dx<1/2 -> division not necessary */ + pushSTACK(Fixnum_0); pushSTACK(x); + } else { + pushSTACK(x); + { var object ln2 =3D ln2_F_float_F(x); /* ln(2) with sufficient accura= cy */ x =3D popSTACK(); - F_F_floor_I_F(x,ln2); # x durch ln(2) dividieren + F_F_floor_I_F(x,ln2); /* x / ln(2) */ }} - # Stackaufbau: originales x, q, r. - {var object temp =3D F_expx_F(STACK_0); # exp(r) - temp =3D F_I_scale_float_F(temp,STACK_1); # mal 2^q - temp =3D F_F_float_F(temp,STACK_2); # (float ... x) als Ergebnis + /* stack layout: original x, q, r. */ + { var object temp =3D F_expx_F(STACK_0); /* exp(r) */ + temp =3D F_I_scale_float_F(temp,STACK_1); /* * 2^q */ + if (end_p !=3D NULL) /* (float ... x) als Ergebnis */ + temp =3D F_R_float_F(temp,*end_p); skipSTACK(3); return temp; - }} + } +} =20 # R_sinh_R(x) liefert zu einer reellen Zahl x die Zahl sinh(x). # can trigger GC @@ -1023,7 +1042,7 @@ { var object temp; pushSTACK(x); pushSTACK(temp =3D F_extend_F(x)); # Rechengenauigkeit erh=C3=B6= hen - temp =3D F_sqrt_F(F_sinhx_F(x)); # Wurzel aus (sinh(x)/x)^2 + temp =3D F_sqrt_F(F_sinhx_F(temp)); # Wurzel aus (sinh(x)/x)^2 temp =3D F_F_mal_F(temp,STACK_0); # mit genauerem x multiplizier= en temp =3D F_F_float_F(temp,STACK_1); # und wieder runden skipSTACK(2); @@ -1032,10 +1051,12 @@ else # e>0 -> verwende exp(x) { var object temp; - pushSTACK(temp =3D R_exp_R(x)); # y:=3Dexp(x) + pushSTACK(x); + pushSTACK(temp =3D R_exp_R(x,true,NULL)); /* y:=3Dexp(x) */ temp =3D F_durch_F(temp); # (/ y) temp =3D F_F_minus_F(popSTACK(),temp); # von y subtrahieren - return F_I_scale_float_F(temp,Fixnum_minus1); # (scale-float ...= -1) + return F_F_float_F(F_I_scale_float_F(temp,Fixnum_minus1), + popSTACK()); /* (scale-float ... -1) */ } } =20 # R_cosh_R(x) liefert zu einer reellen Zahl x die Zahl cosh(x). @@ -1064,10 +1085,12 @@ if (e > 0) # e>0 -> verwende exp(x) { var object temp; - pushSTACK(temp =3D R_exp_R(x)); # y:=3Dexp(x) + pushSTACK(x); + pushSTACK(temp =3D R_exp_R(x,true,NULL)); /* y:=3Dexp(x) */ temp =3D F_durch_F(temp); # (/ y) temp =3D F_F_plus_F(popSTACK(),temp); # zu y addieren - return F_I_scale_float_F(temp,Fixnum_minus1); # (scale-float ..= . -1) + return F_F_float_F(F_I_scale_float_F(temp,Fixnum_minus1), + popSTACK()); /* (scale-float ... -1) */ } else # e<=3D0 @@ -1092,7 +1115,6 @@ =20 # R_cosh_sinh_R_R(x) liefert ((cosh x),(sinh x)), beide Werte auf dem Stac= k. # can trigger GC - local void R_cosh_sinh_R_R (object x); # Methode: # x rational -> bei x=3D0 (1,0) als Ergebnis, sonst x in Float umwandeln. # x Float -> Genauigkeit erh=C3=B6hen, @@ -1108,48 +1130,52 @@ # falls e>0: y:=3Dexp(x) errechnen, # (scale-float (+ y (/ y)) -1) und (scale-float (- y (/ y)) -1) bilden. # Genauigkeit wieder verringern. - local void R_cosh_sinh_R_R(x) - var object x; - { if (R_rationalp(x)) - # x rational - { if (eq(x,Fixnum_0)) { pushSTACK(Fixnum_1); pushSTACK(Fixnum_0); = return; } # x=3D0 -> (1,0) als Ergebnis - x =3D RA_float_F(x); # sonst in Float umwandeln +local void R_cosh_sinh_R_R (object x, bool start_p, object* end_p) +{ + if (R_rationalp(x)) { /* x rational */ + if (eq(x,Fixnum_0)) /* x=3D0 -> return (1,0) */ + { pushSTACK(Fixnum_1); pushSTACK(Fixnum_0); return; } + x =3D RA_float_F(x); /* =3D=3D> Float */ } - # x Float + /* x -- float */ {var sintL e =3D F_exponent_L(x); - if (e > 0) - # e>0 -> verwende exp(x) - { var object temp; - pushSTACK(temp =3D R_exp_R(x)); # y:=3Dexp(x) - pushSTACK(temp =3D F_durch_F(temp)); # (/ y) - # Stackaufbau: exp(x), exp(-x). - temp =3D F_F_minus_F(STACK_1,temp); # von y subtrahieren - temp =3D F_I_scale_float_F(temp,Fixnum_minus1); # (scale-float = ... -1) - {var object temp2 =3D STACK_0; - STACK_0 =3D temp; - temp =3D F_F_plus_F(STACK_1,temp2); # zu y addieren - STACK_1 =3D F_I_scale_float_F(temp,Fixnum_minus1); # (scale-flo= at ... -1) - return; - }} - else - # e<=3D0 - { if (R_zerop(x) - || (e <=3D (sintL)(1-F_float_digits(x))>>1) # e <=3D (1-d)/= 2 <=3D=3D> e <=3D -ceiling((d-1)/2) ? - ) - { pushSTACK(x); pushSTACK(x); STACK_1 =3D I_F_float_F(Fixnum_= 1,x); return; } - {var object temp; + if (e > 0) { /* e>0 -> use exp(x) */ + var object temp; pushSTACK(x); - pushSTACK(temp =3D F_extend_F(x)); # Rechengenauigkeit erh=C3= =B6hen - pushSTACK(F_square_F(temp)); # x*x - pushSTACK(temp =3D F_sinhx_F(STACK_1)); # y:=3D(sinh(x)/x)^2 - # Stackaufbau: originales x, x, x^2, y. - temp =3D F_sqrt_F(temp); # sqrt(y) =3D sinh(x)/x - temp =3D F_F_mal_F(STACK_2,temp); # x*sqrt(y) =3D sinh(x) - STACK_2 =3D F_F_float_F(temp,STACK_3); # und wieder runden - temp =3D F_F_mal_F(STACK_1,STACK_0); # x^2*y - temp =3D F_sqrt_F(R_R_plus_R(Fixnum_1,temp)); # sqrt(1+x^2*y) - STACK_3 =3D F_F_float_F(temp,STACK_3); # und wieder runden + pushSTACK(temp =3D R_exp_R(x,start_p,NULL)); /* y:=3Dexp(x) */ + pushSTACK(temp =3D F_durch_F(temp)); /* (/ y) */ + /* stack layout: x, exp(x), exp(-x). */ + temp =3D F_F_minus_F(STACK_1,temp); /* - y */ + temp =3D F_I_scale_float_F(temp,Fixnum_minus1); /* /2 */ + STACK_2 =3D F_F_float_F(temp,STACK_2); + temp =3D F_F_plus_F(STACK_1,STACK_0); /* + y */ + STACK_1 =3D F_I_scale_float_F(temp,Fixnum_minus1); /* /2 */ + STACK_1 =3D F_F_float_F(STACK_1,STACK_2); + skipSTACK(1); return; + } else { /* e<=3D0 */ + if (R_zerop(x) + || (e <=3D (sintL)(1-F_float_digits(x))>>1)) { /* e <=3D (1-d)/2= <=3D=3D> e <=3D -ceiling((d-1)/2) ? */ + pushSTACK(x); pushSTACK(x); STACK_1 =3D I_F_float_F(Fixnum_1,x); r= eturn; + } + pushSTACK(x); + { var object temp =3D (start_p ? F_extend_F(x) : x); + pushSTACK(temp); + pushSTACK(F_square_F(temp)); /* x*x */ + pushSTACK(temp =3D F_sinhx_F(STACK_1)); /* y:=3D(sinh(x)/x)^2 */ + /* stack layout: original x, x, x^2, y. */ + temp =3D F_sqrt_F(temp); /* sqrt(y) =3D sinh(x)/x */ + temp =3D F_F_mal_F(STACK_2,temp); /* x*sqrt(y) =3D sinh(x) */ + if (end_p !=3D NULL) /* restore the accuracy */ + STACK_2 =3D F_F_float_F(temp,STACK_3); + else STACK_2 =3D temp; + temp =3D F_F_mal_F(STACK_1,STACK_0); /* x^2*y */ + temp =3D F_sqrt_F(R_R_plus_R(Fixnum_1,temp)); /* sqrt(1+x^2*y) */ + if (end_p !=3D NULL) /* restore the accuracy */ + STACK_3 =3D F_F_float_F(temp,STACK_3); + else STACK_3 =3D temp; skipSTACK(2); return; - }} - } } + } + } + } +} =20 Index: src/realelem.d =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D RCS file: /cvsroot/clisp/clisp/src/realelem.d,v retrieving revision 1.15 diff -u -w -b -r1.15 realelem.d --- src/realelem.d 26 Feb 2002 14:18:14 -0000 1.15 +++ src/realelem.d 9 Sep 2002 21:41:43 -0000 @@ -219,6 +219,29 @@ var object x; { return (R_rationalp(x) ? RA_float_F(x) : x); } =20 +/* convert a float(rational) X to the appropriate format for Y + same as F_F_float_F(x,R_float_F(y)) but without consing + an intermediate float + can trigger GC */ +local object F_R_float_F (object x, object y) +{ + defaultfloatcase(S(default_float_format),y, + return F_to_SF(x), + return F_to_FF(x), + return F_to_DF(x), + return F_to_LF(x,I_to_UL(O(LF_digits))), + pushSTACK(x), x =3D popSTACK()); +} +local object RA_R_float_F (object x, object y) +{ + defaultfloatcase(S(default_float_format),y, + return RA_to_SF(x), + return RA_to_FF(x), + return RA_to_DF(x), + return RA_to_LF(x,I_to_UL(O(LF_digits))), + pushSTACK(x), x =3D popSTACK()); +} + # Generiert eine Funktion wie R_floor_I_R #define GEN_R_round(rounding) \ # Liefert ganzzahligen und gebrochenen Anteil einer reellen Zahl. \ Index: src/lisparit.d =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D RCS file: /cvsroot/clisp/clisp/src/lisparit.d,v retrieving revision 1.36 diff -u -w -b -r1.36 lisparit.d --- src/lisparit.d 22 Jul 2002 12:08:56 -0000 1.36 +++ src/lisparit.d 9 Sep 2002 21:41:43 -0000 @@ -978,12 +978,11 @@ mv_count=3D1; set_args_end_pointer(rest_args_pointer); } =20 -LISPFUNN(exp,1) # (EXP number), CLTL S. 203 - { - var object x =3D popSTACK(); - check_number(x); - value1 =3D N_exp_N(x); mv_count=3D1; +LISPFUNN(exp,1) { + check_number(STACK_0); + value1 =3D N_exp_N(STACK_0,true,&STACK_0); mv_count=3D1; + skipSTACK(1); } =20 LISPFUNN(expt,2) @@ -998,17 +997,18 @@ LISPFUN(log,1,1,norest,nokey,0,NIL) # (LOG number [base-number]), CLTL S. 204 { - var object base =3D popSTACK(); - var object arg =3D popSTACK(); + var object base =3D STACK_0; + var object arg =3D STACK_1; check_number(arg); if (eq(base,unbound)) { # LOG mit einem Argument - value1 =3D N_log_N(arg); + value1 =3D N_log_N(arg,true,&STACK_1); } else { # LOG mit zwei Argumenten check_number(base); value1 =3D N_N_log_N(arg,base); } + skipSTACK(2); mv_count=3D1; } =20 @@ -2071,7 +2071,7 @@ # > 1 wachsen lassen, damit es nicht zu h=C3=A4ufig nachberechnet = wird: oldlen +=3D floor(oldlen,2); # oldlen * 3/2 var uintC newlen =3D (d < oldlen ? oldlen : d); - ln_x =3D *objptr =3D R_ln_R(I_to_LF(x,newlen)); # (ln x) als LF mi= t newlen Digits berechnen + ln_x =3D *objptr =3D R_ln_R(I_to_LF(x,newlen),true,objptr); # (ln = x) als LF mit newlen Digits berechnen return (d < newlen ? LF_shorten_LF(ln_x,d) : ln_x); } else { # ein Double-Float reicht Index: src/constsym.d =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D RCS file: /cvsroot/clisp/clisp/src/constsym.d,v retrieving revision 1.161 diff -u -w -b -r1.161 constsym.d --- src/constsym.d 7 Aug 2002 03:28:35 -0000 1.161 +++ src/constsym.d 9 Sep 2002 21:41:43 -0000 @@ -915,6 +915,7 @@ LISPSYM(set_stream_external_format,"SET-STREAM-EXTERNAL-FORMAT",system) LISPSYM(interactive_stream_p,"INTERACTIVE-STREAM-P",lisp) LISPSYM(built_in_stream_close,"BUILT-IN-STREAM-CLOSE",system) +LISPSYM(close,"CLOSE",lisp) /* error reporting in stream.d */ LISPSYM(read_byte,"READ-BYTE",lisp) LISPSYM(read_byte_lookahead,"READ-BYTE-LOOKAHEAD",ext) LISPSYM(read_byte_will_hang_p,"READ-BYTE-WILL-HANG-P",ext) Index: src/comptran.d =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D RCS file: /cvsroot/clisp/clisp/src/comptran.d,v retrieving revision 1.11 diff -u -w -b -r1.11 comptran.d --- src/comptran.d 4 Sep 2002 19:50:27 -0000 1.11 +++ src/comptran.d 9 Sep 2002 21:41:43 -0000 @@ -21,47 +21,65 @@ =20 # N_exp_N(x) liefert (exp x), wo x eine Zahl ist. # can trigger GC - local object N_exp_N (object x); # Methode: # x reell -> klar. # x =3D a+bi -> (exp a) mit (cos b) + i (sin b) multiplizieren: # (complex (* (exp a) (cos b)) (* (exp a) (sin b))) - local object N_exp_N(x) - var object x; +local object N_exp_N (object x, bool start_p, object* end_p) { if (N_realp(x)) { - return R_exp_R(x); - } else { - # x=3Da+bi - pushSTACK(TheComplex(x)->c_real); # a retten - R_cos_sin_R_R(TheComplex(x)->c_imag); # (cos b) und (sin b) errech= nen - # Stackaufbau: a, cos(b), sin(b). - # Da b nicht =3D Fixnum 0 ist, ist auch sin(b) nicht =3D Fixnum 0. - STACK_2 =3D x =3D R_exp_R(STACK_2); # (exp a) - # Stackaufbau: exp(a), cos(b), sin(b). - STACK_0 =3D R_R_mal_R(x,STACK_0); # (* (exp a) (sin b)), nicht =3D= Fixnum 0 - x =3D R_R_mal_R(STACK_2,STACK_1); # (* (exp a) (cos b)) - x =3D R_R_complex_C(x,STACK_0); # (complex ... ...) + return R_exp_R(x,start_p,end_p); + } else { /* x=3Da+bi */ + pushSTACK(TheComplex(x)->c_real); /* save a */ + R_cos_sin_R_R(TheComplex(x)->c_imag,start_p,NULL); /* (cos b), (sin b)= */ + /* stack layout: a, cos(b), sin(b). */ + /* b !=3D Fixnum_0 =3D=3D> sin(b) ~=3D Fixnum_0. */ + STACK_2 =3D R_exp_R(STACK_2,start_p,NULL); /* (exp a) */ + /* stack layout: exp(a), cos(b), sin(b). */ + STACK_0 =3D R_R_mal_R(STACK_2,STACK_0); /* (* (exp a) (sin b)) !=3D Fi= xnum_0 */ + STACK_1 =3D R_R_mal_R(STACK_2,STACK_1); /* (* (exp a) (cos b)) */ + x =3D R_R_complex_C(F_R_float_F(STACK_1,*end_p), + F_R_float_F(STACK_0,*end_p)); /* (complex ... ...) */ skipSTACK(3); return x; } } =20 +local int lfl (object x) { + if (floatp(x)) + floatcase(x,return -1,return -2,return -3,return Lfloat_length(x)); + else return -4; +} + # N_log_N(x) liefert (log x), wo x eine Zahl ist. # can trigger GC - local object N_log_N (object x); # Methode: # (complex (log (abs x)) (phase x)) - local object N_log_N(x) - var object x; +local object N_log_N (object x, bool start_p, object *end_p) { - pushSTACK(x); # x retten - var object r =3D N_abs_R(x); # (abs x) - if (R_zerop(r)) # (abs x) =3D 0 -> Error + pushSTACK(x); /* save x */ + pushSTACK(N_abs_R(x)); /* (abs x) */ + if (R_zerop(STACK_0)) /* (abs x) =3D 0 -> Error */ divide_0(); - r =3D R_ln_R(r); # (log (abs x)) - x =3D STACK_0; STACK_0 =3D r; - x =3D N_phase_R(x); # (phase x) - return R_R_complex_N(popSTACK(),x); # (complex r (phase x)) + STACK_0 =3D R_ln_R(STACK_0,start_p,end_p); /* (log (abs x)) */ + if (start_p) { /* increase precision */ + if (floatp(STACK_1)) + STACK_1 =3D F_extend2_F(STACK_1); + else if (complexp(STACK_1)) { + if (floatp(TheComplex(STACK_1)->c_real)) + TheComplex(STACK_1)->c_real =3D + F_extend2_F(TheComplex(STACK_1)->c_real); + if (floatp(TheComplex(STACK_1)->c_imag)) + TheComplex(STACK_1)->c_imag =3D + F_extend2_F(TheComplex(STACK_1)->c_imag); + }; + } + STACK_1 =3D N_phase_R(STACK_1); /* (phase x) */ + if (end_p !=3D NULL && floatp(STACK_1)) + STACK_1 =3D F_R_float_F(STACK_1,*end_p); + { /* (complex (log (abs x)) (phase x)) */ + var object ret =3D R_R_complex_N(STACK_0,STACK_1); + skipSTACK(2); return ret; + } } =20 # N_N_log_N(a,b) liefert (log a b), wo a und b Zahlen sind. @@ -109,7 +127,7 @@ b =3D STACK_1; if (R_rationalp(b)) b =3D RA_F_float_F(b,angle); - b =3D F_ln_F(b); STACK_0 =3D F_F_durch_F(STACK_0,b); + b =3D F_ln_F(b,true,&STACK_1); STACK_0 =3D F_F_durch_F(STACK_0= ,b); } # Stackaufbau: a, b, Imagin=C3=A4rteil. # Realteil (/ (log (abs a)) (log b)) errechnen: @@ -139,22 +157,23 @@ } } # Keine Chance f=C3=BCr rationalen Realteil - pushSTACK(F_ln_F(N_abs_R(a))); # (log (abs a)), ein Float + pushSTACK(F_ln_F(N_abs_R(a),true,&STACK_3)); /* (log (abs a)), a= float */ # durch (log b) dividieren, liefert den Realteil: b =3D STACK_2; if (R_rationalp(b)) b =3D RA_F_float_F(b,STACK_0); - b =3D F_ln_F(b); STACK_0 =3D F_F_durch_F(STACK_0,b); + b =3D F_ln_F(b,true,&STACK_2); STACK_0 =3D F_F_durch_F(STACK_0,b= ); real_ok: # Stackaufbau: a, b, Imagin=C3=A4rteil, Realteil. var object erg =3D R_R_complex_C(STACK_0,STACK_1); skipSTACK(4); return erg; } - } else { - # normaler komplexer Fall - pushSTACK(b); a =3D N_log_N(a); b =3D STACK_0; # (log a) - STACK_0 =3D a; b =3D N_log_N(b); a =3D popSTACK(); # (log b) - return N_N_durch_N(a,b); # dividieren + } else { /* normal complex case */ + pushSTACK(a); pushSTACK(b); + STACK_1 =3D N_log_N(STACK_1,true,&STACK_1); /* (log a) */ + STACK_0 =3D N_log_N(STACK_0,true,&STACK_0); /* (log b) */ + a =3D N_N_durch_N(a,b); /* divide */ + skipSTACK(2); return a; } } =20 @@ -365,57 +384,63 @@ } } pushSTACK(y); - var object temp =3D N_log_N(x); # (log x) - temp =3D N_N_mal_N(temp,popSTACK()); # mal y - return N_exp_N(temp); # exp davon + { var object temp =3D N_log_N(x,true,NULL); /* (log x) */ + temp =3D N_N_mal_N(temp,F_extend2_F(STACK_0)); /* * y */ + temp =3D N_exp_N(temp,false,&STACK_0); /* exp */ + skipSTACK(1); return temp; + } } =20 # N_sin_N(x) liefert (sin x), wo x eine Zahl ist. # can trigger GC - local object N_sin_N (object x); # Methode: # x reell -> klar # x =3D a+bi -> (complex (* (sin a) (cosh b)) (* (cos a) (sinh b))) - local object N_sin_N(x) - var object x; +local object N_sin_N (object x) { if (N_realp(x)) { return R_sin_R(x); - } else { - # x=3Da+bi - pushSTACK(TheComplex(x)->c_real); # a retten - R_cosh_sinh_R_R(TheComplex(x)->c_imag); # cosh(b), sinh(b) errechn= en - # Stackaufbau: a, cosh(b), sinh(b). - # Da b nicht =3D Fixnum 0 ist, ist auch sinh(b) nicht =3D Fixnum 0. - R_cos_sin_R_R(STACK_2); # cos(a), sin(a) errechnen, cos(a) /=3D Fi= xnum 0 - # Stackaufbau: a, cosh(b), sinh(b), cos(a), sin(a). - STACK_0 =3D R_R_mal_R(STACK_0,STACK_3); # sin(a)*cosh(b) - var object erg =3D R_R_mal_R(STACK_1,STACK_2); # cos(a)*sinh(b), n= icht Fixnum 0 - erg =3D R_R_complex_C(STACK_0,erg); skipSTACK(5); return erg; + } else { /* x=3Da+bi */ + pushSTACK(TheComplex(x)->c_real); /* save a */ + pushSTACK(TheComplex(x)->c_imag); /* save b */ + R_cosh_sinh_R_R(STACK_0,true,NULL); /* cosh(b) sinh(b) */ + /* stack layout: a, b, cosh(b), sinh(b). */ + /* b !=3D Fixnum_0 =3D=3D> sinh(b) !=3D Fixnum_0. */ + R_cos_sin_R_R(STACK_3,true,NULL); /* cos(a)!=3D0, sin(a); */ + /* stack layout: a, b, cosh(b), sinh(b), cos(a), sin(a). */ + STACK_0 =3D R_R_mal_R(STACK_0,STACK_3); /* sin(a)*cosh(b) */ + STACK_2 =3D R_R_contagion_R(STACK_4,STACK_5); + { var object erg =3D R_R_mal_R(STACK_1,STACK_2); /* cos(a)*sinh(b), != =3D Fixnum_0 */ + erg =3D R_R_complex_C(F_R_float_F(STACK_0,STACK_2), + F_R_float_F(erg,STACK_2)); + skipSTACK(6); return erg; + } } } =20 # N_cos_N(x) liefert (cos x), wo x eine Zahl ist. # can trigger GC - local object N_cos_N (object x); # Methode: # x reell -> klar # x =3D a+bi -> (complex (* (cos a) (cosh b)) (- (* (sin a) (sinh b)))) - local object N_cos_N(x) - var object x; +local object N_cos_N (object x) { if (N_realp(x)) { return R_cos_R(x); - } else { - # x=3Da+bi - pushSTACK(TheComplex(x)->c_real); # a retten - R_cosh_sinh_R_R(TheComplex(x)->c_imag); # cosh(b), sinh(b) errechn= en - # Stackaufbau: a, cosh(b), sinh(b). - R_cos_sin_R_R(STACK_2); # cos(a), sin(a) errechnen - # Stackaufbau: a, cosh(b), sinh(b), cos(a), sin(a). - STACK_0 =3D R_minus_R(R_R_mal_R(STACK_0,STACK_2)); # -sin(a)*sinh(= b) - var object erg =3D R_R_mal_R(STACK_1,STACK_3); # cos(a)*cosh(b) - erg =3D R_R_complex_N(erg,STACK_0); skipSTACK(5); return erg; + } else { /* x=3Da+bi */ + pushSTACK(TheComplex(x)->c_real); /* save a */ + pushSTACK(TheComplex(x)->c_imag); /* save b */ + R_cosh_sinh_R_R(STACK_0,true,NULL); /* cosh(b), sinh(b) */ + /* stack layout: a, b, cosh(b), sinh(b). */ + R_cos_sin_R_R(STACK_3,true,NULL); /* cos(a), sin(a) */ + /* stack layout: a, b, cosh(b), sinh(b), cos(a), sin(a). */ + STACK_0 =3D R_minus_R(R_R_mal_R(STACK_0,STACK_2)); /* -sin(a)*sinh(b) = */ + STACK_1 =3D R_R_mal_R(STACK_1,STACK_3); /* cos(a)*cosh(b) */ + STACK_2 =3D R_R_contagion_R(STACK_4,STACK_5); + { var object erg =3D R_R_complex_N(F_R_float_F(STACK_1,STACK_2), + F_R_float_F(STACK_0,STACK_2)); + skipSTACK(6); return erg; + } } } =20 @@ -430,22 +455,24 @@ var object x; { if (N_realp(x)) { - R_cos_sin_R_R(x); - # Stackaufbau: cos(x), sin(x). - var object erg =3D R_R_durch_R(STACK_0,STACK_1); - skipSTACK(2); return erg; - } else { - # x=3Da+bi - pushSTACK(TheComplex(x)->c_real); # a retten - R_cosh_sinh_R_R(TheComplex(x)->c_imag); # cosh(b), sinh(b) errechn= en - # Stackaufbau: a, cosh(b), sinh(b). - R_cos_sin_R_R(STACK_2); # cos(a), sin(a) errechnen - # Stackaufbau: a, cosh(b), sinh(b), cos(a), sin(a). + pushSTACK(x); + R_cos_sin_R_R(x,true,NULL); + { /* stack layout: x, cos(x), sin(x). */ + var object erg =3D F_R_float_F(R_R_durch_R(STACK_0,STACK_1),STAC= K_2); + skipSTACK(3); return erg; + } + } else { /* x=3Da+bi */ + pushSTACK(TheComplex(x)->c_real); /* save a */ + pushSTACK(TheComplex(x)->c_imag); /* save b */ + R_cosh_sinh_R_R(STACK_0,true,NULL); /* cosh(b), sinh(b) */ + /* stack layout: a, cosh(b), sinh(b). */ + R_cos_sin_R_R(STACK_2,true,NULL); /* cos(a), sin(a) */ + # stack layout: a, cosh(b), sinh(b), cos(a), sin(a). var object temp; STACK_4 =3D R_R_mal_R(STACK_0,STACK_3); # sin(a)*cosh(b) temp =3D R_R_mal_R(STACK_1,STACK_2); # cos(a)*sinh(b) /=3D Fixnum 0 STACK_4 =3D R_R_complex_C(STACK_4,temp); # Z=C3=A4hler - # Stackaufbau: Z=C3=A4hler, cosh(b), sinh(b), cos(a), sin(a). + # stack layout: Z=C3=A4hler, cosh(b), sinh(b), cos(a), sin(a). STACK_3 =3D R_R_mal_R(STACK_1,STACK_3); # cos(a)*cosh(b) temp =3D R_minus_R(R_R_mal_R(STACK_0,STACK_2)); # -sin(a)*sinh(b) temp =3D R_R_complex_N(STACK_3,temp); # Nenner @@ -456,113 +483,120 @@ =20 # N_cis_N(x) liefert (cis x), wo x eine Zahl ist. # can trigger GC - local object N_cis_N (object x); # Methode: # x reell -> (complex (cos x) (sin x)) # x =3D a+bi -> (complex (* (exp (- b)) (cos a)) (* (exp (- b)) (sin a))) - local object N_cis_N(x) - var object x; +local object N_cis_N (object x) { if (N_realp(x)) { - R_cos_sin_R_R(x); - # Stackaufbau: cos(x), sin(x). + pushSTACK(x); + R_cos_sin_R_R(x,true,&STACK_0); + /* stack layout: x, cos(x), sin(x). */ var object erg =3D R_R_complex_N(STACK_1,STACK_0); - skipSTACK(2); return erg; - } else { - # x=3Da+bi - pushSTACK(TheComplex(x)->c_imag); # b retten - R_cos_sin_R_R(TheComplex(x)->c_real); # (cos a) und (sin a) errech= nen - # Stackaufbau: b, cos(a), sin(a). - STACK_2 =3D x =3D R_exp_R(R_minus_R(STACK_2)); # (exp (- b)) - # Stackaufbau: exp(-b), cos(a), sin(a). - STACK_0 =3D R_R_mal_R(x,STACK_0); # (* (exp (- b)) (sin a)) - x =3D R_R_mal_R(STACK_2,STACK_1); # (* (exp (- b)) (cos a)) - x =3D R_R_complex_N(x,STACK_0); # (complex ... ...) - skipSTACK(3); return x; + skipSTACK(3); return erg; + } else { /* x=3Da+bi */ + pushSTACK(TheComplex(x)->c_real); /* save a */ + pushSTACK(TheComplex(x)->c_imag); /* save b */ + R_cos_sin_R_R(TheComplex(x)->c_real,true,NULL); /* (cos a), (sin a) */ + /* stack layout: a, b, cos(a), sin(a). */ + STACK_3 =3D R_R_contagion_R(STACK_3,STACK_2); + STACK_2 =3D x =3D R_exp_R(R_minus_R(STACK_2),false,NULL); /* (exp (- b= )) */ + /* stack layout: exp(-b), cos(a), sin(a). */ + STACK_0 =3D R_R_mal_R(x,STACK_0); /* (* (exp (- b)) (sin a)) */ + x =3D R_R_mal_R(STACK_2,STACK_1); /* (* (exp (- b)) (cos a)) */ + x =3D R_R_complex_N(F_R_float_F(x,STACK_3), /* (complex ... ...) */ + F_R_float_F(STACK_0,STACK_3)); + skipSTACK(4); return x; } } =20 # N_sinh_N(x) liefert (sinh x), wo x eine Zahl ist. # can trigger GC - local object N_sinh_N (object x); # Methode: # x reell -> klar # x =3D a+bi -> (complex (* (sinh a) (cos b)) (* (cosh a) (sin b))) - local object N_sinh_N(x) - var object x; +local object N_sinh_N (object x) { if (N_realp(x)) { return R_sinh_R(x); - } else { - # x=3Da+bi - pushSTACK(TheComplex(x)->c_real); # a retten - R_cos_sin_R_R(TheComplex(x)->c_imag); # cos(b), sin(b) errechnen - # Stackaufbau: a, cos(b), sin(b). - # Da b nicht =3D Fixnum 0 ist, ist auch sin(b) nicht =3D Fixnum 0. - R_cosh_sinh_R_R(STACK_2); # cosh(a), sinh(a) errechnen, cosh(a) /= =3D Fixnum 0 - # Stackaufbau: a, cos(b), sin(b), cosh(a), sinh(a). - STACK_0 =3D R_R_mal_R(STACK_0,STACK_3); # sinh(a)*cos(b) - var object erg =3D R_R_mal_R(STACK_1,STACK_2); # cosh(a)*sin(b), n= icht Fixnum 0 - erg =3D R_R_complex_C(STACK_0,erg); skipSTACK(5); return erg; + } else { /* x=3Da+bi */ + pushSTACK(TheComplex(x)->c_real); /* save a */ + pushSTACK(TheComplex(x)->c_imag); /* save b */ + R_cos_sin_R_R(TheComplex(x)->c_imag,true,NULL); /* cos(b), sin(b) */ + /* stack layout: a, b, cos(b), sin(b). */ + /* b !=3D Fixnum_0 =3D=3D> sin(b) !=3D Fixnum_0. */ + R_cosh_sinh_R_R(STACK_3,true,NULL); /* cosh(a), sinh(a); cosh(a) !=3D = Fixnum 0 */ + /* stack layout: a, b, cos(b), sin(b), cosh(a), sinh(a). */ + STACK_0 =3D R_R_mal_R(STACK_0,STACK_3); /* sinh(a)*cos(b) */ + { var object erg =3D R_R_mal_R(STACK_1,STACK_2); /* cosh(a)*sin(b), != =3D Fixnum_0 */ + STACK_5 =3D R_R_contagion_R(STACK_4,STACK_5); + erg =3D R_R_complex_C(F_R_float_F(STACK_0,STACK_5), + F_R_float_F(erg,STACK_5)); + skipSTACK(6); return erg; + } } } =20 # N_cosh_N(x) liefert (cosh x), wo x eine Zahl ist. # can trigger GC - local object N_cosh_N (object x); # Methode: # x reell -> klar # x =3D a+bi -> (complex (* (cosh a) (cos b)) (* (sinh a) (sin b))) - local object N_cosh_N(x) - var object x; +local object N_cosh_N (object x) { if (N_realp(x)) { return R_cosh_R(x); - } else { - # x=3Da+bi - pushSTACK(TheComplex(x)->c_real); # a retten - R_cos_sin_R_R(TheComplex(x)->c_imag); # cos(b), sin(b) errechnen - # Stackaufbau: a, cos(b), sin(b). - R_cosh_sinh_R_R(STACK_2); # cosh(a), sinh(a) errechnen - # Stackaufbau: a, cos(b), sin(b), cosh(a), sinh(a). - STACK_0 =3D R_R_mal_R(STACK_0,STACK_2); # sinh(a)*sin(b) - var object erg =3D R_R_mal_R(STACK_1,STACK_3); # cosh(a)*cos(b) - erg =3D R_R_complex_N(erg,STACK_0); skipSTACK(5); return erg; + } else { /* x=3Da+bi */ + pushSTACK(TheComplex(x)->c_real); /* save a */ + pushSTACK(TheComplex(x)->c_imag); /* save b */ + R_cos_sin_R_R(TheComplex(x)->c_imag,true,NULL); /* cos(b), sin(b) */ + /* stack layout: a, b, cos(b), sin(b). */ + R_cosh_sinh_R_R(STACK_3,true,NULL); /* cosh(a), sinh(a) */ + /* stack layout: a, b,cos(b), sin(b), cosh(a), sinh(a). */ + STACK_0 =3D R_R_mal_R(STACK_0,STACK_2); /* sinh(a)*sin(b) */ + { var object erg =3D R_R_mal_R(STACK_1,STACK_3); /* cosh(a)*cos(b) */ + STACK_5 =3D R_R_contagion_R(STACK_4,STACK_5); + erg =3D R_R_complex_N(F_R_float_F(erg,STACK_4), + F_R_float_F(STACK_0,STACK_4)); + skipSTACK(6); return erg; + } } } =20 # N_tanh_N(x) liefert (tanh x), wo x eine Zahl ist. # can trigger GC - local object N_tanh_N (object x); # Methode: # x reell -> (/ (sinh x) (cosh x)) # x =3D a+bi -> (/ (complex (* (sinh a) (cos b)) (* (cosh a) (sin b))) # (complex (* (cosh a) (cos b)) (* (sinh a) (sin b))) ) - local object N_tanh_N(x) - var object x; +local object N_tanh_N (object x) { if (N_realp(x)) { - R_cosh_sinh_R_R(x); - # Stackaufbau: cosh(x), sinh(x). - var object erg =3D R_R_durch_R(STACK_0,STACK_1); - skipSTACK(2); return erg; - } else { - # x=3Da+bi - pushSTACK(TheComplex(x)->c_real); # a retten - R_cos_sin_R_R(TheComplex(x)->c_imag); # cos(b), sin(b) errechnen - # Stackaufbau: a, cos(b), sin(b). - R_cosh_sinh_R_R(STACK_2); # cosh(a), sinh(a) errechnen - # Stackaufbau: a, cos(b), sin(b), cosh(a), sinh(a). - var object temp; - STACK_4 =3D R_R_mal_R(STACK_0,STACK_3); # sinh(a)*cos(b) - temp =3D R_R_mal_R(STACK_1,STACK_2); # cosh(a)*sin(b) /=3D Fixnum 0 - STACK_4 =3D R_R_complex_C(STACK_4,temp); # Z=C3=A4hler - # Stackaufbau: Z=C3=A4hler, cos(b), sin(b), cosh(a), sinh(a). - STACK_3 =3D R_R_mal_R(STACK_1,STACK_3); # cosh(a)*cos(b) - temp =3D R_R_mal_R(STACK_0,STACK_2); # sinh(a)*sin(b) - temp =3D R_R_complex_N(STACK_3,temp); # Nenner - temp =3D N_N_durch_N(STACK_4,temp); # Z=C3=A4hler/Nenner - skipSTACK(5); return temp; + pushSTACK(x); + R_cosh_sinh_R_R(x,true,NULL); + /* stack layout: x, cosh(x), sinh(x). */ + { var object erg =3D F_R_float_F(R_R_durch_R(STACK_0,STACK_1),STACK_2); + skipSTACK(3); return erg; + } + } else { /* x=3Da+bi */ + pushSTACK(TheComplex(x)->c_real); /* a */ + pushSTACK(TheComplex(x)->c_imag); /* b */ + R_cos_sin_R_R(TheComplex(x)->c_imag,true,NULL); /* cos(b), sin(b) */ + /* stack layout: a, b, cos(b), sin(b). */ + R_cosh_sinh_R_R(STACK_3,true,NULL); /* cosh(a), sinh(a) */ + /* stack layout: a, b, cos(b), sin(b), cosh(a), sinh(a). */ + STACK_5 =3D R_R_contagion_R(STACK_4,STACK_5); + STACK_4 =3D R_R_mal_R(STACK_0,STACK_3); /* sinh(a)*cos(b) */ + { var object temp =3D R_R_mal_R(STACK_1,STACK_2); /* cosh(a)*sin(b) /= =3D Fixnum 0 */ + STACK_4 =3D R_R_complex_C(STACK_4,temp); /* numerator */ + /* stack layout: a, numerator, cos(b), sin(b), cosh(a), sinh(a). */ + STACK_3 =3D R_R_mal_R(STACK_1,STACK_3); /* cosh(a)*cos(b) */ + temp =3D R_R_mal_R(STACK_0,STACK_2); /* sinh(a)*sin(b) */ + temp =3D R_R_complex_N(STACK_3,temp); /* denominator */ + temp =3D N_N_durch_N(STACK_4,temp); /* numerator/denominator */ + temp =3D F_R_float_F(temp,STACK_5); + skipSTACK(6); return temp; + } } } =20 @@ -608,53 +642,49 @@ # rein imagin=C3=A4r ist. =20 # Hilfsfunktion f=C3=BCr beide: u+iv :=3D artanh(x+iy), u,v beide auf den = Stack. - local void R_R_atanh_R_R (object x, object y); - local void R_R_atanh_R_R(x,y) - var object x; - var object y; +local void R_R_atanh_R_R (object x, object y, bool start_p, object *end_p) { - if (eq(x,Fixnum_0)) { # x=3D0 -> u=3D0, v=3Datan(X=3D1,Y=3Dy) (Fall = y=3D0 ist inbegriffen) + if (eq(x,Fixnum_0)) { /* x=3D0 -> u=3D0, v=3Datan(X=3D1,Y=3Dy) (y=3D0 is= included) */ pushSTACK(x); pushSTACK(R_R_atan_R(Fixnum_1,y)); return; } if (eq(y,Fixnum_0)) { if (R_rationalp(x)) - x =3D RA_float_F(x); # x in Float umwandeln - # x Float - if (R_zerop(x)) { # x=3D0.0 -> x als Ergebnis + x =3D RA_float_F(x); /* x --> float */ + /* x -- float */ + if (R_zerop(x)) { /* x=3D0.0 -> return x */ pushSTACK(x); pushSTACK(Fixnum_0); return; } if (F_exponent_L(x) < 0) { - # Exponent e<0, also |x|<1/2 + /* exponent e<0, =3D=3D> |x|<1/2 */ pushSTACK(F_atanhx_F(x)); pushSTACK(Fixnum_0); return; } - # e>=3D0, also |x|>=3D1/2 + /* e>=3D0, =3D=3D> |x|>=3D1/2 */ pushSTACK(x); - pushSTACK(R_R_minus_R(Fixnum_1,x)); # 1-x - # Stackaufbau: x, 1-x. + pushSTACK(R_R_minus_R(Fixnum_1,x)); /* 1-x */ + /* stack layout: x, 1-x. */ var object temp; - temp =3D R_R_plus_R(Fixnum_1,STACK_1); # 1+x - temp =3D F_F_durch_F(temp,STACK_0); # (1+x)/(1-x) + temp =3D R_R_plus_R(Fixnum_1,STACK_1); /* 1+x */ + temp =3D F_F_durch_F(temp,STACK_0); /* (1+x)/(1-x) */ if (!R_minusp(temp)) { - STACK_1 =3D temp; STACK_0 =3D Fixnum_0; # Imagin=C3=A4rteil:=3D0 - if (R_zerop(temp)) # x =3D -1 -> Error + STACK_1 =3D temp; STACK_0 =3D Fixnum_0; /* imag part :=3D0 */ + if (R_zerop(temp)) /* x =3D -1 -> Error */ divide_0(); - } else { - # (1+x)/(1-x) < 0 -> Betrag nehmen, Imagin=C3=A4rteil berechnen: + } else { /* (1+x)/(1-x) < 0 -> negate, compute Im: */ STACK_1 =3D F_minus_F(temp); - temp =3D F_I_scale_float_F(pi(STACK_1),Fixnum_minus1); # (scale-= float pi -1) =3D pi/2 + temp =3D F_I_scale_float_F(pi(STACK_1),Fixnum_minus1); /* (scale-flo= at pi -1) =3D pi/2 */ if (R_minusp(STACK_0)) # 1-x<0 -> dann -pi/2 temp =3D F_minus_F(temp); STACK_0 =3D temp; } - # Stackaufbau: |(1+x)/(1-x)| (>0), Imagin=C3=A4rteil. - STACK_1 =3D F_I_scale_float_F(R_ln_R(STACK_1),Fixnum_minus1); # ln= bilden, durch 2 + /* stack layout: |(1+x)/(1-x)| (>0), Im. */ + STACK_1 =3D F_I_scale_float_F(R_ln_R(STACK_1,true,&STACK_1),Fixnum_min= us1); /* ln / 2 */ return; } pushSTACK(x); pushSTACK(y); - pushSTACK(R_R_plus_R(Fixnum_1,x)); # 1+x - pushSTACK(R_R_minus_R(Fixnum_1,STACK_2)); # 1-x - # Stackaufbau: x, y, 1+x, 1-x. - # x und y in Floats umwandeln: + pushSTACK(R_R_plus_R(Fixnum_1,x)); /* 1+x */ + pushSTACK(R_R_minus_R(Fixnum_1,STACK_2)); /* 1-x */ + /* stack layout: x, y, 1+x, 1-x. */ + /* x , y --> float: */ if (R_rationalp(STACK_3)) { if (R_rationalp(STACK_2)) STACK_2 =3D RA_float_F(STACK_2); @@ -663,13 +693,13 @@ if (R_rationalp(STACK_2)) STACK_2 =3D RA_F_float_F(STACK_2,STACK_3); } - pushSTACK(R_square_R(STACK_2)); # y^2 + pushSTACK(R_square_R(STACK_2)); /* y^2 */ { - var object temp =3D R_square_R(STACK_4); # x^2 - pushSTACK(R_R_plus_R(Fixnum_1,R_R_plus_R(temp,STACK_0))); # 1+x^2+= y^2 + var object temp =3D R_square_R(STACK_4); /* x^2 */ + pushSTACK(R_R_plus_R(Fixnum_1,R_R_plus_R(temp,STACK_0))); /* 1+x^2+y^2= */ } - # Stackaufbau: x, y, 1+x, 1-x, y^2, 1+x^2+y^2. - var object temp =3D # |4x| + /* stack layout: x, y, 1+x, 1-x, y^2, 1+x^2+y^2. */ + var object temp =3D /* |4x| */ F_abs_F(F_I_scale_float_F(STACK_5,fixnum(2))); if (F_F_comp(temp,STACK_0) < 0) { # |4x| < 1+x^2+y^2 ? temp =3D F_I_scale_float_F(STACK_5,Fixnum_1); # 2x @@ -684,12 +714,12 @@ temp =3D F_F_durch_F(STACK_0,temp); # ((1+x)^2+y^2)/((1-x)^2+y^2),= ein Float >=3D0 if (R_zerop(temp)) # sollte >0 sein divide_0(); - temp =3D R_ln_R(temp); # ln davon, ein Float + temp =3D R_ln_R(temp,true,NULL); # ln davon, ein Float temp =3D F_I_scale_float_F(temp,sfixnum(-2)); # .../4 =3D: u } - var signean x_sign =3D R_sign(STACK_5); + { var signean x_sign =3D R_sign(STACK_5); STACK_5 =3D temp; - # Stackaufbau: u, y, 1+x, 1-x, y^2, -. + /* stack layout: u, y, 1+x, 1-x, y^2, -. */ temp =3D R_R_mal_R(STACK_3,STACK_2); # (1+x)(1-x) STACK_0 =3D R_R_minus_R(temp,STACK_1); # (1+x)(1-x)-y^2, ein Float temp =3D F_I_scale_float_F(STACK_4,Fixnum_1); # 2y, ein Float @@ -697,35 +727,40 @@ if (R_minusp(STACK_0) && (x_sign>=3D0) && R_zerop(STACK_4)) # X<0.0 = und x>=3D0.0 und Y=3D0.0 ? temp =3D F_minus_F(temp); # ja -> Vorzeichenwechsel STACK_4 =3D F_I_scale_float_F(temp,Fixnum_minus1); # .../2 =3D: v - # Stackaufbau: u, v, 1+x, 1-x, y^2, -. + STACK_4 =3D F_F_float_F(STACK_4,STACK_3); /* restore the precision */ + STACK_5 =3D F_F_float_F(STACK_5,STACK_3); /* restore the precision */ + /* stack layout: u, v, 1+x, 1-x, y^2, -. */ skipSTACK(4); return; } - # - local object N_atanh_N(z) - var object z; +} + +local object N_atanh_N (object z) { if (N_realp(z)) { - R_R_atanh_R_R(z,Fixnum_0); + pushSTACK(z); + R_R_atanh_R_R(z,Fixnum_0,true,&STACK_0); } else { - R_R_atanh_R_R(TheComplex(z)->c_real,TheComplex(z)->c_imag); + pushSTACK(TheComplex(z)->c_real); + R_R_atanh_R_R(TheComplex(z)->c_real,TheComplex(z)->c_imag,true,&STACK_= 0); } - # Stackaufbau: u, v. - z =3D R_R_complex_N(STACK_1,STACK_0); skipSTACK(2); return z; + /* stack layout: z, u, v. */ + z =3D R_R_complex_N(STACK_1,STACK_0); + skipSTACK(3); return z; } - # - local object N_atan_N(z) - var object z; - { # atanh(iz) errechnen: + +local object N_atan_N (object z) +{ /* atanh(iz) errechnen: */ if (N_realp(z)) { - R_R_atanh_R_R(Fixnum_0,z); + pushSTACK(z); + R_R_atanh_R_R(Fixnum_0,z,true,&STACK_0); } else { pushSTACK(TheComplex(z)->c_real); z =3D R_minus_R(TheComplex(z)->c_imag); - R_R_atanh_R_R(z,popSTACK()); + R_R_atanh_R_R(z,STACK_0,true,&STACK_0); } - # Stackaufbau: u, v. + /* stack layout: z, u, v. */ z =3D R_minus_R(STACK_1); z =3D R_R_complex_N(STACK_0,z); # z :=3D v= -iu - skipSTACK(2); return z; + skipSTACK(3); return z; } =20 # Um f=C3=BCr zwei Zahlen u,v mit u^2-v^2=3D1 und u,v beide in Bild(sqrt) @@ -829,7 +864,7 @@ || (F_exponent_L(y) <=3D (sintL)(-F_float_digits(y))>>1) # e <= =3D -d/2 <=3D=3D> e <=3D -ceiling(d/2) ? ) return; # u=3D0, v=3Dy bereits im Stack - # Stackaufbau: 0, y. + # stack layout: 0, y. var object temp =3D R_R_minus_R(Fixnum_1,F_square_F(y)); # 1-y*y if (!R_minusp(temp)) { # 1-y*y>=3D0, also |y|<=3D1 @@ -844,7 +879,7 @@ else temp =3D F_F_plus_F(temp,y); # temp =3D sqrt(y^2-1)+|y|, ein Float >1 - STACK_1 =3D R_ln_R(temp); # ln(|y|+sqrt(y^2-1)), ein Float >0 + STACK_1 =3D R_ln_R(temp,true,&STACK_0); /* ln(|y|+sqrt(y^2-1)), = Float >0 */ temp =3D F_I_scale_float_F(pi(STACK_1),Fixnum_minus1); # (scale-= float pi -1) =3D pi/2 if (!R_minusp(STACK_0)) { # Vorzeichen von y # y>1 -> v =3D pi/2 @@ -873,10 +908,10 @@ # |x|>=3D1/2 if (!R_minusp(x)) # x>=3D1 - STACK_1 =3D R_ln_R(F_F_plus_F(temp,x)); # u =3D ln(x+sqrt(1+x^= 2)) + STACK_1 =3D R_ln_R(F_F_plus_F(temp,x),true,&STACK_0); /* u =3D= ln(x+sqrt(1+x^2)) */ else # x<=3D-1 - STACK_1 =3D F_minus_F(R_ln_R(F_F_minus_F(temp,x))); # u =3D -l= n(-x+sqrt(1+x^2)) + STACK_1 =3D F_minus_F(R_ln_R(F_F_minus_F(temp,x),true,&STACK_0= )); # u =3D -ln(-x+sqrt(1+x^2)) } return; } @@ -893,7 +928,7 @@ # ist also ein Float, und da Real- und Imagin=C3=A4rteil von z /=3D0= sind, # sind Real- und Imagin=C3=A4rteil von w Floats.) # Daher hat dann atanh(...) Floats als Realteil u und Imagin=C3=A4rt= eil v. - R_R_atanh_R_R(TheComplex(z)->c_real,TheComplex(z)->c_imag); # atanh = nehmen + R_R_atanh_R_R(TheComplex(z)->c_real,TheComplex(z)->c_imag,true,&STAC= K_0); # atanh nehmen # u und v mit 2 multiplizieren: STACK_1 =3D F_I_scale_float_F(STACK_1,Fixnum_1); # u:=3D2*u STACK_0 =3D F_I_scale_float_F(STACK_0,Fixnum_1); # v:=3D2*v @@ -908,7 +943,7 @@ } else { R_R_asinh_R_R(TheComplex(z)->c_real,TheComplex(z)->c_imag); } - # Stackaufbau: u, v. + # stack layout: u, v. z =3D R_R_complex_N(STACK_1,STACK_0); skipSTACK(2); return z; } # @@ -923,7 +958,7 @@ z =3D R_minus_R(TheComplex(z)->c_imag); R_R_asinh_R_R(z,popSTACK()); } - # Stackaufbau: u, v. + # stack layout: u, v. z =3D R_minus_R(STACK_1); z =3D R_R_complex_N(STACK_0,z); # z :=3D v= -iu skipSTACK(2); return z; } @@ -978,8 +1013,9 @@ var object temp =3D STACK_0; # z temp =3D R_R_minus_R(F_square_F(temp),Fixnum_1); # z^2-1, ein Fl= oat >=3D0 temp =3D F_sqrt_F(temp); # sqrt(z^2-1), ein Float >=3D0 - temp =3D F_F_plus_F(popSTACK(),temp); # z+sqrt(z^2-1), ein Float= >1 - temp =3D R_ln_R(temp); # ln(z+sqrt(z^2-1)), ein Float >=3D0 + temp =3D F_F_plus_F(STACK_0,temp); /* z+sqrt(z^2-1), float >1 */ + temp =3D R_ln_R(temp,true,&STACK_0); /* ln(z+sqrt(z^2-1)), float= >=3D0 */ + skipSTACK(1); return R_R_complex_C(Fixnum_0,temp); } R_R_asinh_R_R(Fixnum_0,popSTACK()); @@ -988,7 +1024,7 @@ z =3D R_minus_R(TheComplex(z)->c_imag); R_R_asinh_R_R(z,popSTACK()); } - # Stackaufbau: u, v. + /* stack layout: u, v. */ # Bilde pi/2-v : z =3D STACK_0; z =3D (R_rationalp(z) ? pi(z) : pi_F_float_F(z)); # pi im Float-Form= at von v @@ -1048,7 +1084,7 @@ STACK_0 =3D z =3D RA_float_F(z); # z Float <=3D -1 z =3D F_sqrt_F(R_R_minus_R(F_square_F(z),Fixnum_1)); # sqrt(z^2-= 1), ein Float >=3D0 - STACK_0 =3D R_ln_R(F_F_minus_F(z,STACK_0)); # log(sqrt(z^2-1)-z)= , ein Float >=3D0 + STACK_0 =3D R_ln_R(F_F_minus_F(z,STACK_0),true,&STACK_0); # log(= sqrt(z^2-1)-z), ein Float >=3D0 z =3D pi(STACK_0); # and imaginary part =3D=3D pi return R_R_complex_C(popSTACK(),z); } ```
 Re: Numeric inconsistency From: Raymond Toy - 2002-09-10 14:40:17 ```>>>>> "Sam" == Sam Steingold writes: Sam> the appended patch fixes the expt precision problem which you Sam> reported in another e-mail. I would greatly appreciate it if Sam> you could test it as much as possible. I'll try this when I get a chance. However, let me note that I no longer think it is a problem. I accept Bruno's explanation that clisp returns a complex number for the test case I gave. I might wish it were a bit more accurate, but I don't consider the result given by Clisp to be wrong. It is, after all, a reasonable result. Sam> as far as TAN/COS are concerned, I am at a loss. Sam> I guess using R_cos_sin_R_R to compute just cos would be a solution... This is acceptable to me because I don't run Clisp for its number-crunching speed. :-) I run Clisp because it's small, reasonably fast, and has great bignum arithmetic and the very useful arbitrary precision long-float. If others don't mind the slow down from computing both cos and sin to get just sin or cose, I would find the speed to be acceptable. (I haven't measured the change, though). Thanks for your great work in improving Clisp! Ray ```
 Re: Numeric inconsistency From: Sam Steingold - 2002-09-10 14:56:34 ```> * In message <4n4rcx3o3d.fsf@...> > * On the subject of "Re: Numeric inconsistency" > * Sent on 10 Sep 2002 10:40:06 -0400 > * Honorable Raymond Toy writes: > > >>>>> "Sam" == Sam Steingold writes: > > Sam> the appended patch fixes the expt precision problem which you > Sam> reported in another e-mail. I would greatly appreciate it if > Sam> you could test it as much as possible. > > I'll try this when I get a chance. > > However, let me note that I no longer think it is a problem. I accept > Bruno's explanation that clisp returns a complex number for the test > case I gave. I might wish it were a bit more accurate, but I don't > consider the result given by Clisp to be wrong. It is, after all, a > reasonable result. No, the number returned is, of course, #C() -- Bruno is right. the accuracy is better because all computations are done in the same precision. > Thanks for your great work in improving Clisp! thanks for your testing! -- Sam Steingold (http://www.podval.org/~sds) running RedHat7.3 GNU/Linux ; ; ; ; ; Don't hit a man when he's down -- kick him; it's easier. ```
 Re: Numeric inconsistency From: Sam Steingold - 2002-09-10 23:10:33 ```> * In message <4n4rcx3o3d.fsf@...> > * On the subject of "Re: Numeric inconsistency" > * Sent on 10 Sep 2002 10:40:06 -0400 > * Honorable Raymond Toy writes: > > >>>>> "Sam" == Sam Steingold writes: > > Sam> the appended patch fixes the expt precision problem which you > Sam> reported in another e-mail. I would greatly appreciate it if > Sam> you could test it as much as possible. > > I'll try this when I get a chance. please try the appended patch instead. people, everyone who knows what TAN and TANH are and how they are related, please do \$ cvs up \$ patch < this.diff and play with the transcendental functions. this is important and I mean _you_. please do not think "with an invitation like this lots of people will be doing this, so I don't need to". Yes, you _do_ need to do this - why are you reading clisp-devel otherwise?! consider this the official pre-release test annouoncement. if you will not report bugs in CVS HEAD, it will become CLISP 2.30. -- Sam Steingold (http://www.podval.org/~sds) running RedHat7.3 GNU/Linux ; ; ; ; ; cogito cogito ergo cogito sum Index: src/realtran.d =================================================================== RCS file: /cvsroot/clisp/clisp/src/realtran.d,v retrieving revision 1.8 diff -u -w -b -r1.8 realtran.d --- src/realtran.d 4 Sep 2002 19:48:24 -0000 1.8 +++ src/realtran.d 10 Sep 2002 22:55:45 -0000 @@ -91,7 +91,7 @@ return O(SF_pi), # pi as SF return O(FF_pi), # pi as FF return O(DF_pi), # pi as DF - return O(pi) , # pi as LF with default length + return pi_F_float_F(x), # pi as LF in the same format as x ,); # nothing to save } @@ -494,6 +494,10 @@ }} /* R_cos_sin_R_R(x) places ((cos x),(sin x)), on the Stack. + when start_p, this is a start of a computation, + so the accuracy will be raised + when end_p, this is also the end of the computation, + so the accuracy will be lowed back to this number can trigger GC Method: x rational -> if x=0 ==> (1,0), otherwise x ==> float. @@ -514,7 +518,8 @@ if q = 1 mod 4: ((- (sin r)), (cos r)) if q = 2 mod 4: ((- (cos r)), (- (sin r))) if q = 3 mod 4: ((sin r), (- (cos r))) */ -local void R_cos_sin_R_R (object x) { +local void R_cos_sin_R_R (object x,bool start_p,object *end_p) +{ if (R_rationalp(x)) { if (eq(x,Fixnum_0)) # x=0 -> return (1,0) { pushSTACK(Fixnum_1); pushSTACK(Fixnum_0); return; } @@ -522,14 +527,20 @@ } /* x -- float */ pushSTACK(x); /* save x */ - x = F_extend_F(x); /* increase computational accuracy */ + if (start_p) /* increase computational precision */ + x = F_extend_F(x); F_pi2_round_I_F(x); /* divide by pi/2 */ /* stack layout: Argument, q, r. */ x = STACK_0; if (R_zerop(x) /* r=0.0 -> cos=1.0+O(r^2), sin=r+O(r^3) */ || (F_exponent_L(x) <= (sintL)(-F_float_digits(x))>>1)) { /* e <= -d/2 <==> e <= -ceiling(d/2) ? */ - pushSTACK(I_F_float_F(Fixnum_1,STACK_2)); /* cos=1 */ - pushSTACK(F_F_float_F(STACK_1,STACK_3)); /* sin=r */ + if (end_p != NULL) { + pushSTACK(RA_R_float_F(Fixnum_1,*end_p)); /* cos=1 */ + pushSTACK(F_R_float_F(STACK_1,*end_p)); /* sin=r */ + } else { + pushSTACK(I_F_float_F(Fixnum_1,STACK_0)); /* cos=1 */ + pushSTACK(STACK_1); /* sin=r */ + } } else { pushSTACK(F_I_scale_float_F(STACK_0,Fixnum_minus1)); /* s := r/2 */ pushSTACK(F_sinx_F(STACK_0)); /* y := (sin(s)/s)^2 */ @@ -537,7 +548,9 @@ x = F_F_mal_F(STACK_0,STACK_1); /* y*s */ x = F_F_mal_F(STACK_2,x); /* y*s*r */ x = R_R_minus_R(Fixnum_1,x); /* 1-y*s*r */ - pushSTACK(F_F_float_F(x,STACK_4)); /* scale and save cos(r) */ + if (end_p != NULL) /* scale and save cos(r) */ + pushSTACK(F_R_float_F(x,*end_p)); + else pushSTACK(x); x = F_F_mal_F(STACK_1,STACK_2); /* y*s */ x = F_F_mal_F(x,STACK_2); /* y*s*s */ x = R_R_minus_R(Fixnum_1,x); /* 1-y*s*s = cos(s)^2 */ @@ -545,7 +558,8 @@ x = F_sqrt_F(x); /* cos(s)*sin(s)/s */ x = F_F_mal_F(x,STACK_2); /* cos(s)*sin(s) */ x = F_I_scale_float_F(x,Fixnum_1); /* 2*cos(s)*sin(s) = sin(r) */ - x = F_F_float_F(x,STACK_5); /* scale sin(r) */ + if (end_p != NULL) /* scale sin(r) */ + x = F_R_float_F(x,*end_p); STACK_2 = STACK_0; STACK_1 = x; skipSTACK(1); @@ -567,6 +581,7 @@ } } skipSTACK(2+1); + return; } # F_lnx_F(x) liefert zu einem Float x (>=1/2, <=2) ln(x) als Float. @@ -699,7 +714,6 @@ # R_ln_R(x) liefert zu einer reellen Zahl x>0 die Zahl ln(x). # can trigger GC - local object R_ln_R (object x); # Methode: # x rational -> bei x=1 0 als Ergebnis, sonst x in Float umwandeln. # x Float -> @@ -708,34 +722,35 @@ # (m,e) := (decode-float x), so dass 1/2 <= m < 1. # m<2/3 -> m:=2m, e:=e-1, so dass 2/3 <= m <= 4/3. # ln(m) errechnen, ln(x)=ln(m)+e*ln(2) als Ergebnis. - local object R_ln_R(x) - var object x; - { if (R_rationalp(x)) - { if (eq(x,Fixnum_1)) { return Fixnum_0; } # x=1 -> 0 als Ergebnis - x = RA_float_F(x); # sonst in Float umwandeln +local object R_ln_R (object x, bool start_p, object* end_p) +{ + if (R_rationalp(x)) { + if (eq(x,Fixnum_1)) { return Fixnum_0; } /* x=1 -> return 0 */ + x = RA_float_F(x); /* convert to float */ } - # x Float - pushSTACK(x); # x retten - x = F_extend2_F(x); # Rechengenauigkeit erhöhen - F_decode_float_F_I_F(x); # m,e,s bestimmen - # Stackaufbau: x, m, e, s. + /* x -- float */ + pushSTACK(x); /* save x */ + if (start_p) /* increase accuracy */ + x = F_extend2_F(x); + F_decode_float_F_I_F(x); /* compute m,e,s */ + /* Stack layout: x, m, e, s. */ if (F_F_comp(STACK_2, - make_SF(0,0+SF_exp_mid,floor(bit(SF_mant_len+2),3)) # Short-Float 2/3 - ) - < 0 - ) # m < 2/3 -> - { STACK_2 = F_I_scale_float_F(STACK_2,Fixnum_1); # m verdoppeln - STACK_1 = I_minus1_plus_I(STACK_1); # e decrementieren + make_SF(0,0+SF_exp_mid,floor(bit(SF_mant_len+2),3))) /* short-float 2/3 */ + < 0) { /* m < 2/3 -> */ + STACK_2 = F_I_scale_float_F(STACK_2,Fixnum_1); /* double m */ + STACK_1 = I_minus1_plus_I(STACK_1); /* decrement e */ } - STACK_2 = F_lnx_F(STACK_2); # ln(m) im genaueren Float-Format errechnen + STACK_2 = F_lnx_F(STACK_2); /* ln(m) in the more accurate float format */ {var object temp; - temp = ln2_F_float_F(STACK_0); # ln(2) im genaueren Float-Format - temp = R_R_mal_R(STACK_1,temp); # e*ln(2) - temp = R_R_plus_R(STACK_2,temp); # ln(m)+e*ln(2) - temp = F_F_float_F(temp,STACK_3); # (float ... x) + temp = ln2_F_float_F(STACK_0); /* ln(2) in that float format */ + temp = R_R_mal_R(STACK_1,temp); /* e*ln(2) */ + temp = R_R_plus_R(STACK_2,temp); /* ln(m)+e*ln(2) */ + if (end_p != NULL) /* (float ... x) */ + temp = F_R_float_F(temp,*end_p); skipSTACK(4); return temp; - }} + } +} #define F_ln_F R_ln_R # I_I_log_RA(a,b) liefert zu Integers a>0, b>1 den Logarithmus log(a,b) @@ -902,9 +917,13 @@ STACK_1 = RA_F_float_F(a,b); # a := (float a b) } } # Nun a,b beide Floats. - STACK_1 = R_ln_R(STACK_1); # (ln a) errechnen - {var object lnb = R_ln_R(popSTACK()); # (ln b) errechnen - return F_F_durch_F(popSTACK(),lnb); # (/ (ln a) (ln b)) als Ergebnis + pushSTACK(R_ln_R(STACK_1,true,NULL)); /* (ln a) */ + pushSTACK(R_ln_R(STACK_1,true,NULL)); /* (ln b) */ + STACK_0 = F_F_durch_F(STACK_1,STACK_0); /* (/ (ln a) (ln b)) */ + STACK_1 = R_R_contagion_R(STACK_2,STACK_3); + { var object ret = F_R_float_F(STACK_0,STACK_1); + skipSTACK(4); + return ret; }} # F_expx_F(x) liefert zu einem Float x (betragsmäig <1) exp(x) als Float. @@ -966,7 +985,6 @@ # R_exp_R(x) liefert zu einer reellen Zahl x die Zahl exp(x). # can trigger GC - local object R_exp_R (object x); # Methode: # x rational -> bei x=0 1 als Ergebnis, sonst x in Float umwandeln. # x Float -> @@ -974,32 +992,35 @@ # Genauigkeit um sqrt(d)+max(integer-length(e)) Bits erhöhen, # (q,r) := (floor x ln(2)) # Ergebnis ist exp(q*ln(2)+r) = (scale-float exp(r) q). - local object R_exp_R(x) - var object x; - { if (R_rationalp(x)) - # x rational - { if (eq(x,Fixnum_0)) { return Fixnum_1; } # x=0 -> 1 als Ergebnis - x = RA_float_F(x); # sonst in Float umwandeln +local object R_exp_R (object x, bool start_p, object* end_p) +{ + if (R_rationalp(x)) { /* x rational */ + if (eq(x,Fixnum_0)) { return Fixnum_1; } /* x=0 -> return 1 */ + x = RA_float_F(x); /* ==> float */ } - # x Float - pushSTACK(x); # x retten - x = F_extend2_F(x); # Genauigkeit vergröern - # durch ln(2) dividieren (bei 0<=x<1/2 kann man sofort q:=0 setzen) - if ((!R_minusp(x)) && (F_exponent_L(x)<0)) - { pushSTACK(Fixnum_0); pushSTACK(x); } # x>=0, Exponent <0 -> 0<=x<1/2 -> Division unnötig - else - { pushSTACK(x); - {var object ln2 = ln2_F_float_F(x); # ln(2) mit hinreichender Genauigkeit + /* x -- float */ + pushSTACK(x); /* save x */ + if (start_p) /* increase accuracy */ + x = F_extend2_F(x); + /* divide by ln(2) (if 0<=x<1/2 can immediately set q:=0) */ + if ((!R_minusp(x)) && (F_exponent_L(x)<0)) { + /* x>=0, Exponent <0 -> 0<=x<1/2 -> division not necessary */ + pushSTACK(Fixnum_0); pushSTACK(x); + } else { + pushSTACK(x); + { var object ln2 = ln2_F_float_F(x); /* ln(2) with sufficient accuracy */ x = popSTACK(); - F_F_floor_I_F(x,ln2); # x durch ln(2) dividieren + F_F_floor_I_F(x,ln2); /* x / ln(2) */ }} - # Stackaufbau: originales x, q, r. - {var object temp = F_expx_F(STACK_0); # exp(r) - temp = F_I_scale_float_F(temp,STACK_1); # mal 2^q - temp = F_F_float_F(temp,STACK_2); # (float ... x) als Ergebnis + /* stack layout: original x, q, r. */ + { var object temp = F_expx_F(STACK_0); /* exp(r) */ + temp = F_I_scale_float_F(temp,STACK_1); /* * 2^q */ + if (end_p != NULL) /* (float ... x) als Ergebnis */ + temp = F_R_float_F(temp,*end_p); skipSTACK(3); return temp; - }} + } +} # R_sinh_R(x) liefert zu einer reellen Zahl x die Zahl sinh(x). # can trigger GC @@ -1023,7 +1044,7 @@ { var object temp; pushSTACK(x); pushSTACK(temp = F_extend_F(x)); # Rechengenauigkeit erhöhen - temp = F_sqrt_F(F_sinhx_F(x)); # Wurzel aus (sinh(x)/x)^2 + temp = F_sqrt_F(F_sinhx_F(temp)); # Wurzel aus (sinh(x)/x)^2 temp = F_F_mal_F(temp,STACK_0); # mit genauerem x multiplizieren temp = F_F_float_F(temp,STACK_1); # und wieder runden skipSTACK(2); @@ -1032,10 +1053,12 @@ else # e>0 -> verwende exp(x) { var object temp; - pushSTACK(temp = R_exp_R(x)); # y:=exp(x) + pushSTACK(x); + pushSTACK(temp = R_exp_R(x,true,NULL)); /* y:=exp(x) */ temp = F_durch_F(temp); # (/ y) temp = F_F_minus_F(popSTACK(),temp); # von y subtrahieren - return F_I_scale_float_F(temp,Fixnum_minus1); # (scale-float ... -1) + return F_F_float_F(F_I_scale_float_F(temp,Fixnum_minus1), + popSTACK()); /* (scale-float ... -1) */ } } # R_cosh_R(x) liefert zu einer reellen Zahl x die Zahl cosh(x). @@ -1064,10 +1087,12 @@ if (e > 0) # e>0 -> verwende exp(x) { var object temp; - pushSTACK(temp = R_exp_R(x)); # y:=exp(x) + pushSTACK(x); + pushSTACK(temp = R_exp_R(x,true,NULL)); /* y:=exp(x) */ temp = F_durch_F(temp); # (/ y) temp = F_F_plus_F(popSTACK(),temp); # zu y addieren - return F_I_scale_float_F(temp,Fixnum_minus1); # (scale-float ... -1) + return F_F_float_F(F_I_scale_float_F(temp,Fixnum_minus1), + popSTACK()); /* (scale-float ... -1) */ } else # e<=0 @@ -1092,7 +1117,6 @@ # R_cosh_sinh_R_R(x) liefert ((cosh x),(sinh x)), beide Werte auf dem Stack. # can trigger GC - local void R_cosh_sinh_R_R (object x); # Methode: # x rational -> bei x=0 (1,0) als Ergebnis, sonst x in Float umwandeln. # x Float -> Genauigkeit erhöhen, @@ -1108,48 +1132,55 @@ # falls e>0: y:=exp(x) errechnen, # (scale-float (+ y (/ y)) -1) und (scale-float (- y (/ y)) -1) bilden. # Genauigkeit wieder verringern. - local void R_cosh_sinh_R_R(x) - var object x; - { if (R_rationalp(x)) - # x rational - { if (eq(x,Fixnum_0)) { pushSTACK(Fixnum_1); pushSTACK(Fixnum_0); return; } # x=0 -> (1,0) als Ergebnis - x = RA_float_F(x); # sonst in Float umwandeln +local void R_cosh_sinh_R_R (object x, bool start_p, object* end_p) +{ + if (R_rationalp(x)) { /* x rational */ + if (eq(x,Fixnum_0)) /* x=0 -> return (1,0) */ + { pushSTACK(Fixnum_1); pushSTACK(Fixnum_0); return; } + x = RA_float_F(x); /* ==> Float */ } - # x Float + /* x -- float */ {var sintL e = F_exponent_L(x); - if (e > 0) - # e>0 -> verwende exp(x) - { var object temp; - pushSTACK(temp = R_exp_R(x)); # y:=exp(x) - pushSTACK(temp = F_durch_F(temp)); # (/ y) - # Stackaufbau: exp(x), exp(-x). - temp = F_F_minus_F(STACK_1,temp); # von y subtrahieren - temp = F_I_scale_float_F(temp,Fixnum_minus1); # (scale-float ... -1) - {var object temp2 = STACK_0; - STACK_0 = temp; - temp = F_F_plus_F(STACK_1,temp2); # zu y addieren - STACK_1 = F_I_scale_float_F(temp,Fixnum_minus1); # (scale-float ... -1) + if (e > 0) { /* e>0 -> use exp(x) */ + var object temp; + pushSTACK(x); + pushSTACK(temp = R_exp_R(x,start_p,NULL)); /* y:=exp(x) */ + pushSTACK(temp = F_durch_F(temp)); /* (/ y) */ + /* stack layout: x, exp(x), exp(-x). */ + temp = F_F_plus_F(STACK_1,temp); /* + y */ + temp = F_I_scale_float_F(temp,Fixnum_minus1); /* /2 */ + if (end_p != NULL) /* cosh */ + STACK_2 = F_F_float_F(temp,*end_p); + else STACK_2 = temp; + temp = F_F_minus_F(STACK_1,STACK_0); /* - y */ + STACK_1 = F_I_scale_float_F(temp,Fixnum_minus1); /* /2 */ + if (end_p != NULL) /* sinh */ + STACK_1 = F_F_float_F(STACK_1,*end_p); + skipSTACK(1); return; - }} - else - # e<=0 - { if (R_zerop(x) - || (e <= (sintL)(1-F_float_digits(x))>>1) # e <= (1-d)/2 <==> e <= -ceiling((d-1)/2) ? - ) - { pushSTACK(x); pushSTACK(x); STACK_1 = I_F_float_F(Fixnum_1,x); return; } - {var object temp; + } else { /* e<=0 */ + if (R_zerop(x) + || (e <= (sintL)(1-F_float_digits(x))>>1)) { /* e <= (1-d)/2 <==> e <= -ceiling((d-1)/2) ? */ + pushSTACK(x); pushSTACK(x); STACK_1 = I_F_float_F(Fixnum_1,x); return; + } pushSTACK(x); - pushSTACK(temp = F_extend_F(x)); # Rechengenauigkeit erhöhen - pushSTACK(F_square_F(temp)); # x*x - pushSTACK(temp = F_sinhx_F(STACK_1)); # y:=(sinh(x)/x)^2 - # Stackaufbau: originales x, x, x^2, y. - temp = F_sqrt_F(temp); # sqrt(y) = sinh(x)/x - temp = F_F_mal_F(STACK_2,temp); # x*sqrt(y) = sinh(x) - STACK_2 = F_F_float_F(temp,STACK_3); # und wieder runden - temp = F_F_mal_F(STACK_1,STACK_0); # x^2*y - temp = F_sqrt_F(R_R_plus_R(Fixnum_1,temp)); # sqrt(1+x^2*y) - STACK_3 = F_F_float_F(temp,STACK_3); # und wieder runden + { var object temp = (start_p ? F_extend_F(x) : x); + pushSTACK(temp); + pushSTACK(F_square_F(temp)); /* x*x */ + pushSTACK(temp = F_sinhx_F(STACK_1)); /* y:=(sinh(x)/x)^2 */ + /* stack layout: original x, x, x^2, y. */ + temp = F_sqrt_F(temp); /* sqrt(y) = sinh(x)/x */ + temp = F_F_mal_F(STACK_2,temp); /* x*sqrt(y) = sinh(x) */ + if (end_p != NULL) /* restore the accuracy */ + STACK_2 = F_F_float_F(temp,STACK_3); + else STACK_2 = temp; + temp = F_F_mal_F(STACK_1,STACK_0); /* x^2*y */ + temp = F_sqrt_F(R_R_plus_R(Fixnum_1,temp)); /* sqrt(1+x^2*y) */ + if (end_p != NULL) /* restore the accuracy */ + STACK_3 = F_F_float_F(temp,STACK_3); + else STACK_3 = temp; skipSTACK(2); return; - }} - } } - + } + } + } +} Index: src/realelem.d =================================================================== RCS file: /cvsroot/clisp/clisp/src/realelem.d,v retrieving revision 1.15 diff -u -w -b -r1.15 realelem.d --- src/realelem.d 26 Feb 2002 14:18:14 -0000 1.15 +++ src/realelem.d 10 Sep 2002 22:55:45 -0000 @@ -151,6 +151,23 @@ #undef X } +local object N_N_contagion_R (object x, object y) +{ + if (complexp(x)) + if (complexp(y)) + return R_R_contagion_R(R_R_contagion_R(TheComplex(x)->c_real, + TheComplex(x)->c_imag), + R_R_contagion_R(TheComplex(y)->c_real, + TheComplex(y)->c_imag)); + else return R_R_contagion_R(R_R_contagion_R(TheComplex(x)->c_real, + TheComplex(x)->c_imag),y); + else if (complexp(y)) + return R_R_contagion_R(x,R_R_contagion_R(TheComplex(y)->c_real, + TheComplex(y)->c_imag)); + else + return R_R_contagion_R(x,y); +} + # Macro: verteilt je nach Default-Float-Typ auf 4 Statements. # defaultfloatcase(symbol, SF_statement,FF_statement,DF_statement,LF_statement, save_statement,restore_statement); # symbol sollte ein S(..)-Symbol sein. Dessen Wert sollte SHORT-FLOAT oder @@ -218,6 +235,52 @@ local object R_float_F(x) var object x; { return (R_rationalp(x) ? RA_float_F(x) : x); } + +/* convert a float(rational) X to the appropriate format for Y + same as F_F_float_F(x,R_float_F(y)) but without consing + an intermediate float + can trigger GC */ +local object F_R_float_F (object x, object y) +{ + #if SAFETY>=1 + if (!floatp(x)) abort(); + if_realp(y,,abort();); + #endif + defaultfloatcase(S(default_float_format),y, + return F_to_SF(x), + return F_to_FF(x), + return F_to_DF(x), + return F_to_LF(x,I_to_UL(O(LF_digits))), + pushSTACK(x), x = popSTACK()); +} +local object RA_R_float_F (object x, object y) +{ + defaultfloatcase(S(default_float_format),y, + return RA_to_SF(x), + return RA_to_FF(x), + return RA_to_DF(x), + return RA_to_LF(x,I_to_UL(O(LF_digits))), + pushSTACK(x), x = popSTACK()); +} +local object R_R_float_F (object x, object y) +{ + if (floatp(x)) return F_R_float_F(x,y); + return RA_R_float_F(x,y); +} +local object C_R_float_C (object *c, object y) +{ + TheComplex(*c)->c_real = R_R_float_F(TheComplex(*c)->c_real,y); + TheComplex(*c)->c_imag = R_R_float_F(TheComplex(*c)->c_imag,y); + return *c; +} +local object N_N_float_N (object *n, object y) +{ + var object contagion = !complexp(y) ? y + : R_R_contagion_R(TheComplex(y)->c_real,TheComplex(y)->c_imag); + if (complexp(*n)) + return C_R_float_C(n,contagion); + return R_R_float_F(*n,contagion); +} # Generiert eine Funktion wie R_floor_I_R #define GEN_R_round(rounding) \ Index: src/lisparit.d =================================================================== RCS file: /cvsroot/clisp/clisp/src/lisparit.d,v retrieving revision 1.36 diff -u -w -b -r1.36 lisparit.d --- src/lisparit.d 22 Jul 2002 12:08:56 -0000 1.36 +++ src/lisparit.d 10 Sep 2002 22:55:45 -0000 @@ -978,12 +978,11 @@ mv_count=1; set_args_end_pointer(rest_args_pointer); } -LISPFUNN(exp,1) # (EXP number), CLTL S. 203 - { - var object x = popSTACK(); - check_number(x); - value1 = N_exp_N(x); mv_count=1; +LISPFUNN(exp,1) { + check_number(STACK_0); + value1 = N_exp_N(STACK_0,true,&STACK_0); mv_count=1; + skipSTACK(1); } LISPFUNN(expt,2) @@ -998,17 +997,18 @@ LISPFUN(log,1,1,norest,nokey,0,NIL) # (LOG number [base-number]), CLTL S. 204 { - var object base = popSTACK(); - var object arg = popSTACK(); + var object base = STACK_0; + var object arg = STACK_1; check_number(arg); if (eq(base,unbound)) { # LOG mit einem Argument - value1 = N_log_N(arg); + value1 = N_log_N(arg,true,&STACK_1); } else { # LOG mit zwei Argumenten check_number(base); value1 = N_N_log_N(arg,base); } + skipSTACK(2); mv_count=1; } @@ -2071,7 +2071,7 @@ # > 1 wachsen lassen, damit es nicht zu häufig nachberechnet wird: oldlen += floor(oldlen,2); # oldlen * 3/2 var uintC newlen = (d < oldlen ? oldlen : d); - ln_x = *objptr = R_ln_R(I_to_LF(x,newlen)); # (ln x) als LF mit newlen Digits berechnen + ln_x = *objptr = R_ln_R(I_to_LF(x,newlen),true,objptr); # (ln x) als LF mit newlen Digits berechnen return (d < newlen ? LF_shorten_LF(ln_x,d) : ln_x); } else { # ein Double-Float reicht Index: src/constsym.d =================================================================== RCS file: /cvsroot/clisp/clisp/src/constsym.d,v retrieving revision 1.161 diff -u -w -b -r1.161 constsym.d --- src/constsym.d 7 Aug 2002 03:28:35 -0000 1.161 +++ src/constsym.d 10 Sep 2002 22:55:45 -0000 @@ -915,6 +915,7 @@ LISPSYM(set_stream_external_format,"SET-STREAM-EXTERNAL-FORMAT",system) LISPSYM(interactive_stream_p,"INTERACTIVE-STREAM-P",lisp) LISPSYM(built_in_stream_close,"BUILT-IN-STREAM-CLOSE",system) +LISPSYM(close,"CLOSE",lisp) /* error reporting in stream.d */ LISPSYM(read_byte,"READ-BYTE",lisp) LISPSYM(read_byte_lookahead,"READ-BYTE-LOOKAHEAD",ext) LISPSYM(read_byte_will_hang_p,"READ-BYTE-WILL-HANG-P",ext) Index: src/comptran.d =================================================================== RCS file: /cvsroot/clisp/clisp/src/comptran.d,v retrieving revision 1.11 diff -u -w -b -r1.11 comptran.d --- src/comptran.d 4 Sep 2002 19:50:27 -0000 1.11 +++ src/comptran.d 10 Sep 2002 22:55:45 -0000 @@ -21,47 +21,65 @@ # N_exp_N(x) liefert (exp x), wo x eine Zahl ist. # can trigger GC - local object N_exp_N (object x); # Methode: # x reell -> klar. # x = a+bi -> (exp a) mit (cos b) + i (sin b) multiplizieren: # (complex (* (exp a) (cos b)) (* (exp a) (sin b))) - local object N_exp_N(x) - var object x; +local object N_exp_N (object x, bool start_p, object* end_p) { if (N_realp(x)) { - return R_exp_R(x); - } else { - # x=a+bi - pushSTACK(TheComplex(x)->c_real); # a retten - R_cos_sin_R_R(TheComplex(x)->c_imag); # (cos b) und (sin b) errechnen - # Stackaufbau: a, cos(b), sin(b). - # Da b nicht = Fixnum 0 ist, ist auch sin(b) nicht = Fixnum 0. - STACK_2 = x = R_exp_R(STACK_2); # (exp a) - # Stackaufbau: exp(a), cos(b), sin(b). - STACK_0 = R_R_mal_R(x,STACK_0); # (* (exp a) (sin b)), nicht = Fixnum 0 - x = R_R_mal_R(STACK_2,STACK_1); # (* (exp a) (cos b)) - x = R_R_complex_C(x,STACK_0); # (complex ... ...) + return R_exp_R(x,start_p,end_p); + } else { /* x=a+bi */ + pushSTACK(TheComplex(x)->c_real); /* save a */ + R_cos_sin_R_R(TheComplex(x)->c_imag,start_p,NULL); /* (cos b), (sin b) */ + /* stack layout: a, cos(b), sin(b). */ + /* b != Fixnum_0 ==> sin(b) ~= Fixnum_0. */ + STACK_2 = R_exp_R(STACK_2,start_p,NULL); /* (exp a) */ + /* stack layout: exp(a), cos(b), sin(b). */ + STACK_0 = R_R_mal_R(STACK_2,STACK_0); /* (* (exp a) (sin b)) != Fixnum_0 */ + STACK_1 = R_R_mal_R(STACK_2,STACK_1); /* (* (exp a) (cos b)) */ + x = R_R_complex_C(F_R_float_F(STACK_1,*end_p), + F_R_float_F(STACK_0,*end_p)); /* (complex ... ...) */ skipSTACK(3); return x; } } +local int lfl (object x) { + if (floatp(x)) + floatcase(x,return -1,return -2,return -3,return Lfloat_length(x)); + else return -4; +} + # N_log_N(x) liefert (log x), wo x eine Zahl ist. # can trigger GC - local object N_log_N (object x); # Methode: # (complex (log (abs x)) (phase x)) - local object N_log_N(x) - var object x; +local object N_log_N (object x, bool start_p, object *end_p) { - pushSTACK(x); # x retten - var object r = N_abs_R(x); # (abs x) - if (R_zerop(r)) # (abs x) = 0 -> Error + pushSTACK(x); /* save x */ + pushSTACK(N_abs_R(x)); /* (abs x) */ + if (R_zerop(STACK_0)) /* (abs x) = 0 -> Error */ divide_0(); - r = R_ln_R(r); # (log (abs x)) - x = STACK_0; STACK_0 = r; - x = N_phase_R(x); # (phase x) - return R_R_complex_N(popSTACK(),x); # (complex r (phase x)) + STACK_0 = R_ln_R(STACK_0,start_p,end_p); /* (log (abs x)) */ + if (start_p) { /* increase precision */ + if (floatp(STACK_1)) + STACK_1 = F_extend2_F(STACK_1); + else if (complexp(STACK_1)) { + if (floatp(TheComplex(STACK_1)->c_real)) + TheComplex(STACK_1)->c_real = + F_extend2_F(TheComplex(STACK_1)->c_real); + if (floatp(TheComplex(STACK_1)->c_imag)) + TheComplex(STACK_1)->c_imag = + F_extend2_F(TheComplex(STACK_1)->c_imag); + }; + } + STACK_1 = N_phase_R(STACK_1); /* (phase x) */ + if (end_p != NULL && floatp(STACK_1)) + STACK_1 = F_R_float_F(STACK_1,*end_p); + { /* (complex (log (abs x)) (phase x)) */ + var object ret = R_R_complex_N(STACK_0,STACK_1); + skipSTACK(2); return ret; + } } # N_N_log_N(a,b) liefert (log a b), wo a und b Zahlen sind. @@ -109,7 +127,7 @@ b = STACK_1; if (R_rationalp(b)) b = RA_F_float_F(b,angle); - b = F_ln_F(b); STACK_0 = F_F_durch_F(STACK_0,b); + b = F_ln_F(b,true,&STACK_1); STACK_0 = F_F_durch_F(STACK_0,b); } # Stackaufbau: a, b, Imaginärteil. # Realteil (/ (log (abs a)) (log b)) errechnen: @@ -139,22 +157,23 @@ } } # Keine Chance für rationalen Realteil - pushSTACK(F_ln_F(N_abs_R(a))); # (log (abs a)), ein Float + pushSTACK(F_ln_F(N_abs_R(a),true,&STACK_3)); /* (log (abs a)), a float */ # durch (log b) dividieren, liefert den Realteil: b = STACK_2; if (R_rationalp(b)) b = RA_F_float_F(b,STACK_0); - b = F_ln_F(b); STACK_0 = F_F_durch_F(STACK_0,b); + b = F_ln_F(b,true,&STACK_2); STACK_0 = F_F_durch_F(STACK_0,b); real_ok: # Stackaufbau: a, b, Imaginärteil, Realteil. var object erg = R_R_complex_C(STACK_0,STACK_1); skipSTACK(4); return erg; } - } else { - # normaler komplexer Fall - pushSTACK(b); a = N_log_N(a); b = STACK_0; # (log a) - STACK_0 = a; b = N_log_N(b); a = popSTACK(); # (log b) - return N_N_durch_N(a,b); # dividieren + } else { /* normal complex case */ + pushSTACK(a); pushSTACK(b); + STACK_1 = N_log_N(STACK_1,true,&STACK_1); /* (log a) */ + STACK_0 = N_log_N(STACK_0,true,&STACK_0); /* (log b) */ + a = N_N_durch_N(a,b); /* divide */ + skipSTACK(2); return a; } } @@ -365,204 +384,224 @@ } } pushSTACK(y); - var object temp = N_log_N(x); # (log x) - temp = N_N_mal_N(temp,popSTACK()); # mal y - return N_exp_N(temp); # exp davon + pushSTACK(x); + pushSTACK(N_N_contagion_R(x,y)); + STACK_1 = N_log_N(x,true,NULL); /* (log x) */ + STACK_2 = N_N_float_N(&STACK_2,STACK_1); + x = N_N_mal_N(STACK_2,STACK_1); /* (* (log x) y) */ + x = N_exp_N(x,false,&STACK_0); /* exp */ + skipSTACK(3); return x; } # N_sin_N(x) liefert (sin x), wo x eine Zahl ist. # can trigger GC - local object N_sin_N (object x); # Methode: # x reell -> klar # x = a+bi -> (complex (* (sin a) (cosh b)) (* (cos a) (sinh b))) - local object N_sin_N(x) - var object x; +local object N_sin_N (object x) { if (N_realp(x)) { return R_sin_R(x); - } else { - # x=a+bi - pushSTACK(TheComplex(x)->c_real); # a retten - R_cosh_sinh_R_R(TheComplex(x)->c_imag); # cosh(b), sinh(b) errechnen - # Stackaufbau: a, cosh(b), sinh(b). - # Da b nicht = Fixnum 0 ist, ist auch sinh(b) nicht = Fixnum 0. - R_cos_sin_R_R(STACK_2); # cos(a), sin(a) errechnen, cos(a) /= Fixnum 0 - # Stackaufbau: a, cosh(b), sinh(b), cos(a), sin(a). - STACK_0 = R_R_mal_R(STACK_0,STACK_3); # sin(a)*cosh(b) - var object erg = R_R_mal_R(STACK_1,STACK_2); # cos(a)*sinh(b), nicht Fixnum 0 - erg = R_R_complex_C(STACK_0,erg); skipSTACK(5); return erg; + } else { /* x=a+bi */ + pushSTACK(TheComplex(x)->c_real); /* save a */ + pushSTACK(TheComplex(x)->c_imag); /* save b */ + R_cosh_sinh_R_R(STACK_0,true,NULL); /* cosh(b) sinh(b) */ + /* stack layout: a, b, cosh(b), sinh(b). */ + /* b != Fixnum_0 ==> sinh(b) != Fixnum_0. */ + R_cos_sin_R_R(STACK_3,true,NULL); /* cos(a)!=0, sin(a); */ + /* stack layout: a, b, cosh(b), sinh(b), cos(a), sin(a). */ + STACK_0 = R_R_mal_R(STACK_0,STACK_3); /* sin(a)*cosh(b) */ + STACK_3 = R_R_contagion_R(STACK_4,STACK_5); + { var object erg = R_R_mal_R(STACK_1,STACK_2); /* cos(a)*sinh(b), != Fixnum_0 */ + erg = R_R_complex_C(R_R_float_F(STACK_0,STACK_3), + F_R_float_F(erg,STACK_3)); + skipSTACK(6); return erg; + } } } # N_cos_N(x) liefert (cos x), wo x eine Zahl ist. # can trigger GC - local object N_cos_N (object x); # Methode: # x reell -> klar # x = a+bi -> (complex (* (cos a) (cosh b)) (- (* (sin a) (sinh b)))) - local object N_cos_N(x) - var object x; +local object N_cos_N (object x) { if (N_realp(x)) { return R_cos_R(x); - } else { - # x=a+bi - pushSTACK(TheComplex(x)->c_real); # a retten - R_cosh_sinh_R_R(TheComplex(x)->c_imag); # cosh(b), sinh(b) errechnen - # Stackaufbau: a, cosh(b), sinh(b). - R_cos_sin_R_R(STACK_2); # cos(a), sin(a) errechnen - # Stackaufbau: a, cosh(b), sinh(b), cos(a), sin(a). - STACK_0 = R_minus_R(R_R_mal_R(STACK_0,STACK_2)); # -sin(a)*sinh(b) - var object erg = R_R_mal_R(STACK_1,STACK_3); # cos(a)*cosh(b) - erg = R_R_complex_N(erg,STACK_0); skipSTACK(5); return erg; + } else { /* x=a+bi */ + pushSTACK(TheComplex(x)->c_real); /* save a */ + pushSTACK(TheComplex(x)->c_imag); /* save b */ + R_cosh_sinh_R_R(STACK_0,true,NULL); /* cosh(b), sinh(b) */ + /* stack layout: a, b, cosh(b), sinh(b). */ + R_cos_sin_R_R(STACK_3,true,NULL); /* cos(a), sin(a) */ + /* stack layout: a, b, cosh(b), sinh(b), cos(a), sin(a). */ + STACK_0 = R_minus_R(R_R_mal_R(STACK_0,STACK_2)); /* -sin(a)*sinh(b) */ + STACK_1 = R_R_mal_R(STACK_1,STACK_3); /* cos(a)*cosh(b) */ + STACK_2 = R_R_contagion_R(STACK_4,STACK_5); + { var object erg = (eq(STACK_0,Fixnum_0) ? R_R_float_F(STACK_1,STACK_2) + : R_R_complex_C(R_R_float_F(STACK_1,STACK_2), + F_R_float_F(STACK_0,STACK_2))); + skipSTACK(6); return erg; + } } } # N_tan_N(x) liefert (tan x), wo x eine Zahl ist. # can trigger GC - local object N_tan_N (object x); # Methode: # x reell -> (/ (sin x) (cos x)) # x = a+bi -> (/ (complex (* (sin a) (cosh b)) (* (cos a) (sinh b))) # (complex (* (cos a) (cosh b)) (- (* (sin a) (sinh b)))) ) - local object N_tan_N(x) - var object x; +local object N_tan_N (object x) { if (N_realp(x)) { - R_cos_sin_R_R(x); - # Stackaufbau: cos(x), sin(x). - var object erg = R_R_durch_R(STACK_0,STACK_1); - skipSTACK(2); return erg; - } else { - # x=a+bi - pushSTACK(TheComplex(x)->c_real); # a retten - R_cosh_sinh_R_R(TheComplex(x)->c_imag); # cosh(b), sinh(b) errechnen - # Stackaufbau: a, cosh(b), sinh(b). - R_cos_sin_R_R(STACK_2); # cos(a), sin(a) errechnen - # Stackaufbau: a, cosh(b), sinh(b), cos(a), sin(a). - var object temp; - STACK_4 = R_R_mal_R(STACK_0,STACK_3); # sin(a)*cosh(b) - temp = R_R_mal_R(STACK_1,STACK_2); # cos(a)*sinh(b) /= Fixnum 0 - STACK_4 = R_R_complex_C(STACK_4,temp); # Zähler - # Stackaufbau: Zähler, cosh(b), sinh(b), cos(a), sin(a). - STACK_3 = R_R_mal_R(STACK_1,STACK_3); # cos(a)*cosh(b) - temp = R_minus_R(R_R_mal_R(STACK_0,STACK_2)); # -sin(a)*sinh(b) - temp = R_R_complex_N(STACK_3,temp); # Nenner - temp = N_N_durch_N(STACK_4,temp); # Zähler/Nenner - skipSTACK(5); return temp; + pushSTACK(x); + R_cos_sin_R_R(x,true,NULL); + { /* stack layout: x, cos(x), sin(x). */ + var object erg = F_R_float_F(R_R_durch_R(STACK_0,STACK_1),STACK_2); + skipSTACK(3); return erg; + } + } else { /* x=a+bi */ + pushSTACK(TheComplex(x)->c_real); /* save a */ + pushSTACK(TheComplex(x)->c_imag); /* save b */ + R_cosh_sinh_R_R(STACK_0,true,NULL); /* cosh(b), sinh(b) */ + /* stack layout: a, b, cosh(b), sinh(b). */ + R_cos_sin_R_R(STACK_3,true,NULL); /* cos(a), sin(a) */ + /* stack layout: a, b, cosh(b), sinh(b), cos(a), sin(a). */ + STACK_5 = R_R_contagion_R(STACK_5,STACK_4); + { var object temp; + STACK_4 = R_R_mal_R(STACK_0,STACK_3); /* sin(a)*cosh(b) */ + temp = R_R_mal_R(STACK_1,STACK_2); /* cos(a)*sinh(b) /= Fixnum 0 */ + STACK_4 = R_R_complex_C(STACK_4,temp); /* numerator */ + /* stack layout: contag, numerator, cosh(b), sinh(b), cos(a), sin(a). */ + STACK_3 = R_R_mal_R(STACK_1,STACK_3); /* cos(a)*cosh(b) */ + temp = R_minus_R(R_R_mal_R(STACK_0,STACK_2)); /* -sin(a)*sinh(b) */ + temp = R_R_complex_N(STACK_3,temp); /* denominator */ + STACK_4 = N_N_durch_N(STACK_4,temp); /* numerator/denominator */ + temp = C_R_float_C(&STACK_4,STACK_5); + skipSTACK(6); return temp; + } } } # N_cis_N(x) liefert (cis x), wo x eine Zahl ist. # can trigger GC - local object N_cis_N (object x); # Methode: # x reell -> (complex (cos x) (sin x)) # x = a+bi -> (complex (* (exp (- b)) (cos a)) (* (exp (- b)) (sin a))) - local object N_cis_N(x) - var object x; +local object N_cis_N (object x) { if (N_realp(x)) { - R_cos_sin_R_R(x); - # Stackaufbau: cos(x), sin(x). + pushSTACK(x); + R_cos_sin_R_R(x,true,&STACK_0); + /* stack layout: x, cos(x), sin(x). */ var object erg = R_R_complex_N(STACK_1,STACK_0); - skipSTACK(2); return erg; - } else { - # x=a+bi - pushSTACK(TheComplex(x)->c_imag); # b retten - R_cos_sin_R_R(TheComplex(x)->c_real); # (cos a) und (sin a) errechnen - # Stackaufbau: b, cos(a), sin(a). - STACK_2 = x = R_exp_R(R_minus_R(STACK_2)); # (exp (- b)) - # Stackaufbau: exp(-b), cos(a), sin(a). - STACK_0 = R_R_mal_R(x,STACK_0); # (* (exp (- b)) (sin a)) - x = R_R_mal_R(STACK_2,STACK_1); # (* (exp (- b)) (cos a)) - x = R_R_complex_N(x,STACK_0); # (complex ... ...) - skipSTACK(3); return x; + skipSTACK(3); return erg; + } else { /* x=a+bi */ + pushSTACK(TheComplex(x)->c_real); /* save a */ + pushSTACK(TheComplex(x)->c_imag); /* save b */ + R_cos_sin_R_R(TheComplex(x)->c_real,true,NULL); /* (cos a), (sin a) */ + /* stack layout: a, b, cos(a), sin(a). */ + STACK_3 = R_R_contagion_R(STACK_3,STACK_2); + STACK_2 = x = R_exp_R(R_minus_R(STACK_2),false,NULL); /* (exp (- b)) */ + /* stack layout: exp(-b), cos(a), sin(a). */ + STACK_0 = R_R_mal_R(x,STACK_0); /* (* (exp (- b)) (sin a)) */ + x = R_R_mal_R(STACK_2,STACK_1); /* (* (exp (- b)) (cos a)) */ + x = R_R_complex_N(F_R_float_F(x,STACK_3), /* (complex ... ...) */ + F_R_float_F(STACK_0,STACK_3)); + skipSTACK(4); return x; } } # N_sinh_N(x) liefert (sinh x), wo x eine Zahl ist. # can trigger GC - local object N_sinh_N (object x); # Methode: # x reell -> klar # x = a+bi -> (complex (* (sinh a) (cos b)) (* (cosh a) (sin b))) - local object N_sinh_N(x) - var object x; +local object N_sinh_N (object x) { if (N_realp(x)) { return R_sinh_R(x); - } else { - # x=a+bi - pushSTACK(TheComplex(x)->c_real); # a retten - R_cos_sin_R_R(TheComplex(x)->c_imag); # cos(b), sin(b) errechnen - # Stackaufbau: a, cos(b), sin(b). - # Da b nicht = Fixnum 0 ist, ist auch sin(b) nicht = Fixnum 0. - R_cosh_sinh_R_R(STACK_2); # cosh(a), sinh(a) errechnen, cosh(a) /= Fixnum 0 - # Stackaufbau: a, cos(b), sin(b), cosh(a), sinh(a). - STACK_0 = R_R_mal_R(STACK_0,STACK_3); # sinh(a)*cos(b) - var object erg = R_R_mal_R(STACK_1,STACK_2); # cosh(a)*sin(b), nicht Fixnum 0 - erg = R_R_complex_C(STACK_0,erg); skipSTACK(5); return erg; + } else { /* x=a+bi */ + pushSTACK(TheComplex(x)->c_real); /* save a */ + pushSTACK(TheComplex(x)->c_imag); /* save b */ + R_cos_sin_R_R(TheComplex(x)->c_imag,true,NULL); /* cos(b), sin(b) */ + /* stack layout: a, b, cos(b), sin(b). */ + /* b != Fixnum_0 ==> sin(b) != Fixnum_0. */ + R_cosh_sinh_R_R(STACK_3,true,NULL); /* cosh(a), sinh(a); cosh(a) != Fixnum 0 */ + /* stack layout: a, b, cos(b), sin(b), cosh(a), sinh(a). */ + STACK_0 = R_R_mal_R(STACK_0,STACK_3); /* sinh(a)*cos(b) */ + { var object erg = R_R_mal_R(STACK_1,STACK_2); /* cosh(a)*sin(b), != Fixnum_0 */ + STACK_5 = R_R_contagion_R(STACK_4,STACK_5); + erg = R_R_complex_C(R_R_float_F(STACK_0,STACK_5), + F_R_float_F(erg,STACK_5)); + skipSTACK(6); return erg; + } } } # N_cosh_N(x) liefert (cosh x), wo x eine Zahl ist. # can trigger GC - local object N_cosh_N (object x); # Methode: # x reell -> klar # x = a+bi -> (complex (* (cosh a) (cos b)) (* (sinh a) (sin b))) - local object N_cosh_N(x) - var object x; +local object N_cosh_N (object x) { if (N_realp(x)) { return R_cosh_R(x); - } else { - # x=a+bi - pushSTACK(TheComplex(x)->c_real); # a retten - R_cos_sin_R_R(TheComplex(x)->c_imag); # cos(b), sin(b) errechnen - # Stackaufbau: a, cos(b), sin(b). - R_cosh_sinh_R_R(STACK_2); # cosh(a), sinh(a) errechnen - # Stackaufbau: a, cos(b), sin(b), cosh(a), sinh(a). - STACK_0 = R_R_mal_R(STACK_0,STACK_2); # sinh(a)*sin(b) - var object erg = R_R_mal_R(STACK_1,STACK_3); # cosh(a)*cos(b) - erg = R_R_complex_N(erg,STACK_0); skipSTACK(5); return erg; + } else { /* x=a+bi */ + pushSTACK(TheComplex(x)->c_real); /* save a */ + pushSTACK(TheComplex(x)->c_imag); /* save b */ + R_cos_sin_R_R(TheComplex(x)->c_imag,true,NULL); /* cos(b), sin(b) */ + /* stack layout: a, b, cos(b), sin(b). */ + R_cosh_sinh_R_R(STACK_3,true,NULL); /* cosh(a), sinh(a) */ + /* stack layout: a, b,cos(b), sin(b), cosh(a), sinh(a). */ + STACK_0 = R_R_mal_R(STACK_0,STACK_2); /* sinh(a)*sin(b) */ + { var object erg = R_R_mal_R(STACK_1,STACK_3); /* cosh(a)*cos(b) */ + STACK_5 = R_R_contagion_R(STACK_4,STACK_5); + erg = eq(STACK_0,Fixnum_0) ? F_R_float_F(erg,STACK_4) + : R_R_complex_C(F_R_float_F(erg,STACK_4), + F_R_float_F(STACK_0,STACK_4)); + skipSTACK(6); return erg; + } } } # N_tanh_N(x) liefert (tanh x), wo x eine Zahl ist. # can trigger GC - local object N_tanh_N (object x); # Methode: # x reell -> (/ (sinh x) (cosh x)) # x = a+bi -> (/ (complex (* (sinh a) (cos b)) (* (cosh a) (sin b))) # (complex (* (cosh a) (cos b)) (* (sinh a) (sin b))) ) - local object N_tanh_N(x) - var object x; +local object N_tanh_N (object x) { if (N_realp(x)) { - R_cosh_sinh_R_R(x); - # Stackaufbau: cosh(x), sinh(x). - var object erg = R_R_durch_R(STACK_0,STACK_1); - skipSTACK(2); return erg; - } else { - # x=a+bi - pushSTACK(TheComplex(x)->c_real); # a retten - R_cos_sin_R_R(TheComplex(x)->c_imag); # cos(b), sin(b) errechnen - # Stackaufbau: a, cos(b), sin(b). - R_cosh_sinh_R_R(STACK_2); # cosh(a), sinh(a) errechnen - # Stackaufbau: a, cos(b), sin(b), cosh(a), sinh(a). - var object temp; - STACK_4 = R_R_mal_R(STACK_0,STACK_3); # sinh(a)*cos(b) - temp = R_R_mal_R(STACK_1,STACK_2); # cosh(a)*sin(b) /= Fixnum 0 - STACK_4 = R_R_complex_C(STACK_4,temp); # Zähler - # Stackaufbau: Zähler, cos(b), sin(b), cosh(a), sinh(a). - STACK_3 = R_R_mal_R(STACK_1,STACK_3); # cosh(a)*cos(b) - temp = R_R_mal_R(STACK_0,STACK_2); # sinh(a)*sin(b) - temp = R_R_complex_N(STACK_3,temp); # Nenner - temp = N_N_durch_N(STACK_4,temp); # Zähler/Nenner - skipSTACK(5); return temp; + pushSTACK(x); + R_cosh_sinh_R_R(x,true,NULL); + /* stack layout: x, cosh(x), sinh(x). */ + { var object erg = F_R_float_F(R_R_durch_R(STACK_0,STACK_1),STACK_2); + skipSTACK(3); return erg; + } + } else { /* x=a+bi */ + pushSTACK(TheComplex(x)->c_real); /* a */ + pushSTACK(TheComplex(x)->c_imag); /* b */ + R_cos_sin_R_R(TheComplex(x)->c_imag,true,NULL); /* cos(b), sin(b) */ + /* stack layout: a, b, cos(b), sin(b). */ + R_cosh_sinh_R_R(STACK_3,true,NULL); /* cosh(a), sinh(a) */ + /* stack layout: a, b, cos(b), sin(b), cosh(a), sinh(a). */ + STACK_5 = R_R_contagion_R(STACK_4,STACK_5); + STACK_4 = R_R_mal_R(STACK_0,STACK_3); /* sinh(a)*cos(b) */ + { var object temp = R_R_mal_R(STACK_1,STACK_2); /* cosh(a)*sin(b) /= Fixnum 0 */ + STACK_4 = R_R_complex_C(STACK_4,temp); /* numerator */ + /* stack layout: contag, numerator, cos(b), sin(b), cosh(a), sinh(a). */ + STACK_3 = R_R_mal_R(STACK_1,STACK_3); /* cosh(a)*cos(b) */ + temp = R_R_mal_R(STACK_0,STACK_2); /* sinh(a)*sin(b) */ + temp = R_R_complex_N(STACK_3,temp); /* denominator */ + STACK_4 = N_N_durch_N(STACK_4,temp); /* numerator/denominator */ + temp = C_R_float_C(&STACK_4,STACK_5); + skipSTACK(6); return temp; + } } } @@ -608,114 +647,113 @@ # rein imaginär ist. # Hilfsfunktion für beide: u+iv := artanh(x+iy), u,v beide auf den Stack. - local void R_R_atanh_R_R (object x, object y); - local void R_R_atanh_R_R(x,y) - var object x; - var object y; +local void R_R_atanh_R_R (object x, object y) { - if (eq(x,Fixnum_0)) { # x=0 -> u=0, v=atan(X=1,Y=y) (Fall y=0 ist inbegriffen) + if (eq(x,Fixnum_0)) { /* x=0 -> u=0, v=atan(X=1,Y=y) (y=0 is included) */ pushSTACK(x); pushSTACK(R_R_atan_R(Fixnum_1,y)); return; } if (eq(y,Fixnum_0)) { if (R_rationalp(x)) - x = RA_float_F(x); # x in Float umwandeln - # x Float - if (R_zerop(x)) { # x=0.0 -> x als Ergebnis + x = RA_float_F(x); /* x --> float */ + /* x -- float */ + if (R_zerop(x)) { /* x=0.0 -> return x */ pushSTACK(x); pushSTACK(Fixnum_0); return; } if (F_exponent_L(x) < 0) { - # Exponent e<0, also |x|<1/2 + /* exponent e<0, ==> |x|<1/2 */ pushSTACK(F_atanhx_F(x)); pushSTACK(Fixnum_0); return; } - # e>=0, also |x|>=1/2 + /* e>=0, ==> |x|>=1/2 */ pushSTACK(x); - pushSTACK(R_R_minus_R(Fixnum_1,x)); # 1-x - # Stackaufbau: x, 1-x. + pushSTACK(R_R_minus_R(Fixnum_1,x)); /* 1-x */ + /* stack layout: x, 1-x. */ var object temp; - temp = R_R_plus_R(Fixnum_1,STACK_1); # 1+x - temp = F_F_durch_F(temp,STACK_0); # (1+x)/(1-x) + temp = R_R_plus_R(Fixnum_1,STACK_1); /* 1+x */ + temp = F_F_durch_F(temp,STACK_0); /* (1+x)/(1-x) */ if (!R_minusp(temp)) { - STACK_1 = temp; STACK_0 = Fixnum_0; # Imaginärteil:=0 - if (R_zerop(temp)) # x = -1 -> Error + STACK_1 = temp; STACK_0 = Fixnum_0; /* imag part :=0 */ + if (R_zerop(temp)) /* x = -1 -> Error */ divide_0(); - } else { - # (1+x)/(1-x) < 0 -> Betrag nehmen, Imaginärteil berechnen: + } else { /* (1+x)/(1-x) < 0 -> negate, compute Im: */ STACK_1 = F_minus_F(temp); - temp = F_I_scale_float_F(pi(STACK_1),Fixnum_minus1); # (scale-float pi -1) = pi/2 - if (R_minusp(STACK_0)) # 1-x<0 -> dann -pi/2 + temp = F_I_scale_float_F(pi(STACK_1),Fixnum_minus1); /* (scale-float pi -1) = pi/2 */ + if (R_minusp(STACK_0)) /* 1-x<0 ==> -pi/2 */ temp = F_minus_F(temp); STACK_0 = temp; } - # Stackaufbau: |(1+x)/(1-x)| (>0), Imaginärteil. - STACK_1 = F_I_scale_float_F(R_ln_R(STACK_1),Fixnum_minus1); # ln bilden, durch 2 + /* stack layout: |(1+x)/(1-x)| (>0), Im. */ + STACK_1 = F_I_scale_float_F(R_ln_R(STACK_1,true,&STACK_1),Fixnum_minus1); /* ln / 2 */ return; } pushSTACK(x); pushSTACK(y); - pushSTACK(R_R_plus_R(Fixnum_1,x)); # 1+x - pushSTACK(R_R_minus_R(Fixnum_1,STACK_2)); # 1-x - # Stackaufbau: x, y, 1+x, 1-x. - # x und y in Floats umwandeln: - if (R_rationalp(STACK_3)) { - if (R_rationalp(STACK_2)) - STACK_2 = RA_float_F(STACK_2); - STACK_3 = RA_F_float_F(STACK_3,STACK_2); - } else { - if (R_rationalp(STACK_2)) - STACK_2 = RA_F_float_F(STACK_2,STACK_3); + /* stack layout: x, y */ + /* x , y --> float: */ + if (R_rationalp(STACK_1)) { + if (R_rationalp(STACK_0)) + STACK_0 = RA_float_F(STACK_0); + STACK_1 = RA_F_float_F(STACK_1,STACK_0); + } else { + if (R_rationalp(STACK_0)) + STACK_0 = RA_F_float_F(STACK_0,STACK_1); + } + pushSTACK(R_R_contagion_R(STACK_0,STACK_1)); + STACK_1 = F_extend2_F(STACK_1); /* increase precision y */ + STACK_2 = F_extend2_F(STACK_2); /* increase precision x */ + pushSTACK(R_R_plus_R(Fixnum_1,STACK_2)); /* 1+x */ + pushSTACK(R_R_minus_R(Fixnum_1,STACK_3)); /* 1-x */ + /* stack layout: x, y, contagion, 1+x, 1-x. */ + pushSTACK(R_square_R(STACK_3)); /* y^2 */ + pushSTACK(R_R_plus_R(Fixnum_1,R_R_plus_R(R_square_R(STACK_4), + STACK_0))); /* 1+x^2+y^2 */ + /* stack layout: x, y, contagion, 1+x, 1-x, y^2, 1+x^2+y^2. */ + { var object temp = F_abs_F(F_I_scale_float_F(STACK_6,fixnum(2))); /* |4x| */ + if (F_F_comp(temp,STACK_0) < 0) { /* |4x| < 1+x^2+y^2 ? */ + temp = F_I_scale_float_F(STACK_6,Fixnum_1); /* 2x */ + temp = F_F_durch_F(temp,STACK_0); /* 2x/(1+x^2+y^2) */ + temp = F_atanhx_F(temp); /* atanh */ + temp = F_I_scale_float_F(temp,Fixnum_minus1); /* .../2 =: u */ + } else { + temp = R_square_R(STACK_3); /* (1+x)^2 */ + STACK_0 = R_R_plus_R(temp,STACK_1); /* (1+x)^2+y^2, a float >=0 */ + temp = R_square_R(STACK_2); /* (1-x)^2 */ + temp = R_R_plus_R(temp,STACK_1); /* (1-x)^2+y^2, a float >=0 */ + temp = F_F_durch_F(STACK_0,temp); /* ((1+x)^2+y^2)/((1-x)^2+y^2), a float >=0 */ + if (R_zerop(temp)) /* should be >0 */ + divide_0(); + temp = R_ln_R(temp,true,NULL); /* ln(temp), a float */ + STACK_6 = F_I_scale_float_F(temp,sfixnum(-2)); /* .../4 =: u */ } - pushSTACK(R_square_R(STACK_2)); # y^2 - { - var object temp = R_square_R(STACK_4); # x^2 - pushSTACK(R_R_plus_R(Fixnum_1,R_R_plus_R(temp,STACK_0))); # 1+x^2+y^2 } - # Stackaufbau: x, y, 1+x, 1-x, y^2, 1+x^2+y^2. - var object temp = # |4x| - F_abs_F(F_I_scale_float_F(STACK_5,fixnum(2))); - if (F_F_comp(temp,STACK_0) < 0) { # |4x| < 1+x^2+y^2 ? - temp = F_I_scale_float_F(STACK_5,Fixnum_1); # 2x - temp = F_F_durch_F(temp,STACK_0); # 2x/(1+x^2+y^2) - temp = F_atanhx_F(temp); # atanh davon - temp = F_I_scale_float_F(temp,Fixnum_minus1); # .../2 =: u - } else { - temp = R_square_R(STACK_3); # (1+x)^2 - STACK_0 = R_R_plus_R(temp,STACK_1); # (1+x)^2+y^2, ein Float >=0 - temp = R_square_R(STACK_2); # (1-x)^2 - temp = R_R_plus_R(temp,STACK_1); # (1-x)^2+y^2, ein Float >=0 - temp = F_F_durch_F(STACK_0,temp); # ((1+x)^2+y^2)/((1-x)^2+y^2), ein Float >=0 - if (R_zerop(temp)) # sollte >0 sein - divide_0(); - temp = R_ln_R(temp); # ln davon, ein Float - temp = F_I_scale_float_F(temp,sfixnum(-2)); # .../4 =: u + { var signean x_sign = R_sign(STACK_5); + var object temp = R_R_mal_R(STACK_3,STACK_2); /* (1+x)(1-x) */ + /* stack layout: u, y, contagion, 1+x, 1-x, y^2, -. */ + STACK_0 = R_R_minus_R(temp,STACK_1); /* (1+x)(1-x)-y^2, a float */ + temp = F_I_scale_float_F(STACK_5,Fixnum_1); /* 2y, a float */ + temp = R_R_atan_R(STACK_0,temp); /* atan(X=(1-x)(1+x)-y^2,Y=2y), a float */ + if (R_minusp(STACK_0) && (x_sign>=0) && R_zerop(STACK_4)) /* X<0.0 and x>=0.0 and Y=0.0 ? */ + temp = F_minus_F(temp); /* change sign */ + STACK_5 = F_I_scale_float_F(temp,Fixnum_minus1); /* .../2 =: v */ + STACK_5 = F_F_float_F(STACK_5,STACK_4); /* restore the precision */ + STACK_6 = F_F_float_F(STACK_6,STACK_4); /* restore the precision */ + /* stack layout: u, v, 1+x, 1-x, y^2, -. */ + skipSTACK(5); return; } - var signean x_sign = R_sign(STACK_5); - STACK_5 = temp; - # Stackaufbau: u, y, 1+x, 1-x, y^2, -. - temp = R_R_mal_R(STACK_3,STACK_2); # (1+x)(1-x) - STACK_0 = R_R_minus_R(temp,STACK_1); # (1+x)(1-x)-y^2, ein Float - temp = F_I_scale_float_F(STACK_4,Fixnum_1); # 2y, ein Float - temp = R_R_atan_R(STACK_0,temp); # atan(X=(1-x)(1+x)-y^2,Y=2y), ein Float - if (R_minusp(STACK_0) && (x_sign>=0) && R_zerop(STACK_4)) # X<0.0 und x>=0.0 und Y=0.0 ? - temp = F_minus_F(temp); # ja -> Vorzeichenwechsel - STACK_4 = F_I_scale_float_F(temp,Fixnum_minus1); # .../2 =: v - # Stackaufbau: u, v, 1+x, 1-x, y^2, -. - skipSTACK(4); return; } - # - local object N_atanh_N(z) - var object z; + +local object N_atanh_N (object z) { if (N_realp(z)) { R_R_atanh_R_R(z,Fixnum_0); } else { R_R_atanh_R_R(TheComplex(z)->c_real,TheComplex(z)->c_imag); } - # Stackaufbau: u, v. - z = R_R_complex_N(STACK_1,STACK_0); skipSTACK(2); return z; + /* stack layout: z, u, v. */ + z = R_R_complex_N(STACK_1,STACK_0); + skipSTACK(2); return z; } - # - local object N_atan_N(z) - var object z; - { # atanh(iz) errechnen: + +local object N_atan_N (object z) +{ /* compute atanh(iz): */ if (N_realp(z)) { R_R_atanh_R_R(Fixnum_0,z); } else { @@ -723,7 +761,7 @@ z = R_minus_R(TheComplex(z)->c_imag); R_R_atanh_R_R(z,popSTACK()); } - # Stackaufbau: u, v. + /* stack layout: z, u, v. */ z = R_minus_R(STACK_1); z = R_R_complex_N(STACK_0,z); # z := v-iu skipSTACK(2); return z; } @@ -829,7 +867,7 @@ || (F_exponent_L(y) <= (sintL)(-F_float_digits(y))>>1) # e <= -d/2 <==> e <= -ceiling(d/2) ? ) return; # u=0, v=y bereits im Stack - # Stackaufbau: 0, y. + # stack layout: 0, y. var object temp = R_R_minus_R(Fixnum_1,F_square_F(y)); # 1-y*y if (!R_minusp(temp)) { # 1-y*y>=0, also |y|<=1 @@ -844,7 +882,7 @@ else temp = F_F_plus_F(temp,y); # temp = sqrt(y^2-1)+|y|, ein Float >1 - STACK_1 = R_ln_R(temp); # ln(|y|+sqrt(y^2-1)), ein Float >0 + STACK_1 = R_ln_R(temp,true,&STACK_0); /* ln(|y|+sqrt(y^2-1)), Float >0 */ temp = F_I_scale_float_F(pi(STACK_1),Fixnum_minus1); # (scale-float pi -1) = pi/2 if (!R_minusp(STACK_0)) { # Vorzeichen von y # y>1 -> v = pi/2 @@ -873,10 +911,10 @@ # |x|>=1/2 if (!R_minusp(x)) # x>=1 - STACK_1 = R_ln_R(F_F_plus_F(temp,x)); # u = ln(x+sqrt(1+x^2)) + STACK_1 = R_ln_R(F_F_plus_F(temp,x),true,&STACK_0); /* u = ln(x+sqrt(1+x^2)) */ else # x<=-1 - STACK_1 = F_minus_F(R_ln_R(F_F_minus_F(temp,x))); # u = -ln(-x+sqrt(1+x^2)) + STACK_1 = F_minus_F(R_ln_R(F_F_minus_F(temp,x),true,&STACK_0)); # u = -ln(-x+sqrt(1+x^2)) } return; } @@ -908,7 +946,7 @@ } else { R_R_asinh_R_R(TheComplex(z)->c_real,TheComplex(z)->c_imag); } - # Stackaufbau: u, v. + # stack layout: u, v. z = R_R_complex_N(STACK_1,STACK_0); skipSTACK(2); return z; } # @@ -923,7 +961,7 @@ z = R_minus_R(TheComplex(z)->c_imag); R_R_asinh_R_R(z,popSTACK()); } - # Stackaufbau: u, v. + # stack layout: u, v. z = R_minus_R(STACK_1); z = R_R_complex_N(STACK_0,z); # z := v-iu skipSTACK(2); return z; } @@ -978,8 +1016,9 @@ var object temp = STACK_0; # z temp = R_R_minus_R(F_square_F(temp),Fixnum_1); # z^2-1, ein Float >=0 temp = F_sqrt_F(temp); # sqrt(z^2-1), ein Float >=0 - temp = F_F_plus_F(popSTACK(),temp); # z+sqrt(z^2-1), ein Float >1 - temp = R_ln_R(temp); # ln(z+sqrt(z^2-1)), ein Float >=0 + temp = F_F_plus_F(STACK_0,temp); /* z+sqrt(z^2-1), float >1 */ + temp = R_ln_R(temp,true,&STACK_0); /* ln(z+sqrt(z^2-1)), float >=0 */ + skipSTACK(1); return R_R_complex_C(Fixnum_0,temp); } R_R_asinh_R_R(Fixnum_0,popSTACK()); @@ -988,7 +1027,7 @@ z = R_minus_R(TheComplex(z)->c_imag); R_R_asinh_R_R(z,popSTACK()); } - # Stackaufbau: u, v. + /* stack layout: u, v. */ # Bilde pi/2-v : z = STACK_0; z = (R_rationalp(z) ? pi(z) : pi_F_float_F(z)); # pi im Float-Format von v @@ -1048,7 +1087,7 @@ STACK_0 = z = RA_float_F(z); # z Float <= -1 z = F_sqrt_F(R_R_minus_R(F_square_F(z),Fixnum_1)); # sqrt(z^2-1), ein Float >=0 - STACK_0 = R_ln_R(F_F_minus_F(z,STACK_0)); # log(sqrt(z^2-1)-z), ein Float >=0 + STACK_0 = R_ln_R(F_F_minus_F(z,STACK_0),true,&STACK_0); # log(sqrt(z^2-1)-z), ein Float >=0 z = pi(STACK_0); # and imaginary part == pi return R_R_complex_C(popSTACK(),z); } ```
 Re[2]: Numeric inconsistency From: Arseny Slobodjuck - 2002-09-11 22:28:10 Attachments: hyper.lsp ```Hi, Wednesday, September 11, 2002, 9:10:23 AM, you wrote: > please try the appended patch instead. I was unable to patch with this file - 3 of 14 hunks were failed. But now (when it in the cvs) there is some tanh test: I computed series for sinh and cosh and calculated tanh from it. Other way is using built-in exp (not the best way since I'm sure it used by tanh internally). Series is also not very good for big argument. That is the table for some real arguments. New version marked '+'. argument (tanh x) series exp - 0 0 0 0 + 0 0.0 0 0 3.1415926535897932385L0 0.9962720762207499443L0 0.99627207622074994407L0 0.9962720762207499443L0 - 6.283185307179586477L0 0.9999930253396106107L0 0.9999930253396106108L0 0.9999930253396106107L0 + 6.283185307179586477L0 0.9999930253396106106L0 0.9999930253396106108L0 0.9999930253396106107L0 - 9.424777960769379716L0 0.9999999869751758126L0 0.9999999869751758125L0 0.9999999869751758126L0 + 9.424777960769379716L0 0.99999998697517581266L0 0.9999999869751758125L0 0.9999999869751758126L0 12.566370614359172954L0 0.9999999999756768866L0 0.9999999999756768865L0 0.9999999999756768866L0 - 15.707963267948966192L0 0.99999999999995457795L0 0.9999999999999545782L0 0.99999999999995457795L0 + 15.707963267948966192L0 0.999999999999954578L0 0.9999999999999545782L0 0.99999999999995457795L0 - 18.849555921538759432L0 0.9999999999999999151L0 0.9999999999999999149L0 0.9999999999999999151L0 + 18.849555921538759432L0 0.99999999999999991516L0 0.9999999999999999149L0 0.9999999999999999151L0 - 21.99114857512855267L0 0.9999999999999999999L0 0.9999999999999999996L0 0.9999999999999999999L0 + 21.99114857512855267L0 0.99999999999999999984L0 0.9999999999999999996L0 0.9999999999999999999L0 25.13274122871834591L0 1.0L0 1.0000000000000000002L0 1.0L0 28.274333882308139149L0 1.0L0 1.0000000000000000001L0 1.0L0 31.415926535897932388L0 1.0L0 1.0000000000000000001L0 1.0L0 34.557519189487725626L0 1.0L0 0.99999999999999999995L0 1.0L0 -- Best regards, Arseny ```
 Re: Re[2]: Numeric inconsistency From: Sam Steingold - 2002-09-11 23:18:11 ```> * In message <1829479520.20020912093105@...> > * On the subject of "Re[2]: Numeric inconsistency" > * Sent on Thu, 12 Sep 2002 09:31:05 +1000 > * Honorable Arseny Slobodjuck writes: > > Wednesday, September 11, 2002, 9:10:23 AM, you wrote: > > > please try the appended patch instead. > But now (when it in the cvs) there is some tanh test: thanks! > I computed series for sinh and cosh and calculated tanh from it. Other > way is using built-in exp (not the best way since I'm sure it used by > tanh internally). actually, the EXP _is_ used internally, _but_ with increased precision. the results are not too bad, although you did miss a couple of crashes :-( -- Sam Steingold (http://www.podval.org/~sds) running RedHat7.3 GNU/Linux ; ; ; ; ; Marriage is the sole cause of divorce. ```
 Re[4]: Numeric inconsistency From: Arseny Slobodjuck - 2002-09-12 03:55:39 ```Hi, Sam, Thursday, September 12, 2002, 9:18:04 AM, you wrote: Sam> the results are not too bad, although you did miss a couple of crashes :-( I'm sure I missed something, but I'd not call it crashes. Do you mean that there is precision loss ? Ah, also loop termination for complex argument... Sam> -- Sam> Marriage is the sole cause of divorce. Why don't you add a tagline like 'Uncle Sam wants YOU for clisp development ?'. Sorry if you don't like the idea. -- Best regards, Arseny mailto:ampy@... ```
 Re: Re[4]: Numeric inconsistency From: Sam Steingold - 2002-09-12 13:00:09 ```Hi Arseny, > * In message <1264573946.20020912145549@...> > * On the subject of "Re[4]: Numeric inconsistency" > * Sent on Thu, 12 Sep 2002 14:55:49 +1000 > * Honorable Arseny Slobodjuck writes: > > Thursday, September 12, 2002, 9:18:04 AM, you wrote: > > Sam> the results are not too bad, although you did miss a couple of crashes :-( > I'm sure I missed something, but I'd not call it crashes. Do you mean > that there is precision loss ? I meant seg faults (fixed now). > Ah, also loop termination for complex argument... huh?! > Sam> Marriage is the sole cause of divorce. > Why don't you add a tagline like 'Uncle Sam wants YOU for clisp > development ?'. Sorry if you don't like the idea. :-) -- Sam Steingold (http://www.podval.org/~sds) running RedHat7.3 GNU/Linux ; ; ; ; ; My other CAR is a CDR. ```
 Re[6]: Numeric inconsistency From: Arseny Slobodjuck - 2002-09-13 04:05:16 ```Hello Sam, Thursday, September 12, 2002, 11:00:00 PM, you wrote: >> Sam> the results are not too bad, although you did miss a couple of crashes :-( >> I'm sure I missed something, but I'd not call it crashes. Do you mean >> that there is precision loss ? > I meant seg faults (fixed now). Got it. I thought you meant crashes in my lisp sources. >> Ah, also loop termination for complex argument... > huh?! Never mind. -- Best regards, Arseny ```