Path: csiph.com!v102.xanadu-bbs.net!xanadu-bbs.net!news.glorb.com!border3.nntp.dca.giganews.com!Xl.tags.giganews.com!border1.nntp.dca.giganews.com!nntp.giganews.com!local2.nntp.dca.giganews.com!nntp.earthlink.com!news.earthlink.com.POSTED!not-for-mail
NNTP-Posting-Date: Fri, 06 Apr 2012 18:20:11 -0500
Message-ID: <4F7F7A2A.3303@mindspring.com>
Date: Fri, 06 Apr 2012 19:20:10 -0400
From: pete <pfiland@mindspring.com>
Reply-To: pfiland@mindspring.com
Organization: PF
X-Mailer: Mozilla 3.04Gold (WinNT; I)
MIME-Version: 1.0
Newsgroups: comp.lang.c
Subject: Re: template metaprogramming in C?
References: <jkkcqp$b94$1@speranza.aioe.org> <jklodd$n3r$1@speranza.aioe.org> <fd4f60c4-a447-47fe-94a6-b003918d88c2@fk28g2000vbb.googlegroups.com> <4f7d7070$0$1387$4fafbaef@reader2.news.tin.it> <4f7d7c00$0$1374$4fafbaef@reader1.news.tin.it> <4f7d82b9$0$1381$4fafbaef@reader2.news.tin.it> <4F7F7185.2958@mindspring.com>
Content-Type: text/plain; charset=us-ascii
Content-Transfer-Encoding: 7bit
Lines: 141
X-Usenet-Provider: http://www.giganews.com
NNTP-Posting-Host: 4.154.220.17
X-Trace: sv3-aBd0Ld454+bEhS2LWY9upGhcKzFyQ2FncRkdPjgsgQjvykSt0LREQqlpN94rL+VvxB+0/0G92EKcLMp!oJlRQdQ4OG8zjme/8Y9xBM1/mLZnPiC3BhMJG3u0omVNuJmvFKNRnAKLtYH0vHyHdNb0Hi9YFqH9!rsQw6kDieg==
X-Abuse-and-DMCA-Info: Please be sure to forward a copy of ALL headers
X-Abuse-and-DMCA-Info: Otherwise we will be unable to process your complaint properly
X-Postfilter: 1.3.40
X-Original-Bytes: 5069
Xref: csiph.com comp.lang.c:19418

pete wrote:
> 
> io_x wrote:
> 
> > u32  rsqrt(double*  r, double  p, u32  ncs)
> 
> > what about result? the last number of each line is the sys_sqrt(p)
> 
> If you want to compare the result of rsqrt(&r, p, ncs) to sqrt(p),
> 
> Then you should compare
>     (p / rsqrt(&r, p, ncs) - rsqrt(&r, p, ncs))
> against
>     (p / sqrt(p) - sqrt(p))
> 
> The expression with the smaller magnitude,
> contains the more accurate square root function.

/* BEGIN new.c output*/

square_root(0.000000) == 0.000000
       sqrt(0.000000) == 0.000000

square_root(1.000000) == 1.000000
       sqrt(1.000000) == 1.000000
x / square_root(x) - square_root(x) == 0
x /        sqrt(x) -        sqrt(x) == 0

square_root(1.001000) == 1.000500
       sqrt(1.001000) == 1.000500
x / square_root(x) - square_root(x) == 0
x /        sqrt(x) -        sqrt(x) == 0

square_root(2.000000) == 1.414214
       sqrt(2.000000) == 1.414214
x / square_root(x) - square_root(x) == 2.22045e-016
x /        sqrt(x) -        sqrt(x) == -2.22045e-016

square_root(10.000000) == 3.162278
       sqrt(10.000000) == 3.162278
x / square_root(x) - square_root(x) == 4.44089e-016
x /        sqrt(x) -        sqrt(x) == -4.44089e-016

square_root(100.000000) == 10.000000
       sqrt(100.000000) == 10.000000
x / square_root(x) - square_root(x) == 0
x /        sqrt(x) -        sqrt(x) == 0

square_root(1000.000000) == 31.622777
       sqrt(1000.000000) == 31.622777
x / square_root(x) - square_root(x) == 0
x /        sqrt(x) -        sqrt(x) == 0

square_root(10000.000000) == 100.000000
       sqrt(10000.000000) == 100.000000
x / square_root(x) - square_root(x) == 0
x /        sqrt(x) -        sqrt(x) == 0

square_root(0.100000) == 0.316228
       sqrt(0.100000) == 0.316228
x / square_root(x) - square_root(x) == 0
x /        sqrt(x) -        sqrt(x) == 0

square_root(0.250000) == 0.500000
       sqrt(0.250000) == 0.500000
x / square_root(x) - square_root(x) == 0
x /        sqrt(x) -        sqrt(x) == 0

square_root(0.001000) == 0.031623
       sqrt(0.001000) == 0.031623
x / square_root(x) - square_root(x) == 6.93889e-018
x /        sqrt(x) -        sqrt(x) == 6.93889e-018

square_root(0.000100) == 0.010000
       sqrt(0.000100) == 0.010000
x / square_root(x) - square_root(x) == 0
x /        sqrt(x) -        sqrt(x) == 0

square_root(0.000010) == 0.003162
       sqrt(0.000010) == 0.003162
x / square_root(x) - square_root(x) == 0
x /        sqrt(x) -        sqrt(x) == 0

/* END new.c output*/




/* BEGIN new.c */

#include <stdio.h>
#include <math.h>

double square_root(double x)
{
    double root = x;

    if (root > 0) {
        double next = root / 2 + 0.5;

        do {
            root = next;
            next = (x / root + root) / 2;
        } while (root > next);
    }
    return root;
}

int 
main(void)
{
    double input[] = {
        0, 1, 1.001, 2, 10, 100, 1000, 10000, 
        0.1, 0.25, 0.001, 0.0001, 0.00001
    };
    size_t x;

    puts("/* BEGIN new.c output*/\n");
    for (x = 0; x != sizeof input / sizeof *input; ++x) {
        printf("square_root(%f) == %f\n", input[x], 
            square_root(input[x]));
        printf("       sqrt(%f) == %f\n", input[x], 
            sqrt(input[x]));
        if (input[x] > 0) {
        printf("x / square_root(x) - square_root(x) == %g\n",
            input[x] / square_root(input[x])
            - square_root(input[x]));
        printf("x /        sqrt(x) -        sqrt(x) == %g\n\n",
            input[x] / sqrt(input[x]) - sqrt(input[x]));
        } else {
            putchar('\n');
        }
    }
    puts("/* END new.c output*/");
    return 0;
}

/* END new.c */

-- 
pete