Carl Love

Carl Love

26488 Reputation

25 Badges

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

MaplePrimes Activity


These are replies submitted by Carl Love

I think that in your first line, you tried to enter a formula for h, but it got replaced by characters that look like rectangles with Xs in them.

@zenterix Using string constants as keys is robust and enhances code readabiliy. It seems unlikely to me that anyone with even a tiny knowledge of any standard programming language would think that "chevron2A" could have any value other than the string constant that it manifestly appears to be. And, unlike a variable, it's impossible for you to accidentally assign it a value, which would damage its use as a key.

Gradually over the years, more-and-more Maple library code has been using string constants as keys and keyword values.

In this case, your keys (aka indices) such as chevron2A are the escaped locals, and trying to use them at the worksheet level causes essentially the same problem as described in my 2nd paragraph above. As @acer said, making them module exports ameliorates the problem. (Indeed, module exports are essentially locals designed to be used in escaped form.) However, from the code that you've shown, I see no reason that they should be variables at all, whether local or export. Instead I recommend that they be strings, for example "chevron2A". This would correct all the problems that you've mentioned. 

I love Halloween and pumpkin carving, so I thank you for this, and I vote up. I just have a small numeric nit to pick; we are mathematicians afterall! Halloween/Samhain is the beginning of the darkest quarter of the year (north of the Equator), not the "darker half". The celebration of the end of that quarter is known as Groundhog Day / Imbolc / Candlemass, etc.

@Carl Love I streamlined the above code down to 26 lines and added over 50 lines of comments. Let me know if there's any question about how this code works or about why I chose any particular technique.

(* A161786__2
This procedure returns the primes from ithprime(m) to ithprime(m+n-1), inclusive, that hav exactly D
copies of any "digit" in their base-b representation. Henceforth "digit" always means base-b digit.

Defaults based on  _OEIS_ "A161786":
	The base b defaults to 10, but can be any nonunit positive integer.
	The critical digit count D defaults to 4, but can be any positive integer.

My method is
	1. Use simple combinatorial tools from the Iterator package to construct a sequence of
	digits with the desired count.
	2. Use simple arithmetic to compose the sequence into an integer.
	3. If the integer is in the desired range and prime, save it for return.

This code doesn't use strings or any other characater-based representations in any way.
It doesn't break down integers into digits; rather, it builds large integers from digits.
*)
A161786__2:= proc(m::posint, n::posint, b::And(posint, Not(1)):= 10, D::posint:= 4, $)
option `Author: Carl Love <carl.j.love@gmail.com> 2023-Oct-13`;
uses It= Iterator, NT= NumberTheory;
local  #*1 
	(p0,p1):= op(ithprime~([m, m+n-1])), ds:= ilog[b](p1), D1:= ds+1-D, bs:= {$0..b-1},
	B:= Array(0..ds, k-> b^k), rep1:= add(B), C, BC, rep, d, Rep, k, n1,
	d0:= select(x-> igcd(x,b)=1, bs), mids:= bs$D1-2, d1:= {$map(iquo, p0..p1, B[ds])},
	Prime?:= eval(`if`(p1 < 5*10^9, gmp_isprime, prime)),
	Accept:= p-> if p0 <= p and p <= p1 and Prime?(p) then p fi,
	DigitSets:= rcurry(op@eval, d= ()) @ 
		`if`(D1=1, ()-> [`if`(C[1]=0, d0, bs)], ()-> [`if`(C[1]=0, d0, bs), mids, `if`(C[D1]=ds, d1, bs)])
;
	if D1 <= 0 then `if`(D1 < 0, {}, Accept~(`if`(D>1, {rep1}, bs)))
	elif 0 in d1 then n1:= NT:-pi(B[ds]); thisproc~([m, n1+1], [n1-m, m+n-n1], b, D)
	else  # main loop:
		{for C in It:-Combination(ds+1, D1) do
			rep:= rep1 - add((BC:= [seq](B[[seq](C)])));
			(for d in `if`(C[1]=0, bs, d0) minus `if`(C[D1] = ds, {}, {0}) do
				Rep:= d*rep;
				(for k in It:-CartesianProduct(DigitSets()) do
					Accept(Rep + inner([seq](k), BC))
				od)
			od)
		od}
	fi[]
