Re: erroneous claim at sources.redhat.com/gsl/ref/gsl-ref_4.html# SEC3 2

Re: erroneous claim at sources.redhat.com/gsl/ref/gsl-ref_4.html# SEC3 2

[Achim Gädke]

so, is there anyone outside, knowing similar problems?

I agree, repeated squaring is likely to be the best solution, because 
the binary representation can be obtained easily on a computer.
So the manual should state, that this algorithm is efficient and in most 
cases the fastest.

Achim

keith.briggs@bt.com wrote:

>Thanks for the reply!   No, I do not have an algorithm.
>I think you would have to program special cases up to some maximum n; after
>that repeated squaring is the only simple way even if it is not optimal.
>Keith
>
>	Dr. Keith M. Briggs
>	Senior Mathematician, Complexity Research, BT Exact
>	http://research.btexact.com/teralab/keithbriggs.html
>	phone: +44(0)1473  work: 641 911 home: 610 517  fax: 642 161
>        mail: Keith Briggs, Polaris 134, Adastral Park, Martlesham, Suffolk
>IP5 3RE, UK
>
>
>-----Original Message-----
>From: Achim Gädke [mailto:achim@element.fkp.physik.tu-darmstadt.de]
>Sent: 15 July 2003 10:19
>To: Briggs,KM,Keith,XVR2 R
>Cc: gsl-discuss@sources.redhat.com
>Subject: Re: erroneous claim at
>sources.redhat.com/gsl/ref/gsl-ref_4.html#SEC3 2
>
>
>Hi!
>
>I could not believe that, but it's true and a nice logic training:
>
>x2=x*x;
>x3=x2*x;
>x5=x3*x2;
>x15=x5*x5*x5;
>with 5 * symbols
>
>instead of:
>
>x2=x*x;
>x4=x2*x2;
>x8=x4*x4;
>x15=x8*x4*x2*x;
>here: 6 * symbols.
>
>Keith, do you have a simple approach to the more efficient algorithm, that 
>seems to base on a factorization ( 15=5*3 ).
>
>Achim
>
>On Tue, 15 Jul 2003 keith.briggs@bt.com wrote:
>
>  
>
>>>Function: double gsl_pow_int (double x, int n) 
>>>      
>>>
>>	>This routine computes the power x^n for integer n. The power is
>>computed using the minimum number of multiplications. 
>>A glance at the source code shows that this is wrong.   It uses repeated
>>squaring, so, for example, x^15 is computed
>>with 6 multiplies, whereas it can be done with 5.
>>Keith
>>
>>    
>>
>
>
>  
>

Re: erroneous claim at sources.redhat.com/gsl/ref/gsl-ref_4.html# SEC3 2

[Achim Gädke]

Yes, this is quite true, that I do not really know, what's the optimal way.

1) In order to fix a documentation bug, we should correct the manual.

2) Here we multiply numbers (plain floating point numbers), so the 
effort to do the optimum should be not to large: approx. O(2*log_2(n)) 
is not bad in contrast to O(n).

3) If it is cheap to find a better way or if we deal with matrices (e.g. 
in characteristic polynoms), we may want to find out a better way to 
split this multiplication. Finding prime factors is quite expensive. For 
21 I avoid one multiplication, but what about the costs?

4) I like this problem, because it's simple to explain and becomes 
difficult when solving it. :-)

Achim

jimmy mcguire wrote:

> Breaking my usual rule against speaking proir to knowing for sure, how 
> about use the powers of two approach if the destination exponent is a 
> power of two or if it is prime, otherwise the break the destination 
> exponent into prime addends and construct multiplicative factors 
> corresponding to the prime additive exponents. Repeat (recurse) as 
> necessary.
>
> Achim Gädke wrote:
>
>> so, is there anyone outside, knowing similar problems?
>>
>> I agree, repeated squaring is likely to be the best solution, because 
>> the binary representation can be obtained easily on a computer.
>> So the manual should state, that this algorithm is efficient and in 
>> most cases the fastest.
>>
>> Achim
>>
>> keith.briggs@bt.com wrote:
>

Re: erroneous claim at sources.redhat.com/gsl/ref/gsl-ref_4.html# SEC3 2

[jimmy mcguire]

Breaking my usual rule against speaking proir to knowing for sure, how 
about use the powers of two approach if the destination exponent is a 
power of two or if it is prime, otherwise the break the destination 
exponent into prime addends and construct multiplicative factors 
corresponding to the prime additive exponents. Repeat (recurse) as 
necessary.

Achim Gädke wrote:

