Carl Love

Carl Love

28050 Reputation

25 Badges

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

MaplePrimes Activity


These are answers submitted by Carl Love

I admire your curiosity, but you're asking the wrong questions. When a Maple program crashes merely because of the size of the input, it's usually for a very mundane reason: You've run out of memory. In the case of 10^(10^9)+2, it's computing the number itself that causes the crash before isprime is even called. You can verify that by simply entering 10^(10^9)+2 at a command prompt. So your experiment doesn't reveal anything interesting about isprime.

It's possible to trap errors with the try command. See ?try. Of course, if Maple just crashes, that's not going to help. But that's what's going to happen when you push the memory limits.

You can write the current value to a file many different ways  in Maple. The attached worksheet shows one way at the end.

Now, since you seem curious about these things, I present a much better way to test isprime and much more interesting questions to ask about it. It's all in the worksheet below.

Experimental Computational Number Theory: An experiment to analyze the time complexity of isprime

Written by Carl Love, 2016-Sep-21

restart:


#Initialization:
W:= kernelopts(wordsize);

64


#Adjustable parameters:
EL:= "Exp. #1: Two-word factors": #just a name/label for this experiment

BitsLowerB:= W:  BitsUpperB:= 2*W:  #Experiment with adjusting bounds.

TimeLimit:= 10*60:  #10 minutes; adjust as you wish.
MemLimit:= 2^31:  #2 Gig; adjust to an amount safe for your computer.
  


RunIt:= proc()
local
     N:= 1, st:= time(), RawData:= table(), T,
     RandFactor:= nextprime@rand(map2(Scale2, 1, BitsLowerB..BitsUpperB))
;
     while time() - st < TimeLimit and kernelopts(bytesalloc) < MemLimit do
          N:= N*RandFactor();
          T:= time(isprime(N));
          #Short times (under 8 msec) are reported as 0. This'll only happen near the
          #begining. We just ignore these because they'll mess up the subsequent analysis.
          if T > 0 then RawData[T][ilog2(N)/W]:= () end if
          #Change () to NULL if using 2-D input. ^^
     end do;
     RawData
end proc:

RawData:= RunIt():


#Compensate for time-measurement granularity/resolution: If multiple word counts map to the
#same time, then we take the midrange of those word counts.

FitData:= Matrix(([rhs/~2,lhs]~@op)((add@[min,max]@lhs~@op)~(RawData)), datatype= float[8]);
#                                    ^^^
#Change `add` to `+`@op if using a version older than Maple 2015.

