Carl Love

Carl Love

28055 Reputation

25 Badges

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

MaplePrimes Activity


These are answers submitted by Carl Love

@syhue You should check the help before giving commands that are essentially random garbage. You'll see that read is a command that has nothing to with this, and to_number means nothing at all to Maple.

Here's a complete RSA example, from an original string plaintext, converted to numeric, encrypted, decrypted, and returned to plaintext.

restart:
#I just need any (n,e,d) for this demo; I don't care if they're secure.
(p,q):= 'nextprime(rand())'$2: 
(n,e):= (p*q, 2^16+1):  d:= e^(-1) mod (p-1)*(q-1):
#message in string form:
P:= "Return to headquarters.": 
bby:= 2^8: #base of original bytes 
#base of "chunkified" message:
bch:= bby^iquo(ilog2(n), ilog2(bby)):
#message in numeric "chunk" form:
M:= convert(convert(P, 'bytes'), 'base', bby, bch);
    M := [2055956401727290500434, 2109951468272787744800, 
          199505372532]
C:= `&^`~(M,e) mod n; #Encrypt.
    C := [11266093172382474618748, 59649508547257401648253, 
          57058717428414196529185]
M1:= `&^`~(C,d) mod n; #Decrypt.
    M1 := [2055956401727290500434, 2109951468272787744800, 
           199505372532]
evalb(M1=M); #accuracy check
                              true
#Convert back to string:
P1:= cat(convert~(convert~(M1, 'base', bby), 'bytes')[]);
                P1 := "Return to headquarters."
evalb(P1=P);
                              true

Please let me know if there's anything that you don't understand about the above.

Actually, no loops or seqs are needed, just two lines of code. But there are some problems with your specs! So, I have made the following three assumptions to correct them:

  1. goes from to M-1, not to M-1.
  2. The integral wrt t can be factored out of the integral wrt x.
  3. The ns in the two integrals are different: One corresponds to i and one corresponds to j.

Please confirm that my assumptions are correct! If they are correct, then the matrix can be created by

(k,M):= (2,2):
intpsiomega:= ij-> local t; int(psi[ij](t)*omega[iquo(ij-1,M)+1](t), t= 0..1);
U:= Matrix(2^k/2*M, mul@intpsiomega~@`[]`);

That's coded for Maple 2019. If you're using an earlier Maple, a tiny adjustment will be needed. Let me know.

As I've explained to you before, when doing modular exponentiation you should use &^ instead of ^.  (In 2-D Input mode, enter the operator as &\^ (the backslash won't appear).) If the exponents e or d are as large as would typically be used in cryptography, then you must use &^; using will error out or crash the system.

Compare timings of the two:

C:= CodeTools:-Usage(M^e mod n, iterations= 10^5):
memory used=4.29KiB, alloc change=-3.01MiB, cpu time=17.81us,
real time=17.76us, gc time=8.12us

C:= CodeTools:-Usage(M &^ e mod n, iterations= 10^5):
memory used=104 bytes, alloc change=0 bytes, cpu time=1.09us,
real time=1.01us, gc time=0ns


As the exponent grows, the time of the first method grows exponentially, whereas the second method's time just grows a tiny amount:

C:= CodeTools:-Usage(M^(e^2) mod n, iterations= 10):
memory used=53.89MiB, alloc change=101.21MiB, cpu time=189.10ms,
real time=188.40ms, gc time=0ns

C:= CodeTools:-Usage(M &^ (e^2) mod n, iterations= 10^5):
memory used=144 bytes, alloc change=0 bytes, cpu time=1.25us,
real time=1.27us, gc time=0ns

The time of the first method has increased by a factor of more than 10,000 (note the unit change from microseconds to milliseconds)! But the second method took only 15% longer. And e^2 would still be an extremely small exponent for cryptographic purposes.
 

 

You need to set an artificial finite value for infinity, such as 10.

fsolve needs a bit of help here. I'm not sure why; it may be a precision issue. I didn't investigate because the workaround was obvious to me. Define A = exp(-2*a). This gives the following system of three equations, which fsolve solves immediately:

fsolve({c*A^4+18 = 36, c*A+18 = 55, A= exp(-2*a)}, {A, a, c});
             {A = 0.7864846673, a = 0.1200910258, c = 47.04478236}

Do

Sinus:= evalf@Units:-Standard:-sin:

Then compare:

Sinus(30*Unit(degree));
Sinus(30); #30 radians

Maple is a "symbolic computation system", so those string manipulations are totally unnecessary; just use subs instead.

Your equation can be easily solved with all variables symbolic and using a logarithm to an arbitrary base, B:

