Computing x^y, sort-of quick, somewhat dirty

On Mon, 26 Dec 2005 01:11:39 -0500, Glenn Zelniker wrote:

> Bob Monsen wrote:
>> 
>> How about the generalized binomial theorem? It'll be crappy as x
>> -> 0, but if you can live with that, it is simple to compute. The radius
>> of convergence (as shown) for non-natural n is 0 < x < 2.
>> 
>> double binomial(double x, double n, int terms)
>> {
>>     double y = 1.0;
>>     double z = 1.0;
>>     double i = 1.0;
>>  
>>     x = 1.0 - x;  // It really computes (1 + y)^n
>> 
>>     while (terms-- > 0)
>>     {
>>         z *= -x*n / i;
>>         n -= 1.0;
>>         i += 1.0;
>>         
>>         y += z;
>>     }        
>>         
>>     return y;
>> }
> 
> This is exactly the kind of thing I was looking for. Thank 
> you, thank you. It'll only take a few minutes to code up 
> when I get back to the lab, but I thought I'd jump the gun 
> and ask a leading question. For x not too close to zero, how 
> many iterations is typical to get, say, ten or so bits of 
> accuracy in the result?

I don't know. However, you can easily run a test program on a PC to find
out, at least with double precision on a PC (which is what I did when I
blasted this out the other night). In fact, here is the output from my
trivial test program, massaged to keep error, defined as (x^n -
binomal())/(x^n) = (e) under 0.001. c is the count of iterations required
to do this (at least when compared against the linux pow() function).

0.01^0=1	c=0	e=0
0.01^0.1=0.631582	c=304	e=0.000624887
0.01^0.2=0.398502	c=340	e=0.000395107
0.01^0.3=0.251437	c=352	e=0.00024846
0.01^0.4=0.158647	c=352	e=0.000157706
0.01^0.5=0.1001	c=345	e=9.95747e-05
0.01^0.6=0.0631585	c=332	e=6.27259e-05
0.01^0.7=0.03985	c=313	e=3.92514e-05
0.01^0.8=0.0251439	c=284	e=2.50121e-05
0.01^0.9=0.0158648	c=238	e=1.58481e-05
0.11^0=1	c=0	e=0
0.11^0.1=0.802706	c=26	e=0.00077064
0.11^0.2=0.643734	c=29	e=0.000633494
0.11^0.3=0.516237	c=30	e=0.000512505
0.11^0.4=0.413934	c=31	e=0.000356675
0.11^0.5=0.331968	c=30	e=0.000305913
0.11^0.6=0.266213	c=29	e=0.0002414
0.11^0.7=0.213497	c=27	e=0.00020497
0.11^0.8=0.171199	c=25	e=0.000152972
0.11^0.9=0.137292	c=21	e=0.00012396
0.21^0=1	c=0	e=0
0.21^0.1=0.856235	c=13	e=0.000731476
0.21^0.2=0.73244	c=15	e=0.000552894
0.21^0.3=0.626689	c=15	e=0.000557423
0.21^0.4=0.536144	c=15	e=0.000485893
0.21^0.5=0.458641	c=15	e=0.000383062
0.21^0.6=0.392424	c=14	e=0.000382881
0.21^0.7=0.335643	c=14	e=0.000250614
0.21^0.8=0.287215	c=12	e=0.000285108
0.21^0.9=0.245642	c=11	e=0.00017262
0.31^0=1	c=0	e=0
0.31^0.1=0.890271	c=8	e=0.000790612
0.31^0.2=0.791887	c=9	e=0.000712535
0.31^0.3=0.704201	c=10	e=0.000466796
0.31^0.4=0.626382	c=10	e=0.000424319
0.31^0.5=0.557125	c=10	e=0.000348943
0.31^0.6=0.495678	c=9	e=0.000436995
0.31^0.7=0.440807	c=9	e=0.000299584
0.31^0.8=0.392128	c=8	e=0.000305788
0.31^0.9=0.348756	c=7	e=0.000238144
0.41^0=1	c=0	e=0
0.41^0.1=0.915271	c=6	e=0.00057166
0.41^0.2=0.837114	c=7	e=0.00043884
0.41^0.3=0.765784	c=7	e=0.000478193
0.41^0.4=0.700476	c=7	e=0.000450814
0.41^0.5=0.640697	c=7	e=0.00038464
0.41^0.6=0.585993	c=7	e=0.000300093
0.41^0.7=0.536178	c=6	e=0.000444724
0.41^0.8=0.490309	c=6	e=0.000274357
0.41^0.9=0.448513	c=5	e=0.000278416
0.51^0=1	c=0	e=0
0.51^0.1=0.935651	c=4	e=0.000768879
0.51^0.2=0.874487	c=5	e=0.000481308
0.51^0.3=0.817635	c=5	e=0.000542395
0.51^0.4=0.764414	c=5	e=0.000529082
0.51^0.5=0.71461	c=5	e=0.000467326
0.51^0.6=0.668017	c=5	e=0.000377657
0.51^0.7=0.62444	c=5	e=0.000275912
0.51^0.8=0.584013	c=4	e=0.000492096
0.51^0.9=0.545755	c=4	e=0.000231799
0.61^0=1	c=0	e=0
0.61^0.1=0.952465	c=3	e=0.000692779
0.61^0.2=0.906207	c=4	e=0.000337185
0.61^0.3=0.862571	c=4	e=0.000388886
0.61^0.4=0.820989	c=4	e=0.000388396
0.61^0.5=0.781376	c=4	e=0.000351407
0.61^0.6=0.743649	c=4	e=0.000291023
0.61^0.7=0.707725	c=4	e=0.000217998
0.61^0.8=0.673934	c=3	e=0.000548125
0.61^0.9=0.641177	c=3	e=0.000267026
0.71^0=1	c=0	e=0
0.71^0.1=0.967216	c=2	e=0.000884672
0.71^0.2=0.934101	c=3	e=0.000306058
0.71^0.3=0.902718	c=3	e=0.000363198
0.71^0.4=0.872347	c=3	e=0.000373498
0.71^0.5=0.842963	c=3	e=0.00034821
0.71^0.6=0.814542	c=3	e=0.000297387
0.71^0.7=0.78706	c=3	e=0.00022992
0.71^0.8=0.760492	c=3	e=0.000153582
0.71^0.9=0.735216	c=2	e=0.00047748
0.81^0=1	c=0	e=0
0.81^0.1=0.979375	c=2	e=0.000227138
0.81^0.2=0.959112	c=2	e=0.000380484
0.81^0.3=0.939209	c=2	e=0.000469107
0.81^0.4=0.919668	c=2	e=0.000501881
0.81^0.5=0.900487	c=2	e=0.0004875
0.81^0.6=0.881668	c=2	e=0.000434474
0.81^0.7=0.863209	c=2	e=0.000351136
0.81^0.8=0.845112	c=2	e=0.000245646
0.81^0.9=0.827375	c=2	e=0.000125993
0.91^0=1	c=0	e=0
0.91^0.1=0.991	c=1	e=0.000386735
0.91^0.2=0.982	c=1	e=0.000685359
0.91^0.3=0.973	c=1	e=0.000896699
0.91^0.4=0.963028	c=2	e=4.95754e-05
0.91^0.5=0.953987	c=2	e=4.82986e-05
0.91^0.6=0.945028	c=2	e=4.3173e-05
0.91^0.7=0.937	c=1	e=0.000885495
0.91^0.8=0.928	c=1	e=0.000672554
0.91^0.9=0.919	c=1	e=0.000377131

