You can subscribe to this list here.
2000 
_{Jan}
(81) 
_{Feb}
(55) 
_{Mar}
(459) 
_{Apr}
(159) 
_{May}
(126) 
_{Jun}
(69) 
_{Jul}
(48) 
_{Aug}
(29) 
_{Sep}
(106) 
_{Oct}
(76) 
_{Nov}
(155) 
_{Dec}
(161) 

2001 
_{Jan}
(122) 
_{Feb}
(150) 
_{Mar}
(294) 
_{Apr}
(124) 
_{May}
(197) 
_{Jun}
(266) 
_{Jul}
(111) 
_{Aug}
(259) 
_{Sep}
(163) 
_{Oct}
(142) 
_{Nov}
(101) 
_{Dec}
(86) 
2002 
_{Jan}
(187) 
_{Feb}
(108) 
_{Mar}
(274) 
_{Apr}
(157) 
_{May}
(346) 
_{Jun}
(242) 
_{Jul}
(345) 
_{Aug}
(187) 
_{Sep}
(263) 
_{Oct}
(69) 
_{Nov}
(30) 
_{Dec}
(76) 
2003 
_{Jan}
(125) 
_{Feb}
(191) 
_{Mar}
(87) 
_{Apr}
(69) 
_{May}
(107) 
_{Jun}
(66) 
_{Jul}
(112) 
_{Aug}
(161) 
_{Sep}
(184) 
_{Oct}
(137) 
_{Nov}
(28) 
_{Dec}
(61) 
2004 
_{Jan}
(148) 
_{Feb}
(99) 
_{Mar}
(365) 
_{Apr}
(225) 
_{May}
(311) 
_{Jun}
(204) 
_{Jul}
(95) 
_{Aug}
(214) 
_{Sep}
(256) 
_{Oct}
(290) 
_{Nov}
(239) 
_{Dec}
(152) 
2005 
_{Jan}
(253) 
_{Feb}
(183) 
_{Mar}
(178) 
_{Apr}
(88) 
_{May}
(175) 
_{Jun}
(195) 
_{Jul}
(122) 
_{Aug}
(81) 
_{Sep}
(119) 
_{Oct}
(200) 
_{Nov}
(110) 
_{Dec}
(179) 
2006 
_{Jan}
(154) 
_{Feb}
(64) 
_{Mar}
(55) 
_{Apr}
(69) 
_{May}
(66) 
_{Jun}
(64) 
_{Jul}
(80) 
_{Aug}
(59) 
_{Sep}
(62) 
_{Oct}
(90) 
_{Nov}
(132) 
_{Dec}
(106) 
2007 
_{Jan}
(58) 
_{Feb}
(51) 
_{Mar}
(59) 
_{Apr}
(19) 
_{May}
(33) 
_{Jun}
(52) 
_{Jul}
(15) 
_{Aug}
(50) 
_{Sep}
(41) 
_{Oct}
(259) 
_{Nov}
(323) 
_{Dec}
(136) 
2008 
_{Jan}
(205) 
_{Feb}
(128) 
_{Mar}
(203) 
_{Apr}
(126) 
_{May}
(307) 
_{Jun}
(166) 
_{Jul}
(259) 
_{Aug}
(181) 
_{Sep}
(217) 
_{Oct}
(265) 
_{Nov}
(256) 
_{Dec}
(132) 
2009 
_{Jan}
(104) 
_{Feb}
(81) 
_{Mar}
(27) 
_{Apr}
(21) 
_{May}
(85) 
_{Jun}
(237) 
_{Jul}
(243) 
_{Aug}
(199) 
_{Sep}
(178) 
_{Oct}
(151) 
_{Nov}
(64) 
_{Dec}
(39) 
2010 
_{Jan}
(33) 
_{Feb}
(146) 
_{Mar}
(125) 
_{Apr}
(109) 
_{May}
(52) 
_{Jun}
(135) 
_{Jul}
(103) 
_{Aug}
(68) 
_{Sep}
(99) 
_{Oct}
(88) 
_{Nov}
(45) 
_{Dec}
(56) 
2011 
_{Jan}
(19) 
_{Feb}
(32) 
_{Mar}
(50) 
_{Apr}
(105) 
_{May}
(46) 
_{Jun}
(22) 
_{Jul}
(101) 
_{Aug}
(80) 
_{Sep}
(52) 
_{Oct}
(16) 
_{Nov}
(10) 
_{Dec}
(29) 
2012 
_{Jan}
(8) 
_{Feb}
(22) 
_{Mar}
(17) 
_{Apr}
(68) 
_{May}
(19) 
_{Jun}
(19) 
_{Jul}
(12) 
_{Aug}
(6) 
_{Sep}
(13) 
_{Oct}
(5) 
_{Nov}
(5) 
_{Dec}
(5) 
2013 
_{Jan}
(6) 
_{Feb}
(4) 
_{Mar}
(3) 
_{Apr}
(5) 
_{May}

