ans = 
Results for
If you haven't solved the problem yet, below hints guide how the algorithm should be implemented and clarify subtle rules that are easy to miss.
1. Shield is ONLY defended in HOME matches of the CURRENT holder - Even if a team beats the Shield holder in an away match, that does NOT count as a Shield defense.
2. A team defends the Shield ONLY when:
> They currently hold it.
> They are home team in that match
3. Shield transfer happens ONLY if the HOLDER plays a home match AND loses - A team may lose an away match — no effect.
4. The output ALWAYS includes the initial holder as the first row.
5. Defenses count resets for each new holder. - Every holder accumulates their own count until they lose it at home.
6. Match numbers are 1-indexed in the input, but “0” is used for initial state - The first real match is Match 1, but the output starts with Match 0.
7. Output row is created ONLY WHEN SHIELD CHANGES HANDS - This is an important hidden detail. A new row is appended, When the current holder loses a home match → Shield taken by visitor. If no loss at home occurs after that → no new row until next change.
8. The last holder’s defense count goes until the season ends - Even if they lose away later.
9. If a holder never gets a home match, defenses = 0.
10. In case the holder loses their very first home match → defenses = 0.
11. Shield changes only on HOME LOSS, not on a draw.
I hope above hints will help you in solving the problem.
Thanks and Regards,
Dev
Many MATLAB Cody problems involve solving congruences, modular inverses, Diophantine equations, or simplifying ratios under constraints. A powerful tool for these tasks is the Extended Euclidean Algorithm (EEA), which not only computes the greatest common divisor, gcd(a,b), but also provides integers x and y such that: a*x + b*y = gcd(a,b) - which is Bezout's identity.
Use of the Extended Euclidean Algorithm is very using in solving many different types of MATLAB Cody problems such as:
- Computing modular inverses safely, even for very large numbers
- Solving linear Diophantine equations
- Simplifing fractions or finding nteger coefficients without using symbolic tools
- Avoiding loops (EEA can be implemented recursively)
Below is a recursive implementation of the EEA.
function [g,x,y] = egcd(a,b)
% a*x + b*y = g [gcd(a,b)]
if b == 0
g = a; x = 1; y = 0;
else
[g, x1, y1] = egcd(b, mod(a,b));
x = y1;
y = x1 - floor(a/b)*y1;
end
end
Problem:
Given integers a and m, return the modular inverse of a (mod m).
If the inverse does not exist, return -1.
function inv = modInverse(a,m)
[g,x,~] = egcd(a,m);
if g ~= 1 % inverse doesn't exist
inv = -1;
else
inv = mod(x,m); % Bézout coefficient gives the inverse
end
end
%find the modular inverse of 19 (mod 5)
inv=modInverse(19,5)
Congratulations to all the Relentless Coders who have completed the problem set. I hope you weren't too busy relentlessly solving problems to enjoy the silliness I put into them.
If you've solved the whole problem set, don't forget to help out your teammates with suggestions, tips, tricks, etc. But also, just for fun, I'm curious to see which of my many in-jokes and nerdy references you noticed. Many of the problems were inspired by things in the real world, then ported over into the chaotic fantasy world of Nedland.
I guess I'll start with the obvious real-world reference: @Ned Gulley (I make no comment about his role as insane despot in any universe, real or otherwise.)
Hi Everyone!
As this is the most difficult question in problem group "Cody Contest 2025". To solve this problem, It is very important to understand all the hidden clues in the problem statement. Because everything is not directly visible.
For those who tried the problem, but were not able to solve. You might have missed any of the below hints -
- “The other players do not get to see which card has been shown, but they do know which three cards were asked for and that the player asked had one of them.” - Even when the card identity isn’t revealed (result = 0), you still gain partial knowledge — the asked player must have at least one of those three cards, meaning you can mark other players as not having all three simultaneously.
- "If it is your turn, you know the exact identity of that card" - You only know the exact shown card when result = 1, 2, or 3 — and it must be your turn. If someone else asked (even if you know result = 0), you don’t know which one was shown. So the meaning of result depends on whose turn it was, which is implicit — MATLAB code must assume that turns alternate 1→m→1, so your turn index is determined by (t-1) mod m + 1 == pnum.
- "Any leftover cards are placed face-up so that all players can see them" - These cards (commoncards) are not in anyone’s hand and cannot be in the envelope. So they’re not just visible — they’re logical constraints to eliminate from deduction.
- “It may be possible to determine the solution from less information than is given, but the information given will always be sufficient.”
- "Turn order is implied, not given explicitly" - Players take turns in order (1 to m, and back to 1).
On considering all the clues and constraints in the question, you will definitely be able to card for each category present in envelope.
I hope above clues will be useful for you.
Thank you, wishing you the success!
Regards,
Dev
When solving Cody problems, sometimes your solution takes too long — especially if you’re recomputing large arrays or iterative sequences every time your function is called.
The Cody work area resets between separate runs of your code, but within one Cody test suite, your function may be called multiple times in a single session.
This is where persistent variables come in handy.
A persistent variable keeps its value between function calls, but only while MATLAB is still running your function suite.
This means:
- You can cache results to avoid recomputation.
- You can accumulate data across multiple calls.
- But it resets when Cody or MATLAB restarts.
Suppose you’re asked to find the n-th Fibonacci number efficiently — Cody may time out if you use recursion naively. Here’s how to use persistent to store computed values:
function f = fibPersistent(n)
import java.math.BigInteger
persistent F
if isempty(F)
F=[BigInteger('0'),BigInteger('1')];
for k=3:10000
F(k)=F(k-1).add(F(k-2));
end
end
% Extend the stored sequence only if needed
while length(F) <= n
F(end+1)=F(end).add(F(end-1));
end
f = char(F(n+1).toString); % since F(1) is really F(0)
end
%calling function 100 times
K=arrayfun(@(x)fibPersistent(x),randi(10000,1,100),'UniformOutput',false);
K(100)
The fzero function can handle extremely messy equations — even those mixing exponentials, trigonometric, and logarithmic terms — provided the function is continuous near the root and you give a reasonable starting point or interval.
It’s ideal for cases like:
- Solving energy balance equations
- Finding intersection points of nonlinear models
- Determining parameters from experimental data
Example: Solving for Equilibrium Temperature in a Heat Radiation-Conduction Model
Suppose a spacecraft component exchanges heat via conduction and radiation with its environment. At steady state, the power generated internally equals the heat lost:
Given constants:
= 25 W- k = 0.5 W/K
- ϵ = 0.8
- σ = 5.67e−8 W/m²K⁴
- A = 0.1 m²
= 250 K
Find the steady-state temperature, T.
% Given constants
Qgen = 25;
k = 0.5;
eps = 0.8;
sigma = 5.67e-8;
A = 0.1;
Tinf = 250;
% Define the energy balance equation (set equal to zero)
f = @(T) Qgen - (k*(T - Tinf) + eps*sigma*A*(T.^4 - Tinf^4));
% Plot for a sense of where the root lies before implementing
fplot(f, [250 300]); grid on
xlabel('Temperature (K)'); ylabel('f(T)')
title('Energy Balance: Root corresponds to steady-state temperature')
% Use fzero with an interval that brackets the root
T_eq = fzero(f, [250 300]);
fprintf('Steady-state temperature: %.2f K\n', T_eq);
It’s exciting to dive into a new dataset full of unfamiliar variables but it can also be overwhelming if you’re not sure where to start. Recently, I discovered some new interactive features in MATLAB live scripts that make it much easier to get an overview of your data. With just a few clicks, you can display sparklines and summary statistics using table variables, sort and filter variables, and even have MATLAB generate the corresponding code for reproducibility.
The Graphics and App Building blog published an article that walks through these features showing how to explore, clean, and analyze data—all without writing any code.
If you’re interested in streamlining your exploratory data analysis or want to see what’s new in live scripts, you might find it helpful:
If you’ve tried these features or have your own tips for quick data exploration in MATLAB, I’d love to hear your thoughts!
I set my 3D matrix up with the players in the 3rd dimension. I set up the matrix with: 1) player does not hold the card (-1), player holds the card (1), and unknown holding the card (0). I moved through the turns (-1 and 1) that are fixed first. Then cycled through the conditional turns (0) while checking the cards of each player using the hints provided until it was solved. The key for me in solving several of the tests (11, 17, and 19) was looking at the 1's and 0's being held by each player.
sum(cardState==1,3);%any zeros in this 2D matrix indicate possible cards in the solution
sum(cardState==0,3)>0;%the ones in this 2D matrix indicate the only unknown positions
sum(cardState==1,3)|sum(cardState==0,3)>0;%oring the two together could provide valuable information
Some MATLAB Cody problems prohibit loops (for, while) or conditionals (if, switch, while), forcing creative solutions.
One elegant trick is to use nested functions and recursion to achieve the same logic — while staying within the rules.
Example: Recursive Summation Without Loops or Conditionals
Suppose loops and conditionals are banned, but you need to compute the sum of numbers from 1 to n. This is a simple example and obvisously n*(n+1)/2 would be preferred.
function s = sumRecursive(n)
zero=@(x)0;
s = helper(n); % call nested recursive function
function out = helper(k)
L={zero,@helper};
out = k+L{(k>0)+1}(k-1);
end
end
sumRecursive(10)
- The helper function calls itself until the base case is reached.
- Logical indexing into a cell array (k>0) act as an 'if' replacement.
- MATLAB allows nested functions to share variables and functions (zero), so you can keep state across calls.
Tips:
- Replace 'if' with logical indexing into a cell array.
- Replace for/while with recursion.
- Nested functions are local and can access outer variables, avoiding global state.
Many MATLAB Cody problems involve recognizing integer sequences.
If a sequence looks familiar but you can’t quite place it, the On-Line Encyclopedia of Integer Sequences (OEIS) can be your best friend.
OEIS will often identify the sequence, provide a formula, recurrence relation, or even direct MATLAB-compatible pseudocode.
Example: Recognizing a Cody Sequence
Suppose you encounter this sequence in a Cody problem:
1, 1, 2, 3, 5, 8, 13, 21, ...
Entering it on OEIS yields A000045 – The Fibonacci Numbers, defined by:
F(n) = F(n-1) + F(n-2), with F(1)=1, F(2)=1
You can then directly implement it in MATLAB:
function F = fibSeq(n)
F = zeros(1,n);
F(1:2) = 1;
for k = 3:n
F(k) = F(k-1) + F(k-2);
end
end
fibSeq(15)
When solving MATLAB Cody problems involving very large integers (e.g., factorials, Fibonacci numbers, or modular arithmetic), you might exceed MATLAB’s built-in numeric limits.
To overcome this, you can use Java’s java.math.BigInteger directly within MATLAB — it’s fast, exact, and often accepted by Cody if you convert the final result to a numeric or string form.
Below is an example of using it to find large factorials.
function s = bigFactorial(n)
import java.math.BigInteger
f = BigInteger('1');
for k = 2:n
f = f.multiply(BigInteger(num2str(k)));
end
s = char(f.toString); % Return as string to avoid overflow
end
bigFactorial(100)
From my experience, MATLAB's Deep Learning Toolbox is quite user-friendly, but it still falls short of libraries like PyTorch in many respects. Most users tend to choose PyTorch because of its flexibility, efficiency, and rich support for many mathematical operators. In recent years, the number of dlarray-compatible mathematical functions added to the toolbox has been very limited, which makes it difficult to experiment with many custom networks. For example, svd is currently not supported for dlarray inputs.
This link (List of Functions with dlarray Support - MATLAB & Simulink) lists all functions that support dlarray as of R2026a — only around 200 functions (including toolbox-specific ones). I would like to see support for many more fundamental mathematical functions so that users have greater freedom when building and researching custom models. For context, the core MATLAB mathematics module contains roughly 600 functions, and many application domains build on that foundation.
I hope MathWorks will prioritize and accelerate expanding dlarray support for basic math functions. Doing so would significantly increase the Deep Learning Toolbox's utility and appeal for researchers and practitioners.
Thank you.
Hey Relentless Coders! 😎
Let’s get to know each other. Drop a quick intro below and meet your teammates! This is your chance to meet teammates, find coding buddies, and build connections that make the contest more fun and rewarding!
You can share:
- Your name or nickname
- Where you’re from
- Your favorite coding topic or language
- What you’re most excited about in the contest
Let’s make Team Relentless Coders an awesome community—jump in and say hi! 🚀
Welcome to the Cody Contest 2025 and the Relentless Coders team channel! 🎉
You never give up. When a problem gets tough, you dig in deeper. This is your space to connect with like-minded coders, share insights, and help your team win. To make sure everyone has a great experience, please keep these tips in mind:
- Follow the Community Guidelines: Take a moment to review our community standards. Posts that don’t follow these guidelines may be flagged by moderators or community members.
- Ask Questions About Cody Problems: When asking for help, show your work! Include your code, error messages, and any details needed to reproduce your results. This helps others provide useful, targeted answers.
- Share Tips & Tricks: Knowledge sharing is key to success. When posting tips or solutions, explain how and why your approach works so others can learn your problem-solving methods.
- Provide Feedback: We value your feedback! Use this channel to report issues or share creative ideas to make the contest even better.
Have fun and enjoy the challenge! We hope you’ll learn new MATLAB skills, make great connections, and win amazing prizes! 🚀
I just learned you can access MATLAB Online from the following shortcut in your web browser: https://matlab.new
Thanks @Yann Debray
From his recent blog post: pip & uv in MATLAB Online » Artificial Intelligence - MATLAB & Simulink
I'm working on training neural networks without backpropagation / automatic differentiation, using locally derived analytic forms of update rules. Given that this allows a direct formula to be derived for the update rule, it removes alot of the overhead that is otherwise required from automatic differentiation.
However, matlab's functionalities for neural networks are currently solely based around backpropagation and automatic differentiation, such as the dlgradient function and requiring everything to be dlarrays during training.
I have two main requests, specifically for functions that perform a single operation within a single layer of a neural network, such as "dlconv", "fullyconnect", "maxpool", "avgpool", "relu", etc:
- these functions should also allow normal gpuArray data instead of requiring everything to be dlarrays.
- these functions are currently designed to only perform the forward pass. I request that these also be designed to perform the backward pass if user requests. There can be another input user flag that can be "forward" (default) or "backward", and then the function should have all the necessary inputs to perform that operation (e.g. for "avgpool" forward pass it only needs the avgpool input data and the avgpool parameters, but for the "avgpool" backward pass it needs the deriviative w.r.t. the avgpool output data, the avgpool parameters, and the original data dimensions). I know that there is a maxunpool function that achieves this for maxpool, but it has significant issues when trying to use it this way instead of by backpropagation in a dlgradient type layer, see (https://www.mathworks.com/matlabcentral/answers/2179587-making-a-custom-way-to-train-cnns-and-i-am-noticing-that-avgpool-is-significantly-faster-than-maxpo?s_tid=srchtitle).
I don't know how many people would benefit from this feature, and someone could always spend their time creating these functionalities themselves by matlab scripts, cuDNN mex, etc., but regardless it would be nice for matlab to have this allowable for more customizable neural net training.
Edit 15 Oct 2025: Removed incorrect code. Replaced symmatrix2sym and symfunmatrix2symfun with sym and symfun respectively (latter supported as of 2024b).
The Symbolic Math Toolbox does not have its own dot and and cross functions. That's o.k. (maybe) for garden variety vectors of sym objects where those operations get shipped off to the base Matlab functions
x = sym('x',[3,1]); y = sym('y',[3,1]);
which dot(x,y)
dot(x,y)
which cross(x,y)
cross(x,y)
clearvars
x = symmatrix('x',[3,1]); y = symmatrix('y',[3,1]);
z = symmatrix('z',[1,1]);
sympref('AbbreviateOutput',false);
dot() expands the result, which isn't really desirable for exposition.
eqn = z == dot(x,y)
Also, dot() returns the the result in terms of the conjugate of x, which can't be simplifed away at the symmatrix level
assumeAlso(sym(x),'real')
class(eqn)
try
eqn = z == simplify(dot(x,y))
catch ME
ME.message
end
To get rid of the conjugate, we have to resort to sym
eqn = simplify(sym(eqn))
but again we are in expanded form, which defeats the purpose of symmatrix (et. al.)
But at least we can do this to get a nice equation
eqn = z == x.'*y
dot errors with symfunmatrix inputs
clearvars
syms t real
x = symfunmatrix('x(t)',t,[3,1]); y = symfunmatrix('y(t)',t,[3,1]);
try
dot(x,y)
catch ME
ME.message
end
Cross works (accidentally IMO) with symmatrix, but expands the result, which isn't really desirable for exposition
clearvars
x = symmatrix('x',[3,1]); y = symmatrix('y',[3,1]);
z = symmatrix('z',[3,1]);
eqn = z == cross(x,y)
And it doesn't work at all if an input is a symfunmatrix
syms t
w = symfunmatrix('w(t)',t,[3,1]);
try
eqn = z == cross(x,w);
catch ME
ME.message
end
In the latter case we can expand with
eqn = z == cross(sym(x),symfun(w)) % x has to be converted
But we can't do the same with dot (as shown above, dot doesn't like symfun inputs)
try
eqn = z == dot(sym(x),symfun(w))
catch ME
ME.message
end
Looks like the only choice for dot with symfunmatrix is to write it by hand at the matrix level
x.'*w
or at the sym/symfun level
sym(x).'*symfun(w) % assuming x is real
Ideally, I'd like to see dot and cross implemented for symmatrix and symfunmatrix types where neither function would evaluate, i.e., expand, until both arguments are subs-ed with sym or symfun objects of appropriate dimension.
Also, it would be nice if symmatrix could be assumed to be real. Is there a reason why being able to do so wouldn't make sense?
try
assume(x,'real')
catch ME
ME.message
end
What if you had no isprime utility to rely on in MATLAB? How would you identify a number as prime? An easy answer might be something tricky, like that in simpleIsPrime0.
simpleIsPrime0 = @(N) ismember(N,primes(N));
But I’ll also disallow the use of primes here, as it does not really test to see if a number is prime. As well, it would seem horribly inefficient, generating a possibly huge list of primes, merely to learn something about the last member of the list.
Looking for a more serious test for primality, I’ve already shown how to lighten the load by a bit using roughness, to sometimes identify numbers as composite and therefore not prime.
https://www.mathworks.com/matlabcentral/discussions/tips/879745-primes-and-rough-numbers-basic-ideas
But to actually learn if some number is prime, we must do a little more. Yes, this is a common homework problem assigned to students, something we have seen many times on Answers. It can be approached in many ways too, so it is worth looking at the problem in some depth.
The definition of a prime number is a natural number greater than 1, which has only two factors, thus 1 and itself. That makes a simple test for primality of the number N easy. We just try dividing the number by every integer greater than 1, and not exceeding N-1. If any of those trial divides leaves a zero remainder, then N cannot be prime. And of course we can use mod or rem instead of an explicit divide, so we need not worry about floating point trash, as long as the numbers being tested are not too large.
simpleIsPrime1 = @(N) all(mod(N,2:N-1) ~= 0);
Of course, simpleIsPrime1 is not a good code, in the sense that it fails to check if N is an integer, or if N is less than or equal to 1. It is not vectorized, and it has no documentation at all. But it does the job well enough for one simple line of code. There is some virtue in simplicity after all, and it is certainly easy to read. But sometimes, I wish a function handle could include some help comments too! A feature request might be in the offing.
simpleIsPrime1(9931)
simpleIsPrime1(9932)
simpleIsPrime1 works quite nicely, and seems pretty fast. What could be wrong? At some point, the student is given a more difficult problem, to identify if a significantly larger integer is prime. simpleIsPrime1 will then cause a computer to grind to a distressing halt if given a sufficiently large number to test. Or it might even error out, when too large a vector of numbers was generated to test against. For example, I don't think you want to test a number of the order of 2^64 using simpleIsPrime1, as performing on the order of 2^64 divides will be highly time consuming.
uint64(2)^63-25
Is it prime? I’ve not tested it to learn if it is, and simpleIsPrime1 is not the tool to perform that test anyway.
A student might realize the largest possible integer factors of some number N are the numbers N/2 and N itself. But, if N/2 is a factor, then so is 2, and some thought would suggest it is sufficient to test only for factors that do not exceed sqrt(N). This is because if a is a divisor of N, then so is b=N/a. If one of them is larger than sqrt(N), then the other must be smaller. That could lead us to an improved scheme in simpleIsPrime2.
simpleIsPrime2 = @(N) all(mod(N,2:sqrt(N)));
For an integer of the size 2^64, now you only need to perform roughly 2^32 trial divides. Maybe we might consider the subtle improvement found in simpleIsPrime3, which avoids trial divides by the even integers greater than 2.
simpleIsPrime3 = @(N) (N == 2) || (mod(N,2) && all(mod(N,3:2:sqrt(N))));
simpleIsPrime3 needs only an approximate maximum of 2^31 trial divides even for numbers as large as uint64 can represent. While that is large, it is still generally doable on the computers we have today, even if it might be slow.
Sadly, my goals are higher than even the rather lofty limit given by UINT64 numbers. The problem of course is that a trial divide scheme, despite being 100% accurate in its assessment of primality, is a time hog. Even an O(sqrt(N)) scheme is far too slow for numbers with thousands or millions of digits. And even for a number as “small” as 1e100, a direct set of trial divides by all primes less than sqrt(1e100) would still be practically impossible, as there are roughly n/log(n) primes that do not exceed n. For an integer on the order of 1e50,
1e50/log(1e50)
It is practically impossible to perform that many divides on any computer we can make today. Can we do better? Is there some more efficient test for primality? For example, we could write a simple sieve of Eratosthenes to check each prime found not exceeding sqrt(N).
function [TF,SmallPrime] = simpleIsPrime4(N)
% simpleIsPrime3 - Sieve of Eratosthenes to identify if N is prime
% [TF,SmallPrime] = simpleIsPrime3(N)
%
% Returns true if N is prime, as well as the smallest prime factor
% of N when N is composite. If N is prime, then SmallPrime will be N.
Nroot = ceil(sqrt(N)); % ceil caters for floating point issues with the sqrt
TF = true;
SieveList = true(1,Nroot+1); SieveList(1) = false;
SmallPrime = 2;
while TF
% Find the "next" true element in SieveList
while (SmallPrime <= Nroot+1) && ~SieveList(SmallPrime)
SmallPrime = SmallPrime + 1;
end
% When we drop out of this loop, we have found the next
% small prime to check to see if it divides N, OR, we
% have gone past sqrt(N)
if SmallPrime > Nroot
% this is the case where we have now looked at all
% primes not exceeding sqrt(N), and have found none
% that divide N. This is where we will drop out to
% identify N as prime. TF is already true, so we need
% not set TF.
SmallPrime = N;
return
else
if mod(N,SmallPrime) == 0
% smallPrime does divide N, so we are done
TF = false;
return
end
% update SieveList
SieveList(SmallPrime:SmallPrime:Nroot) = false;
end
end
end
simpleIsPrime4 does indeed work reasonably well, though it is sometimes a little slower than is simpleIsPrime3, and everything is hugely faster than simpleIsPrime1.
timeit(@() simpleIsPrime1(111111111))
timeit(@() simpleIsPrime2(111111111))
timeit(@() simpleIsPrime3(111111111))
timeit(@() simpleIsPrime4(111111111))
All of those times will slow to a crawl for much larger numbers of course. And while I might find a way to subtly improve upon these codes, any improvement will be marginal in the end if I try to use any such direct approach to primality. We must look in a different direction completely to find serious gains.
At this point, I want to distinguish between two distinct classes of tests for primality of some large number. One class of test is what I might call an absolute or infallible test, one that is perfectly reliable. These are tests where if X is identified as prime/composite then we can trust the result absolutely. The tests I showed in the form of simpleIsPrime1, simpleIsPrime2, simpleIsPrime3 and aimpleIsprime4, were all 100% accurate, thus they fall into the class of infallible tests.
The second general class of test for primality is what I will call an evidentiary test. Such a test provides evidence, possibly quite strong evidence, that the given number is prime, but in some cases, it might be mistaken. I've already offered a basic example of a weak evidentiary test for primality in the form of roughness. All primes are maximally rough. And therefore, if you can identify X as being rough to some extent, this provides evidence that X is also prime, and the depth of the roughness test influences the strength of the evidence for primality. While this is generally a fairly weak test, it is a test nevertheless, and a good exclusionary test, a good way to avoid more sophisticated but time consuming tests.
These evidentiary tests all have the property that if they do identify X as being composite, then they are always correct. In the context of roughness, if X is not sufficiently rough, then X is also not prime. On the other side of the coin, if you can show X is at least (sqrt(X)+1)-rough, then it is positively prime. (I say this to suggest that some evidentiary tests for primality can be turned into truth telling tests, but that may take more effort than you can afford.) The problem is of course that is literally impossible to verify that degree of roughness for numbers with many thousands of digits.
In my next post, I'll look at the Fermat test for primality, based on Fermat's little theorem.
Something that I periodically wonder about is whether an integration with the Rubi integration rules package would improve symbolic integration in Matlab's Symbolic Toolbox. The project is open-source with an MIT-licensed, has a Mathematica implementation, and supposedly SymPy is working on an implementation. Much of my intrigue comes from this 2022 report that compared the previous version of Rubi (4.16.1) against various CAS systems, including Matlab 2021a (Mupad):
While not really an official metric for Rubi, this does "feel" similar to my experience computing symbolic integrals in Matlab Symbolic Toolbox vs Maple/Mathematica. What do y'all think?
I saw an interesting problem on a reddit math forum today. The question was to find a number (x) as close as possible to r=3.6, but the requirement is that both x and 1/x be representable in a finite number of decimal places.
The problem of course is that 3.6 = 18/5. And the problem with 18/5 has an inverse 5/18, which will not have a finite representation in decimal form.
In order for a number and its inverse to both be representable in a finite number of decimal places (using base 10) we must have it be of the form 2^p*5^q, where p and q are integer, but may be either positive or negative. If that is not clear to you intuitively, suppose we have a form
2^p*5^-q
where p and q are both positive. All you need do is multiply that number by 10^q. All this does is shift the decimal point since you are just myltiplying by powers of 10. But now the result is
2^(p+q)
and that is clearly an integer, so the original number could be represented using a finite number of digits as a decimal. The same general idea would apply if p was negative, or if both of them were negative exponents.
Now, to return to the problem at hand... We can obviously adjust the number r to be 20/5 = 4, or 16/5 = 3.2. In both cases, since the fraction is now of the desired form, we are happy. But neither of them is really close to 3.6. My goal will be to find a better approximation, but hopefully, I can avoid a horrendous amount of trial and error. It would seem the trick might be to take logs, to get us closer to a solution. That is, suppose I take logs, to the base 2?
log2(3.6)
I used log2 here because that makes the problem a little simpler, since log2(2^p)=p. Therefore we want to find a pair of integers (p,q) such that
log2(3.6) + delta = p + log2(5)*q
where delta is as close to zero as possible. Thus delta is the error in our approximation to 3.6. And since we are working in logs, delta can be viewed as a proportional error term. Again, p and q may be any integers, either positive or negative. The two cases we have seen already have (p,q) = (2,0), and (4,-1).
Do you see the general idea? The line we have is of the form
log2(3.6) = p + log2(5)*q
it represents a line in the (p,q) plane, and we want to find a point on the integer lattice (p,q) where the line passes as closely as possible.
[Xl,Yl] = meshgrid([-10:10]);
plot(Xl,Yl,'k.')
hold on
fimplicit(@(p,q) -log2(3.6) + p + log2(5)*q,[-10,10,-10,10],'g-')
plot([2 4],[0,-1],'ro')
hold off
Now, some might think in terms of orthogonal distance to the line, but really, we want the vertical distance to be minimized. Again, minimize abs(delta) in the equation:
log2(3.6) + delta = p + log2(5)*q
where p and q are integer.
Can we do that using MATLAB? The skill about about mathematics often lies in formulating a word problem, and then turning the word problem into a problem of mathematics that we know how to solve. We are almost there now. I next want to formulate this into a problem that intlinprog can solve. The problem at first is intlinprog cannot handle absolute value constraints. And the trick there is to employ slack variables, a terribly useful tool to emply on this class of problem.
Rewrite delta as:
delta = Dpos - Dneg
where Dpos and Dneg are real variables, but both are constrained to be positive.
prob = optimproblem;
p = optimvar('p',lower = -50,upper = 50,type = 'integer');
q = optimvar('q',lower = -50,upper = 50,type = 'integer');
Dpos = optimvar('Dpos',lower = 0);
Dneg = optimvar('Dneg',lower = 0);
Our goal for the ILP solver will be to minimize Dpos + Dneg now. But since they must both be positive, it solves the min absolute value objective. One of them will always be zero.
r = 3.6;
prob.Constraints = log2(r) + Dpos - Dneg == p + log2(5)*q;
prob.Objective = Dpos + Dneg;
The solve is now a simple one. I'll tell it to use intlinprog, even though it would probably figure that out by itself. (Note: if I do not tell solve which solver to use, it does use intlinprog. But it also finds the correct solution when I told it to use GA offline.)
solve(prob,solver = 'intlinprog')
The solution it finds within the bounds of +/- 50 for both p and q seems pretty good. Note that Dpos and Dneg are pretty close to zero.
2^39*5^-16
and while 3.6028979... seems like nothing special, in fact, it is of the form we want.
R = sym(2)^39*sym(5)^-16
vpa(R,100)
vpa(1/R,100)
both of those numbers are exact. If I wanted to find a better approximation to 3.6, all I need do is extend the bounds on p and q. And we can use the same solution approch for any floating point number.