FitData := Matrix(201, 2, {(1, 1) = 502.53125, (1, 2) = 5.359, (2, 1) = 393.6875, (2, 2) = 2.937, (3, 1) = 375.78125, (3, 2) = 2.64, (4, 1) = 476.71875, (4, 2) = 4.64, (5, 1) = 239.4375, (5, 2) = .859, (6, 1) = 387.71875, (6, 2) = 2.812, (7, 1) = 417.421875, (7, 2) = 3.375, (8, 1) = 494.5625, (8, 2) = 4.984, (9, 1) = 532.25, (9, 2) = 6.093, (10, 1) = 413.484375, (10, 2) = 3.25, (11, 1) = 524.296875, (11, 2) = 5.921, (12, 1) = 484.640625, (12, 2) = 4.843, (13, 1) = 458.96875, (13, 2) = 4.265, (14, 1) = 550.0625, (14, 2) = 6.578, (15, 1) = 460.953125, (15, 2) = 4.296, (16, 1) = 498.546875, (16, 2) = 5.203, (17, 1) = 558.0, (17, 2) = 6.796, (18, 1) = 201.859375, (18, 2) = .546, (19, 1) = 472.75, (19, 2) = 4.578, (20, 1) = 542.09375, (20, 2) = 6.343, (21, 1) = 291.8828125, (21, 2) = 1.406, (22, 1) = 307.6953125, (22, 2) = 1.593, (23, 1) = 268.1015625, (23, 2) = 1.14, (24, 1) = 437.171875, (24, 2) = 3.765, (25, 1) = 204.8359375, (25, 2) = .562, (26, 1) = 451.078125, (26, 2) = 4.078, (27, 1) = 506.5, (27, 2) = 5.453, (28, 1) = 260.25, (28, 2) = 1.062, (29, 1) = 469.7890625, (29, 2) = 4.484, (30, 1) = 510.46875, (30, 2) = 5.515, (31, 1) = 423.3125, (31, 2) = 3.484, (32, 1) = 278.96875, (32, 2) = 1.234, (33, 1) = 411.484375, (33, 2) = 3.234, (34, 1) = 419.375, (34, 2) = 3.437, (35, 1) = 207.828125, (35, 2) = .593, (36, 1) = 389.71875, (36, 2) = 2.875, (37, 1) = 530.265625, (37, 2) = 6.046, (38, 1) = 157.2265625, (38, 2) = .296, (39, 1) = 354.0625, (39, 2) = 2.312, (40, 1) = 288.890625, (40, 2) = 1.375, (41, 1) = 350.109375, (41, 2) = 2.218, (42, 1) = 116.8359375, (42, 2) = .14, (43, 1) = 500.53125, (43, 2) = 5.343, (44, 1) = 544.09375, (44, 2) = 6.468, (45, 1) = 191.90625, (45, 2) = .484, (46, 1) = 296.859375, (46, 2) = 1.484, (47, 1) = 314.65625, (47, 2) = 1.687, (48, 1) = 431.234375, (48, 2) = 3.671, (49, 1) = 559.984375, (49, 2) = 6.843, (50, 1) = 304.703125, (50, 2) = 1.562, (51, 1) = 142.390625, (51, 2) = .234, (52, 1) = 391.703125, (52, 2) = 2.921, (53, 1) = 567.9375, (53, 2) = 7.062, (54, 1) = 555.0, (54, 2) = 6.765, (55, 1) = 405.53125, (55, 2) = 3.14, (56, 1) = 344.171875, (56, 2) = 2.14, (57, 1) = 254.3515625, (57, 2) = .984, (58, 1) = 565.953125, (58, 2) = 7.046, (59, 1) = 518.34375, (59, 2) = 5.734, (60, 1) = 397.625, (60, 2) = 3.0, (61, 1) = 361.9375, (61, 2) = 2.39, (62, 1) = 421.328125, (62, 2) = 3.468, (63, 1) = 491.6015625, (63, 2) = 4.968, (64, 1) = 373.8125, (64, 2) = 2.609, (65, 1) = 383.734375, (65, 2) = 2.765, (66, 1) = 106.0078125, (66, 2) = .109, (67, 1) = 128.5390625, (67, 2) = .171, (68, 1) = 546.078125, (68, 2) = 6.5, (69, 1) = 168.109375, (69, 2) = .359, (70, 1) = 112.9140625, (70, 2) = .125, (71, 1) = 357.0234375, (71, 2) = 2.343, (72, 1) = 215.734375, (72, 2) = .656, (73, 1) = 340.203125, (73, 2) = 2.093, (74, 1) = 520.328125, (74, 2) = 5.828, (75, 1) = 299.7578125, (75, 2) = 1.5, (76, 1) = 189.921875, (76, 2) = .468, (77, 1) = 528.265625, (77, 2) = 6.031, (78, 1) = 371.828125, (78, 2) = 2.593, (79, 1) = 433.203125, (79, 2) = 3.718, (80, 1) = 443.125, (80, 2) = 3.859, (81, 1) = 453.046875, (81, 2) = 4.187, (82, 1) = 180.96875, (82, 2) = .421, (83, 1) = 488.609375, (83, 2) = 4.906, (84, 1) = 486.625, (84, 2) = 4.875, (85, 1) = 163.15625, (85, 2) = .328, (86, 1) = 276.046875, (86, 2) = 1.218, (87, 1) = 464.890625, (87, 2) = 4.406, (88, 1) = 223.65625, (88, 2) = .718, (89, 1) = 213.765625, (89, 2) = .64, (90, 1) = 514.375, (90, 2) = 5.687, (91, 1) = 445.109375, (91, 2) = 3.875, (92, 1) = 282.953125, (92, 2) = 1.281, (93, 1) = 241.4375, (93, 2) = .875, (94, 1) = 280.953125, (94, 2) = 1.265, (95, 1) = 478.71875, (95, 2) = 4.687, (96, 1) = 569.90625, (96, 2) = 7.203, (97, 1) = 399.609375, (97, 2) = 3.015, (98, 1) = 407.515625, (98, 2) = 3.156, (99, 1) = 428.2734375, (99, 2) = 3.593, (100, 1) = 466.859375, (100, 2) = 4.453, (101, 1) = 330.359375, (101, 2) = 1.921, (102, 1) = 516.359375, (102, 2) = 5.718, (103, 1) = 294.875, (103, 2) = 1.437, (104, 1) = 381.734375, (104, 2) = 2.718, (105, 1) = 139.421875, (105, 2) = .218, (106, 1) = 217.703125, (106, 2) = .671, (107, 1) = 395.640625, (107, 2) = 2.984, (108, 1) = 251.375, (108, 2) = .968, (109, 1) = 302.703125, (109, 2) = 1.546, (110, 1) = 326.46875, (110, 2) = 1.828, (111, 1) = 209.78125, (111, 2) = .609, (112, 1) = 249.390625, (112, 2) = .953, (113, 1) = 482.640625, (113, 2) = 4.828, (114, 1) = 150.28125, (114, 2) = .265, (115, 1) = 153.2734375, (115, 2) = .281, (116, 1) = 318.5625, (116, 2) = 1.75, (117, 1) = 328.453125, (117, 2) = 1.875, (118, 1) = 246.40625, (118, 2) = .921, (119, 1) = 310.6875, (119, 2) = 1.656, (120, 1) = 352.09375, (120, 2) = 2.265, (121, 1) = 228.6015625, (121, 2) = .765, (122, 1) = 176.9921875, (122, 2) = .406, (123, 1) = 334.296875, (123, 2) = 1.984, (124, 1) = 134.4609375, (124, 2) = .203, (125, 1) = 526.28125, (125, 2) = 5.937, (126, 1) = 322.5, (126, 2) = 1.796, (127, 1) = 552.03125, (127, 2) = 6.734, (128, 1) = 324.484375, (128, 2) = 1.843, (129, 1) = 447.109375, (129, 2) = 3.921, (130, 1) = 193.890625, (130, 2) = .5, (131, 1) = 243.4375, (131, 2) = .906, (132, 1) = 540.125, (132, 2) = 6.281, (133, 1) = 348.125, (133, 2) = 2.187, (134, 1) = 160.21875, (134, 2) = .312, (135, 1) = 312.6875, (135, 2) = 1.671, (136, 1) = 316.609375, (136, 2) = 1.718, (137, 1) = 359.984375, (137, 2) = 2.375, (138, 1) = 377.75, (138, 2) = 2.671, (139, 1) = 439.15625, (139, 2) = 3.812, (140, 1) = 403.59375, (140, 2) = 3.109, (141, 1) = 455.03125, (141, 2) = 4.156, (142, 1) = 462.921875, (142, 2) = 4.312, (143, 1) = 198.8671875, (143, 2) = .531, (144, 1) = 220.6640625, (144, 2) = .703, (145, 1) = 457.0, (145, 2) = 4.234, (146, 1) = 534.234375, (146, 2) = 6.156, (147, 1) = 265.171875, (147, 2) = 1.109, (148, 1) = 379.75, (148, 2) = 2.703, (149, 1) = 231.578125, (149, 2) = .781, (150, 1) = 480.671875, (150, 2) = 4.781, (151, 1) = 257.328125, (151, 2) = 1.031, (152, 1) = 263.203125, (152, 2) = 1.093, (153, 1) = 365.890625, (153, 2) = 2.453, (154, 1) = 337.28125, (154, 2) = 2.062, (155, 1) = 183.953125, (155, 2) = .437, (156, 1) = 237.5, (156, 2) = .843, (157, 1) = 271.0625, (157, 2) = 1.156, (158, 1) = 121.7109375, (158, 2) = .156, (159, 1) = 320.5625, (159, 2) = 1.781, (160, 1) = 385.71875, (160, 2) = 2.828, (161, 1) = 522.328125, (161, 2) = 5.843, (162, 1) = 538.203125, (162, 2) = 6.265, (163, 1) = 97.15625, (163, 2) = 0.93e-1, (164, 1) = 441.125, (164, 2) = 3.843, (165, 1) = 35.6875, (165, 2) = 0.15e-1, (166, 1) = 346.140625, (166, 2) = 2.156, (167, 1) = 131.4921875, (167, 2) = .187, (168, 1) = 425.28125, (168, 2) = 3.5, (169, 1) = 474.734375, (169, 2) = 4.562, (170, 1) = 235.515625, (170, 2) = .828, (171, 1) = 225.640625, (171, 2) = .734, (172, 1) = 285.9375, (172, 2) = 1.328, (173, 1) = 504.515625, (173, 2) = 5.406, (174, 1) = 88.28125, (174, 2) = 0.78e-1, (175, 1) = 211.765625, (175, 2) = .625, (176, 1) = 512.4375, (176, 2) = 5.64, (177, 1) = 496.5625, (177, 2) = 5.125, (178, 1) = 166.125, (178, 2) = .343, (179, 1) = 435.203125, (179, 2) = 3.781, (180, 1) = 449.09375, (180, 2) = 4.046, (181, 1) = 508.484375, (181, 2) = 5.5, (182, 1) = 548.078125, (182, 2) = 6.515, (183, 1) = 76.390625, (183, 2) = 0.46e-1, (184, 1) = 409.5, (184, 2) = 3.203, (185, 1) = 146.3359375, (185, 2) = .25, (186, 1) = 363.921875, (186, 2) = 2.421, (187, 1) = 369.84375, (187, 2) = 2.562, (188, 1) = 63.46875, (188, 2) = 0.31e-1, (189, 1) = 172.0390625, (189, 2) = .375, (190, 1) = 415.421875, (190, 2) = 3.328, (191, 1) = 186.9375, (191, 2) = .453, (192, 1) = 233.546875, (192, 2) = .796, (193, 1) = 367.875, (193, 2) = 2.515, (194, 1) = 195.875, (194, 2) = .515, (195, 1) = 562.9765625, (195, 2) = 7.0, (196, 1) = 401.609375, (196, 2) = 3.046, (197, 1) = 80.3515625, (197, 2) = 0.62e-1, (198, 1) = 273.046875, (198, 2) = 1.203, (199, 1) = 342.171875, (199, 2) = 2.109, (200, 1) = 332.3125, (200, 2) = 1.953, (201, 1) = 536.21875, (201, 2) = 6.218}, datatype = float[8])


