A nerdy idea:
A C# class that quicly can supply a lot of primes (here limited to about the largest prime about 50G – as arrays with over 2G entries are not currently supported without splitting the arrays)
A struct that supports basic arithmetic for large accurately represented numbers in Q, by factorization into primes.
A test program that tries to find solutions for a!b!=c! with 1<a<b<c
UPDATE: The text (and code) here has beed modified quite a few times, last time 2016-05-09. (I have changed the basics of the ‘motors’ a few times too…)
Currently the Primes can deliver up to 4G primes (in a few minutes).
(I speeded the prime-part up by making a class with a variant of the Sieve of Atkins – slightly improved. I have also made it dynamical expanding (though that could be optimized!). Or you can chose to allocate a bunch initially. adding the first 26 millions primes less than 22222² can be done in a few second. It can do both a bit array of prime yes/no and/or a list of primes.
Flags can change the behaviour. if the isprime (Bit)Array is wanted (adds some limitations) it can be used rather than just for a temporary array for generating the primelist. I made that in two variants. A fast one, and one with a 64-bit entry array that can handle up to 256G, saving odd only, if a prime LIST is also needed it tops around 50G as the largest prime held in an array holding all lesser prime, due to the 2G entry limit.
Currently the LargeQ has been modified so Addition and Subtraction is faster using BigInteger for adding the remaining denominator after removing already found common factors. The BigInt will be factorized when used for multiplications or division. This is a compromise solution, so problems using mainly plus/minus do not need to use time on the factorisation, on the other hand if the problems (like the below) uses multiplication and division only, things are kept normalised. Note that displaying the value, will lead to normalisation (to display “1/2” rather than “3/6”)
Example: 6/35+14/55 = [2*3]/[5*7]+[2*7]/[5*11] = [2]/[5*7*11] [3*11+7*7] = [2]/[5*7*11] (82) and if needed to [2]/[5*7*11] [2*41] = [2*2*41]/[5*7*11]
Here the 33, 49 and 82 are converted to BigInteger, and only factorised if needed. This makes stuff like repeated A +=B a LOT faster
The LargeQ was developed to work with the below prime-problem, but I’m sorry to admit that the overhead in writing this in .Net does that it is actually a bit faster to just use BigInteger for this problem. See it more as an example, than as a place to use it… Using the LargeQ for it, time is wasted on creating and deleting objects, making the memory consumption fluctuate.
The test problem here, is to look for solutions to a!b!=c! where 1<a<b<c. The obvious/banal ones with a!=c=b+1 are thus found also.
Even on a simple laptop, c up to 40 millions (~11!) are done in about 100 seconds. It tests over a billion c‘s per hour, and currently past 100 billions, and only found 6!7!=10! as a ‘real’ solution. Due to the logarithmic ‘filter’. only about 4 b‘s per million c needs to be accurately tested – so the algorithm used for the full factorials is not that important any more (Bigint or largeQ)
Code, Flags and using
//#define UseLargeQ // Use LargeQ OR BigInteger for the test whether an actual product is a Factorial
//#define UsePreallocatedPrimes
//#define DEBUG
//#define UseGlobalIsPrimeArray //Slightly more efficent, but limits the max to short of 2^31 0x7FFEA80F 46340²-1 / next square could not be Int32 addressed
#define UseGlobalIsPrimeUI64Array //Should be able to hold primes (less 2) up to 8G in 1G memory, with a max of 256G in 16GB memory
//#define StorePrimeList //ONLY omit this with "UseGlobalIsPrimeUI64Array" set
#define UseIsPrimeBitArray //Generally more efficent and spacesaving, but for small arrays slightly slower than boolean array
#define gcAllowVeryLargeObjects //Set if the option is set in config, and on 64bit and need above 2GB data-elements
//See this link https://msdn.microsoft.com/en-us/library/hh285054%28v=vs.110%29.aspx
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Numerics;
using System.Collections;
namespace FactorialProducts
{
using Ntyp = UInt64;
using SqNtyp = UInt64; //due to mem-size, limit to 32...
using IdxTyp = Int32; //The only currently allowed
using Factorisation4 = List<Int64>; // N ^ x, xxxxxxxxxNNNNNNNN
Code, Primes
static class Primes {
//public static Ntyp gcd(Ntyp x, Ntyp y) { if (x == 0 || y == 0) return x + y; else return gcd(y, x % y); }
public static SqNtyp gcd(SqNtyp x, SqNtyp y) {
if (x == 0 || y == 0) return x + y;
do {
if ((x %= y) == 0) return y;
if ((y %= x) == 0) return x;
} while (true);
}
public static SqNtyp lcm(Ntyp x, Ntyp y) { if (x == 0 || y == 0) return 0; else return ((SqNtyp)(x / gcd(x, y))) * y; }
private static bool[] IsPrimeOneToNine = new bool[] { false, false, true, true, false, true, false, true, false, false }; //2,3 and 5 are primes to small for the Sieve of Atkins...
private static Ntyp HighestNTested = 9;
public static Ntyp HighestNTestedForPrimes { get { return HighestNTested; } }
#if (StorePrimeList)
private static List<Ntyp> _PrimesList = new List<Ntyp>() { 2, 3, 5, 7 }; //seed with primes less than 3²
public static List<Ntyp> PrimesList { get { return _PrimesList; } }
#endif
#if (UseGlobalIsPrimeUI64Array)
private static List<UInt64> GlobalOddIsPrimeU64Bits=new List<UInt64>{0xE};// =1110b (that is 7,5,3 and not 1) only set for 1-7 EXCEPT 2 initially...
public static Ntyp? NextPrimeLargerThan(Ntyp N) {
if (N<=1) return 2;
Ntyp M = (N+ (N & 1)) / 2; // This is (N'-1)/2, with N' the next odd N, that is odd add two, even add one
while ((M >> 6) < (UInt64)GlobalOddIsPrimeU64Bits.Count) {
if (0 != (GlobalOddIsPrimeU64Bits[(IdxTyp)(M >> 6)] & (((UInt64)1) << (Byte)(M & 63)))) return M * 2 + 1;
M++;
}
return null; //######??????
}
public static Ntyp? NextPrimeSmallerThan(Ntyp N) {
if (N<=2 || (UInt64)GlobalOddIsPrimeU64Bits.Count< (N >> 7)) return null;
if (N==3) return 2;
Ntyp M = (N-2- (N & 1)) / 2; // This is (N'-1)/2, with N' the previous odd N, that is odd subtract two, even substract one
while (0<M) { //Condition to keep compiler happy
if (0 != (GlobalOddIsPrimeU64Bits[(IdxTyp)(M >> 6)] & (((UInt64)1) << (Byte)(M & 63)))) return M * 2 + 1;
M--;
}
return null; //Keep compiler happy
}
public static UInt64 BitCount(UInt64 x) {
//string xx = x.ToString("X");
//+-------------------------------+
//| 1 | 1 | 1 | 0 | 0 | 1 | 0 | 0 | <-x
//| 1 0 | 0 1 | 0 1 | 0 0 | <-first time merge
//| 0 0 1 1 | 0 0 0 1 | <-second time merge
//| 0 0 0 0 0 1 0 0 | <-third time Half bytes to bytes
//+-------------------------------+ The rest is plain add of bytes
x -= (x >> 1) & 0x5555555555555555; ; //Here lies the 'magic' counting in quarter-bytes
x = (x & 0x3333333333333333) +
( (x >> 2) & 0x3333333333333333);
x = (x + (x >> 4)) & 0x0F0F0F0F0F0F0F0F; // Half Bytes to bytes .n.n.n.n.n.n.n.n
x += (x >> 8); //Every even byte -- ..nn..nn..nn..nn
x += (x >> 16);//Every even byte in every even word -- ......nn......nn
x += (x >> 32);//Every even byte in every even word in every even dword -- ............nn
return x & 0xFF; //leaving only lowest byte
}
public static UInt64 PrimesCount() {
#if (!StorePrimeList)
UInt64 Cnt = 0;
foreach (UInt64 U in GlobalOddIsPrimeU64Bits) Cnt+=BitCount(U) ;
return Cnt;
#else
return _PrimesList.Count;
#endif
}
#else
#if (UseGlobalIsPrimeArray)
#if (UseIsPrimeBitArray)
private static BitArray IsPrimeAtkinDyn = new BitArray (IsPrimeOneToNine);
#else
private static List<bool> IsPrimeAtkinDyn = new List<bool>(IsPrimeOneToNine);
#endif
#endif
#endif
#if (!StorePrimeList)
public static Ntyp IsNull(Ntyp? X, Ntyp Default) {
return X != null ? (Ntyp)X : Default;
}
#else
/// <summary>
/// Prime by number. Note that a dynamic list is maintained so a request for a high number might take a while, populating the list with the intermediate ones
/// </summary>
/// <param name="no">index in primes list, 0 for "2"</param>
public static Ntyp PrimesByNo(IdxTyp no) {
//Ntyp HighestNTested = _PrimesList.Last();
while (_PrimesList.Count<=no && FillPrimesAtkinDynAddAbout25PctPlus1000(HighestNTested)) {;};
while (_PrimesList.Count<=no) {
HighestNTested += 2;
bool skip = false;
foreach (Ntyp p in _PrimesList) { if (HighestNTested < p * p || (skip = gcd(p, HighestNTested) != 1)) break; }
if (!skip) _PrimesList.Add(HighestNTested);
}
return _PrimesList[(Int32)no]; //This will give the needed exception if the callers int-increments wrap around to negative...
}
#endif
/// <summary>
/// Is the number a prime?
/// Behind the scene a dynamic list of primes is maintained, so this could take a while if large primes has to be found first
/// If NOT within the capacity of the list, numbers up the square of the highest cachable prime could be checked by simple check on all primes less √N
/// </summary>
/// <param name="CachableOnly">Only test if within the reach of cache, otherwise return false</param>
public static bool isprime(Ntyp N, bool CachableOnly=false) {
if (HighestNTested<N) if (!FillPrimesAtkinDynAddAbout25PctPlus1000((Ntyp)N)) if (CachableOnly) return false;
#if (UseGlobalIsPrimeUI64Array)
if (N <= HighestNTested) return (N==2 || ((N&1)==1 && 0!=(GlobalOddIsPrimeU64Bits[(IdxTyp)(N >> 7)] & (((UInt64)1) << ((Byte)(N & 126) / 2)))));
#else
if (N <= HighestNTested) return 0 <= _PrimesList.BinarySearch((UInt32)N); //In some contexts, this could speed up things.
#endif
bool skip = false; //Try the simple way, this should work up to the square of the largest stored prime...
IdxTyp ppoi = 0;
Ntyp P=1; //SMALLER than first prime
while (true) {
#if (!StorePrimeList)
P = (Ntyp)NextPrimeLargerThan(P);
#else
P = PrimesByNo(ppoi++);
#endif
if (N < (SqNtyp)P * P || (skip = gcd(P, N) != 1)) break;
}
return !skip;
}
private static UInt64 MaxCacheMB { get {
UInt64 TmpOut =
#if (gcAllowVeryLargeObjects)
2048* sizeof(Ntyp);
#else
2047;
#endif
UInt64.TryParse(System.Configuration.ConfigurationManager.AppSettings.Get("MaxCacheMB"), out TmpOut);
return TmpOut;
} }
#if (StorePrimeList)
private static IdxTyp MaxListLenNTyp //This is a 'constant', but COULD be modified on OM-error
= (IdxTyp) Math.Min( (UInt64)1048576* MaxCacheMB /sizeof(Ntyp), 0x7FFFFFC7);
//0x7FEFFFFF=2G-1M-1, maxvalue allowed for any type array. (Byte-type arrays might be 0x7FFFFFC7=2G-56-1 though)
//see also https://msdn.microsoft.com/en-us/library/hh285054%28v=vs.110%29.aspx
#endif
public static bool FillPrimesAtkinDynAddAbout25PctPlus1000(Ntyp SeedN) {return FillPrimesAtkinDyn(SeedN + (SeedN >> 2)+1000);}
public static bool FillPrimesAtkinDyn(Ntyp NextMax) {
#if (!UseGlobalIsPrimeArray)
bool res=false;
while (HighestNTested < NextMax) {
Ntyp LastHighestNTested = HighestNTested;
//Console.WriteLine("{0:X16} {1:X16} {2:X16} {3:X16}", NextMax, (UInt64)HighestNTested + UInt32.MaxValue, (UInt64)HighestNTested , Int32.MaxValue);
Ntyp NextMaxCalc = (Ntyp)Math.Min(NextMax, (UInt64)HighestNTested+Int32.MaxValue);
res |= FillPrimesAtkinDynInner(NextMaxCalc);
if (LastHighestNTested== HighestNTested) break;
#if (DEBUG && StorePrimeList)
Console.WriteLine(" ## {0}, max {1}", Primes.PrimesList.Count, Primes.PrimesList.Last()); // ###########
#endif
}
return res;
}
public static bool FillPrimesAtkinDynInner(Ntyp limit) {
#endif
// en.wikipedia.org/wiki/Sieve_of_Atkin
// for (i,j) in NxN
//3.1) 4x²+y², y odd : flip when mod 60 in {1,13,17,29,37,41,49,53}
//3.2) 3x²+y², x odd, y even : flip when mod 60 in {7,19,31,43}
//3.3) 3x²-y², 2≤x, 1≤y<x x-y odd : flip when mod 60 in {11,23,47,59}
// for (i,j) in NxN
// Note that:
// (y² % 60) are {0,1,4,9,16,21,24,25,36,40,45,49}
// (y² % 12) are {0,1,4,9} since 0,..,5 gives 0,1,4,9,4,1
/// The fact that (x+6n)² % 12 = (x²+12nx+36n²) % 12 = x² % 12 are use to ignre adds of 6
//on 3.1)
// 4x²+y² mod 60 possible values : 1 5 9 13 17 21 25 29 33 37 41 45 49 53 57
// (4x² % 12) are {0,4} since 0,..,5 gives 0,4,4,0,4,4
// (4x²+y² % 12) are {0,1,4,5,8,9}
// mod 60 {1,13,17,29,37,41,49,53}={1,9,13,17,21,29,33,37,41,45,49,53,57}\{9+12i}={1+4j}\{5+20k}\{9+12i}
//on 3.2)
// 3x²+y² mod 60 possible values : 1 5 9 13 17 21 25 29 33 37 41 45 49 53 57
// (3x² % 12) are {0,3} (multipli 0,1,4,9 by 3, and mod 12)
// mod 60 (7,19,31,43} = {7+12i}\{55+60j}, and as 55 not possible no need to exclude
// So only when (y² % 12) = 4 we can get 7, so that is the start limit
//on 3.3)
// 3x²-y² mod 60 possible values : 3 11 15 23 27 35 39 47 51 59
// 3x²-y², lowest value for fixed y when x=y+1 is 2y²+6y+3, so we start at {0,1,4,9}*2+3+n*6 % 12 ={3,5,11,9}
// (3x² % 12) are {0,3} (multiply 0,1,4,9 by 3, and mod 12)
// (11,23,47,59) mod 60 -> (11-12i)\(33+*60j), and as 33 not possible no need to exclude
// So only when (y² % 12) = 11 we can get 11, so that is the start limit, and as +3 is not wanted only every second x is needed, adding 0 mod 12 every time
#if (StorePrimeList)
if (MaxListLenNTyp <= _PrimesList.Count) return false; //List already full....
#endif
UInt64 LastMax = HighestNTested;
#if (UseGlobalIsPrimeUI64Array)
limit = limit|127; //Round op to one below closest 127 bit border, that is "...x1111111" Note the last bit is NOT used as only Odd primes are saved
if (0x7FFFFFC7 <= (limit-LastMax)) limit -= 128; // Too many added in the line above, back down...
#endif
#if (UseGlobalIsPrimeArray)
UInt64 ArrOfs = 0;
limit = Math.Min(limit, 0x7FFFFFC7);
#else
UInt64 ArrOfs = LastMax + 1;
#if (UseIsPrimeBitArray)
BitArray IsPrimeAtkinDyn = new BitArray(0);
#else
List<bool> IsPrimeAtkinDyn = new List<bool>(){};
#endif
#endif
if (limit <= HighestNTested) return false; //Already got enough....
Ntyp sqrtNmax = (Ntyp)Math.Ceiling(Math.Sqrt(limit)); //Round up
#if (UseIsPrimeBitArray)
IsPrimeAtkinDyn.Length += (int)(limit - LastMax); for (SqNtyp i = LastMax + 1; i <= limit; i++) IsPrimeAtkinDyn[(int)(i - ArrOfs)] = false;
#else
IsPrimeAtkinDyn.AddRange(Enumerable.Repeat(false, (int)(limit - LastMax)));
#endif
UInt64 y2;
Ntyp x;
// NOTE: (y2 % 12) are {0,1,4,9} since 0,..,5 gives 0,1,4,9,4,1 and (x+6n)² % 12 = (x²+12nx+36n²) % 12 = x² % 12
//// in the increments below: Delta x² from x-1 to x is (x-1)²-x² = 2x-1 , used with 4x or 3x so 8x-4 or 6x-3
UInt64 nstart;
for (Ntyp y = 1; y <= sqrtNmax; y++) { //First value is y²+4, so Floor(Sqrt(LastMax-4)) should be a good start
y2 = (UInt64)y * y;
x = 1; //4x²+y² , 3.1) see above
nstart = 4 + y2;
if ((y2 % 4) == 1) // 1,5 or 9 mod 12, since nn only can adjust start with 0 or 8 mod 12
//for (SqNtyp nn = 4 + y2; nn <= limit; nn += (SqNtyp)8 * ++x - 4) if ((nn % 12) != 9 && (nn % 20) != 5) IsPrime[nn] ^= true;
for (UInt64 nn = nstart; nn <= limit; nn += (UInt64)8 * ++x - 4)
if (LastMax < nn) if ((nn % 12) != 9) IsPrimeAtkinDyn[(int)(nn-ArrOfs)] ^= true;
}
for (Ntyp y = 1; y <= sqrtNmax; y++) { //First value is y²+3, so Floor(Sqrt(LastMax-3)) should be a good start
y2 = (UInt64)y * y;
x = 1; //3x²+y² , 3.2) see above
nstart = 3 + y2;
if ((y2 % 12) == 4) // 7 +/- 3 as nn only can adjust start with 0 or 3 mod 12, so only four from {0,1,4,9} can give 7
for (UInt64 nn = nstart; nn <= limit; nn += (UInt64)12 * ++x, x++)
if (LastMax < nn) IsPrimeAtkinDyn[(int)(nn - ArrOfs)] ^= true; //as only the +3 not the +0 can give 7, only every second is relevant!!
}
for (Ntyp y = 1; y <= sqrtNmax; y++) { //First value is 2y²+6y+3=2(y+1½)²-1½, so Floor(Sqrt((LastMax+1.5)/2)-1.5) should be a good start
y2 = (UInt64)y * y;
x = y + 1; //3x²-y², so first is 2y²+6y+3 , 3.3) see above
nstart = 2 * y2 + 6 * y + 3; //{0,1,4,9}*2+3+n*6 % 12 ={3,5,11,9}
if ((nstart % 12) == 11) // 11 +/- 3 as nn only can adjust start with 0 or 3 mod 12, so only 11 from {3,5,11,9} can give 11
for (UInt64 nn = nstart; nn <= limit; nn += (UInt64)12 * ++x, x++)
if (LastMax < nn) IsPrimeAtkinDyn[(int)(nn-ArrOfs)] ^= true; //as only the +3 not the +0 can give 11, only every second is relevant!!
}
UInt64 n2;
#if (!StorePrimeList)
Ntyp P = 3; //Start by testing with 5 (=Next prime from 3)
while (P <= sqrtNmax) {
if (LastMax<P) P+=2; //Simple guess, next odd
else P = IsNull(NextPrimeLargerThan(P) , (LastMax + 1)|1 );
if (LastMax<P) //NO! NOT part of previous IF !!!
if (!IsPrimeAtkinDyn[(int)(P - ArrOfs)]) continue;
n2 = (SqNtyp)P * P;
for (UInt64 q = (((UInt64)LastMax) / n2 + 1) * (n2); q <= limit; q += n2) IsPrimeAtkinDyn[(int)(q - ArrOfs)] = false;
}
#else
int ppoi = 2; //5 is index 2 on the zerobased primelist
for (Ntyp n = 5; n <= sqrtNmax; n++) {
if (n <= LastMax) {
if (_PrimesList.Last()<n || n < _PrimesList[ppoi]) continue; //not on primelist
ppoi++;
} else if (!IsPrimeAtkinDyn[(int)(n - ArrOfs)]) continue;
n2 = (SqNtyp)n * n;
for (UInt64 q = (((UInt64)LastMax) / n2 + 1) * (n2); q <= limit; q += n2) IsPrimeAtkinDyn[(int)(q - ArrOfs)] = false;
}
int PrimesAddCou = 0;
//Pass 1
for (SqNtyp n = HighestNTested + 1; n <= limit; n++)
if (IsPrimeAtkinDyn[(int)(n - ArrOfs)]) {
if (MaxListLenNTyp <= _PrimesList.Count+ PrimesAddCou) { limit = n-1; break; }
PrimesAddCou++;
}
//Extend if needed and possible
int OM_Reduce = Math.Max(0,PrimesList.Count - MaxListLenNTyp + PrimesAddCou);
#if (UseGlobalIsPrimeUI64Array)
OM_Reduce= (OM_Reduce+127)&(-128); //Round up to multipla of 128
#endif
while ( _PrimesList.Capacity < PrimesList.Count + PrimesAddCou - OM_Reduce)
try { _PrimesList.Capacity = PrimesList.Count + PrimesAddCou - OM_Reduce; }
catch (OutOfMemoryException) { OM_Reduce += 1000000;} //Try to backout a few MB
if (0 < OM_Reduce) {
MaxListLenNTyp = PrimesList.Count - OM_Reduce + PrimesAddCou;
Console.WriteLine("\r Out of cache space, further caching is stopped a few MB short");
}
//Pass 2
for (SqNtyp n = HighestNTested + 1; n <= limit; n++) {
if (IsPrimeAtkinDyn[(int)(n - ArrOfs)]) {
if (MaxListLenNTyp <= _PrimesList.Count ) { limit = n-1; break; }
_PrimesList.Add((Ntyp)n);
}
}
if (MaxListLenNTyp <= _PrimesList.Count)
Console.WriteLine("\r (Prime-cache is full {0} primes, 2-{1})", _PrimesList.Count, _PrimesList.Last());
#endif
#if (UseGlobalIsPrimeUI64Array)
GlobalOddIsPrimeU64Bits.AddRange(Enumerable.Repeat((UInt64)0, (int)((limit - HighestNTested)>>7))) ;
for (SqNtyp n = HighestNTested + 1; n <= limit; n++)
if (IsPrimeAtkinDyn[(int)(n - ArrOfs)])
GlobalOddIsPrimeU64Bits[(IdxTyp)(n >> 7)] |= (((UInt64)1) << ((Byte)(n & 126)/2));
#endif
bool res = HighestNTested != (SqNtyp)limit;
HighestNTested = (Ntyp)limit;
return res;
}
}
////
// http://webprimes.com/?index=1073741824
// prime[1073741824] 24563311309
Code, LargeQ
public struct LargeQ
{
public BigInteger SignMM { get; private set; }
public int Sign{get {return SignMM.Sign;}} //-1, 0, +1
private Factorisation4 _factorisation4;// null: Zero , empty list One....
static public LargeQ Zero() { return new LargeQ((Int64)0); }
private void NormaliseToFactors() {
if (SignMM==0) return; //0
BigInteger temp = BigInteger.Abs(SignMM);
if (temp==1) return; // +/- 1: already factorised
//Factorise...
Factorisation4 Fact = new Factorisation4 { };
int ppoi = 0;
Ntyp P=1; //SMALLER than first prime
BigInteger DivRes, rem;
while (1 < temp) {
#if (!StorePrimeList)
P = (Ntyp)Primes.NextPrimeLargerThan(P);
#else
P = Primes.PrimesByNo(ppoi);
#endif
int pow = 0;
while (1 < temp) {
DivRes=BigInteger.DivRem(temp,P, out rem);
if (rem == 0) {
pow++;
temp = DivRes;
} else break;
};
if(pow!=0)Fact.Add(((Int64)pow<<32)|(UInt32)P);
ppoi++;
}
SignMM = Sign; // +/- 1
FactListMultiplyPow(_factorisation4, Fact,+1);
}
/// <summary>
/// Create a number from an BigInteger
/// </summary>
public LargeQ(BigInteger x) : this(){
_factorisation4 = new Factorisation4 { }; //SignMM * 1
SignMM = x;
}
//public LargeQ(Int64 x) : this((BigInteger)x){}
public LargeQ(Int64 x) : this(){
_factorisation4 = new Factorisation4 { }; //Sign * 1
SignMM = Math.Sign(x);
if (x == 0) {return;}
//Factorise...
UInt64 temp = (UInt64)(x * Sign);
int ppoi = 0;
Ntyp P=1; //SMALLER than first prime
while (1 < temp) {
#if (!StorePrimeList)
P = (Ntyp)Primes.NextPrimeLargerThan(P);
#else
P = Primes.PrimesByNo(ppoi++);
#endif
int pow = 0;
while (1 < temp && temp % P == 0) {
pow++;
temp /= P;
};
if (pow != 0) _factorisation4.Add(((Int64)pow<<32)|(UInt32)P);
}
}
//public LargeQ(Int32 x) : this((BigInteger)x){}
public LargeQ(Int32 x) : this((Int64)x) { }
static public implicit operator LargeQ(Int32 value) { return new LargeQ(value); }
static public implicit operator Int32(LargeQ value) { return (Int32)value.deFactorizeI(); } //Might loose some...
static public implicit operator LargeQ(Int64 value) { return new LargeQ(value); }
static public implicit operator Int64(LargeQ value) { return value.deFactorizeI(); }
static public implicit operator LargeQ(BigInteger value) { return new LargeQ(value); }
static public implicit operator BigInteger(LargeQ value) { return value.deFactorizeBi(); }
static public implicit operator double(LargeQ value) { return value.deFactorizeDbl(); }
public bool isone() {
CleanTail();
return _factorisation4.Count==0 && SignMM==1;
}
public bool isint(){
CleanTail();
foreach (Int64 FNP in _factorisation4) if (FNP < 0) return false; //x in Q\N
return true;
}
/// <summary>
/// Clone the number
/// </summary>
public LargeQ(LargeQ x) : this() {
_factorisation4 = new Factorisation4(x._factorisation4);
SignMM = x.SignMM;
}
/// <summary>
/// Attempt to create the Int64 representing the number, might give overflow, and will be inaccurate for non integral numbers.
/// </summary>
/// <returns></returns>
public Int64 deFactorizeI() {
Int64 res = 1;
for (int ppoi = 0; ppoi < _factorisation4.Count; ppoi++)
for (int i = 1; i <= (Int32)(_factorisation4[ppoi]>>32); i++) res *= (_factorisation4[ppoi] & 0xFFFFFFFF);
for (int ppoi = 0; ppoi < _factorisation4.Count; ppoi++)
for (int i = (Int32)(_factorisation4[ppoi]>>32); i < 0; i++) res /= (_factorisation4[ppoi] & 0xFFFFFFFF);
return res * ((1 < BigInteger.Abs(SignMM)) ? (Int64)SignMM : Sign); //Might overflow
}
/// <summary>
/// Create the BigInteger representing the number, will be inaccurate for non integral numbers.
/// </summary>
/// <returns></returns>
public BigInteger deFactorizeBi() {
BigInteger res = 1;
for (int ppoi = 0; ppoi < _factorisation4.Count; ppoi++)
for (int i = 1; i <= (Int32)(_factorisation4[ppoi]>>32); i++) res *= (_factorisation4[ppoi] & 0xFFFFFFFF);
for (int ppoi = 0; ppoi < _factorisation4.Count; ppoi++)
for (int i = (Int32)(_factorisation4[ppoi]>>32); i < 0; i++) res /= (_factorisation4[ppoi] & 0xFFFFFFFF);
return res * ((1<BigInteger.Abs(SignMM))? SignMM:Sign);
}
/// <summary>
/// Log10(number) as double. This is usualy possible to calculate (for positive numbers of course) without owerflow
/// </summary>
/// <returns></returns>
public double Log10deFactorizeDbl() {
if (SignMM<=0) throw new ArithmeticException("Logarithm to none positive number");
double res = 0.0;
if (1<SignMM) res = BigInteger.Log10(SignMM); //We have checked that it is positive.
for (int ppoi = 0; ppoi < _factorisation4.Count; ppoi++) res += (Int32)(_factorisation4[ppoi]>>32) * Math.Log10(_factorisation4[ppoi] & 0xFFFFFFFF);
return res;
}
/// <summary>
/// Convert the number to double, often possible, but may well loose accuracy, especially for large numbers
/// </summary>
/// <returns></returns>
public double deFactorizeDbl() {
if (Sign==0) return 0;
double res = 0.0;
for (int ppoi = 0; ppoi < _factorisation4.Count; ppoi++) res += (Int32)(_factorisation4[ppoi]>>32) * Math.Log10(_factorisation4[ppoi] & 0xFFFFFFFF);
res=Math.Pow(10,res); //May throw exception
return res * ((1 < BigInteger.Abs(SignMM)) ? (double)SignMM : Sign);
}
private static void FactListMultiplyPow(Factorisation4 XF, Factorisation4 YF, Int32 k) {
int ComLng = Math.Min(XF.Count, YF.Count);
if (YF.Count == 0) return;
UInt32 PX; Int32 VX = 0;
UInt32 PY; Int32 VY = 0;
int poiX = XF.Count - 1; int poiY = YF.Count - 1;
while (0<=poiX || 0<=poiY ) {
if (poiX < 0) PX = 0; else { VX = (Int32)(XF[poiX]>>32); PX=(UInt32)XF[poiX]; }
if (poiY < 0) PY = 0; else { VY = (Int32)(YF[poiY]>>32); PY=(UInt32)YF[poiY]; }
if (PY < PX) { poiX--; }
else if (PX < PY) { XF.Insert(poiX+1,(((Int64)(VY*k))<<32)+PY); poiY--;}
else if (PX < PY) { XF.Insert(poiX+1,YF[poiY]); poiY--;}
else {//(PX== PY)
Int32 XY = VX + VY*k;
if (XY==0) XF.RemoveAt(poiX); else XF[poiX]=((Int64)XY<<32)+PX;
poiX--; poiY--;
}
}
}
//private static void FactListMultiply(Factorisation4 XF, Factorisation4 YF) {
// FactListMultiplyPow(XF, YF, +1);
//}
//private static void FactListDivide(Factorisation4 XF, Factorisation4 YF) {
// FactListMultiplyPow(XF, YF, -1);
//}
/// <summary>
/// Multiply the element with another raised to power. X = X * (Y^k)
/// </summary>
/// <param name="x"></param>
/// <param name="y">The Y that X is multipied/divided with a suitable number of times</param>
/// <param name="k">Power, eg +1 and -1 for plain multiplikation and division</param>
/// <returns></returns>
public static LargeQ MultiplyPow(LargeQ x, LargeQ y, Int32 k) {
if (k<0 && y.Sign==0) throw new DivideByZeroException();
if (x.Sign==0|| y.Sign==0) return Zero();
x.NormaliseToFactors();
y.NormaliseToFactors();
LargeQ res = new LargeQ(x);
FactListMultiplyPow(res._factorisation4, y._factorisation4, k);
res.SignMM = x.Sign*(k%1==1?y.Sign:1); //if the sign of x swapped only if k is odd and y<0
return res;
}
public static LargeQ operator *(LargeQ x, LargeQ y) {
return MultiplyPow(x, y, +1);
//if (x.Sign==0|| y.Sign==0) return Zero();
//x.NormaliseToFactors();
//y.NormaliseToFactors();
//LargeQ res = new LargeQ(x);
//FactListMultiply(res._factorisation4, y._factorisation4);
//res.SignMM = x.Sign*y.Sign;
//return res;
}
public static LargeQ operator /(LargeQ x, LargeQ y) {
return MultiplyPow(x, y, -1);
//if (x.Sign==0) return Zero();
//if (y.Sign==0) throw new DivideByZeroException();
//x.NormaliseToFactors();
//y.NormaliseToFactors();
//LargeQ res = new LargeQ(x);
//FactListDivide(res._factorisation4, y._factorisation4);
//res.SignMM = x.Sign * y.Sign;
//return res;
}
/// <summary>
/// Remove zero-entries from the tail of the primefactor lists
/// </summary>
private void CleanTail(){
if (Sign==0 && 0 < _factorisation4.Count) _factorisation4.Clear();
NormaliseToFactors();
}
public bool isdivisor(LargeQ y){
//return ((this/y).isint());
NormaliseToFactors();
y.NormaliseToFactors();
if (y._factorisation4.Count == 0) return true;
Factorisation4 XF = _factorisation4;
Factorisation4 YF = y._factorisation4;
UInt32 PX; Int32 VX = 0;
UInt32 PY; Int32 VY = 0;
int poiX = XF.Count - 1; int poiY = YF.Count - 1;
while (0<=poiX || 0<=poiY ) {
if (poiX < 0) PX = 0; else { VX = (Int32)(XF[poiX]>>32); PX=(UInt32)XF[poiX]; }
if (poiY < 0) PY = 0; else { VY = (Int32)(YF[poiY]>>32); PY=(UInt32)YF[poiY]; }
if (PY < PX) { poiX--; }
else if (PX < PY) {if( 0<VY)return false; poiY--; }
else {//(PX== PY)
if(VX<VY)return false;
poiX--; poiY--;
}
}
return true;
}
public static bool operator ==(LargeQ x, LargeQ y) {
if (x.Sign != y.Sign) return false;
x.NormaliseToFactors();
y.NormaliseToFactors();
if (x._factorisation4.Count != y._factorisation4.Count) return false;
for (int ppoi = 0; ppoi < x._factorisation4.Count; ppoi++)
if (x._factorisation4[ppoi] != y._factorisation4[ppoi]) return false;
return true;
}
public static bool operator !=(LargeQ x, LargeQ y) {return !(x == y);}
public override bool Equals(object obj) { return obj==null ? false : this == (LargeQ)obj; }
public override int GetHashCode() { return this._factorisation4.GetHashCode(); }
public static LargeQ Factorial(Int64 x) {
LargeQ Res = x;
while (1 < x) Res *= (--x);
return Res;
}
/// <summary>
/// Check if the element could be a!, start without calculating a!, by simply checking for a minimal occurence of each prime p≤a
/// every 1 in w has w as a divisor, so for each w the occurrence nw of w in a! will abide a+1≤w(nw+1), this also go for every prime.
/// </summary>
/// <param name="a"></param>
/// <returns></returns>
public bool IsFactorial(Int64 a) {
Int64 P =2; IdxTyp ppoi =0;
Int32 pow;
while (P<=a) {
if (_factorisation4.Count<=ppoi || ((UInt32)_factorisation4[ppoi] & 0xFFFFFFFF) !=P) return false; //all primes must be present
pow = (Int32)(_factorisation4[ppoi]>>32);
if (P*(pow+1)<a+1) return false; //every 1 in p number has p as a divisor, at the least once, this goes for all p.
ppoi++;
#if (!StorePrimeList)
P = (Int64)Primes.NextPrimeLargerThan((Ntyp)P);
#else
P = (Int64)Primes.PrimesByNo(ppoi);
#endif
}
return (Factorial(a) == this);
}
public static LargeQ Factorial(UInt32 x)
{
return Factorial((Int64)x);
//LargeQ Res = new LargeQ((Int64)x);
//while (1<x) Res*=((Int64)(--x));
//return Res;
}
public static LargeQ operator -(LargeQ x) {
LargeQ res = new LargeQ(x);
res.SignMM = -res.SignMM;
return res;
}
public static LargeQ operator +(LargeQ x) { return x; }
/// <summary>
/// Split a number into its integral denominator and nominator
/// </summary>
/// <param name="denominator"></param>
/// <param name="nominator"></param>
public void splitFraction(out LargeQ denominator, out LargeQ nominator) {
denominator = new LargeQ(this);
if (this.isint()) { nominator = 1; return; }
denominator.NormaliseToFactors();
nominator = new LargeQ(denominator);
nominator.SignMM=+1;
for(int i=denominator._factorisation4.Count-1;0<=i;i--) if ( denominator._factorisation4[i]<0) denominator._factorisation4.RemoveAt(i);
for(int i= nominator._factorisation4.Count-1;0<=i;i--) if (0<=nominator._factorisation4[i] ) nominator._factorisation4.RemoveAt(i);
else nominator._factorisation4[i] = (nominator._factorisation4[i] & 0xFFFFFFFF)*2
-nominator._factorisation4[i]; //Swap to positive! X+P -> 2P-(X+P) = -X + P for the prime P and power X
denominator.CleanTail();
nominator.CleanTail();
}
public void splitFraction(out BigInteger denominatorBi, out BigInteger nominatorBi) {
LargeQ denominator;
LargeQ nominator;
splitFraction(out denominator, out nominator);
denominatorBi = denominator;
nominatorBi = nominator;
}
public override string ToString() {
return ToString("R");
}
public string ToString(string format) {
//return "####DEBUG_DISABLED####"; //####DEBUG#### to avoid the normalisition to be done implicit
if (isint()) return (deFactorizeBi().ToString());
BigInteger denominator;
BigInteger nominator;
BigInteger IntegralPart=0;
splitFraction(out denominator, out nominator);
bool AbsOverOne=BigInteger.Abs(nominator)<BigInteger.Abs(denominator);
switch (format){
case "G": format=(AbsOverOne?"{2:R};":"")+"{0:R}/{1:R}"; break;
case "R": format= "{0:R}/{1:R}"; break;
}
if (AbsOverOne && 0<=format.IndexOf("{2")) {
IntegralPart=BigInteger.DivRem(denominator, nominator, out denominator);
denominator =BigInteger.Abs(denominator);
}
return String.Format(format, denominator, nominator, IntegralPart);
}
/// <summary>
/// Largest common divisor. including fractions. 6/77 and 26/35 returns 2/(5*7*11), BUT ignoring UNNORMALISED factors for speed!!
/// </summary>
public static LargeQ gcd_NORMALIZEDONLY(LargeQ x, LargeQ y) {
if (x.Sign==0 && y.Sign==0) return 1; //A matter of definition.... Exception?
if (x.Sign==0) return new LargeQ(y);
if (y.Sign==0) return new LargeQ(x);
//#### OBS OBS Ignore the value part of SignMM;
LargeQ res = 1;
Factorisation4 XF = x._factorisation4;
Factorisation4 YF = y._factorisation4;
int poiX = 0; int poiY = 0;
UInt32 PX; Int32 VX = 0;
UInt32 PY; Int32 VY = 0;
while (poiX<XF.Count || poiY<YF.Count) {
if (XF.Count<=poiX) PX = UInt32.MaxValue; else { VX = (Int32)(XF[poiX] >> 32); PX = (UInt32)XF[poiX]; }
if (YF.Count<=poiY) PY = UInt32.MaxValue; else { VY = (Int32)(YF[poiY] >> 32); PY = (UInt32)YF[poiY]; }
if (PX < PY) { {if((Int32)VX<0) res._factorisation4.Add(XF[poiX]); }; poiX++; }
else if (PY < PX) { {if((Int32)VY<0) res._factorisation4.Add(YF[poiY]); }; poiY++; }
else {VX=Math.Min(VX,VY); if(VX!=0){ res._factorisation4.Add((((Int64)VX)<<32)+PX); }; poiX++; poiY++; }
}
return res;
}
/// <summary>
/// Largest common divisor. including fractions. 6/77 and 26/35 returns 2/(5*7*11), dividing with this leaves integers (3*5 and 11*13)
/// </summary>
public static LargeQ gcd(LargeQ x, LargeQ y) {
x.NormaliseToFactors();
y.NormaliseToFactors();
return gcd_NORMALIZEDONLY(x, y);
}
public static LargeQ operator +(LargeQ x, LargeQ y) {
//Simplified by dividing out largest common divisor before the add, as less needs to be expanded for the sum
// Method example: 6/35+14/55 = [2*3]/[5*7]+[2*7]/[5*11] = [2/5] [3/7+7/11] = [2]/[5*7*11] [3*11+7*7] = [2]/[5*7*11] (82) = [2]/[5*7*11] [2*41] = [2*2*41]/[5*7*11]
LargeQ Res = gcd_NORMALIZEDONLY(x, y); //Note this INCLUDES both fractional parts for fractions
LargeQ x_r = new LargeQ(x); FactListMultiplyPow(x_r._factorisation4, Res._factorisation4,-1); //The result is a (potentially large) integer
LargeQ y_r = new LargeQ(y); FactListMultiplyPow(y_r._factorisation4, Res._factorisation4,-1); //do
//LargeQ x_r = new LargeQ(x); FactListDivide(x_r._factorisation4, Res._factorisation4); //The result is a (potentially large) integer
//LargeQ y_r = new LargeQ(y); FactListDivide(y_r._factorisation4, Res._factorisation4); //do
Res.SignMM = ((BigInteger)x_r) + ((BigInteger)y_r); //If large BigInteger is needed, the decoding back to LargeQ is going to be slow.
return Res;
}
public static LargeQ operator -(LargeQ x, LargeQ y) {
return x+(-y);
}
}
Code, Main with test
class Program
{
static List<double> SortedLog10nFac = new List<double> { 0 };
static void FillSortedLog10nFacUpToLog10X(double Log10X)
{
while (SortedLog10nFac.Last() < Log10X) SortedLog10nFac.Add(SortedLog10nFac.Last() + Math.Log10(SortedLog10nFac.Count));
}
public static BigInteger FactorialBi(Ntyp x)
{
BigInteger Res = 1;
while (1 < x) Res *= (x--);
return Res;
}
static void Main(string[] args)
{
#if (DEBUG)
//Console.WriteLine(Environment.Is64BitProcess);
for (int q1 = 3; q1 <= 1003; q1 += 10) Primes.FillPrimesAtkinDyn((Ntyp)(q1 * q1));
Primes.FillPrimesAtkinDyn((Ntyp)(1003 * 1003));
#if (StorePrimeList)
Int32 XX = Primes.PrimesList.Count; //78941
Ntyp YY = Primes.PrimesList.Last(); //1006007
#else
Ntyp? YY = Primes.NextPrimeSmallerThan(8); //1006007
YY = Primes.NextPrimeSmallerThan((Ntyp)YY); //5
YY = Primes.NextPrimeSmallerThan((Ntyp)YY); //3
YY = Primes.NextPrimeSmallerThan((Ntyp)YY); //2
YY = Primes.NextPrimeSmallerThan((Ntyp)YY); //Null
YY = Primes.HighestNTestedForPrimes; //1006079 10003^2=1006009 rounded up to nearest 128 minus 1
YY = Primes.NextPrimeSmallerThan(Primes.HighestNTestedForPrimes); //1006063
YY = Primes.PrimesCount();//78943
#endif
bool ZZ = Primes.isprime(9973); //true; IsPrimeAtkinDyn.Count=1006010 = 1003*1003+1
//for (UInt32 i = 23171; i < 999999;i+=1 ) Primes.FillPrimesAtkinDyn(i * i);
LargeQ Xq = ((new LargeQ(77) * 99));
LargeQ Yq = Xq + Xq + Xq + Xq;
Yq -= Xq;
Console.WriteLine(Yq.isdivisor(11));
Yq /= 11;
Console.WriteLine(Yq.isdivisor(11));
Yq /= 11;
Console.WriteLine(Yq.isdivisor(11));
Yq /= 11;
Console.WriteLine(Yq.ToString());
Yq *= 11;
Xq = (LargeQ)1 / 2;
Console.WriteLine(Xq.ToString());
Console.WriteLine(Xq.ToString());
Console.WriteLine(Xq.ToString());
Yq = (LargeQ)1 / 6;
Console.WriteLine(Yq.ToString());
Console.WriteLine(Yq.ToString());
Console.WriteLine(Yq.ToString());
Xq -= Yq;
Xq -= Yq;
Xq -= 7;
Console.WriteLine(Xq.ToString());
Console.WriteLine(Xq.ToString("R"));
Console.WriteLine(Xq.ToString("G"));
Xq -= Yq;
Xq += Yq;
Xq += 7;
Xq += Yq;
Xq += Yq;
Xq += Yq;
// Method example: 6/35+14/55 = [2*3]/[5*7]+[2*7]/[5*11] = [2]/[5*7*11] [3*11+7*7] = [2]/[5*7*11] (82) = [2]/[5*7*11] [2*41] = [2*2*41]/[5*7*11]
LargeQ A = (LargeQ)6 / 35;
LargeQ B = (LargeQ)14 / 55;
LargeQ AB = A * B;
LargeQ AdB = A / B;
LargeQ C = A + B;
LargeQ D = C * (LargeQ)(5 * 7 * 11);
Int64 Di = D;
LargeQ E = A - B;
LargeQ F = E * (LargeQ)(5 * 7 * 11);
Int64 Fi = F;
#endif
//Console.WriteLine((Int64)LargeQ.Factorial(5));
//Console.WriteLine((Int64)LargeQ.Factorial(20));
//Console.WriteLine((Int64)(LargeQ.Factorial(123456) / LargeQ.Factorial(123455)));
////Check that the inverting-algorithm works, even with 10 times narrower limit
//double XXX = 0;// Math.Log10(2);
//for (int chk=2;chk<33000000;chk++) {
// XXX += Math.Log10(chk);
// FillSortedLog10nFacUpToLog10X(XXX + 0.01);
// int Ix = SortedLog10nFac.BinarySearch(XXX - 0.01); //-0.01 to ensure that we do not get above the real number due to numerical errors
// if (Ix < 0) Ix = ~Ix; //First element larger, minus one=last element smaller
// if (Math.Abs(SortedLog10nFac[Ix] - XXX) < 0.000001 || Ix!=chk)
// { int x = chk; }
// else
// Console.WriteLine("Reverse algorithm failed for {0}!, found IX", chk, Ix);
// if(chk%1000==0)Console.Write(chk.ToString() + "\r");
//};
//Console.Write(" \r");
Console.WriteLine("Solutions for a!b!=c! where 1<a<b<c");
#if (UsePreallocatedPrimes)
Console.Write("Primes preload");
Primes.FillPrimesAtkinDyn(0x7FFEA80F); //46340²-1 adressing within 2Gb
while (true)
{
Console.Write(".");
Ntyp LastHighestNTested = Primes.PrimesList.Last();
if (!Primes.FillPrimesAtkinDynAddAbout25PctPlus1000(LastHighestNTested)) break;
if (LastHighestNTested == Primes.PrimesList.Last()) break;
}
Console.WriteLine("\rPrimes preloaded {0}, max {1}", Primes.PrimesList.Count, Primes.PrimesList.Last());
#endif
Ntyp c = 1;
if (0 < args.Length) { if (!(Ntyp.TryParse(args[0], out c))) c = 1; }
Ntyp CouTested = 0;
do
{
if (!Primes.isprime((Ntyp)c) //&& !Primes.isprime((Ntyp)c-1) //Adding this term, could the simple a!=c solutions
)
{
double Log10CFdivBF = 0;
for (Ntyp b = c - 1; 3 < b; b--)
{
Log10CFdivBF += Math.Log10(b + 1);
FillSortedLog10nFacUpToLog10X(Log10CFdivBF + 0.01);
int Ix = SortedLog10nFac.BinarySearch(Log10CFdivBF - 0.01); //-0.01 to ensure that we do not get above the real number due to numerical errors
if (Ix < 0) Ix = ~Ix; //First element larger than search=last element smaller
Ntyp a = 0; //Only do the real test if the estimate is close
//Limit changed to changed, way beyond estimated calculation error, based on (1e6!) where the error where less than 1:1.000000000001 (1:1e-12)
if (Math.Abs(SortedLog10nFac[Ix] - Log10CFdivBF) < 0.000000001 * Log10CFdivBF ) a = (Ntyp)Ix; // was fixed limit <0.00001 changed to factor
if (0 < a && a < b)
{
#if (UseLargeQ)
LargeQ CFdivBF = 1; for (Ntyp bb = b + 1; bb <= c; bb++) CFdivBF *= (Int64)bb;
LargeQ aF = LargeQ.Factorial((Int64)a);
if (CFdivBF.IsFactorial((Int64)a))
#else
BigInteger CFdivBF = 1; for (Ntyp bb = b + 1; bb <= c; bb++) CFdivBF *= bb;
BigInteger aF = FactorialBi(a);
if (aF == CFdivBF)
#endif
Console.WriteLine("{0}! * {1}! = {2}!" + (b + 1 == c ? " - simple solution: c=a!" : " - ******* REAL SOLUTION *******"), a, b, c);
CouTested++;
}
if (Primes.isprime((Ntyp)b)) break; // next b=p-1: a<b<p<c p divisor in c! but not in b! and certainly not in a! So no need to check further!
}
}
if (c % 10000 == 0) Console.Write("{0} {1} {2} \r",c, SortedLog10nFac.Count, CouTested);
c++;
} while (c != 0);
}
}
}