> so, is there anyone outside, knowing similar problems?
>
> I agree, repeated squaring is likely to be the best solution, because 
> the binary representation can be obtained easily on a computer.
> So the manual should state, that this algorithm is efficient and in 
> most cases the fastest.
>
> Achim
>
> keith.briggs@bt.com wrote:
>
>> Thanks for the reply!   No, I do not have an algorithm.
>> I think you would have to program special cases up to some maximum n; 
>> after
>> that repeated squaring is the only simple way even if it is not optimal.
>> Keith
>>
>>     Dr. Keith M. Briggs
>>     Senior Mathematician, Complexity Research, BT Exact
>>     http://research.btexact.com/teralab/keithbriggs.html
>>     phone: +44(0)1473  work: 641 911 home: 610 517  fax: 642 161
>>        mail: Keith Briggs, Polaris 134, Adastral Park, Martlesham, 
>> Suffolk
>> IP5 3RE, UK
>>
>>
>> -----Original Message-----
>> From: Achim Gädke [mailto:achim@element.fkp.physik.tu-darmstadt.de]
>> Sent: 15 July 2003 10:19
>> To: Briggs,KM,Keith,XVR2 R
>> Cc: gsl-discuss@sources.redhat.com
>> Subject: Re: erroneous claim at
>> sources.redhat.com/gsl/ref/gsl-ref_4.html#SEC3 2
>>
>>
>> Hi!
>>
>> I could not believe that, but it's true and a nice logic training:
>>
>> x2=x*x;
>> x3=x2*x;
>> x5=x3*x2;
>> x15=x5*x5*x5;
>> with 5 * symbols
>>
>> instead of:
>>
>> x2=x*x;
>> x4=x2*x2;
>> x8=x4*x4;
>> x15=x8*x4*x2*x;
>> here: 6 * symbols.
>>
>> Keith, do you have a simple approach to the more efficient algorithm, 
>> that seems to base on a factorization ( 15=5*3 ).
>>
>> Achim
>>
>> On Tue, 15 Jul 2003 keith.briggs@bt.com wrote:
>>
>>  
>>
>>>> Function: double gsl_pow_int (double x, int n)     
>>>
>>>     >This routine computes the power x^n for integer n. The power is
>>> computed using the minimum number of multiplications. A glance at 
>>> the source code shows that this is wrong.   It uses repeated
>>> squaring, so, for example, x^15 is computed
>>> with 6 multiplies, whereas it can be done with 5.
>>> Keith
>>>
>>>   
>>
>>
>>
>>  
>>
>
>
>
>

Re: erroneous claim at sources.redhat.com/gsl/ref/gsl-ref_4.html#SEC3 2

[Brian Gough]

keith.briggs@bt.com writes:
 > > Function: double gsl_pow_int (double x, int n) 
 > 	>This routine computes the power x^n for integer n. The power is
 > computed using the minimum number of multiplications. 
 > A glance at the source code shows that this is wrong.   It uses repeated
 > squaring, so, for example, x^15 is computed
 > with 6 multiplies, whereas it can be done with 5.

Thanks for pointing that out.  I've changed the manual to say the
algorithm is "efficient" rather than using the minimum number of
multiplications.

-- 
Brian

RE: erroneous claim at sources.redhat.com/gsl/ref/gsl-ref_4.html#SEC3 2

[Billinghurst, David (CRTS)]

> From: Achim Gädke [mailto:achim@element.fkp.physik.tu-darmstadt.de]
> Sent: Tuesday, 15 July 2003 7:19 PM
> To: keith.briggs@bt.com
> Cc: gsl-discuss@sources.redhat.com
>Subject: Re: erroneous claim at sources.redhat.com/gsl/ref/gsl-ref_4.html#SEC3 2
>
> Keith, do you have a simple approach to the more efficient algorithm, that 
> seems to base on a factorization ( 15=5*3 ).
>
> Achim

