# Find min sum for indices to match equation

I want to reduce time complexity of my code:

I have an array of integers arr say [1, 2, 3, 6, 67]

I have an equation :

a*x+b*y=z, I can use the array values in this equation to replace x and y, by some array positions say arr[i], arr[j] i and j can be different or same.

so it looks like this:

a*arr[i] + b*arr[j] = z;


Now I have another array Zarr = [1, 4, 30, 7, 88]

The task is to find values for a and b such that a*arr[i] + b*arr[j] = Zarr[k], where k is the index position of Zarr.

once a and b are calculated I need to get the sum (a + b) and store in result, the condition is the sum (a+b) should be minimum.

I have 1 more restriction, if the sum a+b is more than input variable Z_MAX, then set result for that index to be -1;

Example:

so for arr =  [1, 2, 3, 6, 67]
Zarr = [1, 4, 30, 7, 88]
Z_MAX = 5
Result is [1,2,5,2,-1]


Explanation:

Zarr[0] = 1
i=0, j=0
a=1, b=0
sum = a+b = 1
Eq = 1*1 + 0*1 = 1

Zarr[1] = 4
i=1, j=1
a=1, b=1
sum = a+b = 2
Eq = 1*2 + 1*2 = 4

Zarr[2] = 30
i=3, j=3
a=5, b=0
sum = a+b = 5
Eq = 5*6 + 0*6 = 30

Zarr[3] = 7
i=0, j=3
a=1, b=1
sum = a+b = 2
Eq = 1*1 + 1*6 = 7

Zarr[4] = 88
It's not possible to find a and b whose sum to be less than or equal to Z_MAX, so -1

Finally result = [1,2,5,2,-1]


Here is my code, using brute force approach:

public static void main(String[] args) {
int[] a = solve(new int[]{1, 2, 3, 6, 67}, new int[]{1, 4, 30, 7, 88}, 5);
System.out.println(Arrays.toString(a));//1,2,5,2,-1
}

public static int[] solve(int[] arr, int[] zArr, int Z_MAX) {
int n = arr.length;
int[] result = new int[zArr.length];
for(int i=0; i<zArr.length; i++) {
result[i] = Integer.MAX_VALUE;
}
for (int p = 0; p < zArr.length; p++) {
int z = zArr[p];

for(int i=0; i<n; i++) {
for(int j=0; j<n;j++) {

for (int a = 0; a <= Z_MAX; a++) {
for (int b = 0; b <= Z_MAX; b++) {
if (a * arr[i] + b * arr[j] == z) {
result[p] = Math.min(result[p], a + b);
break;
}
}
}
}
}
if (result[p] > Z_MAX) {
result[p] = -1;
}
}
return result;
}


Constraints:

array size is of range 1 to 1000
array element values in range 1 to 10^7
Z_MAX is 1 to 20
zArr size is 1 to 20
zArr element values in range 1 to 10^9


How can the time complexity be decreased?

• Where is k defined? Commented Feb 29 at 21:51
• @Reinderien, its a typo error, I corrected it now. Commented Feb 29 at 23:09
• The problem looks ill-defined. There could be more than one pair of a,b satisfying the constrains. Which one is the result?
– vnp
Commented Mar 1 at 5:05
• @vnp, I updated now, the sum of a+b should be minimum to get result. Commented Mar 1 at 14:04
• Looks like you're supposed to use Bézout's identity, but I haven't worked it out. Commented Mar 2 at 2:23

# reference

Cite your references, please. Or at least tell us where this problem comes from. What motivated it?

The details matter. For example, perhaps it makes sense to turn discretized integers into continuous reals, and minimize a loss function
loss = a*arr[i] + b*arr[j] - Zarr[k]

Also, the relevant Business Domain might give us additional clues, like telling us that certain inputs cannot be zero or cannot be negative.

# generic names

arr is a terrible identifier. It doesn't tell us anything helpful. (Its type already revealed it's an integer array.) Choose a better name.

zArr could be a good name, if the math in your cited reference used a z variable.

result is an OK name, but a name like numWidgets or dollars would be much better, if that's the meaning of the result figures.

# sort()

I recommend sorting those three vectors upon entry, or requiring that caller supply ordered inputs. We should be able to exploit that via binary search. Or at the very least, it would allow abandoning some loops early.

# GCD

