Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]

Groups > comp.programming > #337

Re: Fourth root of integer

Show all headers | View raw

On Mon, 16 May 2011 13:28:29 -0700, Dann Corbit wrote:
>  oldwolf@inspire.net.nz says...
>> 
>> Looking for an algorithm to find the fourth root of a large integer (I
>> have a multi-precision integer library already). If there is no exact
>> fourth root then the floor or ceiling is OK, I don't need a 'remainder'
>> or anything.
...
> Assuming you can get any sort of decent estimate of the root, you can
> implement this root polisher in multiple precision:
> 
> /*
> xprimeB :=
> 
>             16        12          8  2        3  4       4    4
>         x (x   + 120 x   a + 546 x  a  + 336 a  x  + 21 a ) (x  - a)
>     x - -------------------------------------------------------------
>          20        16           12  2         8  3        4  4      5
>         x   + 203 x   a + 1554 x   a  + 1918 x  a  + 413 x  a  + 7 a
> 
> */
> double xprimeB(double x, double a)

[snip code and sample output]

If this is a 7th-order approximation (as suggested in Christopher 
Henrich's followup) it won't require many iterations, even for numbers
with many digits.  However, in multiple precision arithmetic, even a 
single iteration of the above may be slower than just finding the 
square root twice.  Typical MP packages use efficient square root 
computations that take slightly longer than a single divide.

For example, on my computer using mpfr (gmp) and with default 
precision 256 bits, one iteration of xprimeB takes ~ 3.8 times as 
long as 2 mpfr_sqrt calls.  Here are some times, in seconds per 
half-million passes.  DeltaI is t_I-t_1.

I   Item          t_I
1. Empty loop:  0.247415   Passes: 500000
2. 2-Sqrt:      1.015866   Delta2:  0.768451
3. xprimeB:     3.172006   Delta3:  2.924591
4. divide:      0.500479   Delta4:  0.253064
(t_3-t_1)/(t_2-t_1):    3.805826

Code for loop 2 that took about .77 seconds longer than the 
empty loop consists of following, where RM is MPFR_RNDU rounding 
mode, RSEED=13343, ir is a 32-bit unsigned long int, and ttime() 
is as shown below sig:

  srandom(RSEED);
  t0 = ttime(0);
  for (i=0; i<REPS; ++i) {
    ir = random ();
    mpfr_set_ui (q, ir, RM);
    mpfr_sqr  (q2, q, RM);
    mpfr_sqr  (q4, q2, RM);
    mpfr_sqrt (r2, q4, RM);
    mpfr_sqrt (r,  r2, RM);
  }
  t2 = ttime(t0);

The "empty loop" is like loop 2 but without mpfr_sqrt() 
calls.

Loop 3 is like loop 2 but with 
    xprimeB (q, q4, pkb);
in place of the two mpfr_sqrt() calls.   xprimeB() was
coded with mpfr_mul(), mpfr_add(), mpfr_div() etc, using
global mpfr_t variables initialized when main() began.

Loop 4 is like loop 2 but has
    mpfr_div (pkb, q4, q2, RM);
in place of the two mpfr_sqrt() calls.  For numbers as in 
this test, 2 sqrts take about as long as 3 divides.

-- 
jiw

double ttime(double base) {
  struct timeval tod;
  gettimeofday(&tod, NULL);
  return tod.tv_sec + tod.tv_usec/1e6 - base;
}

Back to comp.programming | Previous | Next — Previous in thread | Next in thread | Find similar | Unroll thread

Thread

Fourth root of integer Old Wolf <oldwolf@inspire.net.nz> - 2011-05-15 20:30 -0700
  Re: Fourth root of integer kym@kymhorsell.com - 2011-05-16 04:10 +0000
    Re: Fourth root of integer Jacob <asdf@hsu.edu> - 2011-05-16 02:49 -0400
  Re: Fourth root of integer pete <pfiland@mindspring.com> - 2011-05-16 09:46 -0400
    Re: Fourth root of integer Charles Richmond <frizzle@tx.rr.com> - 2011-05-16 22:39 -0500
      Re: Fourth root of integer pete <pfiland@mindspring.com> - 2011-05-17 09:17 -0400
  Re: Fourth root of integer rossum <rossum48@coldmail.com> - 2011-05-16 16:41 +0100
  Re: Fourth root of integer Ralph Boland <rpboland@gmail.com> - 2011-05-16 11:06 -0700
  Re: Fourth root of integer Dann Corbit <dcorbit@connx.com> - 2011-05-16 13:28 -0700
    Re: Fourth root of integer Christopher Henrich <chenrich@monmouth.com> - 2011-05-18 23:41 -0400
      Re: Fourth root of integer rossum <rossum48@coldmail.com> - 2011-05-19 09:54 +0100
        Re: Fourth root of integer Christopher Henrich <chenrich@monmouth.com> - 2011-05-20 15:16 -0400
      Re: Fourth root of integer - fourth-root.mws [1/1] Dann Corbit <dcorbit@connx.com> - 2011-05-20 12:42 -0700
        Re: Fourth root of integer - fourth-root.mws [1/1] Christopher Henrich <chenrich@monmouth.com> - 2011-05-20 16:10 -0400
    Re: Fourth root of integer James Waldby <not@valid.invalid> - 2011-05-19 07:01 +0000
      Re: Fourth root of integer James Waldby <not@valid.invalid> - 2011-05-21 00:37 +0000

csiph-web