@2cUniverse Sorry, I mistakenly thought that you said that you were using Maple 2023. What version are you using? I modified my uses of $ so that they'll work for earlier Maple. Here is the worksheet link:

Carmichael.mw

@MalakMMK Thanks. You did that so quickly. You've made quite-sophisticated changes to the code. You know a lot more Maple than I originally thought. 

Any Maple computation that can be done in 2D Input can be done in 1D. (But the converse of that is very far from true.) Indeed, any 2D Input is internally converted to 1D as the input is being parsed.

Since you've made a bona fide effort to work with me on this, I'll show you how to do for loops and whatnot in 1D. It might take me a day to get to it.

@MalakMMK Okay, I read the paper. They're using Mathematica, not Maple. I don't think that the difference in the black-point count is significant. The only difference that might be significant is the comparison with the counts for the other algorithms.

It's possible to rewrite the code so that every arithmetic step is done in 64-bit hardware floating point "double precision", and then perhaps it would function exactly like the authors' version. Since (as far as I can see), the authors didn't provide actual code, that's hard to say for sure.

I'm not willing to continue working with code in 2D Input. If you want further advice on this from me, retype it in 1D. The first and last sections of your worksheet are already in 1D input anyway.

> 1D input is the stuff in red boldface monospaced characters, like this.

@MalakMMK I don't see where in the code N occurs. Is it the loop upper bound, which is 20 in your original worksheet? The difference is likely due to a different model of floating-point computation with different rounding errors. Were the authors using Maple?

@MalakMMK Your error is due to using 2D Input (the default input mode). It doesn't recognize the ,= or ++ operators (as well as a host of other modern niceties).

Perhaps going against my better judgement, I'll show you how to do this in 2D Input. (But, fair warning, I'm not about to go down the nearly bottomless rabbit hole of correcting errors caused by 2D Input.)

Make the line in the procedure
    :-BlackCt:= :-BlackCt + 1;  :-BlackZ(:-BlackCt):= x0 + I*y0;
and make the line outside the procedure
    BlackCt:= 0:  BlackZ:= Array(1..0, datatype= complex[8]):
The indexing of BlackZ inside the procedure must be done with parentheses ( ), as shown, not with the usual square brackets [ ].

If you run this on the code exactly as you had it, there are indeed 0 black points. But if you change the function F to z^4 - 1, then there will be 84.

Please use the green uparrow on the toolbar of the MaplePrimes editor to upload a link to the file that you're trying to upload to Maple. Or, if the file is publically accessible on the Internet, you may just give the URL.

@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." 
);