How can I make use of Saramäki and MATLAB implementation of the Remez multiple exchange algorithm for FIR filter design

Question

Seeking assistance with an optimized and well-documented implementation of the Parks-McClellan algorithm, I explored beyond the original FORTRAN code provided in A Computer Program for Designing Optimum FIR Linear Phase Digital Filter.

Instead, I discovered other available options such as an open-source C implementation, which was essentially a direct translation of the FORTRAN version. Further research led me to uncover texts like Two Novel Implementations of the Remez Multiple Exchange Algorithm for Optimum FIR Filter Design and A MATLAB based optimum multiband FIR filters design program following the original idea of the Remez multiple exchange algorithm These programs assumes to be an enhanced version of original code, resulting in significantly faster execution times.

I'm encountering two main issues:

1. I've struggled to utilize the algorithm effectively to generate an FIR filter of an order higher than approximately 700.
2. I'm uncertain whether I'm correctly implementing the code.

Note: I understand that the Remez implementation doesn't directly provide the impulse response. Moreover, I am aware of the necessity to adjust the desired magnitude response and weight function based on the filter type (symmetry and length). Any guidance on this aspect would be greatly appreciated.

Currently, my code usage resembles the following:

% the number of indices in the alternation set (keep even?)
nz = 1000;
% resolution of the frequency grid
ngrid = nz * 20;
% frequency grid (length ngrid)
grid = (0:ngrid-1).*(0.5/ngrid);
% indices of extremal frequencies (length nz)
iext = round(linspace(1, ngrid, nz));
% desired amplitude response
nyq = 0.5;
fpass = 0.45*nyq; fstop = 0.55*nyq; fcut = (fpass + fstop)/2;
des = (grid <= fcut);
% weight function
Kp = 1.5; Ks = 3.5;
wt = (grid<=fpass).*Kp +(grid>=fstop).*Ks;
wt = wt + 0.3;
[x,y,ad,dev]=remez_imp1(nz, iext, ngrid, grid, des, wt);

Answer 1

700 coefficients seems a high number for remez.

If I try using Matlab's design tools, it seems get numerical problems between 300 and 350 coefficients.