mersenneforum.org  

Go Back   mersenneforum.org > Great Internet Mersenne Prime Search > Math > Number Theory Discussion Group

Reply
 
Thread Tools
Old 2021-11-13, 20:46   #23
bhelmes
 
bhelmes's Avatar
 
Mar 2016

5678 Posts
Default

Quote:
Originally Posted by bhelmes View Post
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.
bhelmes is offline   Reply With Quote
Old 2021-11-13, 21:50   #24
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

2×32×13×17 Posts
Default

Here is my solution written in C

Code:
#include <gmp.h>
#include <stdio.h>
#include <stdlib.h>

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
paulunderwood is offline   Reply With Quote
Old 2021-11-13, 22:24   #25
bhelmes
 
bhelmes's Avatar
 
Mar 2016

3×53 Posts
Default

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.
bhelmes is offline   Reply With Quote
Old 2021-11-13, 22:47   #26
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

2×32×13×17 Posts
Default

Quote:
Originally Posted by bhelmes View Post
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);
        mpz_add(r, r, a);
    } while (mpz_sizeinbase(r, 2) > p);
    mpz_clear(a);
}

Last fiddled with by paulunderwood on 2021-11-13 at 23:08
paulunderwood is offline   Reply With Quote
Old 2021-11-15, 17:33   #27
bhelmes
 
bhelmes's Avatar
 
Mar 2016

1011101112 Posts
Default

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.


bhelmes is offline   Reply With Quote
Old 2021-11-15, 21:28   #28
R. Gerbicz
 
R. Gerbicz's Avatar
 
"Robert Gerbicz"
Oct 2005
Hungary

2×32×5×17 Posts
Default

Quote:
Originally Posted by bhelmes View Post
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);
        mpz_add(r, r, tmp);
    } 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.
R. Gerbicz is offline   Reply With Quote
Old 2021-11-26, 04:03   #29
Happy5214
 
Happy5214's Avatar
 
"Alexander"
Nov 2008
The Alamo City

30816 Posts
Default

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.
Happy5214 is offline   Reply With Quote
Old 2021-11-26, 05:05   #30
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

2·32·13·17 Posts
Default

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)<abs(d).
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
paulunderwood is offline   Reply With Quote
Old 2021-11-26, 19:47   #31
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

2×32×13×17 Posts
Default

Quote:
Originally Posted by paulunderwood View Post
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
File Type: c bernhard.c (1.7 KB, 20 views)

Last fiddled with by paulunderwood on 2021-11-26 at 19:55
paulunderwood is offline   Reply With Quote
Reply

Thread Tools


Similar Threads
Thread Thread Starter Forum Replies Last Post
Factorial modulo a prime axn Computer Science & Computational Number Theory 66 2011-09-01 21:55
modulo operation for polynomials? smslca Math 3 2011-04-18 17:18
Order of 3 modulo a Mersenne prime T.Rex Math 7 2009-03-13 10:46
N! modulo for large N Cyclamen Persicum Math 2 2003-12-10 10:52
The modulo operation, how is it computed? eepiccolo Math 7 2003-01-08 03:07

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


Wed Jan 19 04:05:42 UTC 2022 up 179 days, 22:34, 0 users, load averages: 1.95, 1.68, 1.36

Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2022, Jelsoft Enterprises Ltd.

This forum has received and complied with 0 (zero) government requests for information.

Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.2 or any later version published by the Free Software Foundation.
A copy of the license is included in the FAQ.

≠ ± ∓ ÷ × · − √ ‰ ⊗ ⊕ ⊖ ⊘ ⊙ ≤ ≥ ≦ ≧ ≨ ≩ ≺ ≻ ≼ ≽ ⊏ ⊐ ⊑ ⊒ ² ³ °
∠ ∟ ° ≅ ~ ‖ ⟂ ⫛
≡ ≜ ≈ ∝ ∞ ≪ ≫ ⌊⌋ ⌈⌉ ∘ ∏ ∐ ∑ ∧ ∨ ∩ ∪ ⨀ ⊕ ⊗ 𝖕 𝖖 𝖗 ⊲ ⊳
∅ ∖ ∁ ↦ ↣ ∩ ∪ ⊆ ⊂ ⊄ ⊊ ⊇ ⊃ ⊅ ⊋ ⊖ ∈ ∉ ∋ ∌ ℕ ℤ ℚ ℝ ℂ ℵ ℶ ℷ ℸ 𝓟
¬ ∨ ∧ ⊕ → ← ⇒ ⇐ ⇔ ∀ ∃ ∄ ∴ ∵ ⊤ ⊥ ⊢ ⊨ ⫤ ⊣ … ⋯ ⋮ ⋰ ⋱
∫ ∬ ∭ ∮ ∯ ∰ ∇ ∆ δ ∂ ℱ ℒ ℓ
𝛢𝛼 𝛣𝛽 𝛤𝛾 𝛥𝛿 𝛦𝜀𝜖 𝛧𝜁 𝛨𝜂 𝛩𝜃𝜗 𝛪𝜄 𝛫𝜅 𝛬𝜆 𝛭𝜇 𝛮𝜈 𝛯𝜉 𝛰𝜊 𝛱𝜋 𝛲𝜌 𝛴𝜎 𝛵𝜏 𝛶𝜐 𝛷𝜙𝜑 𝛸𝜒 𝛹𝜓 𝛺𝜔