It seems like extracting the prime factors of each input number should be an aid to your Diophantine analysis.

Once you have fixed a, arr[i], b, and z, a division and a binary search would give you arr[j]. Or it might be worth your while to build a Map from integers to j index, if you get to reuse it for a bunch of p (z) values.

# mixed-integer

There are several OR libraries available for solving mixed integer problems. You might benefit from recasting this problem in terms of {0, 1} indicator decision variables multiplied by the arr candidate values.

As J_H says, one of the ways to frame this problem is as a MIP/MILP. It might beat your (quintuple!) nested loop, but your mileage may vary based on input size and chosen solver, so you need to benchmark. There's a long list of gotchas:

• you'll either need to write a (maybe MPS) problem file generator and send off to a stand-alone solver, or pull in a Java FFI library that exposes one of the solver APIs for you.
• depending on which API you pick, you might have to reduce the list of expressions I demonstrate to matrix form.
• the problem cannot be expressed directly as ar_i + br_j = z because that's non-linear. You need to linearise via big-M constraints.
• the value of M should generally not be hard-coded, and leverage knowledge of the maximum value of arr.

## Python equivalent

In PuLP/Python this looks like

import pulp

arr = (1, 2, 3, 6, 67)
Zarr = (1, 4, 30, 7, 88)
Z_MAX = 5

'''
a.arr[i] + b.arr[j] = z[k] for all k
what is a + b? a + b <= ZMAX
a >= 0, a int
b >= 0, b int

aarri + barrj = z
aarri >= a*1 - M(1 - seli)
'''

def solve(z: int) -> str | None:
isel = pulp.LpVariable.matrix(name='i', indices=range(len(arr)), cat=pulp.LpBinary)
jsel = pulp.LpVariable.matrix(name='j', indices=range(len(arr)), cat=pulp.LpBinary)
a = pulp.LpVariable(name='a', cat=pulp.LpInteger, lowBound=0, upBound=Z_MAX)
b = pulp.LpVariable(name='b', cat=pulp.LpInteger, lowBound=0, upBound=Z_MAX)
ari = pulp.LpVariable(name='ari', cat=pulp.LpContinuous)
brj = pulp.LpVariable(name='brj', cat=pulp.LpContinuous)

prob = pulp.LpProblem(name='diophantine', sense=pulp.LpMinimize)
prob.setObjective(a + b)
prob.addConstraint(name='z', constraint=ari + brj == z)

M = 2*Z_MAX*max(arr)
for i, (arr_value, sel) in enumerate(zip(arr, isel)):
av = a*arr_value
relax = M*(1 - sel)
prob.addConstraint(name=f'lower_ari_{i}', constraint=ari >= av - relax)
prob.addConstraint(name=f'upper_ari_{i}', constraint=ari <= av + relax)
for j, (arr_value, sel) in enumerate(zip(arr, jsel)):
bv = b*arr_value
relax = M*(1 - sel)
prob.addConstraint(name=f'lower_brj_{j}', constraint=brj >= bv - relax)
prob.addConstraint(name=f'upper_brj_{j}', constraint=brj <= bv + relax)

# print(prob)
prob.solve()
if prob.status != pulp.LpStatusOptimal:
return None

ri = next(v for sel, v in zip(isel, arr) if sel.value() > 0.5)
rj = next(v for sel, v in zip(jsel, arr) if sel.value() > 0.5)
return f'{a.value()}*{ri} + {b.value()}*{rj} = {z}'

for z in Zarr:
print(solve(z))


## LP construction

The above has the following problem description output, which shows the affine expressions necessary to constrain it:

diophantine:
MINIMIZE
1*a + 1*b + 0
SUBJECT TO
z: ari + brj = 1