_{Jun}

_{Jul}

_{Aug}

_{Sep}

_{Oct}

_{Nov}
(6) 
_{Dec}

2014 
_{Jan}

_{Feb}

_{Mar}
(16) 
_{Apr}
(1) 
_{May}
(8) 
_{Jun}

_{Jul}
(1) 
_{Aug}
(1) 
_{Sep}

_{Oct}

_{Nov}

_{Dec}

2015 
_{Jan}

_{Feb}
(8) 
_{Mar}
(23) 
_{Apr}
(5) 
_{May}

_{Jun}

_{Jul}

_{Aug}
(7) 
_{Sep}
(1) 
_{Oct}

_{Nov}

_{Dec}
(5) 
S  M  T  W  T  F  S 



1
(14) 
2
(16) 
3
(7) 
4
(4) 
5
(11) 
6
(9) 
7
(14) 
8
(4) 
9
(1) 
10
(1) 
11
(2) 
12
(2) 
13
(1) 
14
(4) 
15
(1) 
16
(1) 
17
(2) 
18
(9) 
19

20

21
(3) 
22
(3) 
23
(3) 
24

25
(1) 
26
(2) 
27

28
(4) 
29
(2) 
30
(5) 



From: SourceForge.net <noreply@so...>  20080408 23:00:59

