1

I have tried to implement the algorithm described in here to find primitive roots for a prime number.

It works for small prime numbers, however as I try big numbers, it doesn't return correct answers anymore.

I then notice that a^(p-1)/pi tends to be a big number, it returns inf in MATLAB, so I thought factorizing (p-1) could help, but I am failing to see how.

I wrote a small piece of code in MATLABand here it is.

clear all
clc

%works with prime =23,31,37,etc.

prime=761; %doesn't work for this value

F=factor(prime-1); % the factors of prime-1
 for i = 2: prime-1
        a=i;
        tag =1;
        for j= 1 :prime-1
            if (isprime(j))
                 p_i = j;
                 if(mod(a^((prime-1)/p_i),prime)== 1)
                     tag=0;
                     break 
                else
                    tag = tag +1;
                 end
            end
        end
        if (tag > 1 )
            a %it should only print the primitive root
            break
        end
end

Any input is welcome. Thanks

Adriaan
  • 17,741
  • 7
  • 42
  • 75
dede9714
  • 73
  • 2
  • 13
  • 2
    You're getting `INF`... probably because of overflow. 761 is probably too high of a number which makes the input into `mod` very large. You may be better off using MATLAB's [`vpa`](http://www.mathworks.com/help/symbolic/vpa.html) module instead. – rayryeng Aug 26 '15 at 06:06
  • 2
    The largest number MATLAB can operate with is [`1.7977e+308`](http://se.mathworks.com/help/matlab/ref/realmax.html), unless you have the symbolic toolbox and use vpa as @rayryeng suggests. Your numbers are much higher than the max value. Often, the solution is to work with the logarithm of the numbers you're interested in. Unfortunately I can't immediately see a way you can use it though... – Stewie Griffin Aug 26 '15 at 06:15
  • Thanks for the suggestions I will try it and come back to you. very much appreciated. – dede9714 Aug 26 '15 at 10:48

1 Answers1

1

What Matlab does in this case is it calculates a^((p-1)/p) before taking the modulus. As a^((p-1)/p) quite quickly becomes too large to handle, Matlab seems to resolve this by turning it into a floating point number, losing some resolution and yielding the wrong result when you take the modulus.

As mentioned by @rayreng, you could use an arbitrary precision toolbox to resolve this.

Alternatively, you could split the exponentiation into parts, taking the modulus at each stage. This should be faster, as it is less memory intensive. You could dump this in a function and just call that.

% Calculates a^b mod c
i = 0;
result = 1;

while i < b
    result = mod(result*a, c);
    i = i + 1;
end
RPM
  • 1,704
  • 12
  • 15