Model:= Statistics:-PowerFit(FitData, words, output= solutionmodule):


T:= Model:-Results("leastsquaresfunction");

HFloat(1.2921379147178833e-6)*words^HFloat(2.447839189637252)


plots:-display([
          plot(FitData, style= point, symbolsize= 3, symbol= solidcircle),
          plot(T, words= (min..max)(FitData[..,1]), color= blue, thickness= 0)
     ],        
     labels= ["size of N (words)", "time for isprime(N) (seconds)"],
     labeldirections= [HORIZONTAL, VERTICAL],
     labelfont= [TIMES, ROMAN, 12],
     title= EL,
     caption= sprintf(
          "\n              Order = %f, R-squared = %f\n",
          Model:-Results("parametervalues")[2], Model:-Results("rsquared")
     )
);


ModelPlot:= %:

 

The most interesting thing to me about this whole experiment is the exponent on words in T (which I call the Order in the caption of the plot). Notice the extrememly high correlation (called R-squared in the caption). This exponent is a key to understanding the time complexity of the isprime algorithm. (The coefficient in T is totally irrelevant. (Why?)) How does changing the way that the experiment is conducted change this exponent? What if N is prime? What if N is the product of two large primes rather than many two-word primes? Come up with some more interesting modifications of the experiment that may change the exponent. And perhaps the correlation may not be high, in which case something other than Statistics:-PowerFit will be needed.

 