Bugs item #1934968, was opened at 20080404 22:51 Message generated for change (Comment added) made by sds You can respond by visiting: https://sourceforge.net/tracker/?func=detail&atid=101355&aid=1934968&group_id=1355 Please note that this message will contain a full copy of the comment thread, including the initial issue submission, for this request, not just the latest update. >Category: clisp >Group: lisp error >Status: Closed >Resolution: Fixed Priority: 5 Private: No Submitted By: Raymond Toy (rtoy) >Assigned to: Sam Steingold (sds) Summary: log(complex) sometimes inaccurate Initial Comment: (log #c(1.000000001d0 1d5)) > #C(1.0500000863261394d9 9.999999989666666d6) But consider: (setf (longfloatdigits) 256) (log (coerce #c(1.000000001d0 1d5) '(complex longfloat))) > #C(1.05000008213787091674976747623722615763692379335589165541504453172940877702943L9 9.9999999896666666683134966416887953589434016059498897741294881863894719391568L6) so the real part is a bit off.  Comment By: Sam Steingold (sds) Date: 20080408 18:42 Message: Logged In: YES user_id=5735 Originator: NO thank you for your bug report. the bug has been fixed in the CVS tree. you can either wait for the next release (recommended) or check out the current CVS tree (see http://clisp.cons.org) and build CLISP from the sources (be advised that between releases the CVS tree is very unstable and may not even build on your platform).  Comment By: Raymond Toy (rtoy) Date: 20080408 09:09 Message: Logged In: YES user_id=28849 Originator: YES Of course there's room for improvement. CMUCL computes a better result using just doublefloats. No extended precision. The basic algorithm was given previously. All you need is a log1p(x) = log(1+x) routine. You already have lnx. I think it uses a series, so a slight modification of it would work as a log1p function. Then you would have an accurate answer. And for this application you only need log1p for small x, so it doesn't have to be fully general. (But log1p is useful in it's own right so having one for all x > 1 would be beneficial. An expm1(x) = exp(x)  1 would be good too.)  Comment By: Raymond Toy (rtoy) Date: 20080408 09:09 Message: Logged In: YES user_id=28849 Originator: YES Of course there's room for improvement. CMUCL computes a better result using just doublefloats. No extended precision. The basic algorithm was given previously. All you need is a log1p(x) = log(1+x) routine. You already have lnx. I think it uses a series, so a slight modification of it would work as a log1p function. Then you would have an accurate answer. And for this application you only need log1p for small x, so it doesn't have to be fully general. (But log1p is useful in it's own right so having one for all x > 1 would be beneficial. An expm1(x) = exp(x)  1 would be good too.)  Comment By: Sam Steingold (sds) Date: 20080407 17:15 Message: Logged In: YES user_id=5735 Originator: NO > (checklogabs 1d0 1d8) 3.657306783374019d4 > (realpart (log (coerce #c(1d0 1d8) '(complex longfloat)))) 4.9999999999997260093L17 > (realpart (log #c(1d0 1d8))) 4.9981720151581754d17  Comment By: Sam Steingold (sds) Date: 20080407 17:03 Message: Logged In: YES user_id=5735 Originator: NO (defun checklogabs (x y) (let ((r (realpart (log (complex x y))))) (/ ( r (float (realpart (log (complex (float x 0l0) (float y 0l0)))) x)) (abs r)))) (z 1.000000001d0 1d5) 4.027036309576161d9 ; before 4.7058695262211d12 ; after so, the improvement is 3 decimal places, but we are still 4 decimal places behind doublefloatepsilon  Comment By: Sam Steingold (sds) Date: 20080407 16:49 Message: Logged In: YES user_id=5735 Originator: NO I committed a patch which improves accuracy by 2.5 decimal places. [5]> (log #c(1.000000001d0 1d5)) #C(1.050000082142812d9 9.999999989666666d6) [6]> (log (coerce #c(1.000000001d0 1d5) '(complex longfloat))) #C(1.0500000821378709167L9 9.9999999896666666686L6) still no cigar, but I don't see any space for improvement here anymore.  Comment By: Raymond Toy (rtoy) Date: 20080407 14:56 Message: Logged In: YES user_id=28849 Originator: YES That's because sbcl uses cmucl's implementation (I think), which, in turn, comes from an article by Kahan. Nothing says you should use naive implementations of the special functions. :) Clisp's sinh(x) isn't consistent either: (= (sinh x) (* 1/2 ( (exp x) (exp ( x))))) returns NIL for small values of x. :)  Comment By: Sam Steingold (sds) Date: 20080407 13:46 Message: Logged In: YES user_id=5735 Originator: NO at least clisp is consistent: (= (realpart (log #c(1.000000001d0 1d5))) (log (abs #c(1.000000001d0 1d5)))) ==> T as opposed to sbcl which returns NIL. :)  Comment By: Raymond Toy (rtoy) Date: 20080407 13:36 Message: Logged In: YES user_id=28849 Originator: YES I think this is only a problem if the realpart is near 1 and the imagpart is small (or vice versa). The computation of log(abs(z)) loses precision. It might be better computed as 1/2*log(x^2+y^2) = 1/2*log(1+((x1)*(x+1)+y^2), using a special log(1+x) function for small x.  You can respond by visiting: https://sourceforge.net/tracker/?func=detail&atid=101355&aid=1934968&group_id=1355 
From: <clispcvsrequest@li...>  20080408 19:10:36