end proc

(* Code comments:

Locals:
	Static:
		p0:	smallest acceptable prime;
		p1:	largest acceptable prime;
		ds:	digit count of p1, counting the units digits as "zeroeth";
		D1:	number of digit positions NOT being specially counted;
		bs:	set of possible single digits;
		B:	array of powers of b such that B[k] = b^k;
		rep1: integer with desired number of digits, all of which are 1;
		d0:	set of allowed units digits for a multidigits prime;
		d1:	set of allowed leading digits (possibly including 0) in p1..p2;
		mids: middle digits: sequence of max(D1-2, 0) copies of bs;
		
	Variable:
		C:	iterates over all D1-combinations of digit positions (0-relative)
			(C is returned by Iterator:-Combination as a sorted vector);
		BC:	sorted list of powers of b corresponding to the postions in C;
		rep: integer with D 1-digits and D1 0-digits with the positions of the 0s given by C;
		d:	iterates over the allowed values for the specially counted digits, taking into
			account whether their positions are first and/or last;
		Rep:	d*rep, thus the digits become all D d's and D1 0s in some order;
		k:	iterates over all digit sequences that can replace the 0s in Rep
			(k is returned by Iterator:-CartesianProduct as a vector in digit-position order);
		n1:	new value of n potentially used for recursive call;

	Procedures:
		Prime?:
			Replaces isprime with the more-efficient undocumented builtin gmp_isprime if all
			its inputs will be < 5 billion;
		Accept:
			If its input is in the desired range and prime, return it;
		DigitSets:
			Constructs the sequence of allowed-digit sets to pass to Iterator:-CartesianProduct.
			C and d are unassigned at the time DigitSets is assigned, but they'll be assigned
			when it's called.
			
Body:
	The 1st branch of the if-statement is to handle cases where there are no digits other than those
	being specially counted. It's easier to handle these as separate cases.

	The 2nd branch handles the case where there are fewer digits in p0 than in p1. Since the main loop
	requires that their lengths be equal, this branch uses recursion to split the cases by length.

	Main loop:
		Consider all possible D1-combinations C of digit positions:
			BC is the powers of b corresponding to those positions;
			rep has all digits 1 or 0.
			Consider all possible digit values d for the specially counted value:
				Rep is rep with 1s replaced by d;
				Consider all possible digit combinations k to replace the 0s in Rep:
					`inner` is an undocumented builtin command to compute the inner- (aka dot-)
					product of lists as if they were vectors;
					so Rep + inner(..) is Rep with its 0s replaced by the digits in k. 
*)
:		


 

restart:

