> p:=nextprime(876632563576546313879901876677637463831870864);
> p mod 4;
> q:=nextprime(3793673643546353637363698378338376974637453365);
> q mod 4;
> N:=p*q;
> convert(N, binary);
> evalf(log(N)/log(2));
> s:=7547496497676383275849409470847478469763763513631635376638574575940
> 7843163735531763768383578537;
> w:=s^2 mod N;
> 
> convert(w, binary);evalf(log(w)/log(2));
> line:=array([0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]):
> for i from 1 to 100 do 
>   for j from 1 to 25 do
>     w:=w^2 mod N: 
>     line[j]:=w mod 2:
>   end:
>   print(line);
> end:
> 
> u:=16355645439:
> v:=3544584873717:
> P:=0:
> while not isprime(P) do
>  u:=u+2; v:=v+2;
>  t:=2^225*u*v:
>  P:=t+1:
>  print(isprime(P));
> end do:
> u;v;
> isprime(P);
> P;
> ifactor(P-1);
> convert(P, binary);evalf(log(P)/log(2));
> g:=numtheory[primroot](1,P);
> line:=array([0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]):
> bound:=(P+1)/2:
> for i from 1 to 100 do 
>   for j from 1 to 25 do
>     w:=g&^w mod P: 
>     if w< bound then line[j]:=0 else line[j]:=1 fi:
>   end:
>   print(line);
> end:
>