Send clispcvs mailing list submissions to clispcvs@... To subscribe or unsubscribe via the World Wide Web, visit https://lists.sourceforge.net/lists/listinfo/clispcvs or, via email, send a message with subject or body 'help' to clispcvsrequest@... You can reach the person managing the list at clispcvsowner@... When replying, please edit your Subject line so it is more specific than "Re: Contents of clispcvs digest..." CLISP CVS commits for today Today's Topics: 1. clisp/tests number2.tst,1.44,1.45 ChangeLog,1.539,1.540 (Sam Steingold) 2. clisp/src comptran.d,1.47,1.48 ChangeLog,1.6115,1.6116 (Sam Steingold)  Message: 1 Date: Mon, 07 Apr 2008 20:46:34 +0000 From: Sam Steingold <sds@...> Subject: clisp/tests number2.tst,1.44,1.45 ChangeLog,1.539,1.540 To: clispcvs@... MessageID: <E1JiyEc0002vrUj@...> Update of /cvsroot/clisp/clisp/tests In directory sc8prcvs4.sourceforge.net:/tmp/cvsserv15101/tests Modified Files: number2.tst ChangeLog Log Message: partially fix bug #[ 1934968 ]: log(complex) sometimes inaccurate (N_log_abs_R): new function: compute (log (abs z)) without using hypot, thus avoiding a square root (N_log_N, N_N_log_N): use it Index: number2.tst =================================================================== RCS file: /cvsroot/clisp/clisp/tests/number2.tst,v retrieving revision 1.44 retrieving revision 1.45 diff u d r1.44 r1.45  number2.tst 28 Mar 2008 03:02:54 0000 1.44 +++ number2.tst 7 Apr 2008 20:46:32 0000 1.45 @@ 71,6 +71,8 @@ (log 8d0 2d0) #C(1.0928406470908163d0 0.4207872484158604d0) (log z) #C(2.3025850929940455d0 0d0) z #C(1d1 0d0) +(log #c(1 0)) 0 +(log #c(0 1)) #C(0 1.5707964) (cis 10) #c(0.8390715 0.5440211) (cis 123) #c(0.8879689 0.45990348) @@ 393,9 +395,9 @@ (ACOSH #C(1d0 2d0)) #C(1.5285709194809982d0 1.1437177404024204d0) (ACOSH #C(1L0 2L0)) #C(1.528570919480998161L0 1.1437177404024204937L0) (LOG #C(1s0 2s0)) #C(0.80471s0 1.10715s0) +(LOG #C(1s0 2s0)) #C(0.80472s0 1.10715s0) (LOG #C(1f0 2f0)) #C(0.804719f0 1.1071488f0) (LOG #C(1d0 2d0)) #C(0.8047189562170503d0 1.1071487177940904d0) +(LOG #C(1d0 2d0)) #C(0.8047189562170501d0 1.1071487177940904d0) (LOG #C(1L0 2L0)) #C(0.8047189562170501873L0 1.107148717794090503L0) (ATANH #C(1s0 2s0)) #C(0.173286s0 1.1781s0) Index: ChangeLog =================================================================== RCS file: /cvsroot/clisp/clisp/tests/ChangeLog,v retrieving revision 1.539 retrieving revision 1.540 diff u d r1.539 r1.540  ChangeLog 6 Apr 2008 21:11:32 0000 1.539 +++ ChangeLog 7 Apr 2008 20:46:32 0000 1.540 @@ 1,3 +1,7 @@ +20080407 Sam Steingold <sds@...> + + * number2.tst: updated for improved accuracy of log(complex) + 20080406 Sam Steingold <sds@...> * setf.tst: check that macro doc strings are preserved by compile  Message: 2 Date: Mon, 07 Apr 2008 20:46:37 +0000 From: Sam Steingold <sds@...> Subject: clisp/src comptran.d,1.47,1.48 ChangeLog,1.6115,1.6116 To: clispcvs@... MessageID: <E1JiyEh0007qVTt@...> Update of /cvsroot/clisp/clisp/src In directory sc8prcvs4.sourceforge.net:/tmp/cvsserv15101/src Modified Files: comptran.d ChangeLog Log Message: partially fix bug #[ 1934968 ]: log(complex) sometimes inaccurate (N_log_abs_R): new function: compute (log (abs z)) without using hypot, thus avoiding a square root (N_log_N, N_N_log_N): use it Index: comptran.d =================================================================== RCS file: /cvsroot/clisp/clisp/src/comptran.d,v retrieving revision 1.47 retrieving revision 1.48 diff u d r1.47 r1.48  comptran.d 27 Jan 2008 20:05:05 0000 1.47 +++ comptran.d 7 Apr 2008 20:46:32 0000 1.48 @@ 84,33 +84,55 @@ } } /* N_log_N(x,&end_precision) returns (log x), being x a number. +/* N_log_abs_R(z,&end_precision) returns (log (abs z)), z being a number + can trigger GC + Method: (/ (log (+ x^2 y^2)) 2) */ +local maygc object N_log_abs_R (object z, gcv_object_t *end_p) { + if (N_realp(z)) { + z = N_abs_R(z); + if (R_zerop(z)) + divide_0(); + return R_ln_R(z,end_p); + } else { + pushSTACK(z); + var object tmp = R_square_R(TheComplex(z)>c_real); pushSTACK(tmp); /*x^2*/ + tmp = R_square_R(TheComplex(STACK_1)>c_imag); /* y^2 */ + tmp = R_R_plus_R(tmp,STACK_0); /* x^2 + y^2 */ + if (R_zerop(tmp)) + divide_0(); + tmp = R_ln_R(tmp,end_p); /* log(x^2 + y^2) */ + tmp = floatp(tmp) /* log(x^2 + y^2)/2 */ + ? F_I_scale_float_F(tmp,Fixnum_minus1) + : RA_RA_div_RA(tmp,fixnum(2)); + skipSTACK(2); + return tmp; + } +} + +/* N_log_N(x,&end_precision) returns (log x), x being a number. can trigger GC Method: (complex (log (abs x)) (phase x)) */ local maygc object N_log_N (object x, gcv_object_t *end_p) { pushSTACK(x); /* save x */  pushSTACK(N_abs_R(x)); /* (abs x) */  if (R_zerop(STACK_0)) /* (abs x) = 0 > Error */  divide_0();  STACK_0 = R_ln_R(STACK_0,end_p); /* (log (abs x)) */ /* Increase precision: */  if (floatp(STACK_1))  STACK_1 = F_extend_F(STACK_1);  else if (complexp(STACK_1)  && (floatp(TheComplex(STACK_1)>c_real)   floatp(TheComplex(STACK_1)>c_imag))) {  var object realpart = TheComplex(STACK_1)>c_real; + if (floatp(STACK_0)) + STACK_0 = F_extend_F(STACK_0); + else if (complexp(STACK_0) + && (floatp(TheComplex(STACK_0)>c_real) +  floatp(TheComplex(STACK_0)>c_imag))) { + var object realpart = TheComplex(STACK_0)>c_real; if (floatp(realpart)) realpart = F_extend_F(realpart); pushSTACK(realpart);  var object imagpart = TheComplex(STACK_(1+1))>c_imag; + var object imagpart = TheComplex(STACK_(0+1))>c_imag; if (floatp(imagpart)) imagpart = F_extend_F(imagpart); realpart = popSTACK();  STACK_1 = R_R_complex_C(realpart,imagpart); + STACK_0 = R_R_complex_C(realpart,imagpart); } + pushSTACK(N_log_abs_R(STACK_0,end_p)); /* (log (abs x)) */ STACK_1 = N_phase_R(STACK_1,true); /* (phase x) */ if (end_p != NULL && floatp(STACK_1)) STACK_1 = F_R_float_F(STACK_1,*end_p); @@ 190,7 +212,7 @@ } } /* no chance for rational realpart */  pushSTACK(F_ln_F(N_abs_R(a),&STACK_3)); /* (log (abs a)), a float */ + pushSTACK(N_log_abs_R(a,&STACK_3)); /* (log (abs a)), a float */ /* divide by (log b), return the realpart: */ b = STACK_2; if (R_rationalp(b)) Index: ChangeLog =================================================================== RCS file: /cvsroot/clisp/clisp/src/ChangeLog,v retrieving revision 1.6115 retrieving revision 1.6116 diff u d r1.6115 r1.6116  ChangeLog 6 Apr 2008 21:11:31 0000 1.6115 +++ ChangeLog 7 Apr 2008 20:46:32 0000 1.6116 @@ 1,3 +1,10 @@ +20080407 Sam Steingold <sds@...> + + partially fix bug #[ 1934968 ]: log(complex) sometimes inaccurate + * comptran.d (N_log_abs_R): new function: compute (log (abs z)) + without using hypot, thus avoiding a square root + (N_log_N, N_N_log_N): use it + 20080406 Sam Steingold <sds@...> fix bug #[ 1936255 ]: COMPILE discards MACRO doc strings   This SF.net email is sponsored by the 2008 JavaOne(SM) Conference Register now and save $200. Hurry, offer ends at 11:59 p.m., Monday, April 7! Use priority code J8TLD2. http://ad.doubleclick.net/clk;198757673;13503038;p?http://java.sun.com/javaone  _______________________________________________ clispcvs mailing list clispcvs@... https://lists.sourceforge.net/lists/listinfo/clispcvs End of clispcvs Digest, Vol 24, Issue 7 **************************************** 
From: SourceForge.net <noreply@so...>  20080408 13:09:36