iexcl: i_0 + i_1 + i_2 + i_3 + i_4 = 1
jexcl: j_0 + j_1 + j_2 + j_3 + j_4 = 1
lower_ari_0: - a + ari - 670 i_0 >= -670
upper_ari_0: - a + ari + 670 i_0 <= 670
lower_ari_1: - 2 a + ari - 670 i_1 >= -670
upper_ari_1: - 2 a + ari + 670 i_1 <= 670
lower_ari_2: - 3 a + ari - 670 i_2 >= -670
upper_ari_2: - 3 a + ari + 670 i_2 <= 670
lower_ari_3: - 6 a + ari - 670 i_3 >= -670
upper_ari_3: - 6 a + ari + 670 i_3 <= 670
lower_ari_4: - 67 a + ari - 670 i_4 >= -670
upper_ari_4: - 67 a + ari + 670 i_4 <= 670
lower_brj_0: - b + brj - 670 j_0 >= -670
upper_brj_0: - b + brj + 670 j_0 <= 670
lower_brj_1: - 2 b + brj - 670 j_1 >= -670
upper_brj_1: - 2 b + brj + 670 j_1 <= 670
lower_brj_2: - 3 b + brj - 670 j_2 >= -670
upper_brj_2: - 3 b + brj + 670 j_2 <= 670
lower_brj_3: - 6 b + brj - 670 j_3 >= -670
upper_brj_3: - 6 b + brj + 670 j_3 <= 670
lower_brj_4: - 67 b + brj - 670 j_4 >= -670
upper_brj_4: - 67 b + brj + 670 j_4 <= 670
VARIABLES
__dummy = 0 Continuous
0 <= a <= 5 Integer
ari free Continuous
0 <= b <= 5 Integer
brj free Continuous
0 <= i_0 <= 1 Integer
0 <= i_1 <= 1 Integer
0 <= i_2 <= 1 Integer
0 <= i_3 <= 1 Integer
0 <= i_4 <= 1 Integer
0 <= j_0 <= 1 Integer
0 <= j_1 <= 1 Integer
0 <= j_2 <= 1 Integer
0 <= j_3 <= 1 Integer
0 <= j_4 <= 1 Integer


This produces solutions

0.0*67 + 1.0*1 = 1
2.0*2 + 0.0*2 = 4
5.0*6 + 0.0*67 = 30
1.0*1 + 1.0*6 = 7
None


## MPS problem format

If you can generate MPS (from Java or otherwise) looking like the following, then a few open-source solvers including CBC will accept it. This output is generated by PuLP.

