Looking for a Primitive Roots program

[Joe Horn]

Has anybody written an HP Prime CAS program that inputs an integer (not necessarily a prime) and outputs the complete list of its primitive roots? For example, an input of 311 should yield a list of the 120 integers which are primitive roots of 311. Info regarding "primitive roots" (a concept used in Number Theory) can be found here: http://mathworld.wolfram.com/PrimitiveRoot.html

A brute-force search algorithm would be too slow for large inputs. It is much faster to perform an intelligent search for any one primitive root, then use that one to find all the rest.

For testing purposes, here are all 120 primitive roots of 311:

{17 19 22 23 29 31 33 34 37 38 43 44 55 57 58 59 62 66 69
71 74 76 82 85 88 92 93 97 99 101 102 103 110 111 114 115
118 119 122 123 124 129 131 132 133 136 138 145 148 149
151 152 153 154 155 161 164 167 170 172 174 176 177 181
183 184 186 191 194 199 202 203 204 205 207 211 213 215
217 221 227 230 231 232 233 236 238 239 241 244 246 247
248 251 255 257 258 261 263 266 269 271 272 276 281 283
284 285 286 290 295 297 299 301 302 303 306 307 308 309}

Thanks in advance!

[douganc]

Look at this website for a C program then translate to Prime.
http://e-maxx-eng.github.io/algebra/primitive-root.html
Colm

[Joe Horn]

(08-27-2015 05:48 AM)douganc Wrote:  Look at this website for a C program then translate to Prime.
http://e-maxx-eng.github.io/algebra/primitive-root.html
Colm

Thanks. Almost makes me wish I knew both C and PPL well enough to perform that translation.

[roadrunner]

Below is a primitive roots program, for prime numbers only, that I was working on a while ago but had to put it aside because real work got in the way. I recall that somewhere above around 1000 it started giving incorrect answers; I'm not sure why.

Code:

#cas
n_to_z_mod_p(n,z,p):=
BEGIN
 return n^z MOD p;
END;
#end

EXPORT smallest_p_root(p)
BEGIN
//debug;
 if p>1000 then return "that hurts my brain";end;
 local a,b,c,j,k;
 a:=p-1;
 b:=ifactors(a); 
 c:=length(b);
 for j from 2 to a do
  for k from 1 to c step 2 do 
   if n_to_z_mod_p(j,a/b(k),p) == 1 then
    break;
   end;
   if k == c-1 then return j;end;
  end;
 end; 
END;

EXPORT p_root(p)
BEGIN
//debug
 if p>1000 then return "that hurts my brain";end;
 local a,b,c,d,j,k;
 a:=p-1;
 b:=ifactors(a); 
 c:=length(b);
 d:={};
 for j from 2 to a do
  print();
  print("checking " + j);
  for k from 1 to c step 2 do 
   if n_to_z_mod_p(j,a/b(k),p) == 1 then
    break;
   end;
   if k == c-1 then
    d(0):=j;
   end;
  end;
 end;
 return d; 
END;

[Gerald H]

Below a programme for the 50G.

Stack

541
0

returns the prime & first primitive root > 0

541
2

running the programme again returns

541
10

as 10 is the next primitive root > 2 & so on.

Code:


::
  CK2&Dispatch
  # FFFF
  ::
    SWAP
    FPTR2 ^ZABS
    FPTR2 ^DupZIsTwo?
    case
    ::
      SWAPDROP
      Z1_
    ;
    DUP
    ::
      FPTR2 ^ISPRIME
      %0<>
      ?SEMI
      # DE1E
      ERROROUT
    ;
    DUP
    Z1_
    FPTR2 ^QSub
    DUP
    FPTR2 ^MZSQFF
    NULL{}
    SWAP
    #2/
    DUPUNROT
    ZERO_DO
    2SWAP
    DROP
    >TCOMP
    LOOP
    SWAPROT
    FPTR2 ^2LAMBIND
    ROT
    BEGIN
    Z1_
    FPTR2 ^QAdd
    ::
      DUP
      FPTR2 ^ZSQRT
      SWAPDROP
      caseFALSE
      3PICK
      FPTR2 ^ZMod
      FPTR2 ^DupQIsZero?
      caseFALSE
      TRUE
      4UNROLL
      2GETLAM
      #1+_ONE_DO
      DUP
      1GETLAM
      4PICK
      INDEX@
      NTHCOMPDROP
      FPTR2 ^ZQUOText
      5PICK
      FPTR2 ^ModPow
      4PICK
      FPTR2 ^QMod
      FPTR2 ^ZIsOne?
      IT
      ::
        4ROLLDROP
        FALSE
        4UNROLL
        ExitAtLOOP
      ;
      LOOP
      4ROLL
    ;
    UNTIL
    ABND
    SWAPDROP
  ;
;

[Joe Horn]

(08-31-2015 02:19 PM)roadrunner Wrote:  

Code:

...
n_to_z_mod_p(n,z,p):=
BEGIN
 return n^z MOD p;
END;
...

Wouldn't the built-in powmod function be much faster?

[roadrunner]

(08-31-2015 03:29 PM)Joe Horn Wrote:  Wouldn't the built-in powmod function be much faster?

There's a powmod function??? Go figure.

Thanks for your help Joe. I made that change and it seems to have ended the issues I was having with large values of p.

Below is revision A; it seems to be working correctly with large values of p as long as the emulator doesn't run out of memory at p=20000 or so.

Code:

#pragma mode( separator(.,;) integer(h32))

EXPORT smallest_p_root(p)
BEGIN
 local a,b,c,j,k;
 a:=euler(p);
 b:=ifactors(a);
 c:=length(b);
 for j from 2 to a do
  for k from 1 to c step 2 do
   if powmod(j,a/b(k),p) == 1 then
    break;
   end;
   if k == c-1 then return j;end;
  end;
 end;
END;

EXPORT p_root(p)
BEGIN
 local a,b,c,d,j,k;
 a:=euler(p);
 b:=ifactors(a);
 c:=length(b);
 d:={};
 for j from 2 to a do
  print();
  print("checking " + j);
  for k from 1 to c step 2 do
   if powmod(j,a/b(k),p) == 1 then
    break;
   end;
   if k == c-1 then
    d(0):=j;
   end;
  end;
 end;
 return d;
END;

Thanks again!!

road

[Joe Horn]

(08-31-2015 06:02 PM)roadrunner Wrote:  Below is revision A; it seems to be working correctly with large values of p as long as the emulator doesn't run out of memory at p=20000 or so.

Awesome! Many thanks!

[Arno K]

Well, as I yesterday had started to write that program as a CAS-program, I finished it today and want to give it as another example which is much shorter and needs one loop less, as it uses more of the built-in list capabilities, so it should be much faster, but I haven't checked that

Code:


#cas
PrimW(n):=
BEGIN
local l,i,phin,l1;
phin:=euler(n);
l:=ifactors(phin);
l1:=phin/revlist(MAKELIST(l(i),i,1,size(l),2));
l:={};
FOR i FROM 2 TO n-1 DO
IF pos(powmod(i,l1,n),1)==0 THEN
l:=append(l,i);
END;
END;
 return l;
END;
#end

Perhaps it is even possible to omit the second loop, as powmod({2,3,...n-1},{...},n) works, but that provides a list of lists, and here the first objects of all lists belong to 2, the second objects to 3 etc, which makes checking the results difficult.
Arno