Bugs item #1934968, was opened at 20080404 22:51 Message generated for change (Comment added) made by rtoy You can respond by visiting: https://sourceforge.net/tracker/?func=detail&atid=101355&aid=1934968&group_id=1355 Please note that this message will contain a full copy of the comment thread, including the initial issue submission, for this request, not just the latest update. Category: None Group: None Status: Open Resolution: None Priority: 5 Private: No Submitted By: Raymond Toy (rtoy) Assigned to: Nobody/Anonymous (nobody) Summary: log(complex) sometimes inaccurate Initial Comment: (log #c(1.000000001d0 1d5)) > #C(1.0500000863261394d9 9.999999989666666d6) But consider: (setf (longfloatdigits) 256) (log (coerce #c(1.000000001d0 1d5) '(complex longfloat))) > #C(1.05000008213787091674976747623722615763692379335589165541504453172940877702943L9 9.9999999896666666683134966416887953589434016059498897741294881863894719391568L6) so the real part is a bit off.  >Comment By: Raymond Toy (rtoy) Date: 20080408 09:09 Message: Logged In: YES user_id=28849 Originator: YES Of course there's room for improvement. CMUCL computes a better result using just doublefloats. No extended precision. The basic algorithm was given previously. All you need is a log1p(x) = log(1+x) routine. You already have lnx. I think it uses a series, so a slight modification of it would work as a log1p function. Then you would have an accurate answer. And for this application you only need log1p for small x, so it doesn't have to be fully general. (But log1p is useful in it's own right so having one for all x > 1 would be beneficial. An expm1(x) = exp(x)  1 would be good too.)  Comment By: Raymond Toy (rtoy) Date: 20080408 09:09 Message: Logged In: YES user_id=28849 Originator: YES Of course there's room for improvement. CMUCL computes a better result using just doublefloats. No extended precision. The basic algorithm was given previously. All you need is a log1p(x) = log(1+x) routine. You already have lnx. I think it uses a series, so a slight modification of it would work as a log1p function. Then you would have an accurate answer. And for this application you only need log1p for small x, so it doesn't have to be fully general. (But log1p is useful in it's own right so having one for all x > 1 would be beneficial. An expm1(x) = exp(x)  1 would be good too.)  Comment By: Sam Steingold (sds) Date: 20080407 17:15 Message: Logged In: YES user_id=5735 Originator: NO > (checklogabs 1d0 1d8) 3.657306783374019d4 > (realpart (log (coerce #c(1d0 1d8) '(complex longfloat)))) 4.9999999999997260093L17 > (realpart (log #c(1d0 1d8))) 4.9981720151581754d17  Comment By: Sam Steingold (sds) Date: 20080407 17:03 Message: Logged In: YES user_id=5735 Originator: NO (defun checklogabs (x y) (let ((r (realpart (log (complex x y))))) (/ ( r (float (realpart (log (complex (float x 0l0) (float y 0l0)))) x)) (abs r)))) (z 1.000000001d0 1d5) 4.027036309576161d9 ; before 4.7058695262211d12 ; after so, the improvement is 3 decimal places, but we are still 4 decimal places behind doublefloatepsilon  Comment By: Sam Steingold (sds) Date: 20080407 16:49 Message: Logged In: YES user_id=5735 Originator: NO I committed a patch which improves accuracy by 2.5 decimal places. [5]> (log #c(1.000000001d0 1d5)) #C(1.050000082142812d9 9.999999989666666d6) [6]> (log (coerce #c(1.000000001d0 1d5) '(complex longfloat))) #C(1.0500000821378709167L9 9.9999999896666666686L6) still no cigar, but I don't see any space for improvement here anymore.  Comment By: Raymond Toy (rtoy) Date: 20080407 14:56 Message: Logged In: YES user_id=28849 Originator: YES That's because sbcl uses cmucl's implementation (I think), which, in turn, comes from an article by Kahan. Nothing says you should use naive implementations of the special functions. :) Clisp's sinh(x) isn't consistent either: (= (sinh x) (* 1/2 ( (exp x) (exp ( x))))) returns NIL for small values of x. :)  Comment By: Sam Steingold (sds) Date: 20080407 13:46 Message: Logged In: YES user_id=5735 Originator: NO at least clisp is consistent: (= (realpart (log #c(1.000000001d0 1d5))) (log (abs #c(1.000000001d0 1d5)))) ==> T as opposed to sbcl which returns NIL. :)  Comment By: Raymond Toy (rtoy) Date: 20080407 13:36 Message: Logged In: YES user_id=28849 Originator: YES I think this is only a problem if the realpart is near 1 and the imagpart is small (or vice versa). The computation of log(abs(z)) loses precision. It might be better computed as 1/2*log(x^2+y^2) = 1/2*log(1+((x1)*(x+1)+y^2), using a special log(1+x) function for small x.  You can respond by visiting: https://sourceforge.net/tracker/?func=detail&atid=101355&aid=1934968&group_id=1355 
From: SourceForge.net <noreply@so...>  20080408 13:09:31

