arm_mult_q31: saturate to 0x7FFFFFFF, as the other two mult kernels do

[DanielMBouyou]

There are two separate things in this PR and I would rather split them myself than
have them read as one. Part A is the defect: a saturation the documentation promises
and that does not happen. Part B is the full effect of the fix, which is larger than
the defect and is a trade-off that belongs to you, not to me.

---

Part A: the defect

A1. The test is already there, and the reference already holds the right value

Testing/Source/Tests/BasicTestsQ31.cpp, TEST_MULT_Q31_29 ("Test mul sat") loads
MAXNEG2_Q31_ID, which is 0x80000000, into both inputs, and compares against
REF_POSSAT_16_Q31_ID. That reference comes from Testing/PatternGeneration/BasicMaths.py:

d3  = 1.0*(data2) / 2**format      # data2 = -2**31   ->  d3 = -1.0
ref = d3*d3                        # 1.0
config.writeReference(nb+5, ref, "PosSat")

and Tools.to_q31 clamps on the way out:

r = int(round(v * 2**31))          # 1.0 -> 0x80000000
if (r > 0x07FFFFFFF): r = 0x07FFFFFFF

so the pattern committed to the repository holds 0x7FFFFFFF. The kernel returns
0x7FFFFFFE. The test passes because the assertion is

ASSERT_NEAR_EQ(output,ref,ABS_ERROR_Q31);   /* ABS_ERROR_Q31 = 4 */

A four-LSB tolerance on a saturation case absorbs a one-LSB miss. Nothing about the
expected value is in dispute: it is already in your pattern.

This PR tightens that one assertion to exact equality, for the saturation case only.
ABS_ERROR_Q31 is untouched and the other arm_mult_q31 tests keep it.

A2. All three kernels of the family carry the same sentence

Verbatim, in arm_mult_q31.c, arm_mult_q15.c and arm_mult_q7.c:

The function uses saturating arithmetic.
Results outside of the allowable Q31 range[0x80000000 0x7FFFFFFF] are saturated.

A3. The two siblings honour it

arm_mult_q15 and arm_mult_q7 shift by the full fractional width, then saturate:

/* arm_mult_q15.c */  *pDst++ = (q15_t) __SSAT((((q31_t) a * b) >> 15), 16);
/* arm_mult_q7.c  */  *pDst++ = (q7_t)  __SSAT((((q15_t) a * b) >>  7),  8);

Checked against an integer model derived from the documented Q format, compared with ==:

kernel	input pairs compared	mismatches
arm_mult_q7	65 536 (the complete 256 x 256 domain)	0
arm_mult_q15	27 921 (every container edge x edge, plus random)	0

A4. One input pair in the whole domain is concerned

Q1.31 x Q1.31 leaves the binary point at bit 62. Getting back to Q1.31 is a right shift
of 31. The committed code shifts by 32, which puts the binary point at bit 30, saturates
to 31 bits, which is 2^30 - 1, and then doubles:

out = ((q63_t) *pSrcA++ * *pSrcB++) >> 32;
out = __SSAT(out, 31);
*pDst++ = out << 1U;

so the reachable ceiling is 2^31 - 2. The only product that leaves the Q31 range is
-1.0 x -1.0 = +1.0, that is 0x80000000 x 0x80000000; every other pair fits. Measured:

arm_mult_q31(0x80000000, 0x80000000)  ->  0x7FFFFFFE      documented: 0x7FFFFFFF

I built twelve out-of-range cases across eleven q31 kernels that document saturating
arithmetic, and arm_scale_q31 carries two of them. Before this patch, nine of the
twelve returned 0x7FFFFFFF, one did not, this one, and two could not be driven out of
range at all, so they are recorded as untested rather than as passing: 9 + 1 + 2 = 12.
After this patch: ten saturate, none fail, the same two remain untested.

The replacement is the idiom already used elsewhere in the tree, including in the
non-ARM_MATH_DSP branch of arm_cmplx_mult_real_q31.c:

