org.apfloat.internal.LongElementaryModMath Maven / Gradle / Ivy
Show all versions of apfloat Show documentation
/*
* MIT License
*
* Copyright (c) 2002-2023 Mikko Tommila
*
* Permission is hereby granted, free of charge, to any person obtaining a copy
* of this software and associated documentation files (the "Software"), to deal
* in the Software without restriction, including without limitation the rights
* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
* copies of the Software, and to permit persons to whom the Software is
* furnished to do so, subject to the following conditions:
*
* The above copyright notice and this permission notice shall be included in all
* copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
* SOFTWARE.
*/
package org.apfloat.internal;
/**
* Elementary modulo arithmetic functions for long data.
*
* Modular addition and subtraction are trivial, when the modulus is less
* than 263 and overflow can be detected easily.
*
* Modular multiplication is more complicated, and since it is usually
* the single most time consuming operation in the whole program execution,
* the very core of the Number Theoretic Transform (NTT), it should be
* carefully optimized.
*
* The algorithm for multiplying two longs and taking the
* remainder is not entirely obvious. The basic problem is to get the
* full 128-bit result of multiplying two 64-bit integers. It would be
* possible to do this by splitting the arguments to high and low 32-bit
* words and performing four multiplications. The performance of this
* solution would be not very good.
*
* Another approach is to use longs only for getting the lowest
* 64 bits of the result. Casting the operands to double and
* multiplying as floating-point numbers, we can get the highest (roughly) 52
* bits of the result. However since only 116 bits can be acquired this
* way, it would be possible to only use 58 bits in each of the multiplication
* operands (not the full 64 or 63 bits). Furthermore, round-off errors in
* the floating-point multiplications, as allowed by the IEEE specification,
* actually prevent getting even 52 of the top bits accurately, and actually
* only 57 bits can be used in the multiplication operands. This is the
* approach chosen in this implementation.
*
* The first observation is that since the modulus is practically
* constant, it should be more efficient to calculate (once) the inverse
* of the modulus, and then subsequently multiply by the inverse modulus
* instead of dividing by the modulus.
*
* The second observation is that to get the remainder of the division,
* we don't necessarily need the actual result of the division (we just
* want the remainder). So, we should discard the topmost 50 bits of the
* full 114-bit result whenever possible, to save a few operations.
*
* The basic approach is to get an approximation of a * b / modulus
* (using floating-point operands, that is doubles). The approximation
* should be within +1 or -1 of the correct result. We first calculate
* a * b - approximateDivision * modulus to get the initial remainder.
* This calculation can use the lowest 64 bits only and is done using longs.
* It is enough to use a double to do the approximate division, as it eliminates
* at least 51 bits from the top of the 114-bit multiplication result, leaving at
* most 63 bits in the remainder. The calculation result - approximateDivision * modulus
* must then be done once more to reduce the remainder since the original multiplication operands
* are only 57-bit numbers. The second reduction reduces the results to the correct value ±modulus.
* It is then easy to detect the case when the approximate division was off by one (and the
* remainder is ±modulus off) as the final step of the algorithm.
*
* @version 1.0
* @author Mikko Tommila
*/
public class LongElementaryModMath
{
/**
* Default constructor.
*/
public LongElementaryModMath()
{
}
/**
* Modular multiplication.
*
* @param a First operand.
* @param b Second operand.
*
* @return a * b % modulus
*/
public final long modMultiply(long a, long b)
{
long r = a * b - this.modulus * (long) ((double) a * (double) b * this.inverseModulus);
r -= this.modulus * (int) ((double) r * this.inverseModulus);
r = (r >= this.modulus ? r - this.modulus : r);
r = (r < 0 ? r + this.modulus : r);
return r;
}
/**
* Modular addition.
*
* @param a First operand.
* @param b Second operand.
*
* @return (a + b) % modulus
*/
public final long modAdd(long a, long b)
{
long r = a + b;
return (r >= this.modulus ? r - this.modulus : r);
}
/**
* Modular subtraction. The result is always >= 0.
*
* @param a First operand.
* @param b Second operand.
*
* @return (a - b + modulus) % modulus
*/
public final long modSubtract(long a, long b)
{
long r = a - b;
return (r < 0 ? r + this.modulus : r);
}
/**
* Get the modulus.
*
* @return The modulus.
*/
public final long getModulus()
{
return this.modulus;
}
/**
* Set the modulus.
*
* @param modulus The modulus.
*/
public final void setModulus(long modulus)
{
this.inverseModulus = 1.0 / modulus;
this.modulus = modulus;
}
private long modulus;
private double inverseModulus;
}