*SENSE:Minimize
NAME          MODEL
ROWS
N  OBJ
E  C0000000
E  C0000001
E  C0000002
G  C0000003
L  C0000004
G  C0000005
L  C0000006
G  C0000007
L  C0000008
G  C0000009
L  C0000010
G  C0000011
L  C0000012
G  C0000013
L  C0000014
G  C0000015
L  C0000016
G  C0000017
L  C0000018
G  C0000019
L  C0000020
G  C0000021
L  C0000022
COLUMNS
X0000000  OBJ        1.000000000000e+00
MARK      'MARKER'                 'INTORG'
X0000001  C0000003  -1.000000000000e+00
X0000001  C0000004  -1.000000000000e+00
X0000001  C0000005  -2.000000000000e+00
X0000001  C0000006  -2.000000000000e+00
X0000001  C0000007  -3.000000000000e+00
X0000001  C0000008  -3.000000000000e+00
X0000001  C0000009  -6.000000000000e+00
X0000001  C0000010  -6.000000000000e+00
X0000001  C0000011  -6.700000000000e+01
X0000001  C0000012  -6.700000000000e+01
MARK      'MARKER'                 'INTEND'
X0000002  C0000000   1.000000000000e+00
X0000002  C0000003   1.000000000000e+00
X0000002  C0000004   1.000000000000e+00
X0000002  C0000005   1.000000000000e+00
X0000002  C0000006   1.000000000000e+00
X0000002  C0000007   1.000000000000e+00
X0000002  C0000008   1.000000000000e+00
X0000002  C0000009   1.000000000000e+00
X0000002  C0000010   1.000000000000e+00
X0000002  C0000011   1.000000000000e+00
X0000002  C0000012   1.000000000000e+00
MARK      'MARKER'                 'INTORG'
X0000003  C0000013  -1.000000000000e+00
X0000003  C0000014  -1.000000000000e+00
X0000003  C0000015  -2.000000000000e+00
X0000003  C0000016  -2.000000000000e+00
X0000003  C0000017  -3.000000000000e+00
X0000003  C0000018  -3.000000000000e+00
X0000003  C0000019  -6.000000000000e+00
X0000003  C0000020  -6.000000000000e+00
X0000003  C0000021  -6.700000000000e+01
X0000003  C0000022  -6.700000000000e+01
MARK      'MARKER'                 'INTEND'
X0000004  C0000000   1.000000000000e+00
X0000004  C0000013   1.000000000000e+00
X0000004  C0000014   1.000000000000e+00
X0000004  C0000015   1.000000000000e+00
X0000004  C0000016   1.000000000000e+00
X0000004  C0000017   1.000000000000e+00
X0000004  C0000018   1.000000000000e+00
X0000004  C0000019   1.000000000000e+00
X0000004  C0000020   1.000000000000e+00
X0000004  C0000021   1.000000000000e+00
X0000004  C0000022   1.000000000000e+00
MARK      'MARKER'                 'INTORG'
X0000005  C0000001   1.000000000000e+00
X0000005  C0000003  -6.700000000000e+02
X0000005  C0000004   6.700000000000e+02
MARK      'MARKER'                 'INTEND'
MARK      'MARKER'                 'INTORG'
X0000006  C0000001   1.000000000000e+00
X0000006  C0000005  -6.700000000000e+02
X0000006  C0000006   6.700000000000e+02
MARK      'MARKER'                 'INTEND'
MARK      'MARKER'                 'INTORG'
X0000007  C0000001   1.000000000000e+00
X0000007  C0000007  -6.700000000000e+02
X0000007  C0000008   6.700000000000e+02
MARK      'MARKER'                 'INTEND'
MARK      'MARKER'                 'INTORG'
X0000008  C0000001   1.000000000000e+00
X0000008  C0000009  -6.700000000000e+02
X0000008  C0000010   6.700000000000e+02
MARK      'MARKER'                 'INTEND'
MARK      'MARKER'                 'INTORG'
X0000009  C0000001   1.000000000000e+00
X0000009  C0000011  -6.700000000000e+02
X0000009  C0000012   6.700000000000e+02
MARK      'MARKER'                 'INTEND'
MARK      'MARKER'                 'INTORG'
X0000010  C0000002   1.000000000000e+00
X0000010  C0000013  -6.700000000000e+02
X0000010  C0000014   6.700000000000e+02
MARK      'MARKER'                 'INTEND'
MARK      'MARKER'                 'INTORG'
X0000011  C0000002   1.000000000000e+00
X0000011  C0000015  -6.700000000000e+02
X0000011  C0000016   6.700000000000e+02
MARK      'MARKER'                 'INTEND'
MARK      'MARKER'                 'INTORG'
X0000012  C0000002   1.000000000000e+00
X0000012  C0000017  -6.700000000000e+02
X0000012  C0000018   6.700000000000e+02
MARK      'MARKER'                 'INTEND'
MARK      'MARKER'                 'INTORG'
X0000013  C0000002   1.000000000000e+00
X0000013  C0000019  -6.700000000000e+02
X0000013  C0000020   6.700000000000e+02
MARK      'MARKER'                 'INTEND'
MARK      'MARKER'                 'INTORG'
X0000014  C0000002   1.000000000000e+00
X0000014  C0000021  -6.700000000000e+02
X0000014  C0000022   6.700000000000e+02
MARK      'MARKER'                 'INTEND'
RHS
RHS       C0000000   1.000000000000e+00
RHS       C0000001   1.000000000000e+00
RHS       C0000002   1.000000000000e+00
RHS       C0000003  -6.700000000000e+02
RHS       C0000004   6.700000000000e+02
RHS       C0000005  -6.700000000000e+02
RHS       C0000006   6.700000000000e+02
RHS       C0000007  -6.700000000000e+02
RHS       C0000008   6.700000000000e+02
RHS       C0000009  -6.700000000000e+02
RHS       C0000010   6.700000000000e+02
RHS       C0000011  -6.700000000000e+02
RHS       C0000012   6.700000000000e+02
RHS       C0000013  -6.700000000000e+02
RHS       C0000014   6.700000000000e+02
RHS       C0000015  -6.700000000000e+02
RHS       C0000016   6.700000000000e+02
RHS       C0000017  -6.700000000000e+02
RHS       C0000018   6.700000000000e+02
RHS       C0000019  -6.700000000000e+02
RHS       C0000020   6.700000000000e+02
RHS       C0000021  -6.700000000000e+02
RHS       C0000022   6.700000000000e+02
BOUNDS
FX BND       X0000000   0.000000000000e+00
UP BND       X0000001   5.000000000000e+00
FR BND       X0000002
UP BND       X0000003   5.000000000000e+00
FR BND       X0000004
BV BND       X0000005
BV BND       X0000006
BV BND       X0000007
BV BND       X0000008
BV BND       X0000009
BV BND       X0000010
BV BND       X0000011
BV BND       X0000012
BV BND       X0000013
BV BND       X0000014
ENDATA