*pDst++ = (q31_t) clip_q63_to_q31(((q63_t) *pSrcA++ * *pSrcB++) >> 31);

---

Part B: what the fix actually changes, which is much more than that one pair

Part A is about one input pair. The patch is not. I am stating this myself because
the diff would otherwise understate it.

out = (p >> 32); out = __SSAT(out, 31); out << 1 produces 2 * floor(p / 2^32).
That value is always even. The correctly truncated Q31 result is floor(p / 2^31),
which is odd for half the domain. So the committed kernel loses the least significant
bit by construction, and the patch gives it back.

Measured by running the committed kernel and the patched kernel side by side in one
binary, over the same inputs. Both numbers come from compiled CMSIS-DSP code, and nothing
is modelled:

uniform random pairs                 : 1000000
  outputs that change                : 500079  (50.01 %)
  difference range (new - old)       : [0, 1]
  odd outputs, committed kernel      : 0
  odd outputs, patched kernel        : 500079  (50.01 %)

So: about half of all outputs change, each by exactly one LSB, never more, never in the
other direction. On an edge-weighted input set (every container edge against every other,
plus random) the figure was 26 251 changed out of 54 225, and after the patch that set
matches the exact integer model at every one of its 54 225 points, where before it matched
at 27 974.

This is accuracy traded for speed, and the trade is very probably deliberate.
(a*b) >> 32 followed by << 1 is the shape that falls onto a single most-significant-word
multiply on Cortex-M; the replacement needs the full 64-bit product and a shift across the
register pair.

I have measured no cycle cost, on any core, and I claim none. I have no way to weigh
one LSB of accuracy in a vector multiply against instructions in a hot loop on the parts
you care about. That balance is yours.

If you want the current form kept, then the correct change is the opposite of this PR:
document the exception in the three files, the way #327 documented the missing saturation
in arm_pid_q31, and close this. I will follow either way, and I would rather you told me
than merged this because it was quiet.

---

Spotted, not fixed

arm_cmplx_mult_real_q31.c carries both idioms inside the same function, selected by
ARM_MATH_DSP, and so returns two different numbers for the same inputs:

0x80000000 x 0x80000000,  ARM_MATH_DSP not defined  ->  0x7FFFFFFF
0x80000000 x 0x80000000,  ARM_MATH_DSP defined      ->  0x7FFFFFFE

Same documented sentence as above. Deliberately left untouched here: it is a separate
decision, and this change is meant to stay small.

---

Where this was run

An x86-64 Linux container: gcc 12.2.0, cmake -S . -B build-host -DHOST=ON -DDISABLEFLOAT16=ON, and the Testing/ framework built natively.

Nothing was run on an Arm target: no board, no FVP, no simulator of any kind.
Someone with hardware should confirm before this is merged. That applies to the numbers in
Part B as much as to Part A.

One note on running Testing/ off Arm: Testing/CMakeLists.txt hardcodes
-march=armv8.2-a+fp16fml at three sites, so as committed it only builds on an AArch64
host. I dropped that flag locally to build it on x86-64. That is a local build-side change
only; it is not part of this PR and it touches no library or test source.

Results after the change:

• Testing/, 19 suites: 502 tests passed, 1 failed, namely test_mat_qr_f32, SNR
  87.0749 < 90. That one fails identically with and without this patch; I re-ran the suite
  with arm_mult_q31.c restored to its committed content to check, rather than assume.
• Python wrapper, PythonWrapper/examples/testall.py: 309 tests, 4 failures, all
  pre-existing and all floating point: test_cholesky_f32, test_nuttall4,
  test_nuttall4a, test_nuttall4b.
• MFCCQ31 passes; arm_mfcc_q31.c:141 is the only in-tree caller of arm_mult_q31, so
  it is the one place where the Part B change propagates inside the library.
• No new compiler warning on the patched translation unit under
  -Wall -Wextra -Wsign-compare -Wdouble-promotion -Wfloat-conversion -Wmissing-prototypes.