mersenneforum.org > Math calculation of modulo Mp
 Register FAQ Search Today's Posts Mark Forums Read

2021-11-13, 20:46   #23
bhelmes

Mar 2016

1011101112 Posts

Quote:
 Originally Posted by bhelmes for a=0 the loop with mpz_sizeinbase is not entered.
That was a logical error, a<>0, the do loop is entered and by chance a=0 for mpz_fdiv_q_2exp (a, a, mp);

f=174224571863520493293247799005065324265471
and p=137 by the way.

We are approaching the error.

 2021-11-13, 21:50 #24 paulunderwood     Sep 2002 Database er0rr F8D16 Posts Here is my solution written in C Code: #include #include #include void mpz_mod_mp (mpz_t r, mpz_t a, mp_bitcnt_t p) { // assumed that r is initialized and a is expendable do { mpz_fdiv_r_2exp(r, a, p); mpz_fdiv_q_2exp(a, a, p); mpz_add(r, r, a); } while (mpz_sizeinbase(r, 2) > p); } int main() { mpz_t r, a, f; mp_bitcnt_t p; mpz_init(r); mpz_init(a); mpz_init(f); p=137; printf("%d is exponent.\n",p); mpz_set_ui(f, 1); mpz_mul_2exp(f, f, p); mpz_sub_ui(f, f, 1); printf("Mp is "); mpz_out_str(stdout, 10, f); printf(".\n"); mpz_set_str(a,"10584739298444444444444444444409185657387898722659423897260",10); printf("a is "); mpz_out_str(stdout, 10, a); printf(".\n"); mpz_mod(r, a, f); printf("Result in-built is "); mpz_out_str(stdout, 10, r); printf(".\n"); mpz_mod_mp(r, a, p); // a was changed by mpz_mod_mp!! printf("Result written is "); mpz_out_str(stdout, 10, r); printf(".\n"); } A general rule: Write your tests outside the function being tested, not inside! Last fiddled with by paulunderwood on 2021-11-13 at 22:17
 2021-11-13, 22:24 #25 bhelmes     Mar 2016 3×53 Posts In gmp I could use the same mpz_t variable A like mpz_mod (A, A, f); I tried to write a function like this mpz_mod_mp (A, A, Mp), therefore I initialized a help variable. Thanks for your patience.
2021-11-13, 22:47   #26
paulunderwood

Sep 2002
Database er0rr

3·1,327 Posts

Quote:
 Originally Posted by bhelmes In gmp I could use the same mpz_t variable A like mpz_mod (A, A, f); I tried to write a function like this mpz_mod_mp (A, A, Mp), therefore I initialized a help variable. Thanks for your patience.
In that case you need:

Code:
void mpz_mod_mp (mpz_t r, mpz_t A, mp_bitcnt_t p)
{
mpz_t a;
mpz_init_set(a, A);
do
{
mpz_fdiv_r_2exp(r, a, p);
mpz_fdiv_q_2exp(a, a, p);
} while (mpz_sizeinbase(r, 2) > p);
mpz_clear(a);
}

Last fiddled with by paulunderwood on 2021-11-13 at 23:08

 2021-11-15, 17:33 #27 bhelmes     Mar 2016 17716 Posts Something is boogie: I used my variable a and put it in your solution: Result is wrong: 137 is exponent. Mp is 174224571863520493293247799005065324265471. a is 15580960185386716076816921462913615019618152413172348145118933359820136234923109492. Result in-built is 10584739298409185657387898722659423897260. Result written is 89430325577680987480245132687307139815991.
2021-11-15, 21:28   #28
R. Gerbicz

"Robert Gerbicz"
Oct 2005
Hungary

1,531 Posts

Quote:
 Originally Posted by bhelmes Something is boogie: I used my variable a and put it in your solution: Result is wrong: 137 is exponent. Mp is 174224571863520493293247799005065324265471. a is 15580960185386716076816921462913615019618152413172348145118933359820136234923109492. Result in-built is 10584739298409185657387898722659423897260. Result written is 89430325577680987480245132687307139815991.
Yeah, you're really messed up with a trivial algorithm, a corrected one: (here you don't need to use ptr, so pointers, but give a clean code)
Code:
void mpz_mod_mp (mpz_ptr r, mpz_ptr a, mp_bitcnt_t p)
{// the result is r=a%(2^p-1)
mpz_t tmp;
mpz_init(tmp);

mpz_set(r,a);
do
{
mpz_fdiv_r_2exp(tmp, r, p);
mpz_fdiv_q_2exp(r, r, p);
} while (mpz_sizeinbase(r, 2) > p);

mpz_clear(tmp);
}

Last fiddled with by R. Gerbicz on 2021-11-15 at 21:37 Reason: corrected the code.

 2021-11-26, 04:03 #29 Happy5214     "Alexander" Nov 2008 The Alamo City 23×97 Posts mpz_t is formally defined in gmp.h as a typedef for an array of exactly one __mpz_struct internal struct, so it's implicitly a pointer. By the way, note to mods, this is a programming thread in the Math forum.
2021-11-26, 05:05   #30
paulunderwood

Sep 2002
Database er0rr

F8D16 Posts

Quoting the excellent GMP website documentation:

Quote:
 Divide n by d, forming a quotient q and/or remainder r. For the 2exp functions, d=2^b. The rounding is in three styles, each suiting different applications. cdiv rounds q up towards +infinity, and r will have the opposite sign to d. The c stands for “ceil”. fdiv rounds q down towards -infinity, and r will have the same sign as d. The f stands for “floor”. tdiv rounds q towards zero, and r will have the same sign as n. The t stands for “truncate”. In all cases q and r will satisfy n=q*d+r, and r will satisfy 0<=abs(r)
I changed "fdiv" to "cdiv" and all is well with Bernhard's problematic example.

Bernhard should think of these artifacts as a "feature". For example, if he is taking GCDs then there is no problem.

For extra speed, I would shift Mp left to the word boundary and calculate at the word level, before any necessary drop-downs to the bit level. A sample program using the "word level" logic is attached. The "adjustment" might need some tweaking to account for the sign bit.

Last fiddled with by paulunderwood on 2021-11-26 at 19:43

2021-11-26, 19:47   #31
paulunderwood

Sep 2002
Database er0rr

3×1,327 Posts

Quote:
 Originally Posted by paulunderwood For extra speed, I would shift Mp left to the word boundary and calculate at the word level, before any necessary drop-downs to the bit level. A sample program using the "word level" logic is attached. The "adjustment" might need some tweaking to account for the sign bit.
I have done the program with a while loop rather than a do-while loop. You can strip out the shifted functionality as it makes no difference to timings -- for a random number it does make sense but for Mersenne numbers not.

Here is the updated code.
Attached Files
 bernhard.c (1.7 KB, 21 views)

Last fiddled with by paulunderwood on 2021-11-26 at 19:55

 Similar Threads Thread Thread Starter Forum Replies Last Post axn Computer Science & Computational Number Theory 66 2011-09-01 21:55 smslca Math 3 2011-04-18 17:18 T.Rex Math 7 2009-03-13 10:46 Cyclamen Persicum Math 2 2003-12-10 10:52 eepiccolo Math 7 2003-01-08 03:07

All times are UTC. The time now is 00:04.

Sat Jan 22 00:04:18 UTC 2022 up 182 days, 18:33, 0 users, load averages: 2.24, 1.46, 1.31