## Matrix form

From Scipy or a Java equivalent, this is the sparse matrix form that you can use to define the problem constraints. It is trickier but more efficient to set up, particularly if you share common portions between your outer z loop iteration.

from functools import partial
from typing import Callable

import numpy as np
from scipy.optimize import milp, Bounds, LinearConstraint, OptimizeResult
import scipy.sparse as sp

def setup(arr: np.ndarray, ab_max: int) -> tuple[
Callable[[LinearConstraint], OptimizeResult],  # binding to milp
sp.csc_matrix,  # A matrix to constraint
np.ndarray,     # lower and upper constraint bounds, z unpopulated
]:
n = len(arr)
M = 2*ab_max*arr.max()

'''
Columns:
a         first term variable, int
b         second term variable, int
ari       = a*r_i, continuous
brj       = b*r_j, continuous
isel (n)  a binary assignment vector
jsel (n)  b binary assignment vector
'''

# Common terms to big-M linearisation constraint rows
ab_rij = (
sp.kron(sp.eye_array(2), -arr[:, np.newaxis]),
sp.kron(sp.eye_array(2), np.ones((n, 1))),
)
msel = sp.diags_array(np.full(shape=2*n, fill_value=M))

A = sp.block_array((
(   # ari + brj == z
[(0, 0)], [(1, 1)], None,
),
(   # Kronecker: each term has one assignment
# sum(isel) == 1, sum(jsel) == 1
None, None, sp.kron(sp.eye_array(2), np.ones((1, n))),
),
(   # Kronecker: lower bound for a*r_i, b*r_j
# ari >= a*arr_value - M*(1 - isel)       -a*arr_value + ari - M*isel >= -M
# brj >= b*arr_value - M*(1 - jsel)       -b*arr_value + brj - M*jsel >= -M
*ab_rij, -msel,
),
(   # Kronecker: upper bound for a*r_i, b*r_j
# ari <= a*arr_value + M*(1 - isel)       -a*arr_value + ari + M*isel <= M
# brj <= b*arr_value + M*(1 - jsel)       -b*arr_value + brj + M*jsel <= M
*ab_rij, msel,
),
), format='csc')

# nan to be replaced with z
constraint_bounds = np.block([
[
np.array((np.nan, 1, 1)),
np.full(shape=2*n, fill_value=-M),
np.full(shape=2*n, fill_value=-np.inf),
],
[
np.array((np.nan, 1, 1)),
np.full(shape=2*n, fill_value=+np.inf),
np.full(shape=2*n, fill_value=+M)
],
])

c = np.zeros(2 + 2 + n + n)
c[:2] = 1  # minimize a + b
integrality = np.ones(shape=(2 + 2 + n + n), dtype=np.uint8)
integrality[2:4] = 0  # only ari/brj are continuous
bounds = Bounds(
#        a       b      ari      brj    sel_ij
lb=(     0,      0, -np.inf, -np.inf) + (0,)*(2*n),
ub=(ab_max, ab_max, +np.inf, +np.inf) + (1,)*(2*n),
)
milp_bind = partial(
milp, c=c, integrality=integrality, bounds=bounds,
)
return milp_bind, A, constraint_bounds

def solve(
arr: np.ndarray,
milp_bind: Callable[[LinearConstraint], OptimizeResult],
A: sp.csc_matrix,
constraint_bounds: np.ndarray,
z: int,
) -> None | tuple[int, int, int, int]:
cbounds = constraint_bounds.copy()
cbounds[:, 0] = z
lba, uba = cbounds
result = milp_bind(
constraints=LinearConstraint(A=A, lb=lba, ub=uba),
)
if not result.success:
return None

(a, b, ari, brj), asel, bsel = np.split(
np.round(result.x).astype(int),
(4, 4 + len(arr)),
)
aval, = arr[asel.nonzero()]
bval, = arr[bsel.nonzero()]
return a, aval, b, bval