Bugs item #1934968, was opened at 20080404 22:51 Message generated for change (Comment added) made by rtoy You can respond by visiting: https://sourceforge.net/tracker/?func=detail&atid=101355&aid=1934968&group_id=1355 Please note that this message will contain a full copy of the comment thread, including the initial issue submission, for this request, not just the latest update. Category: None Group: None Status: Open Resolution: None Priority: 5 Private: No Submitted By: Raymond Toy (rtoy) Assigned to: Nobody/Anonymous (nobody) Summary: log(complex) sometimes inaccurate Initial Comment: (log #c(1.000000001d0 1d5)) > #C(1.0500000863261394d9 9.999999989666666d6) But consider: (setf (longfloatdigits) 256) (log (coerce #c(1.000000001d0 1d5) '(complex longfloat))) > #C(1.05000008213787091674976747623722615763692379335589165541504453172940877702943L9 9.9999999896666666683134966416887953589434016059498897741294881863894719391568L6) so the real part is a bit off.  >Comment By: Raymond Toy (rtoy) Date: 20080408 09:09 Message: Logged In: YES user_id=28849 Originator: YES Of course there's room for improvement. CMUCL computes a better result using just doublefloats. No extended precision. The basic algorithm was given previously. All you need is a log1p(x) = log(1+x) routine. You already have lnx. I think it uses a series, so a slight modification of it would work as a log1p function. Then you would have an accurate answer. And for this application you only need log1p for small x, so it doesn't have to be fully general. (But log1p is useful in it's own right so having one for all x > 1 would be beneficial. An expm1(x) = exp(x)  1 would be good too.)  Comment By: Sam Steingold (sds) Date: 20080407 17:15 Message: Logged In: YES user_id=5735 Originator: NO > (checklogabs 1d0 1d8) 3.657306783374019d4 > (realpart (log (coerce #c(1d0 1d8) '(complex longfloat)))) 4.9999999999997260093L17 > (realpart (log #c(1d0 1d8))) 4.9981720151581754d17  Comment By: Sam Steingold (sds) Date: 20080407 17:03 Message: Logged In: YES user_id=5735 Originator: NO (defun checklogabs (x y) (let ((r (realpart (log (complex x y))))) (/ ( r (float (realpart (log (complex (float x 0l0) (float y 0l0)))) x)) (abs r)))) (z 1.000000001d0 1d5) 4.027036309576161d9 ; before 4.7058695262211d12 ; after so, the improvement is 3 decimal places, but we are still 4 decimal places behind doublefloatepsilon  Comment By: Sam Steingold (sds) Date: 20080407 16:49 Message: Logged In: YES user_id=5735 Originator: NO I committed a patch which improves accuracy by 2.5 decimal places. [5]> (log #c(1.000000001d0 1d5)) #C(1.050000082142812d9 9.999999989666666d6) [6]> (log (coerce #c(1.000000001d0 1d5) '(complex longfloat))) #C(1.0500000821378709167L9 9.9999999896666666686L6) still no cigar, but I don't see any space for improvement here anymore.  Comment By: Raymond Toy (rtoy) Date: 20080407 14:56 Message: Logged In: YES user_id=28849 Originator: YES That's because sbcl uses cmucl's implementation (I think), which, in turn, comes from an article by Kahan. Nothing says you should use naive implementations of the special functions. :) Clisp's sinh(x) isn't consistent either: (= (sinh x) (* 1/2 ( (exp x) (exp ( x))))) returns NIL for small values of x. :)  Comment By: Sam Steingold (sds) Date: 20080407 13:46 Message: Logged In: YES user_id=5735 Originator: NO at least clisp is consistent: (= (realpart (log #c(1.000000001d0 1d5))) (log (abs #c(1.000000001d0 1d5)))) ==> T as opposed to sbcl which returns NIL. :)  Comment By: Raymond Toy (rtoy) Date: 20080407 13:36 Message: Logged In: YES user_id=28849 Originator: YES I think this is only a problem if the realpart is near 1 and the imagpart is small (or vice versa). The computation of log(abs(z)) loses precision. It might be better computed as 1/2*log(x^2+y^2) = 1/2*log(1+((x1)*(x+1)+y^2), using a special log(1+x) function for small x.  You can respond by visiting: https://sourceforge.net/tracker/?func=detail&atid=101355&aid=1934968&group_id=1355 