Issue 55735 - Improve the accuracy of ERF/ERFC for large x value
Summary: Improve the accuracy of ERF/ERFC for large x value
Status: CLOSED FIXED
Alias: None
Product: Calc
Classification: Application
Component: programming (show other issues)
Version: OOo 2.0
Hardware: All All
: P3 Trivial (vote)
Target Milestone: ---
Assignee: oc
QA Contact: issues@sc
URL:
Keywords:
Depends on:
Blocks:
 
Reported: 2005-10-10 18:52 UTC by kyoshida
Modified: 2013-08-07 15:13 UTC (History)
4 users (show)

See Also:
Issue Type: PATCH
Latest Confirmation in: ---
Developer Difficulty: ---


Attachments
implementation of erf(x) and erfc(x) in BASIC (26.97 KB, application/vnd.oasis.opendocument.spreadsheet)
2005-10-19 23:03 UTC, stair
no flags Details
Erf Erfc comparison for -26.5 < x 26.5 (45.20 KB, application/vnd.oasis.opendocument.spreadsheet)
2005-10-30 03:06 UTC, kyoshida
no flags Details
Proposed patch (8.71 KB, patch)
2005-10-30 03:15 UTC, kyoshida
no flags Details | Diff
Make pn/qn static variables (7.91 KB, patch)
2005-11-11 04:59 UTC, kyoshida
no flags Details | Diff
better patch (8.06 KB, patch)
2005-11-11 12:03 UTC, kyoshida
no flags Details | Diff
The variant I used in ooo-build (8.07 KB, patch)
2005-11-11 15:52 UTC, jodygoldberg
no flags Details | Diff

Note You need to log in before you can comment on or make changes to this issue.
Description kyoshida 2005-10-10 18:52:06 UTC
This issue is a spin-off of Issue 19991.

The patch in Issue 19991 fixes incorrect ERF/ERFC values for large x values
where the original algorithm produced erroneous results.  However, the algorithm
should be further improved to match the values computed by octave/gnumeric.  The
revised algorithm in Issue 19991, while it does a better job than the original
algorithm, converges to 1 too quickly.
Comment 1 stair 2005-10-19 23:02:19 UTC
I have a calculation method in the paper 'Algorithms with guarenteed error
bounds for the Error Function and Complementary Error Function - Walter Kramer
and Frithjof Blomquist', read it here
http://www.math.uni-wuppertal.de/org/WRST/preprints/prep_00_2.pdf

Assuming that the results given by Gnumeric are correct then these algorithms
match to better than 8 significant figures for -26.5<x<26.5

Note that there appears to be an error in the published algorithm for x>6 with a
factor of 1/sqrt(pi) missing, in fact even with this correction extending the
algorithm given for 2.2<x<6 to 2.2<x<26.5 gives a much better result at all
points than the published algoritm

The attached spreadsheet includes user functions myerf(x) and myerfc(x)based on
the above and comparison with the Gnumeric results
Comment 2 stair 2005-10-19 23:03:46 UTC
Created attachment 30648 [details]
implementation of erf(x) and erfc(x) in BASIC
Comment 3 kyoshida 2005-10-22 00:06:40 UTC
Oh yeah, this is right on!

I will work on integrating this algorithm into Calc as soon as I get a chance. 
Thanks for finding the algorithm.  I tried previously, but I wasn't successful.

Kohei
Comment 4 kyoshida 2005-10-30 03:06:10 UTC
Created attachment 30998 [details]
Erf Erfc comparison for -26.5 < x 26.5
Comment 5 kyoshida 2005-10-30 03:09:53 UTC
This file contains values of erf & erfc calculated by Calc's new algorithm
(pasted as values) in comparison to Gnumeric's.  They match pretty darn well. 
Many thanks to stair for finding the algorithm reference!

A proposed patch will follow.
Comment 6 kyoshida 2005-10-30 03:15:51 UTC
Created attachment 31000 [details]
Proposed patch
Comment 7 kyoshida 2005-10-30 03:22:49 UTC
Setting target milestone to 2.0.1 and priority set to P3.
Comment 8 kyoshida 2005-10-30 03:23:51 UTC
re-assigning it to Niklas.
Comment 9 kyoshida 2005-10-30 03:25:16 UTC
kohei->nn: Niklas, can we possibly integrate this patch for 2.0.1, or is it too
late?
Comment 10 kyoshida 2005-11-04 14:39:45 UTC
Just to be clear, my patch is a clean implementation from the publication by
Kramer and Blomquist.  The discrepancies in the resulting values as stair
mentioned were not present in the original algorithm. -Kohei
Comment 11 kyoshida 2005-11-09 01:44:36 UTC
Per Niklas' comment on dev@sc.ooo, this issue is re-targeted for 2.0.2 (too late
for 2.0.1), provided that the patch be good enough for integration.  -Kohei
Comment 12 jodygoldberg 2005-11-09 14:52:05 UTC
Nice work I'll put this in ooo-build temporarily.
One small style/performance comment on the patch.
Why use stl vectors for constants ?  You would be better off using
static const double qn[] = { ..... } rather than setting up the constants each
time it's called.  For the erfc0600 case jsut replicate the loop or use pointers
to select which arrays to use.
Comment 13 kyoshida 2005-11-09 16:04:14 UTC
Hi Jody, thanks for your commit and comment.  You're right.  It makes no sense
to define constant for each call.  I'll look into it shortly.  -Kohei
Comment 14 kyoshida 2005-11-11 04:59:16 UTC
Created attachment 31358 [details]
Make pn/qn static variables
Comment 15 kyoshida 2005-11-11 05:01:49 UTC
This patch reflects minor performance improvement per Jody's suggestion.
Comment 16 kyoshida 2005-11-11 12:03:20 UTC
Created attachment 31370 [details]
better patch
Comment 17 kyoshida 2005-11-11 12:07:32 UTC
Perhaps this patch is better?  This one declares the static variables outside
those Erf(c) functions.  The last patch was actually no better than the
original, because the varialbes were still declared in each call.
Comment 18 jodygoldberg 2005-11-11 15:52:24 UTC
Created attachment 31375 [details]
The variant I used in ooo-build
Comment 19 kyoshida 2005-11-11 15:58:03 UTC
Ah.  So those static array declarations are not executed in the subsequent
calls.  I wasn't sure about that.  Thanks!
Comment 20 kyoshida 2005-11-11 17:56:45 UTC
declarations => initializations
I got the terms mixed up.  :P
Comment 21 niklas.nebel 2005-11-15 19:35:50 UTC
I committed the last version (with Jody's changes) into the CWS "dr43".
It will go into OOo 2.0.2.
Comment 22 niklas.nebel 2006-01-06 11:17:06 UTC
back to QA for verification.

re-open issue and reassign to oc@openoffice.org
Comment 23 niklas.nebel 2006-01-06 11:17:16 UTC
reassign to oc@openoffice.org
Comment 24 niklas.nebel 2006-01-06 11:17:21 UTC
reset resolution to FIXED
Comment 25 oc 2006-01-11 12:35:23 UTC
verified in internal build cws_dr43
Comment 26 oc 2006-02-13 11:15:59 UTC
closed because fix available in OOo2.0.2m156