@@ -8714,10 +8714,10 @@ mul_var(const NumericVar *var1, const NumericVar *var2, NumericVar *result,
8714
8714
}
8715
8715
8716
8716
/*
8717
- * If var1 has 1-4 digits and the exact result was requested, delegate to
8717
+ * If var1 has 1-6 digits and the exact result was requested, delegate to
8718
8718
* mul_var_short() which uses a faster direct multiplication algorithm.
8719
8719
*/
8720
- if (var1ndigits <= 4 && rscale == var1 -> dscale + var2 -> dscale )
8720
+ if (var1ndigits <= 6 && rscale == var1 -> dscale + var2 -> dscale )
8721
8721
{
8722
8722
mul_var_short (var1 , var2 , result );
8723
8723
return ;
@@ -8876,7 +8876,7 @@ mul_var(const NumericVar *var1, const NumericVar *var2, NumericVar *result,
8876
8876
/*
8877
8877
* mul_var_short() -
8878
8878
*
8879
- * Special-case multiplication function used when var1 has 1-4 digits, var2
8879
+ * Special-case multiplication function used when var1 has 1-6 digits, var2
8880
8880
* has at least as many digits as var1, and the exact product var1 * var2 is
8881
8881
* requested.
8882
8882
*/
@@ -8898,7 +8898,7 @@ mul_var_short(const NumericVar *var1, const NumericVar *var2,
8898
8898
8899
8899
/* Check preconditions */
8900
8900
Assert (var1ndigits >= 1 );
8901
- Assert (var1ndigits <= 4 );
8901
+ Assert (var1ndigits <= 6 );
8902
8902
Assert (var2ndigits >= var1ndigits );
8903
8903
8904
8904
/*
@@ -8925,6 +8925,13 @@ mul_var_short(const NumericVar *var1, const NumericVar *var2,
8925
8925
* carry up as we go. The i'th result digit consists of the sum of the
8926
8926
* products var1digits[i1] * var2digits[i2] for which i = i1 + i2 + 1.
8927
8927
*/
8928
+ #define PRODSUM1 (v1 ,i1 ,v2 ,i2 ) ((v1)[(i1)] * (v2)[(i2)])
8929
+ #define PRODSUM2 (v1 ,i1 ,v2 ,i2 ) (PRODSUM1(v1,i1,v2,i2) + (v1)[(i1)+1] * (v2)[(i2)-1])
8930
+ #define PRODSUM3 (v1 ,i1 ,v2 ,i2 ) (PRODSUM2(v1,i1,v2,i2) + (v1)[(i1)+2] * (v2)[(i2)-2])
8931
+ #define PRODSUM4 (v1 ,i1 ,v2 ,i2 ) (PRODSUM3(v1,i1,v2,i2) + (v1)[(i1)+3] * (v2)[(i2)-3])
8932
+ #define PRODSUM5 (v1 ,i1 ,v2 ,i2 ) (PRODSUM4(v1,i1,v2,i2) + (v1)[(i1)+4] * (v2)[(i2)-4])
8933
+ #define PRODSUM6 (v1 ,i1 ,v2 ,i2 ) (PRODSUM5(v1,i1,v2,i2) + (v1)[(i1)+5] * (v2)[(i2)-5])
8934
+
8928
8935
switch (var1ndigits )
8929
8936
{
8930
8937
case 1 :
@@ -8936,9 +8943,9 @@ mul_var_short(const NumericVar *var1, const NumericVar *var2,
8936
8943
* ----------
8937
8944
*/
8938
8945
carry = 0 ;
8939
- for (int i = res_ndigits - 2 ; i >= 0 ; i -- )
8946
+ for (int i = var2ndigits - 1 ; i >= 0 ; i -- )
8940
8947
{
8941
- term = ( uint32 ) var1digits [ 0 ] * var2digits [ i ] + carry ;
8948
+ term = PRODSUM1 ( var1digits , 0 , var2digits , i ) + carry ;
8942
8949
res_digits [i + 1 ] = (NumericDigit ) (term % NBASE );
8943
8950
carry = term / NBASE ;
8944
8951
}
@@ -8954,23 +8961,17 @@ mul_var_short(const NumericVar *var1, const NumericVar *var2,
8954
8961
* ----------
8955
8962
*/
8956
8963
/* last result digit and carry */
8957
- term = ( uint32 ) var1digits [ 1 ] * var2digits [ res_ndigits - 3 ] ;
8964
+ term = PRODSUM1 ( var1digits , 1 , var2digits , var2ndigits - 1 ) ;
8958
8965
res_digits [res_ndigits - 1 ] = (NumericDigit ) (term % NBASE );
8959
8966
carry = term / NBASE ;
8960
8967
8961
8968
/* remaining digits, except for the first two */
8962
- for (int i = res_ndigits - 3 ; i >= 1 ; i -- )
8969
+ for (int i = var2ndigits - 1 ; i >= 1 ; i -- )
8963
8970
{
8964
- term = (uint32 ) var1digits [0 ] * var2digits [i ] +
8965
- (uint32 ) var1digits [1 ] * var2digits [i - 1 ] + carry ;
8971
+ term = PRODSUM2 (var1digits , 0 , var2digits , i ) + carry ;
8966
8972
res_digits [i + 1 ] = (NumericDigit ) (term % NBASE );
8967
8973
carry = term / NBASE ;
8968
8974
}
8969
-
8970
- /* first two digits */
8971
- term = (uint32 ) var1digits [0 ] * var2digits [0 ] + carry ;
8972
- res_digits [1 ] = (NumericDigit ) (term % NBASE );
8973
- res_digits [0 ] = (NumericDigit ) (term / NBASE );
8974
8975
break ;
8975
8976
8976
8977
case 3 :
@@ -8982,34 +8983,21 @@ mul_var_short(const NumericVar *var1, const NumericVar *var2,
8982
8983
* ----------
8983
8984
*/
8984
8985
/* last two result digits */
8985
- term = ( uint32 ) var1digits [ 2 ] * var2digits [ res_ndigits - 4 ] ;
8986
+ term = PRODSUM1 ( var1digits , 2 , var2digits , var2ndigits - 1 ) ;
8986
8987
res_digits [res_ndigits - 1 ] = (NumericDigit ) (term % NBASE );
8987
8988
carry = term / NBASE ;
8988
8989
8989
- term = (uint32 ) var1digits [1 ] * var2digits [res_ndigits - 4 ] +
8990
- (uint32 ) var1digits [2 ] * var2digits [res_ndigits - 5 ] + carry ;
8990
+ term = PRODSUM2 (var1digits , 1 , var2digits , var2ndigits - 1 ) + carry ;
8991
8991
res_digits [res_ndigits - 2 ] = (NumericDigit ) (term % NBASE );
8992
8992
carry = term / NBASE ;
8993
8993
8994
8994
/* remaining digits, except for the first three */
8995
- for (int i = res_ndigits - 4 ; i >= 2 ; i -- )
8995
+ for (int i = var2ndigits - 1 ; i >= 2 ; i -- )
8996
8996
{
8997
- term = (uint32 ) var1digits [0 ] * var2digits [i ] +
8998
- (uint32 ) var1digits [1 ] * var2digits [i - 1 ] +
8999
- (uint32 ) var1digits [2 ] * var2digits [i - 2 ] + carry ;
8997
+ term = PRODSUM3 (var1digits , 0 , var2digits , i ) + carry ;
9000
8998
res_digits [i + 1 ] = (NumericDigit ) (term % NBASE );
9001
8999
carry = term / NBASE ;
9002
9000
}
9003
-
9004
- /* first three digits */
9005
- term = (uint32 ) var1digits [0 ] * var2digits [1 ] +
9006
- (uint32 ) var1digits [1 ] * var2digits [0 ] + carry ;
9007
- res_digits [2 ] = (NumericDigit ) (term % NBASE );
9008
- carry = term / NBASE ;
9009
-
9010
- term = (uint32 ) var1digits [0 ] * var2digits [0 ] + carry ;
9011
- res_digits [1 ] = (NumericDigit ) (term % NBASE );
9012
- res_digits [0 ] = (NumericDigit ) (term / NBASE );
9013
9001
break ;
9014
9002
9015
9003
case 4 :
@@ -9021,45 +9009,128 @@ mul_var_short(const NumericVar *var1, const NumericVar *var2,
9021
9009
* ----------
9022
9010
*/
9023
9011
/* last three result digits */
9024
- term = ( uint32 ) var1digits [ 3 ] * var2digits [ res_ndigits - 5 ] ;
9012
+ term = PRODSUM1 ( var1digits , 3 , var2digits , var2ndigits - 1 ) ;
9025
9013
res_digits [res_ndigits - 1 ] = (NumericDigit ) (term % NBASE );
9026
9014
carry = term / NBASE ;
9027
9015
9028
- term = (uint32 ) var1digits [2 ] * var2digits [res_ndigits - 5 ] +
9029
- (uint32 ) var1digits [3 ] * var2digits [res_ndigits - 6 ] + carry ;
9016
+ term = PRODSUM2 (var1digits , 2 , var2digits , var2ndigits - 1 ) + carry ;
9030
9017
res_digits [res_ndigits - 2 ] = (NumericDigit ) (term % NBASE );
9031
9018
carry = term / NBASE ;
9032
9019
9033
- term = (uint32 ) var1digits [1 ] * var2digits [res_ndigits - 5 ] +
9034
- (uint32 ) var1digits [2 ] * var2digits [res_ndigits - 6 ] +
9035
- (uint32 ) var1digits [3 ] * var2digits [res_ndigits - 7 ] + carry ;
9020
+ term = PRODSUM3 (var1digits , 1 , var2digits , var2ndigits - 1 ) + carry ;
9036
9021
res_digits [res_ndigits - 3 ] = (NumericDigit ) (term % NBASE );
9037
9022
carry = term / NBASE ;
9038
9023
9039
9024
/* remaining digits, except for the first four */
9040
- for (int i = res_ndigits - 5 ; i >= 3 ; i -- )
9025
+ for (int i = var2ndigits - 1 ; i >= 3 ; i -- )
9041
9026
{
9042
- term = (uint32 ) var1digits [0 ] * var2digits [i ] +
9043
- (uint32 ) var1digits [1 ] * var2digits [i - 1 ] +
9044
- (uint32 ) var1digits [2 ] * var2digits [i - 2 ] +
9045
- (uint32 ) var1digits [3 ] * var2digits [i - 3 ] + carry ;
9027
+ term = PRODSUM4 (var1digits , 0 , var2digits , i ) + carry ;
9046
9028
res_digits [i + 1 ] = (NumericDigit ) (term % NBASE );
9047
9029
carry = term / NBASE ;
9048
9030
}
9031
+ break ;
9049
9032
9050
- /* first four digits */
9051
- term = (uint32 ) var1digits [0 ] * var2digits [2 ] +
9052
- (uint32 ) var1digits [1 ] * var2digits [1 ] +
9053
- (uint32 ) var1digits [2 ] * var2digits [0 ] + carry ;
9054
- res_digits [3 ] = (NumericDigit ) (term % NBASE );
9033
+ case 5 :
9034
+ /* ---------
9035
+ * 5-digit case:
9036
+ * var1ndigits = 5
9037
+ * var2ndigits >= 5
9038
+ * res_ndigits = var2ndigits + 5
9039
+ * ----------
9040
+ */
9041
+ /* last four result digits */
9042
+ term = PRODSUM1 (var1digits , 4 , var2digits , var2ndigits - 1 );
9043
+ res_digits [res_ndigits - 1 ] = (NumericDigit ) (term % NBASE );
9055
9044
carry = term / NBASE ;
9056
9045
9057
- term = (uint32 ) var1digits [0 ] * var2digits [1 ] +
9058
- (uint32 ) var1digits [1 ] * var2digits [0 ] + carry ;
9059
- res_digits [2 ] = (NumericDigit ) (term % NBASE );
9046
+ term = PRODSUM2 (var1digits , 3 , var2digits , var2ndigits - 1 ) + carry ;
9047
+ res_digits [res_ndigits - 2 ] = (NumericDigit ) (term % NBASE );
9048
+ carry = term / NBASE ;
9049
+
9050
+ term = PRODSUM3 (var1digits , 2 , var2digits , var2ndigits - 1 ) + carry ;
9051
+ res_digits [res_ndigits - 3 ] = (NumericDigit ) (term % NBASE );
9060
9052
carry = term / NBASE ;
9061
9053
9062
- term = (uint32 ) var1digits [0 ] * var2digits [0 ] + carry ;
9054
+ term = PRODSUM4 (var1digits , 1 , var2digits , var2ndigits - 1 ) + carry ;
9055
+ res_digits [res_ndigits - 4 ] = (NumericDigit ) (term % NBASE );
9056
+ carry = term / NBASE ;
9057
+
9058
+ /* remaining digits, except for the first five */
9059
+ for (int i = var2ndigits - 1 ; i >= 4 ; i -- )
9060
+ {
9061
+ term = PRODSUM5 (var1digits , 0 , var2digits , i ) + carry ;
9062
+ res_digits [i + 1 ] = (NumericDigit ) (term % NBASE );
9063
+ carry = term / NBASE ;
9064
+ }
9065
+ break ;
9066
+
9067
+ case 6 :
9068
+ /* ---------
9069
+ * 6-digit case:
9070
+ * var1ndigits = 6
9071
+ * var2ndigits >= 6
9072
+ * res_ndigits = var2ndigits + 6
9073
+ * ----------
9074
+ */
9075
+ /* last five result digits */
9076
+ term = PRODSUM1 (var1digits , 5 , var2digits , var2ndigits - 1 );
9077
+ res_digits [res_ndigits - 1 ] = (NumericDigit ) (term % NBASE );
9078
+ carry = term / NBASE ;
9079
+
9080
+ term = PRODSUM2 (var1digits , 4 , var2digits , var2ndigits - 1 ) + carry ;
9081
+ res_digits [res_ndigits - 2 ] = (NumericDigit ) (term % NBASE );
9082
+ carry = term / NBASE ;
9083
+
9084
+ term = PRODSUM3 (var1digits , 3 , var2digits , var2ndigits - 1 ) + carry ;
9085
+ res_digits [res_ndigits - 3 ] = (NumericDigit ) (term % NBASE );
9086
+ carry = term / NBASE ;
9087
+
9088
+ term = PRODSUM4 (var1digits , 2 , var2digits , var2ndigits - 1 ) + carry ;
9089
+ res_digits [res_ndigits - 4 ] = (NumericDigit ) (term % NBASE );
9090
+ carry = term / NBASE ;
9091
+
9092
+ term = PRODSUM5 (var1digits , 1 , var2digits , var2ndigits - 1 ) + carry ;
9093
+ res_digits [res_ndigits - 5 ] = (NumericDigit ) (term % NBASE );
9094
+ carry = term / NBASE ;
9095
+
9096
+ /* remaining digits, except for the first six */
9097
+ for (int i = var2ndigits - 1 ; i >= 5 ; i -- )
9098
+ {
9099
+ term = PRODSUM6 (var1digits , 0 , var2digits , i ) + carry ;
9100
+ res_digits [i + 1 ] = (NumericDigit ) (term % NBASE );
9101
+ carry = term / NBASE ;
9102
+ }
9103
+ break ;
9104
+ }
9105
+
9106
+ /*
9107
+ * Finally, for var1ndigits > 1, compute the remaining var1ndigits most
9108
+ * significant result digits.
9109
+ */
9110
+ switch (var1ndigits )
9111
+ {
9112
+ case 6 :
9113
+ term = PRODSUM5 (var1digits , 0 , var2digits , 4 ) + carry ;
9114
+ res_digits [5 ] = (NumericDigit ) (term % NBASE );
9115
+ carry = term / NBASE ;
9116
+ /* FALLTHROUGH */
9117
+ case 5 :
9118
+ term = PRODSUM4 (var1digits , 0 , var2digits , 3 ) + carry ;
9119
+ res_digits [4 ] = (NumericDigit ) (term % NBASE );
9120
+ carry = term / NBASE ;
9121
+ /* FALLTHROUGH */
9122
+ case 4 :
9123
+ term = PRODSUM3 (var1digits , 0 , var2digits , 2 ) + carry ;
9124
+ res_digits [3 ] = (NumericDigit ) (term % NBASE );
9125
+ carry = term / NBASE ;
9126
+ /* FALLTHROUGH */
9127
+ case 3 :
9128
+ term = PRODSUM2 (var1digits , 0 , var2digits , 1 ) + carry ;
9129
+ res_digits [2 ] = (NumericDigit ) (term % NBASE );
9130
+ carry = term / NBASE ;
9131
+ /* FALLTHROUGH */
9132
+ case 2 :
9133
+ term = PRODSUM1 (var1digits , 0 , var2digits , 0 ) + carry ;
9063
9134
res_digits [1 ] = (NumericDigit ) (term % NBASE );
9064
9135
res_digits [0 ] = (NumericDigit ) (term / NBASE );
9065
9136
break ;
0 commit comments