#Save everything as one record in a file. Thus we can accumulate a database of timing
#data.
Description:=
     "The times required to run isprime(N), where N has unique strictly two-word prime \n"
     "factors, were collected for succesively larger values of N. The times were \n"
     "modelled by a power function of the word length of N."
:
fd:= FileTools:-Text:-Open(cat(kernelopts(homedir), "/isprime timing data.txt"), append):
fprintf(
     fd, "%a\n",
     Record(
          #Change "garbage" to "good" if it's one that you really want to save.
          'RecordType'::string= "garbage",
          'Label'::string= EL,
          'Description'::string= Description,
          'TimeStamp'::string= StringTools:-FormatTime("%Y-%m-%d %T"),
          'Proc'::procedure= eval(RunIt),
          'RawData'= eval(RawData),
          'FitData'::Matrix= FitData,
          'Model'= eval(Model),
          'ModelPlot'= ModelPlot
     )
):    
fclose(fd):

 

Download isprime_complexity.mw

Now, if you want to continue this line of research, here's what to do:

1. By using Maple's help system and by asking Questions here if necessary, figure what the above program does. Write a brief paragraph explaining it.

2. By doing some research on the Internet, figure out why what the above program does is interesting, and write a brief paragraph about it. I'd start by looking up "time complexity" and "primality test" in Wikepedia.

