@2cUniverse I ask again Why do you use square-free numbers? The techniques that we're discussing here will work for any positive integer (other than 1).

Use the green uparrow on the toolbar of the MaplePrimes editor to upload your worksheet. Even if it can't display your worksheet inline (which often happens), it will still put in a link such that your worksheet can be downloaded.

For the example number that you show (den = 11842585), we have Totient(den)/CarmichaelLambda(den) = 704, which is fairly high. I can probably speed it up some, but when that ratio is large, it becomes more difficult. It can also be speeded up significantly if you'd be content with some bases of the specified period rather than all bases.

Please define "divergence point". Is it any point every neighborhood of which intersects all basins?

@jonsimoniowa I confirm your report that there is a bug with respect to .graphml files in Maple 2023.1 (the current version). Futhermore, the bug is specifically in the Import. I inspected the result of your Export with a plain-text editor, and the file contained the correct weights. (The erroneous weight matrix is <0, 2, 2; 3, 0, 3; 4, 4, 0>.) Using ExportGraph and/or ImportGraph produce identical erroneous results.

The erroneous values aren't necessarily lower than the originals. It seems like the first n nonzero values encountered (reading left-to-right, then top-to-bottom) are duplicated as constant values for the nonzero entries in each of the n rows.

@2cUniverse You wrote:

• I have checked numbers up to 10^8 with your code. With MultiplicativeOrder its possible up t 10^6.

  I use this for visualizing small periods (< 40)  of a period-base combination (see attached file) .

  My code produces a lineplot which has a lot of polygones. Then I fill the Polygones random with colors.

I don't see any attached file in your Reply. If you attach the file, I can advise further on speeding up the code. For example, if you want to compute the multiplicative order of all members of the multiplicative group modulo n, there are ways to do that that are much much faster than the sum of the times to compute them individually.

@2cUniverse I think that this is what you're looking for: Carmichael's lambda function returns the largest divisor of the totient that can be a period, and all other periods will be divisors of this number. So

NumberTheory:-CarmichaelLambda(78);
                               12

NumberTheory:-CarmichaelLambda(123);
                               40

Like the totient, this function can be easily computed from the prime factorization of its argument. Details are given in this Wikipedia article https://en.wikipedia.org/wiki/Carmichael_function

CarmichaelLambda(n) is the largest possible value of MultiplicativeOrder(a,n), and any value of MultiplicativeOrder(a,n) is a divisor of CarmichaelLambda(n).

@RezaZanjirani You can't split a procedure definition, or any other statement, over multiple execution groups. In other words, you can't have a "> " input prompt in the middle of a procedure, or in the middle of any other statement.

@2cUniverse You asked:

• This (Totient, fundamental period) is only for primes or can I use it also for squarefree numbers in general.

I think for now you should ignore the distinction that I was trying to make between period and fundamental period. What you've been calling period all along is what I'm calling fundamental period, and I'll just say period for now.

The totient is meaningful for any positive integer. As far as I know, there's nothing special about square-free integers with respect to the period of the inverse, or anything else we've discussed in this thread. So, I don't know why you've focused on them.

• Example:  n = 123    ( = 3 x 41)
  Torient(123) = 80   (the number of positive integers coprime to n and not greater  than n)
  Divisors(80) = {1, 2, 4, 5, 8, 10, 16, 20, 40, 80}
  possible periods calculated with MultiplicaticeOrder : { 2, 4, 5, 8, 10, 16, 20, 40}
  Question: why is the 80 not in the period-list ?
  Has this to do with 2, 4, p^k, 2 p^k ?

Yes, since 123 is not of the form p^k or 2*p^k, its totient is not a period. Nonetheless all periods are divisors of the totient.

Note that totient(123) = (3-1)*(41-1). This is easier than explicitly counting the coprimes.

• Is this primitive Root calculating only for primes or may I use it for square free numbers (combined prime factors that have no exponent) ?