def demo() -> None:
arr = np.array((1, 2, 3, 6, 67))
milp_bind, A, constraint_bounds = setup(arr=arr, ab_max=5)

np.set_printoptions(linewidth=200)
print(A.astype(int).toarray())

for z in (1, 4, 30, 7, 88):
result = solve(
arr=arr, milp_bind=milp_bind, A=A, constraint_bounds=constraint_bounds, z=z,
)
if result is None:
print(f'z={z}: no solution')
else:
a, aval, b, bval = result
print(f'{a}*{aval} + {b}*{bval} = {z}')

if __name__ == '__main__':
demo()

[[   0    0    1    1    0    0    0    0    0    0    0    0    0    0]
[   0    0    0    0    1    1    1    1    1    0    0    0    0    0]
[   0    0    0    0    0    0    0    0    0    1    1    1    1    1]
[  -1    0    1    0 -670    0    0    0    0    0    0    0    0    0]
[  -2    0    1    0    0 -670    0    0    0    0    0    0    0    0]
[  -3    0    1    0    0    0 -670    0    0    0    0    0    0    0]
[  -6    0    1    0    0    0    0 -670    0    0    0    0    0    0]
[ -67    0    1    0    0    0    0    0 -670    0    0    0    0    0]
[   0   -1    0    1    0    0    0    0    0 -670    0    0    0    0]
[   0   -2    0    1    0    0    0    0    0    0 -670    0    0    0]
[   0   -3    0    1    0    0    0    0    0    0    0 -670    0    0]
[   0   -6    0    1    0    0    0    0    0    0    0    0 -670    0]
[   0  -67    0    1    0    0    0    0    0    0    0    0    0 -670]
[  -1    0    1    0  670    0    0    0    0    0    0    0    0    0]
[  -2    0    1    0    0  670    0    0    0    0    0    0    0    0]
[  -3    0    1    0    0    0  670    0    0    0    0    0    0    0]
[  -6    0    1    0    0    0    0  670    0    0    0    0    0    0]
[ -67    0    1    0    0    0    0    0  670    0    0    0    0    0]
[   0   -1    0    1    0    0    0    0    0  670    0    0    0    0]
[   0   -2    0    1    0    0    0    0    0    0  670    0    0    0]
[   0   -3    0    1    0    0    0    0    0    0    0  670    0    0]
[   0   -6    0    1    0    0    0    0    0    0    0    0  670    0]
[   0  -67    0    1    0    0    0    0    0    0    0    0    0  670]]
1*1 + 0*1 = 1
0*67 + 2*2 = 4
0*1 + 5*6 = 30
1*6 + 1*1 = 7
z=88: no solution


There is lots of symmetry and short-circuiting you can exploit here. Every loop-iteration is independent and duplicates are pointless to test. E.g a = 1 and b = 2 gives the same result as a = 1 and b = 2. Moreover, we can order the iterations by the sum of a and b and exit early if we have found a result. Like this:

for (int sum = 0; sum < 2 * Z_MAX + 1; sum++) {
int lo = max(sum - Z_MAX, 0);
int hi = min(Z_MAX + 1, sum / 2 + 1);
for (int a = lo; a < hi; a++) {
int b = s - a;
if (a * arr[i] + b * arr[j] == z) {
...
}
}
}


If you work it out on paper, you'll see that the assignments for a and b will be:

(0, 0), (0, 1), (0, 2), (1, 1), (0, 3), (1, 2), (0, 4), ...


That is, they will be in ascending order on their sum. This shaves of a lot of redundant calculation so you want to reorder your loops so that a and b are outermost:

for (int p = 0; p < zArr.length; p++) {
int z = zArr[p];
for (int sum = 0; sum < 2 * Z_MAX + 1; sum++) {
int lo = max(sum - Z_MAX, 0);
int hi = min(Z_MAX + 1, sum / 2 + 1);
for (int a = lo; a < hi; a++) {
int b = s - a;
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
if (a * arr[i] + b * arr[j] == z) {
result[p] = a + b;
break 4 for-loops here... cause you won't
find a better result
}
}
}
}
}
}


You also want to consult the literature on Diophantine equations. Immediately when you have a and b check:

(gcd(a, b) % z) == 0


If it is false, there are no solutions and you can immediately set result[p] = -1. In my testing these two changes improved performance by about 100.