3. Answer the questions that I pose at the end of the worksheet.

4. Come up with some similar questions to ask, and answer them.

If you're unsure how to modify the program to address any of the questions, feel free to ask.

You don't show the erroneous command in the worksheet. Thus I am guessing that it is

coeffs(%, [u[t], u[x], u[x, t], u[x, x]]);

In the future, please explicitly show the error message and the command that generates it.

The problem with the above is that the expression contains g(u[x]). Since it also contains g[u[x]], I'm guessing that you intended the former to be the latter.

Change the command solve(sys) to solve(sys, explicit). Then the results will be expressed as square-root expressions rather than RootOf.

The easiest change that you can make to correct the code is to change vector to Vector (uppercase V). The vector, with a lowercase v, is from a much older version of Maple. There's no reason to use the old form.

Yes, it's very easy. Wherever 10 appears in the format string, replace it with an asterisk *. Then in the sequence of values that fill formats, use decimals, value1, decimals, value2. This is IMO more robust than Acer's method and also consistent with standard C-language usage.  However, it does produce more-cryptic code. 

The observed phenomenon has nothing to do with tilde or map. The difference is between `convert/string`(...and convert(..., string) when those are applied an object of type symbol that contains nonstandard characters and is thus enclosed in backquotes. The conversion from symbol to string is a primitive, kernel operation when done as convert(..., string); it doesn't involve a call to the Maple-level procedure `convert/string`. When that procedure is called, the extra quotes are used. 

The proper use of map is map(convert, L, string). The proper use of tilde is convert~(L, string). They produce identical results. Clearly, the tilde form is more compact and produces easier-to-read code.

It's not ~y; the relevant strings are *~ and /~. These are elementwise operators. See ?elementwise. For example, suppose the A and B are sets or lists of the same size or rtables of the same size and shape. Then A *~ B is equivalent to zip(`*`, A, B), i.e., it produces a set, list, or rtable whose entries are the products of the corresponding elements of A and B.

Set interface(prettyprint= 2, labelling= true) and it'll work.

Here's an overload of eval that makes it smarter. My intention is that it's smart enough to convert any bound variable used inside a procedure to a local.

restart:
local eval;

Warning....
eval:= proc(e::uneval, V::{posint, name= anything, {list,set}(name= anything)})
local r, v;
     if nargs=1 then :-eval(e)
     elif V::'posint' then :-eval(e,V)
     else
          v:= {`if`(V::`=`, V, V[])};
          r:= :-eval(
               evalindets(
                    e,
                    function,
                    f-> subs(map(x-> lhs(x)= convert(lhs(x), `local`), v), :-eval(op(0,f)))(op(f))
               ),
               v
          );         
          if procname::'indexed' and op(procname) = 'recurse' then
              r:= `if`(:-eval(r = e), :-eval(r), eval['recurse'](:-eval(r),v))
          else r:= :-eval(r)
          end if;
          `if`(op(0,r) = :-eval, 'eval'(op(r)), r)   
     end if
end proc:         

Now, your examples:

g:= a-> int(x+a, x= a..2*a):
:-eval(g(x)=g(z), [x,z]=~ 1);  #without Carl's correction

     3 = 5/2
eval(g(x)=g(z), [x,z]=~ 1);  #with Carl's correction
     5/2 = 5/2
(:-eval = eval)(g(x), x= 2*x);
     12*x^2 = 10*x^2

g:= a-> int(sin(sin(x+a)), x= a..2*a):
evalf(:-eval(g(x)=g(z), [x,z]=~ 1));  #without Carl's correction
   
     .105207050696364 = .529440545264378
evalf(eval(g(x)=g(z), [x,z]=~ 1));  #with Carl's correction
     .529440545264378 = .529440545264378

And here are some examples for Edgardo to consider:

g:= a-> sum(f(x+a), x= a..2*a):
:-eval(g(x)=g(z), [x,z]=~ 1);  #without Carl's correction

     f(2)+f(4) = f(2)+f(3)
eval(g(x)=g(z), [x,z]=~ 1);  #with Carl's correction
     f(2)+f(3) = f(2)+f(3)

Physics:-Setup(redefinesum= true):
:-eval(g(x)=g(z), [x,z]=~ 1);  #without Carl's correction

     sum(f(2*x), x = x .. 2*x) = f(2)+f(3)
eval(g(x)=g(z), [x,z]=~ 1);  #with Carl's correction
    
f(2)+f(3) = f(2)+f(3)

[Edit 1: Code corrected and example added due to VV's comment.]
[Edit 2: Code updated to work with eval's recurse option and to return unevaluated when appropriate.]

What you're calling a chimera is technically known to logicians and computer scientists as a "bound variable". The command depends and the types dependent and freeof are used to detect them or their absence. However, those commands are irrelevant to your dsolve question because there's an easy way to deal with that: option known. Here's an example:

alpha:= proc(x)
local t;
     if not x::numeric then return 'procname'(args) end if;
     evalf(Int(sin(t^3)*cos(t^2), t= 0..x))
end proc:

Sol:= dsolve({diff(y(x), x) = alpha(x)*y(x), y(0)=1}, numeric, known= [alpha]);

Note that the procedure alpha must contain the line that makes it return unevaluated when it gets non-numeric input.

Please don't post Questions as Posts. I already moved this one.

The letter e isn't a predefined value in Maple. To get e^5, you need to use exp(5).

The problem is that (unfortunately) Heaviside(a::algebraic) evaluates to something other than itself: It evaluates to Heaviside(a)! My preference would be that it return an error message. Perhaps this should be considered a bug of Heaviside. We can quickly look through the code showstat(Heaviside) and see where things go wrong. The argument a::algebraic passes through the type check in the proc line because an expression of the form a::b is of type algebraic! (I think that this should be considered the real bug.) Then between lines 23 and 24, the argument passes the type check linear(a). (This should also be considered a bug.) Then in line 30, the return value is set to Heaviside(a). Many other standard procedures, such as sin, return unevaluated when passed a::b.

All of the bugs mentioned above would be fixed if type(a::b, algebraic) returned false. According to ?type,algebraic, it should return false because `::` is not one of the 13 types listed, and a::b typechecks as false against all 13.

Now, turning to your patmatch command: patmatch(Heaviside(x), Heaviside(a)) returns false because the arguments aren't identical. There's no real pattern in the second argument because there's no ::. A perfect identity match is expected for any parts of the pattern for which there's no ::. This also explains why it returns true when you change a to x.

The difference between patmatch and typematch is that patmatch accepts as its second argument patterns which are far more complex than those allowed by typematch. The latter basically only allows type expressions prepended with name::. Examples:

patmatch(x+2, a::name + n::posint) returns true.
typematch(x+2, a::name + n::posint) returns an error message because name+posint isn't a valid type.
typematch(Heaviside(x), Heaviside(x)) returns an error message because x isn't a valid type.

I find Christian's explanation of the difference quite confusing and misleading, although it's technically correct, I guess.

 

If you use lprint on the results, you'll see that they are essentially the same. A 1-D Array whose indices start at 1 prettyprints as if it were a Vector, and 2-D Array both of whose indices start at 1 prettyprints as if it were a Matrix. Using lprint reveals the true nature of things.

Basically, prettyprinting allows almost anything to be displayed as if it were something else without changing the underlying value. If I wanted sin(1) to display as 42, I could write a one-line procedure to do that. Of course, prettyprinting is only useful when it shows some useful mathematical representation of the underlying object. Since a 1-D Array with indices starting at 1 is essentially a Vector, that's appropriate here.

Your code above the fsolve does not define, or even use at all, any variables named Eq1 or Eq2, yet you refer to them in the fsolve command.

See ?Statistics,Lowess. One of the examples shows doing an integration.

First 210 211 212 213 214 215 216 Last Page 212 of 395