> with(numtheory);
Warning, new definition for order

[B, F, GIgcd, J, L, M, bernoulli, bigomega, cfrac, cfracpol,

    cyclotomic, divisors, euler, factorEQ, factorset, fermat, ifactor,

    ifactors, imagunit, index, invcfrac, invphi, isolve, isprime,

    issqrfree, ithprime, jacobi, kronecker, lambda, legendre,

    mcombine, mersenne, minkowski, mipolys, mlog, mobius, mroot,

    msqrt, nearestp, nextprime, nthconver, nthdenom, nthnumer, nthpow,

    order, pdexpand, phi, pprimroot, prevprime, primroot, quadres,

    rootsunity, safeprime, sigma, sq2factor, sum2sqr, tau, thue]

> p:=nextprime(6152423417838736357537136764543212343431121);
> q:=nextprime(342525163874787309109318363313893003849749817);

           p := 6152423417838736357537136764543212343431397


          q := 342525163874787309109318363313893003849750057

> n:=p*q;

n := 2107359839422292204771737681451028490472663954664841603555314\
    313876849796609124176339629

> x:=8372615243536338449443838276252611617282394948372726 mod n;

      x := 8372615243536338449443838276252611617282394948372726

> BBS2:=proc(a,b)
> local m,i,s,t;
> m:=a;
> s:=0;
> for i from 1 to floor(evalf(log[2](n))) do
> s:=2*s+(m mod 2);
> m:=m^2 mod b;
> od;
> t:=0;
> for i from 1 to floor(evalf(log[2](n))) do
> t:=2*t+(m mod 2);
> m:=m^2 mod b;
> od;
> RETURN(s,t);
> end;
> 

BBS2 := proc(a, b)
local m, i, s, t;
    m := a;
    s := 0;
    for i to floor(evalf(log[2](n))) do
        s := 2*s + (m mod 2); m := m^2 mod b
    od;
    t := 0;
    for i to floor(evalf(log[2](n))) do
        t := 2*t + (m mod 2); m := m^2 mod b
    od;
    RETURN(s, t)
end

> BBS := proc(a, b)
> local m, i, t;
>     t := 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, 0, 0, 0, 0, 0, 0, 0]);
>     m := a;
>     while 1 < m do
>         for i to 32 do t[i] := m mod 2; m := m^2 mod b od;
>         print(t)
>     od;
> end;

BBS := proc(a, b)
local m, i, t;
    t := 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, 0, 0, 0, 0, 0, 0, 0]);
    m := a;
    while 1 < m do
        for i to 32 do t[i] := m mod 2; m := m^2 mod b od;
        print(t)
    od
end

> BBS2(x,n);

200265730973064923233013530898033801499285124246418507293313412114\

    331286774217787910597, 4833649762141637826649400700627065014109\
    05309268936550566180751665595073753719195138678

> PHI := proc(x, y)
> local m, i, j, s, t;
>     m := y;
>     s := x;
>     for j from 1 to 160 do
>         i := BBS(s, n);
>         if m mod 2 = 0 then s := i[1]; m := 1/2*m
>         else s := i[2]; m := 1/2*m - 1/2
>         fi;
>         #print(s);
>     od;
>     if m mod 2 = 0 then RETURN(i[1]) else RETURN(i[2]) fi;
> end;

PHI := proc(x, y)
local m, i, j, s, t;
    m := y;
    s := x;
    for j to 160 do
        i := BBS(s, n);
        if m mod 2 = 0 then s := i[1]; m := 1/2*m
        else s := i[2]; m := 1/2*m - 1/2
        fi
    od;
    if m = 0 then RETURN(i[1]) else RETURN(i[2]) fi
end

> PHI(x,0);

135746924991083367666122404235855108061410208181231446712723704743\
    1286956132944238889493

> cha:=9994309205154329926376315592441594389238946449388542547506334647615771507206022569348266793813603644792118986108097649580413414156959932616705066279005364723047157382432633314821460193215501784480779;

cha := 99943092051543299263763155924415943892389464493885425475063\
    34647615771507206022569348266793813603644792118986108097649580\
    41341415695993261670506627900536472304715738243263331482146019\
    3215501784480779

>