Code for this was put into gcc recently
(http://gcc.gnu.org/ml/gcc-patches/2003-06/msg02472.html)

Comments in the code given below.  Check the follow-up posts (or grab the current
code from CVS).  I recall there was a bug somewhere.

+ /* To evaluate powi(x,n), the floating point value x raised to the
+    constant integer exponent n, we use a hybrid algorithm that
+    combines the "window method" with look-up tables.  For an
+    introduction to exponentiation algorithms and "addition chains",
+    see section 4.6.3, "Evaluation of Powers" of Donald E. Knuth,
+    "Seminumerical Algorithms", Vol. 2, "The Art of Computer Programming",
+    3rd Edition, 1998, and Daniel M. Gordon, "A Survey of Fast Exponentiation
+    Methods", Journal of Algorithms, Vol. 27, pp. 129-146, 1998.  */
+
+ /* Provide a default value for POWI_MAX_MULTS, the maximum number of
+    multiplications to inline before calling the system library's pow
+    function.  powi(x,n) requires at worst 2*bits(n)-2 multiplications,
+    so this default never requires calling pow, powf or powl.  */
+
+ #ifndef POWI_MAX_MULTS
+ #define POWI_MAX_MULTS  (2*HOST_BITS_PER_WIDE_INT-2)
+ #endif
+
+ /* The size of the "optimal power tree" lookup table.  All
+    exponents less than this value are simply looked up in the
+    powi_table below.  This threshold is also used to size the
+    cache of pseudo registers that hold intermediate results.  */
+ #define POWI_TABLE_SIZE 256
+
+ /* The size, in bits of the window, used in the "window method"
+    exponentiation algorithm.  This is equivalent to a radix of
+    (1<<POWI_WINDOW_SIZE) in the corresponding "m-ary method".  */
+ #define POWI_WINDOW_SIZE 3
+
+ /* The following table is an efficient representation of an
+    "optimal power tree".  For each value, i, the corresponding
+    value, j, in the table states than an optimal evaluation
+    sequence for calculating pow(x,i) can be found by evaluating
+    pow(x,j)*pow(x,i-j).  An optimal power tree for the first
+    100 integers is given in Knuth's "Seminumerical algorithms".  */
+

RE: erroneous claim at sources.redhat.com/gsl/ref/gsl-ref_4.html# SEC3 2

[keith.briggs]

There is a section discussing this problem in quite a lot of detail in
Knuth's Seminumerical Algorithms.
Keith

	Dr. Keith M. Briggs
	Senior Mathematician, Complexity Research, BT Exact
	http://research.btexact.com/teralab/keithbriggs.html
	phone: +44(0)1473  work: 641 911 home: 610 517  fax: 642 161
      mail: Keith Briggs, Polaris 134, Adastral Park, Martlesham, Suffolk
IP5 3RE, UK


-----Original Message-----
From: Achim Gädke [mailto:Achim.Gaedke@physik.tu-darmstadt.de]
Sent: 15 July 2003 13:22
To: jimmy mcguire; gsl-discuss@sources.redhat.com
Subject: Re: erroneous claim at
sources.redhat.com/gsl/ref/gsl-ref_4.html# SEC3 2


Yes, this is quite true, that I do not really know, what's the optimal way.

1) In order to fix a documentation bug, we should correct the manual.

2) Here we multiply numbers (plain floating point numbers), so the 
effort to do the optimum should be not to large: approx. O(2*log_2(n)) 
is not bad in contrast to O(n).

3) If it is cheap to find a better way or if we deal with matrices (e.g. 
in characteristic polynoms), we may want to find out a better way to 
split this multiplication. Finding prime factors is quite expensive. For 
21 I avoid one multiplication, but what about the costs?

4) I like this problem, because it's simple to explain and becomes 
difficult when solving it. :-)

Achim

jimmy mcguire wrote:

> Breaking my usual rule against speaking proir to knowing for sure, how 
> about use the powers of two approach if the destination exponent is a 
> power of two or if it is prime, otherwise the break the destination 
> exponent into prime addends and construct multiplicative factors 
> corresponding to the prime additive exponents. Repeat (recurse) as 
> necessary.
>
> Achim Gädke wrote:
>
>> so, is there anyone outside, knowing similar problems?
>>
>> I agree, repeated squaring is likely to be the best solution, because 
>> the binary representation can be obtained easily on a computer.
>> So the manual should state, that this algorithm is efficient and in 
>> most cases the fastest.
>>
>> Achim
>>
>> keith.briggs@bt.com wrote:
>

Re: erroneous claim at sources.redhat.com/gsl/ref/gsl-ref_4.html#SEC3 2

[Achim Gädke]

Hi!

I could not believe that, but it's true and a nice logic training:

x2=x*x;
x3=x2*x;
x5=x3*x2;
x15=x5*x5*x5;
with 5 * symbols

instead of:

x2=x*x;
x4=x2*x2;
x8=x4*x4;
x15=x8*x4*x2*x;
here: 6 * symbols.

Keith, do you have a simple approach to the more efficient algorithm, that 
seems to base on a factorization ( 15=5*3 ).

Achim

On Tue, 15 Jul 2003 keith.briggs@bt.com wrote:

> > Function: double gsl_pow_int (double x, int n) 
> 	>This routine computes the power x^n for integer n. The power is
> computed using the minimum number of multiplications. 
> A glance at the source code shows that this is wrong.   It uses repeated
> squaring, so, for example, x^15 is computed
> with 6 multiplies, whereas it can be done with 5.
> Keith
> 

erroneous claim at sources.redhat.com/gsl/ref/gsl-ref_4.html#SEC3 2

[keith.briggs]

> Function: double gsl_pow_int (double x, int n) 
	>This routine computes the power x^n for integer n. The power is
computed using the minimum number of multiplications. 
A glance at the source code shows that this is wrong.   It uses repeated
squaring, so, for example, x^15 is computed
with 6 multiplies, whereas it can be done with 5.
Keith

	Dr. Keith M. Briggs
	Senior Mathematician, Complexity Research, BT Exact
	http://research.btexact.com/teralab/keithbriggs.html
	phone: +44(0)1473  work: 641 911 home: 610 517  fax: 642 161
        mail: Keith Briggs, Polaris 134, Adastral Park, Martlesham, Suffolk
IP5 3RE, UK