Carl Love

Carl Love

28045 Reputation

25 Badges

12 years, 333 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are replies submitted by Carl Love

@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 (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.

@raj2018 You have several for statements of the form 

for x in a..b do

I guess that you expect it to choose some numbers between and b? But it doesn't: It only uses x=a and x=b. How many numbers do you want it to use? Do you want them evenly spaced?

@sursumCorda In native Maple, initializing arrays with the functional form of initializer is faster than using a loop. I used the loop because you're not allowed to create arrays in compiled code.

The reason that your code is slower is the repeated formation and garbage collection of expression sequences for multiple-assignment statements such as

x,y:= a,b;

(This is disappointing to me because I often prefer the neater presentation of the code with the multiple assignments.)

Here I've done the case of "positive" years, which was sufficient to do your 5.7-million-year example. The compiled code finished in under 1.5 seconds without using multi-threading. I'll leave the "negative"-year case to you.

restart
:
EasterC_Mpl:= proc(
    BEGIN::integer[4], END::integer[4], R::Array(datatype= integer[4])
)
local 
    Year::integer[4], a::integer[4], b::integer[4], c::integer[4], g::integer[4],
    j::integer[4], gjk::integer[4], Month::integer[4], k::integer[4]
;
    for k to END-BEGIN+1 do
        Year:= BEGIN+k-1;
        a:= irem(Year, 19);
        b:= iquo(Year, 100); c:= irem(Year, 100);   
        g:= irem(19*a + b - iquo(b,4) - iquo(8*b+13, 25) + 15, 30); 
        j:= iquo(a + 11*g, 319); 
        gjk:= g - j + irem(2*irem(b,4) + 2*iquo(c,4) - irem(c,4) - g + j + 32, 7);
        Month:= iquo(gjk + 90, 25);
        R[k, 1]:= Month; 
        R[k, 2]:= irem(gjk + Month + 19, 32)
    od;
    return 
end proc
:
EasterC:= Compiler:-Compile(EasterC_Mpl)
:
Easter:= proc(BEGIN::integer, END::integer)
local N:= END-BEGIN+1, E:= Array(1..N, 1..2, datatype= integer[4]);
    EasterC(BEGIN, END, E);
    ArrayTools:-Alias(E, [BEGIN..END, 1..2])
end proc
:
E:= CodeTools:-Usage(Easter(1, 57*10^5)):
memory used=43.49MiB, alloc change=42.48MiB, 
cpu time=1.44s, real time=1.43s, gc time=15.62ms

MD:= md-> local r; [iquo(trunc(lhs(md)),100,r), r] = rhs(md):
T:= (MD~@sort)(Statistics:-Tally(E[..,1]*100+~E[..,2]), key= lhs);
   T := [
     [3, 22] = 27550, [3, 23] = 54150, [3, 24] = 81225, 
     [3, 25] = 110200, [3, 26] = 133000, [3, 27] = 165300, 
     [3, 28] = 186200, [3, 29] = 192850, [3, 30] = 189525, 
     [3, 31] = 189525, [4, 1] = 192850, [4, 2] = 186200, 
     [4, 3] = 192850, [4, 4] = 186200, [4, 5] = 192850, 
     [4, 6] = 189525, [4, 7] = 189525, [4, 8] = 192850, 
     [4, 9] = 186200, [4, 10] = 192850, [4, 11] = 186200, 
     [4, 12] = 192850, [4, 13] = 189525, [4, 14] = 189525, 
     [4, 15] = 192850, [4, 16] = 186200, [4, 17] = 192850, 
     [4, 18] = 197400, [4, 19] = 220400, [4, 20] = 189525, 
     [4, 21] = 162450, [4, 22] = 137750, [4, 23] = 106400, 
     [4, 24] = 82650, [4, 25] = 42000
  ]

 

@RezaZanjirani I don't see your attached code.

@2cUniverse I guess that I need to rephrase my last Reply above because it seems like you either didn't read it or didn't understand it: All fundamental periods are divisors of the totient. The totient is very easy to compute from the prime factorization. Prime factorizations are easy to compute for numbers on the order that you're considering.

The totient itself will only be a fundamental period for x = 2, 4, p^k, 2*p^k for odd prime p. In these cases, the bases for which the totient is a fundamental period are called primitive roots. You might benefit from this Wikipedia article: https://en.wikipedia.org/wiki/Primitive_root_modulo_n

@Carl Love The possibilities for the period are the divisors (NumberTheory:-Divisors) of the totient (NumberTheory:-Totient) of 78; but not every divisor is necessarily a fundamental period. In this case, the totient is 24, but and 24 are not fundamental periods.

First 45 46 47 48 49 50 51 Last Page 47 of 709