restart:
equation1:= filler(a*x+b) - filler(c*x+d) = rs:
equationAct:= subs(filler= log[B], equation1):
simplify(solve(equationAct, x));
             (b - B^rs*d)/(c*B^rs - a)

The condition for both log arguments to be positive is (b*c - a*d)/(c*B^rs - a) > 0.

If you ever again get the urge to use string manipulations to do math, ask here for an easier way.

Your procedure has highly symbolic components such as isolate; it could never beat pure integer arithmetic.

Here's my version of repeated squaring. It's like VV's except that I unrolled his recursion into a do loop. It works for any associative operator (such as square matrix multiplication or modular multiplication), with the default being ordinary `*` multiplication.

RepSq:= proc(X, N::nonnegint, `*`:= :-`*`, Identity:= 1)
local x:= X, R:= Identity, n:= N;   
     do
          if n::odd then R:= `*`(x,R); n:= n-1 fi;
          if n=0 then return R else n:= n/2 fi;
          x:= `*`(x,x)
     od
end proc
: 

 

Some pointers on algorithm timing, referring to mistakes that you made in yours:

1. Don't display output and make a timing at the same time. Suppress the output and display it later if it's needed. For example, note that the number needs to be converted to base 10 to be displayed, and you don't want to be timing that.

2. You need much, much larger examples and many, many more of them to get meaningful timings. It's inconceivable that your puny numbers took 0.04 seconds to exponentiate; all you're measuring is overhead. Try a list of a thousand numbers with exponents in the hundreds of thousands. (Make absolutely sure that you suppress the output in this case!! When the GUI gets its output buffer overloaded, there's nothing that you can do except either wait hours for the output to display or kill the GUI (which unfortunately kills all your open worksheets).)

3. CPU time is a more meaningful comparison than real time for raw algorithms, such as this. Real time is more meaningful for comparing parallelized code or code that does disk I/O operations.

No Maple assumptions are needed for this:

ysol:= piecewise(x < 1/2, 2, x >= 1/2, 0); #1/2, NOT 0.5
​​​​​​diff(ysol, x);

In Maple, you should never use a decimal number such as 0.5 as an abbreviation for an exact fraction such as 1/2

Your code had some minor syntax errors, as corrected by Tom. Your algorithm (as opposed to your code) has two major flaws (not errors): It uses matrix inversion, and it unnecessarily recomputes the symbolic Jacobian on each iteration. Here is my rewrite that avoids those things. Inversion is replaced with solving a linear system, and the Jacobian is computed once before the loop. I also added argument-error checking and made TOL and itmax optional with default values. The solution is returned if convergence (as measured by TOL) is achieved; otherwise an error message is given.

restart
:

NewtonRaphsonSYS:= proc(
   X::list(And(name, Not(constant))),
   F::And(
         list(algebraic),
         satisfies(F->
            nops(F)=nops(X) and
            indets(F, And(name, Not(constant))) subset {X[]}
      )
   ),
   X0::And(list(complexcons), satisfies(X0-> nops(X0)=nops(X))),
   #optional keyword arguments and their default values:
   {  
      TOL::And(float,positive):= 10^(1+nops(X)-Digits),
      itmax::nonnegint:= 99
   }
)
local
   n, f:= <F>, Xn:= <X0>, H:= <seq(TOL, n= 1..nops(X))>,
   J:= VectorCalculus:-Jacobian(f,X)
;
   for n to itmax while LinearAlgebra:-Norm(H,2) >= TOL do
      H:= LinearAlgebra:-LinearSolve(
         rtable(evalf(eval(J, X=~ seq(Xn))), datatype= float),
         rtable(evalf(eval(f, X=~ seq(Xn))), datatype= float)
      );
      Xn:= Xn - H;
      userinfo(1, 'Newton', 'NoName', 'n'=n, 'Xn'=evalf[8]([seq(Xn)]));
      userinfo(1, 'Newton', 'NoName', `      `, 'H'=evalf[3]([seq(H)]))
   end do;
   if n > itmax then error "did not converge" end if;
   userinfo(1, 'Newton', 'NoName', `\n`); #blank line
   [seq(Xn)]
end proc
:   

#Usage example:
F:= [x*y-z^2-1, x*y*z+y^2-x^2-2, exp(x)+z-exp(y)-3];
X:= [x,y,z]:
X0:= [1,1,1]:

[x*y-z^2-1, x*y*z-x^2+y^2-2, exp(x)+z-exp(y)-3]

infolevel[Newton]:= 1: #Turn on userinfo statements.
Sol:= NewtonRaphsonSYS(X, F, X0, TOL= 1e-8, itmax= 10);

n = 1 Xn = [2.1893261 1.5984752 1.3939006]

       H = [-1.19 -.598 -.394]
n = 2 Xn = [1.8505896 1.4442514 1.2782240]
       H = [.339 .154 .116]
n = 3 Xn = [1.7801612 1.4244360 1.2392924]
       H = [0.704e-1 0.198e-1 0.389e-1]
n = 4 Xn = [1.7776747 1.4239609 1.2374738]
       H = [0.249e-2 0.475e-3 0.182e-2]
n = 5 Xn = [1.7776719 1.4239606 1.2374711]
       H = [0.279e-5 0.328e-6 0.270e-5]
n = 6 Xn = [1.7776719 1.4239606 1.2374711]
       H = [0.314e-11 0.422e-13 0.441e-11]

[HFloat(1.7776719180107403), HFloat(1.423960597888489), HFloat(1.2374711177317033)]

#Check residuals:
eval(F, X=~ Sol);

[HFloat(0.0), HFloat(-1.8388068845354155e-16), HFloat(-1.7763568394002505e-15)]

 


 

Download Newton.mw

The ​​​​​​solve command, used by Rouben, says that's there's no unconditional or "free" solutions, which is not suprising given that there are 12 equations for 6 unknowns. The eliminate command is an alternative that gives conditions under which solutions exist. I think that there's a good chance that that's what you're interested in. The following worksheet is meant to be put at the end of Rouben's. 
 

Sols:= eliminate(sys, vars):

These are the expressions that must be equal to 0 in order for the following solution to be valid:

Conds:= <Sols[2][]>;

Vector(6, {(1) = (p[1, 2]-p[1, 3])*alpha, (2) = -(p[2, 2]-p[2, 3])*alpha, (3) = p[1, 1]-p[1, 3], (4) = -p[2, 1]+p[2, 3], (5) = alpha*p[1, 2]-alpha*p[1, 3]-p[1, 1]+p[1, 3], (6) = -alpha*p[2, 2]+alpha*p[2, 3]+p[2, 1]-p[2, 3]})

(1)

And these are the solutions for the q's:

Sol:= < (<lhs~([Sols[1][]])> =~ ``) | <rhs~([Sols[1][]])>>;

"[[[q[1,1]=,-(beta p[2,3]-beta-p[1,3]+1)/((beta^2-1) (alpha+1))],[q[1,2]=,-(alpha beta p[2,2]-alpha beta p[2,3]-alpha p[1,2]+alpha p[1,3]+beta p[2,2]-beta-p[1,2]+1)/((beta^2-1) (alpha+1))],[q[1,3]=,-(beta p[2,3]-beta-p[1,3]+1)/((beta^2-1) (alpha+1))],[q[2,1]=,-(beta p[1,3]-beta-p[2,3]+1)/(alpha beta^2+beta^2-alpha-1)],[q[2,2]=,-(alpha beta p[1,2]-alpha beta p[1,3]-alpha p[2,2]+alpha p[2,3]+beta p[1,2]-beta-p[2,2]+1)/((beta^2-1) (alpha+1))],[q[2,3]=,-(beta p[1,3]-beta-p[2,3]+1)/(alpha beta^2+beta^2-alpha-1)]]]"

(2)

 


 

Download Conditional.mw

Numerous applications, particularly in number theory, require only the integer part of the base-2 logarithm. There is a special command for that, ilog2(6), which is much more efficient than floor(log[2](6)).

This is done with the identity option to solve.

solve(
   identity(
      cos(t)^6+a*cos(t)^4*sin(t)^2+b*cos(t)^2*sin(t)^4+c*sin(t)^6 = cos(6*t), 
      t
   ), 
   {a,b,c}
);
                   {a = -15, b = 15, c = -1}

 

These very small imaginary parts are due to round-off error in float computations. I call them spurious imaginary parts, because the true answers are real. They can be removed by

simplify(fnormal([%]), zero);

In some exceptional cases, you may need to specify a tolerance (a number of digits) as a second argument to fnormal. See help page ?fnormal.

The RealDomain package is not appropriate for this. It's better to generate the answers with the imaginary parts and then remove them than it is to try to generate them without the imaginary parts.

The command to turn an expression into a function is unapply.

Mat:= module()
option package;
uses CF= CurveFitting;
export 
     LinModel:= proc(xliste::list, yliste::list, $)
     local a, b, x;
          unapply(
               CF:-LeastSquares(xliste, yliste, x, 'curve'= a*x+b),
               x
          )
     end proc; # LinModel 
end module; #Mat

If you are returning a function, there's no need to pass x to LinModel. Thus, I've made it a local.

First 126 127 128 129 130 131 132 Last Page 128 of 395