Carl Love

Carl Love

28015 Reputation

25 Badges

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

MaplePrimes Activity


These are replies submitted by Carl Love

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.

@dharr Using null indices also works for copying:

prev[]:= curr[];

This works for any rtables, regardless of the number of dimensions (as long as it's the same on both sides) and regardless of the starting indices.

@SitecoreQedge This AI-generated Answer is much too superficial (although mostly correct), but it does have an outright error towards the end.

@Traruh Synred For me to do anything to help you with this, you need to post a complete and executed worksheet showing the occurrence of the error, and containing all the code needed to execute the worksheet.

Was it truly your goal to combine your original three 3-vectors into a single 9-vector, as shown? Or did you intend to combine them into a 3x3 matrix? Although there's nothing incorrect about the single-vector goal, the matrix goal seems a more-normal option, and it would also simplify the syntax needed to retrieve the original vectors.

@KIRAN SAJJAN 

The appliable module (essentially the same thing as a procedureSolve is the heart of the code that you're using---the thing that accomplishes my titular goal of "Numerically solving BVPs with many parameters", with an ease-of-use similar to dsolve's parameters option for IVPs, albeit not with the same efficiency. Also, Solve passively accumulates a table matching the numeric values of the input parameters with the numeric values of the user's desired output parameters. Output parameters could be anything computed using values from the solution curves, such as the Nusselt number or the skin-friction coefficient.

There is a small "compile-time" syntax error in the version of Solve that you have, which means that Solve isn't being used at all! I don't know when it became an error. In Solve, look for the line
    ModuleLoad:= eval(Init)
and change it to
    ModuleLoad:= ()-> Init()

Now, start with smaller goals: Try to make some plots of solution curves for specific parameter values rather than full-blown parameter vs. parameter plots. These problems are extremely complicated from the perspective of numeric computation, and there are almost always parameter combinations that have convergence problems. Sometimes those convergence problems are relatively easy to fix, and sometimes not.

@vv I vote up.

Would you please explain some of the theory behind your algorithm? In particular, I note that using your default R = 1000, there are about 4 million loop iterations (a reasonably small number), but it seems to reach the correct minimum even when the search space is far larger than 4 million.

I don't need elaborate details about LLL, just why it works well in this case.

Would it be possible to heuristically automatically adjust R within your algorithm based on the magnitude of LCM(m)?

First 35 36 37 38 39 40 41 Last Page 37 of 708