-- 
Regards,
  Bob Monsen

Our minds are finite, and yet even in those circumstances of finitude, we 
are surrounded by possibilities that are infinite, and the purpose of human 
life is to grasp as much as we can out of that infinitude.
- Alfred North Whitehead

Bob,

Do you have a reference for the iterative gen. binomial algorithm? I
need to look into the convergence properties a bit more carefully. I
couldn't find anything by googling. Also, could you post or email a
snippet of the code you used for the table you generated? I'm not sure
about the pseudocode you posted the other day.

GZ

On Tue, 27 Dec 2005 15:01:13 -0800, Glenn Zelniker wrote:

> Bob,
> 
> Do you have a reference for the iterative gen. binomial algorithm? I
> need to look into the convergence properties a bit more carefully. I
> couldn't find anything by googling. Also, could you post or email a
> snippet of the code you used for the table you generated? I'm not sure
> about the pseudocode you posted the other day.
> 
> GZ

The binomial theorem was Fermat and Pascal, I think. Newton was the first
to use it with negative exponents.

http://mathworld.wolfram.com/BinomialTheorem.html

Here is the example:
----------------------------
#include <stdio.h>
#include <math.h>

double binomial_power(double x, double n, int * cnt)
{
    double y = 1.0;                     /* result */
    double z = 1.0;                     /* running term */
    double i = 1.0;                     /* for n! */
    double p = pow(x,n);

    x = 1.0 - x;                        /* computes (1+z)^n */

    while (fabs((p-y)/p) > 0.001)
    {
        z *= -x*n / i;
        n -= 1.0;
        i += 1.0;
        
        y += z;

        *cnt = *cnt + 1;
    }        
        
    return y;
}

test_binomial()
{
    double x = 0.1;
    double n = 0.1;
    int i, j;

    x = .01;
    for (i = 0; i < 10; i++)
    {
        n = 0.0;
        for (j = 0; j < 10; j++)
        {
            int cnt = 0;
            double b;

            b = binomial_power(x,n,&cnt);

            printf("%lg^%lg=%lg\tc=%d\te=%lg\n", x, n, b, cnt, fabs(pow(x,n)-b));
            n += 0.1;
        }
        x += 0.1;
    }
}

int main()
{
    test_binomial();
}


-- 
Regards,
  Bob Monsen

One cannot inquire into the foundations and nature of mathematics without 
delving into the question of the operations by which the mathematical 
activity of the mind is conducted.  If one failed to take that into account, 
then one would be left studying only the language in which mathematics is 
represented rather than the essence of mathematics.
- Luitzen Brouwer