Only numbers of the form {2, 4, p^k, 2*p^k} (for odd prime p) have primitive roots. Having a primitive root is equivalent to the totient itself being a period. However, for any positive integer, all possible periods are divisors of the totient. The reason that I mentioned primitive roots is that the process of computing a primitive root is pretty much the same as computing the bases corresponding to any period.

• For our example n=78:
  numbers that have a primitive root  < Totient (78) are 
  {2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24} 

??? I don't know the significance of that set! It's certainly not "numbers that have a primitive root < totient(78)" because 8, 12, 15, 16, 20, 21, and 24 do not have primitive roots.

• Divisors(Totient(78)) are  {2, 3, 4, 6, 8, 12, 24}
  Intersection of above Sets gives:   {2, 3, 4, 6}

??? I don't know what sets you're intersecting! Certainly the 1st set you showed contains {8, 12, 24}.

• Now 8 and 12 are removed as you mentioned.

??? I didn't say that 12 was removed! I said 8 and 24.

@minhthien2016 Like this:

L := map2(`~`[`=`], [a,b], [[2,3], [3,7], [9,10]]):
r1 := Display~(eval~(temp1=temp2, L),inert=false):
r2 := eval~(p,L):
mrow~(Typeset~(EV~(r1)),mo("&equals;"),Typeset~(EV~(r2)));
 