(* A161786__2

 

The vast majority of the the time used in this first call is for compiling the Iterators. This only needs to be done once per restart.

CodeTools:-Usage(A161786__2(9, 9, 2, 1));

memory used=27.91MiB, alloc change=12.01MiB, cpu time=265.00ms, real time=6.14s, gc time=0ns

23, 29, 47, 59, 61

CodeTools:-Usage(A161786__2(9, 99, 10, 2));

memory used=0.80MiB, alloc change=0 bytes, cpu time=0ns, real time=9.00ms, gc time=0ns

101, 113, 131, 151, 181, 191, 199, 211, 223, 227, 229, 233, 277, 311, 313, 331, 337, 353, 373, 383, 433, 443, 449, 499, 557, 577

P2:= CodeTools:-Usage(A161786__2(10^6, 10^6)):

memory used=1.33GiB, alloc change=0 bytes, cpu time=12.95s, real time=12.17s, gc time=1.64s

nops([P2]), P2[-1];

42720, 32449441

The computational effort is largely determined by the size of the upper bound (p1 in the code); the size of p0 has relattively little effect. Here we compute the entire sequence up to the two-millionth prime using only slightly more time than above:

P2a:= CodeTools:-Usage(A161786__2(1, 2*10^6)):

memory used=1.48GiB, alloc change=0 bytes, cpu time=14.75s, real time=13.78s, gc time=1.97s

nops([P2a]), P2a[-1];

80331, 32449441

 

Download DigitCountPrimes.mw

@Jamie128 Yes, it can be animated by just making a few small changes to the code above:

R:= 2:  #max r of featured surface 
spokes:= 24:  #number of radial gridlines in base

plots:-animate(
    plot3d, #the feature
    [
        [r, theta, (1+k)*(csc(theta)/r)^2], r= 0..R, theta= -Pi..Pi,
        coords= cylindrical, style= surface, shading= zhue, transparency= 0.15
    ],
    k= -1..1,
    background= [
        plot3d(  #the base
            [r, theta, 0], r= 0..3/2*R, theta= -Pi..Pi, coords= cylindrical,
            grid= [7, 2*spokes+1], color= COLOR(RGB, 0.97$3), thickness= 2.5, glossiness= 0
        ),
        plots:-textplot3d(  #the spoke labels
            [seq](
                [7/4*R, i, 0, i],
                i= (2*Pi/spokes)*~([$0..spokes-1] -~ iquo(spokes,2))
            ), 
            coords= cylindrical, font= ["times", 10]
        )
    ],
    view= [map(`*`, -1..1, 2*R)$2, 0..20], axes= framed, 
    orientation= [40,40,0], labels= [r$2, z], labelfont= ["times", 16], 
    axis[1,2]= [tickmarks= [-R, R]]
);

 

We need to see the code that you entered that produced that error. You can use the green uparrow on the toolbar of the MaplePrimes editor to upload a worksheet.

@sursumCorda Okay, that's a reasonable reason.

@sursumCorda I'm curious why you usually do module references as A[':-B'] rather than A:-B?

@sursumCorda 

[Edit: This comment was directed at your file-deletion version of editing one's libraries.]

You're right, and there's an easier way to confirm it: Just edit the global sequence libname:

restart:
sav_libname:= libname;
  sav_libname := 
    "C:\Users\carlj\maple\toolbox\2023\Physics Updates\lib", 
    "C:\Program Files\Maple 2023\lib", 
    "C:\Users\carlj\maple\toolbox\DirectSearch\lib", 
    "C:\Users\carlj\maple\toolbox\Glyph2\lib", 
    "C:\Users\carlj\maple\toolbox\OrthogonalExpansions\lib"

libname:= select(L-> SearchText("Physics", L) = 0, [libname])[]:
inttrans:-invlaplace(1/(s+a), s, t);
                           exp(-a t)

libname:= sav_libname:
forget(inttrans:-invlaplace);
inttrans:-invlaplace(1/(s+a), s, t);
                            exp(a t)

 

@acer Based on my seat-of-the-pants attempt to intepret that Matlab code, the OP might have this in mind, a slight variation of your densityplot:

plots:-densityplot(
    argument(-exp(1/(x+I*y))), (x,y)=~ -0.5..0.5, grid= [401$2],
    axes= none, colorbar= false, scaling= constrained,               
    colorscheme= ["zgradient", ["Red","Violet"], colorspace= "HSV"]
);

 

@Preben Alsholm I have the bug, and I confirmed that I have the latest Physics Update:

Physics:-Version();
The "Physics Updates" version in the MapleCloud is 1561 and is 
the same as the version installed in this computer, created
2023, October 20, 2:37 hours Pacific Time.


inttrans:-invlaplace(1/(s+a), s, t);
                            exp(a t)

This is probably unrelated, but I'll note that the given time and date is nearly 3 hours in the future.

@sursumCorda I decided that the D=0 case (in other words, Is there any digit that doesn't appear?) was too distinct from the the D > 0 case to handle with this same code, so I simply forbade it.

By making some minor improvements, I got the time to under a quarter minute, indeed, under 12 seconds.

(* A161786__2
This procedure returns the primes from ithprime(m) to ithprime(m+n-1), inclusive, that have
exactly D copies of any "digit" in their base-b representation.
The base b defaults to 10 but can be any nonunit positive integer.
The critical digit count defaults to 4.

My method is
	1. Use simple combinatorial tools from the Iterator package to construct a sequence of
	base-b digits with the desired count.
	2. Use simple arithmetic to compose the sequence into an integer.
	3. If the integer is in the desired range and prime, save it.

This code doesn't use strings or any other characater-based representations in any way.
It doesn't break down integers into digits; rather, it builds large integers from digits.

Comments of the form #*n refer to footnotes ar the end of the procedure.
*)
A161786__2:= proc(m::posint, n::posint, b::And(posint, Not(1)):= 10, D::posint:= 4, $)
option `Author: Carl Love <carl.j.love@gmail.com> 2023-Oct-13`;
uses It= Iterator, NT= NumberTheory;
local 
	(p0,p1):= op(ithprime~([m, m+n-1])), ds:= ilog[b](p1), D1:= ds+1-D, bs:= {$0..b-1},
	B:= Array(0..ds, k-> b^k), rep1:= add(B), rep, Rep, BC, k, d, C, p, mids:= bs$D1-2,
	d0:= select(x-> igcd(x,b)=1, bs), d1:= {$map(iquo, p0..p1, B[ds])},
	FactorSelector:=  #*1
		if D1 = 1 then ()-> `if`(C[1]=0, d0, bs) minus {d}
		else ()-> map(`minus`, [`if`(C[1]=0, d0, bs), mids, `if`(C[D1]=ds, d1, bs)], {d})[]
		fi,
	Isprime:= eval(`if`(p1 < 5*10^9, gmp_isprime, isprime))  #*2 
;
	# Short-circuit returns for corner cases:
	if D1 <= 0 then
		return
		     if D1 < 0 then {}
	    	 	elif D>1 then `if`(p0 <= rep1 and rep1 <= p1 and Isprime(rep1), {rep1}, {})
			else select(p-> p0 <= p and p <= p1 and Isprime(p), bs)
			fi
	fi; 
	if p0 < B[ds] then 
		return thisproc(m, (p:= NT:-pi(B[ds]))-m, b, D) union thisproc(p+1, n+m-p, b, D)
	fi;

	# Main code:
	{
		for C in It:-Combination(ds+1, D1) do
			rep:= rep1 - add((BC:= [seq](B[[seq](C)])));
			(
				for d in `if`(C[1]<>0, d0, bs) minus `if`(C[D1] = ds, {}, {0}) do
					Rep:= d*rep;				
					(
						for k in It:-CartesianProduct(FactorSelector()) do   
							p:= Rep + inner([seq](k), BC);  #*3
							if p0 <= p and p <= p1 and Isprime(p) then p fi
						od
					)
				od
			)
		od
	}
end proc

(* Footnotes:
	*1: FactorSelector builds the factor sets of digits for the CartesianProduct Iterators.
	Variables C and d are not defined at the time this procedure is decalred, but they'll be defined
	before it's called.
	
	*2: gmp_isprime is an undocumented builtin procedure that is faster 
	than isprime for most small cases.

	*3: 'inner' is an undocumented builtin procedure that computes the
	inner or "dot" product of 2 lists as if they were vectors.
*)
:

CodeTools:-Usage(A161786__2(10,10,10,1));
memory used=462.70KiB, alloc change=0 bytes, cpu time=16.00ms, real time=6.00ms, gc time=0ns

            {29, 31, 37, 41, 43, 47, 53, 59, 61, 67}

P2:= CodeTools:-Usage(A161786__2(10^6, 10^6)): 
memory used=1.31GiB, alloc change=0 bytes, cpu time=12.64s, real time=11.44s, gc time=1.88s

nops(P2);
                             42720

Download worksheet: DigitCountPrimes.mw

@mmcdara For the 6-digit numbers, you've selected those that have 5 occurences of a digit rather than the requested "exactly four" occurences.

@Scot Gould The command prev[]:= curr[] copies the elements of curr to the already-existing rtable represented by prev. If prev is not already an rtable, then the command does something else (very likely undesirable). The meaning of the [] when appended to an rtable is "all indices in all dimensions". The command could be replaced by prev[]:= rtable(curr).  (rtable could be replaced by Vector or copy, but these are relatively inefficient.) On the left side of the assignment, the [] is essential, just like [1] would be essential to assign to only the first element of prev.

First 15 16 17 18 19 20 21 Last Page 17 of 688