Here is a plot of the Tally list (T from @acer 's or my code). Although this is tangential to the main point of this thread (speeding up some code), it's quite satisfying for a calendar nerd such as myself.

The steep incline for the first 7 days and steep decline for the last 7 days is expected given that the date is always a Sunday. The sharp peak on April 19th is unexpected though.

#This code depends on the list T produced at the end of the previous code.
dataplot(
    rhs~(T), style= point, view= [default, 0..max(rhs~(T))], size= [6,3]*99,
    color= purple, symbolsize= 9, symbol= soliddiamond,
    axis[1]= [
        tickmarks= [seq](
            n= nprintf(`#mfrac(mn("%d"),mn("%d"));`, lhs(T[n])[]),
            n= 1..nops(T)
        )
    ],
    axesfont= [TIMES, 9],
    labels= ["\nDate", "Frequency\n"], labeldirections= ["horizontal", "vertical"],
    labelfont= [TIMES, ITALIC, 10],
    title= "Easter date frequency\n", titlefont= [HELVETICA, 16],
    caption= 
        "\nThe Gregorian date of Easter is the Sunday after the full Moon after"
        " the Vernal Equinox." 
);

@Carl Love I was able to make an additional factor-of-2 time improvement simply by replacing trunc(n/d) with iquo(n,d,r), which means that the iquo must be several times faster than integer division followed by trunc. That's surprising. With this change, the k=6 is done in 4.2 s, and I also did the k=7 case in 8.2 minutes.

restart
:
(* Maple's library 'ceil' has a highly symbolic aspect that slows it significantly.
(You can confirm this by showstat(ceil).) But 'iquo' is kernel builtin and doesn't
have this problem. Thus, I wrote this version of 'ceil' that only uses 'iquo'.
It returns ceil(n/d) for integer n and positive integer d: *)
Ceil:= (n,d)-> local r;
    if n::integer then (thisproc(n,d):= iquo(n,d,r) + `if`(r=0 or n<0, 0, 1))
    else 'procname'(n,d)
    fi
:
nequaln:= (n::And(posint, Not(1)))->
local 
    i, %add, n1:= n-1, 

    (* This trick creates an arbitrarily long sequence of catenated *local* names.
    (I've tried to figure out for decades how to do this; I finally got it!)
    Catenated names (such as v7) work significantly faster than indexed names
    (such as v[7]) as the index variables in the loops. *) 
    v:= (v-> subs(_x= v, ()-> local _x; _x)())~([_v||(1..n1)]),
    V:= (n-2)*24 -~ (0, seq['scan'= `+`](v[i], i= 2..n1))  #seq of partial sums 
;
    -(eval@subs)(
        v[1]= V[1] - n1, %add= add,
        foldl(
            %add,
            min(0, Ceil(V[-1], 2) - min(v[-1] + 1, V[-1])),
            seq(v[i]= Ceil(V[i-1], n-i+2)..v[i-1], i= n1..2, -1)
        )        
    )           
:
[seq](CodeTools:-Usage(nequaln(k)), k= 2..7);
memory used=1.40MiB, alloc change=0 bytes, 
cpu time=62.00ms, real time=66.00ms, gc time=0ns

memory used=0.70MiB, alloc change=0 bytes, 
cpu time=16.00ms, real time=23.00ms, gc time=0ns

memory used=499.52KiB, alloc change=0 bytes, 
cpu time=0ns, real time=7.00ms, gc time=0ns

memory used=2.09MiB, alloc change=0 bytes,
cpu time=62.00ms, real time=49.00ms, gc time=0ns

memory used=137.33MiB, alloc change=32.00MiB,
cpu time=4.16s, real time=4.17s, gc time=0ns

memory used=13.33GiB, alloc change=0 bytes,
cpu time=8.16m, real time=8.17m, gc time=343.75ms

              [0, 48, 816, 10642, 117788, 1143743]

@minhthien2016 To join the two commands as you propose, you just need to type them together, like

Display(temp, inert= false) = value(temp);

@acer Okay, it was the option threadsafe that I was missing earlier. If I include that, and use Threads:-Seq['tasksize'= 1] and Iterator:-SplitRanks, then I consistently get accurate results with a time just slightly longer (4%) than yours, such as 311 ms.

Can you give an example of Maple code where this behavior happens? There's no need at this point to give any elaborate explanation of it or to verify that indeed the results are wrong..

@acer You wrote:

• I only mentioned the bookkeeping because I made a slight adjustment to EasterC_Mpl's loop (k from BEGIN to END), and I'd had just the one coffee. That seemed simpler than fiddling with the Array indexing. It seemed reasonably straightforward because EasterC_Mpl already had the useful BEGIN & END parameters.

Ah. The adjustment that I made to the indexing is needed for BEGIN <> 1.

While you were writing the above, I was adding an assessment of multithreading time overhead to my last Reply. Check it out.

@acer Thank you. It's quite satisfying for me to make such a performance improvement on an OP's code with so little programming effort spent on refactoring (and I did initially predict that that refactoring would be "trivial") .

Your multi-threaded performance is impressive. My only attempt at multi-threading this last night produced accurate resuts but using slightly longer time (2.3 s vs. the 1.5 s for single-threaded). I used something like (recalling this from memory -- I didn't save the code):

Threads:-Seq['tasksize'= 1](EasterC(r[1], add(r)-1, E), r= Iterator:-SplitRanks(57*10^5));

Here is the timing of your code run on my computer, Intel Core i7-7700HQ @ 2.80 GHz x 4 cpus x 2 threads/cpu:

E := CodeTools:-Usage(Easter(1, 57*10^5)):
memory used=43.49MiB, alloc change=42.48MiB, 
cpu time=1.39s, real time=1.39s, gc time=15.62ms

EasterTask(1, 2): # first time overhead
F := CodeTools:-Usage(EasterTask(1, 57*10^5)):
memory used=43.69MiB, alloc change=43.49MiB, 
cpu time=2.27s, real time=298.00ms, gc time=0ns

#multi-threading overhead percentage:
(1 - 1.39 / .298 / kernelopts(numcpus))*100;
                          41.69463088

That's quite impressive performance! That overhead of 42% is lower than I usually get on this computer, though my experience of that is based on non-compiled Maple. Also, it's surprising to me that your multi-threading incurs almost zero overhead on "memory used", 0.2 MiB.

I don't know what you mean by "Did I get the bookkeeping right?" Your results are obviously accurate---identical to the single-threaded.