-
Notifications
You must be signed in to change notification settings - Fork 2
/
Mlucas.c
6693 lines (6233 loc) · 374 KB
/
Mlucas.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
/*******************************************************************************
* *
* (C) 1997-2021 by Ernst W. Mayer. *
* *
* This program is free software; you can redistribute it and/or modify it *
* under the terms of the GNU General Public License as published by the *
* Free Software Foundation; either version 2 of the License, or (at your *
* option) any later version. *
* *
* This program is distributed in the hope that it will be useful, but WITHOUT *
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or *
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for *
* more details. *
* *
* You should have received a copy of the GNU General Public License along *
* with this program; see the file GPL.txt. If not, you may view one at *
* http://www.fsf.org/licenses/licenses.html, or obtain one by writing to the *
* Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA *
* 02111-1307, USA. *
* *
*******************************************************************************/
/******************************************************************************************
* COMPILING AND RUNNING THE PROGRAM: see https://www.mersenneforum.org/mayer/README.html *
******************************************************************************************/
#include "Mlucas.h"
#ifndef imul_macro_h_included
#error imul_macro.h file not included in build!
#endif
// Oct 2021: Fixed non-threadsafe bug in v20-added FFT-length reversion code, but add preprocessor
// flag to allow builders to disable it at compile time if they encounter any further bugs or performance issues:
#ifndef USE_FFTLEN_REVERSION
#define USE_FFTLEN_REVERSION 1 // Builder would invoke -DUSE_FFTLEN_REVERSION=0 at compile time to override the default
#endif
/* Make sure none of the factoring-module-only flags are active: */
#if(defined(P4WORD) || defined(P3WORD) || defined(P2WORD))
#error multiword exponents only allowed for factor.c built in standalone mode!
#endif
#ifdef FACTOR_STANDALONE
#error FACTOR_STANDALONE flag only allowed for factor.c built in standalone mode!
#endif
/******************************************************************************************************/
/* Allocate storage for Globals (externs). Unless specified otherwise, these are declared in Mdata.h: */
/******************************************************************************************************/
#if INCLUDE_HWLOC
hwloc_topology_t hw_topology;
int HWLOC_AFFINITY = 0; // Is per-thread LPU-binding (affinity) supported?
#endif
// System-related globals:
uint32 SYSTEM_RAM = 0; // Total usable main memory size in MB, and max. % of that to use per instance.
#ifdef OS_TYPE_LINUX // Linux: this is based on the value of the sysinfo "freeram" field; default is 90%.
uint32 MAX_RAM_USE = 90;
#else // MacOS: I've not found a reliable way to obtain free-RAM numbers, so default is 50%
uint32 MAX_RAM_USE = 50; // of available RAM, based on the sysctl "hw.memsize" value.
#endif
// Used to force local-data-tables-reinits in cases of suspected table-data corruption:
int REINIT_LOCAL_DATA_TABLES = 0;
// Normally = True; set = False on quit-signal-received to allow desired code sections to and take appropriate action:
int MLUCAS_KEEP_RUNNING = 1;
// v18: Enable savefile-on-interrupt-signal, access to argc/argv outside main():
char **global_argv;
FILE *dbg_file = 0x0;
double*ADDR0 = 0x0; // Allows for easy debug on address-read-or-write than setting a watchpoint
// Define FFT-related globals (declared in Mdata.h):
uint32 N2,NRT,NRT_BITS,NRTM1;
int PFETCH_BLOCK_IDX[MAX_RADIX];// Need this for prefetch-block-index arrays
uint32 NRADICES, RADIX_VEC[10]; // NRADICES, RADIX_VEC[] store number & set of complex FFT radices used.
#ifdef MULTITHREAD
uint64 CORE_SET[MAX_CORES>>6]; // Bitmap for user-controlled affinity setting, as specified via the -cpu flag
#endif
int ROE_ITER = 0; // Iteration of any dangerously high ROE encountered during the current iteration interval.
uint32 NERR_ROE = 0; // v20: Add counter for dangerously high ROEs encountered during test
// This must be > 0, but make signed to allow sign-flip encoding of retry-fail.
double ROE_VAL = 0.0; // Value (must be in (0, 0.5)) of dangerously high ROE encountered during the current iteration interval
int USE_SHORT_CY_CHAIN = 0;
int ITERS_BETWEEN_CHECKPOINTS; /* number of iterations between checkpoints */
int DO_GCHECK = FALSE; // If Mersenne/PRP or Fermat/Peoin test, Toggle to TRUE at runtime
uint32 NERR_GCHECK = 0; // v20: Add counter for Gerbicz-check errors encountered during test
int ITERS_BETWEEN_GCHECK_UPDATES = 1000; // iterations between Gerbicz-checkproduct updates
int ITERS_BETWEEN_GCHECKS = 1000000; // #iterations between Gerbicz-checksum residue-integrity checks
char ESTRING[STR_MAX_LEN]; // Mersenne exponent or Fermat-number index in string form - for M(p) this == p, for F(m) this == m
char BIN_EXP[STR_MAX_LEN]; // Binary exponent in string form - for M(p) this == p, for F(m) this == 2^m
char PSTRING[STR_MAX_LEN]; // Modulus being used in string form, e.g. "M110916929" and "F33".
// The index following 'mask' here = log2(#doubles in SIMD register) = log2(#bits in SIMD register) - 6.
// The corrsponding mask masks off one my bit than this, since we are dealing with complex FFT data which
// occupy pairs of such registers:
#ifdef USE_AVX512
const uint32 mask03 = 0xfffffff0,
br16[16] = {0,8,1,9,2,10,3,11,4,12,5,13,6,14,7,15}, // length-16 index-scramble array for mapping from scalar-complex to AVX512 (8 x re,8 x im)
brinv16[16] = {0,2,4,6,8,10,12,14,1,3,5,7,9,11,13,15}; // length-16 index-unscramble array: br[brinv[i]] = brinv[br[i]] = i .
// EZ way to enumerate (i)th element: 'in which slot of br16[] is i?'
#endif
#ifdef USE_AVX // AVX and AVX2 both use 256-bit registers
const uint32 mask02 = 0xfffffff8,
br8[8] = {0,4,1,5,2,6,3,7}, // length-8 index-scramble array for mapping from scalar-complex to AVX (re,re,re,re,im,im,im,im)
brinv8[8] = {0,2,4,6,1,3,5,7}; // length-8 index-unscramble array: br[brinv[i]] = brinv[br[i]] = i .
#endif
#ifdef USE_SSE2
const uint32 mask01 = 0xfffffffc,
br4[4] = {0,2,1,3}; // length-4 index-scramble array for mapping from scalar-complex (re,im,re,im) to SSE2 (re,re,im,im)
// For length-4 this is its own inverse.
#endif
const int hex_chars[16] = {'0','1','2','3','4','5','6','7','8','9','A','B','C','D','E','F'};
char cbuf[STR_MAX_LEN],cstr[STR_MAX_LEN];
char in_line[STR_MAX_LEN];
char *char_addr;
FILE *fp = 0x0, *fq = 0x0; // Our convention is to always set these = 0x0 on fclose, so nullity correlates with "no open file attached"
/* File access mode is a 3-character string, with the last of these being a mandatory null terminator: */
char FILE_ACCESS_MODE[3] = {'x','y','\0'};
/* Matrix of supported moduli/test-types - init'ed in Mlucas_init */
int ASSIGNMENT_TYPE_MATRIX[MODULUS_TYPE_DIM][TEST_TYPE_DIM];
/* These need to be kept updated to match the #defines in Mdata.h: */
const char *err_code[ERR_MAX] = {
"ERR_INCORRECT_RES64",
"ERR_RADIX0_UNAVAILABLE",
"ERR_RADIXSET_UNAVAILABLE",
"ERR_TESTTYPE_UNSUPPORTED",
"ERR_EXPONENT_ILLEGAL",
"ERR_FFTLENGTH_ILLEGAL",
"ERR_ECHECK_NOTBOOL",
"ERR_TESTITERS_OUTOFRANGE",
"ERR_ROUNDOFF",
"ERR_CARRY",
"ERR_RUN_SELFTEST_FORLENGTH",
"ERR_ASSERT",
"ERR_UNKNOWN_FATAL",
"ERR_SKIP_RADIX_SET",
"ERR_INTERRUPT",
"ERR_GERBICZ_CHECK"
};
// Shift count and auxiliary arrays used to support rotated-residue computations:
uint64 RES_SHIFT = 0xFFFFFFFFFFFFFFFFull; // 0 is a valid value here, so init to UINT64_MAX, which value is treated as "uninited"
uint64 GCHECK_SHIFT = 0ull;
uint32 RES_SIGN = 0; // Feb 2020: uint32 to keep track of shifted-residue sign flips, needed for rotated residue Fermat-mod arithmetic.
uint64 *BIGWORD_BITMAP = 0x0;
uint32 *BIGWORD_NBITS = 0x0;
// For PRP tests, the base. For Pépin tests via the "FermatTest" worktype, the base defaults to 3; to do a
// Pépin test to another base, the more-general PRP-worktype must be specified with appropriate parameters.
uint32 PRP_BASE = 0;
uint64 *BASE_MULTIPLIER_BITS = 0x0; // Runtime-allocated bitwise multiply-by-base array
// Nov 2020: p-1 stuff:
uint64 *PM1_S1_PRODUCT = 0x0, PM1_S1_PROD_RES64 = 0ull; // Vector to hold Stage 1 prime-powers product, and (mod 2^64) checksum on same
uint32 PM1_S1_PROD_B1 = 0, PM1_S1_PROD_BITS = 0; // Stage 1 bound to which the current value of PM1_S1_PRODUCT corresponds, and its #bits
uint32 PM1_S2_NBUF = 0; // # of floating-double residue-length memblocks available for Stage 2
// Allow Stage 2 bounds to be > 2^32; B2_start defaults to B1, but can be set > B1 to allow for arbitrary Stage 2 prime intervals:
uint32 B1 = 0;
uint64 B2 = 0ull, B2_start = 0ull;
// Bit-depth of TF done on a given exponent. This is currently only used for auto-setting p-1 bounds:
uint32 TF_BITS = 0;
/* These should all be set to a valid (nonzero) value at the time the appropriate test is begun */
uint32 TEST_TYPE = 0;
uint32 MODULUS_TYPE = 0;
uint32 TRANSFORM_TYPE = 0;
const char HOMEPAGE [] = "https://www.mersenneforum.org/mayer/README.html";
/* Program version with patch suffix:
For the foreseeable future, version numbers will be of form x.yz, with x a 1-digit integer,
y an int having 3 digits or less, and z an optional alphabetic patch suffix. We choose these such that
retiming is only necessary between versions that differ in x or in the leading digit of y,
i.e. we'd need to regenerate .cfg files in going from version 2.9 to 3.0 or from 3.0 to 3.141,
but not between versions 3.0x and 3.0g or between 3.1 and 3.141.
This reduces the "retime?" decision in the cfgNeedsUpdating() function to a simple comparison
of the leading 3 characters of the two version strings in question.
*/
/* Dec 2014: Implement version numbering according to the scheme:
Major index = year - 2000
Minor index = release # of that year, zero-indexed.
A version suffix of x, y, or z following the above numeric index indicates an [alpha,beta,gamma] (experimental,unstable) code.
A third index following release # indicates a patch number relative to that release. No 3rd index can be read as "patch number 0".
*/
const char VERSION [] = "21.0.1";
const char OFILE [] = "results.txt"; /* ASCII logfile containing FINAL RESULT ONLY for each
assignment - detailed intermediate results for each assignment
are written to the exponent-specific STATFILE (see below). */
const char WORKFILE [] = "worktodo.txt"; /* File containing exponents to be tested. New exponents
may be appended at will, while the program is running. */
const char MLUCAS_INI_FILE[] = "mlucas.ini"; /* File containing user-customizable configuration settings [currently unused] */
char CONFIGFILE[15]; /* Configuration File: contains allowed FFT lengths
and allows user to control (at runtime, and in
a modifiable way) which of a predefined set
of FFT radices gets used for each runlength.
This file is re-read after each checkpoint.
File is named mlucas.cfg or fermat.cfg, depending
on whether Mersenne or Fermat number test is being done. */
/* These are set at runtime, based on the exponent being processed. */
char STATFILE [STR_MAX_LEN]; /* ASCII logfile for the current exponent */
char RESTARTFILE[STR_MAX_LEN]; /* Restart file name(s) */
uint64 KNOWN_FACTORS[40]; // Known prime-factors input to p-1 runs ... for now limit to 10 factors, each < 2^256
int INTERACT;
double AME,MME; /* Avg and Max per-iteration fractional error for a given iteration interval */
uint32 AME_ITER_START; /* Iteration # at which to start collecting RO Err data for AME & MME computation: */
// These externs declared in platform.h:
int MAX_THREADS = 0; /* Max. allowable No. of threads. */
int NTHREADS = 0; /* actual No. of threads. If multithreading disabled, set = 1. */
uint64 PMIN; /* minimum exponent allowed */
uint64 PMAX; /* maximum exponent allowed depends on max. FFT length allowed
and will be determined at runtime, via call to given_N_get_maxP(). */
/****** END(Allocate storage for Globals (externs)). ******/
#ifndef NO_USE_SIGNALS
void sig_handler(int signo)
{
if (signo == SIGINT) {
fprintf(stderr,"received SIGINT signal.\n"); sprintf(cbuf,"received SIGINT signal.\n");
} else if(signo == SIGTERM) {
fprintf(stderr,"received SIGTERM signal.\n"); sprintf(cbuf,"received SIGTERM signal.\n");
#ifndef __MINGW32__
} else if(signo == SIGHUP) {
fprintf(stderr,"received SIGHUP signal.\n"); sprintf(cbuf,"received SIGHUP signal.\n");
} else if(signo == SIGALRM) {
fprintf(stderr,"received SIGALRM signal.\n"); sprintf(cbuf,"received SIGALRM signal.\n");
} else if(signo == SIGUSR1) {
fprintf(stderr,"received SIGUSR1 signal.\n"); sprintf(cbuf,"received SIGUSR1 signal.\n");
} else if(signo == SIGUSR2) {
fprintf(stderr,"received SIGUSR2 signal.\n"); sprintf(cbuf,"received SIGUSR2 signal.\n");
#endif
}
// Dec 2021: Until resolve run-to-run inconsistencies in signal handling, kill it with fire:
exit(1);
// Toggle a global to allow desired code sections to detect signal-received and take appropriate action:
MLUCAS_KEEP_RUNNING = 0;
}
#endif
/****************/
/*
!...Code to test primality of Mersenne numbers, using arbitrary-precision (array-integer) arithmetic.
! Author: Ernst W. Mayer.
!
! Accomplishes the Mersenne-mod squaring via the weighted discrete Fourier transform technique
! of Crandall and Fagin (Mathematics of Computation 62 (205), pp.305-324, January 1994; also
! available online at http://www.faginfamily.net/barry/Papers/Discrete%20Weighted%20Transforms.pdf ).
!
! For each binary exponent p, generates first p-2 terms of the Lucas-Lehmer sequence,
!
! s(n+1) = s(n)**2 - 2, s(0) = 4
!
! modulo N = 2^p - 1. If N divides s(p-2) (specifically, the (p-2)th residue == 0 modulo N), N is prime.
!
! See https://www.mersenneforum.org/mayer/README.html for build instructions, recent revision history
! of the code and a summary of what has been changed/added/deleted in the latest release.
!
!***TO DO LIST:
!
! Performance-related:
! (0) Continue to improve cache/TLB behavior of code. This should include the following enhancements at a minimum:
!
! - Extend block-FFT strategy to include more than just the initial FFT pass, thus allowing the blocks
! to be smaller than n/radix(pass 1). This should help at very large runlengths and on systems with
! small L1 caches.
!
! (1) Continue to optimize SIMD assembler; add support for next-gen AVX512 vector instructions.
!
! (2) Implement hybrid complex floating-point FFT with a modular complex (i.e. Gaussian integer) transform
! over a suitable Galois field GF(Mp^2), where Mp = 2^31 - 1. Based on vector integer-math capabilities of
! current bleeding-edge x86 processors, p = 31 is the prime candidate (pun intended).
! For background on the mathematics, see Richard Crandall's preprint, "Integer convolution via split-radix fast
! Galois transform", freely available at: http://academic.reed.edu/physics/faculty/crandall/papers/confgt.pdf .
!
! Functionality-related:
! (*) P-1 factoring module (together with hand-rolled subquadratic, memory-efficient GCD)
! (*) An (optional) GUI would be nice...
!
!
!***ACKNOWLEDGEMENTS: special thanks are due to the following people:
!
! * Richard Crandall and Barry Fagin - for the discrete weighted transform.
!
! * Peter Montgomery - for help with theoretical and implementation questions,
! loan of his fast factoring code, and for testing early versions of the program on his SGI
! (and uncovering several compiler bugs along the way).
!
! * Luke Welsh - for historical perspective, und besonders fu"r das Max und Moritz Weissbier Glas.
!
! * George Woltman - for organizing the Great Internet Mersenne Prime Search,
! and for freely sharing his time and keen computational insights.
!
! * Scott Kurowski - for the GIMPS PrimeNet server.
!
! * Jason Papadopoulos - for his valuable perspectives regarding FFT machine
! implementations and algorithmic optimization for various machine architectures.
!
! * Guillermo Ballester Valor - for many interesting discussions on optimization
! and useful suggestions for improving the code.
!
! * John Pierce - for hosting the ftp archive on his hogranch.com server for many years.
!
! * David Stanfill - for generous access to high-end Xeon server machinery and for hosting the "GIMPS KNL".
!
! ...as well as the numerous people who have been kind enough to try the code out
! and send timings, compiler options for various architectures and suggestions
! for improvement.
!
!***For UPDATES, PATCHES and USAGE INSTRUCTIONS, see
!
! https://www.mersenneforum.org/mayer/README.html
!
!***Send QUESTIONS, COMMENTS or SUGGESTIONS to me at <[email protected]>.
!
! Alright, let's go hunt for some primes!
!
*/
// Mod-2^64 sum of elements of uint64 array a[] ... in the context of the G-checkproduct integrity
// checking we send it a double-float array treated as a uint64 one via type-punning cast of a[]:
uint64 sum64(uint64 a[], uint32 n) {
uint32 i;
uint64 sum = 0ull;
for(i = 0; i < n; i++)
sum += a[i];
return sum;
}
// Simple majority-vote consensus:
uint64 consensus_checksum(uint64 s1, uint64 s2, uint64 s3) {
if(s1 == s2) return s1;
if(s1 == s3) return s1;
if(s2 == s3) return s2;
return 0;
}
/* The return value here is considered as consisting of bytewise subfields. The lowest byte contains
one of the error codes declared in Mdata.h (or 0 if successful exit); the upper bytes should only be nonzero
for certain specific value of the lowest-order byte, as indicated in Mdata.h:
*/
uint32 ernstMain
(
int mod_type,
int test_type,
uint64 exponent,
uint32 fft_length,
int radix_set,
uint32 maxFFT,
uint32 iterations, /* Use to store log2[max factor depth] in TF mode; ignore in p-1 mode */
uint64 *sh0, /* Reference/Computed mod-(2^64) residue */
uint64 *sh1, /* Reference/Computed Selfridge/Hurwitz Residue mod 2^35-1 */
uint64 *sh2, /* Reference/Computed Selfridge/Hurwitz Residue mod 2^36-1 */
int scrnFlag,
double *runtime
)
{
int resFlag = 0;
/*...Various general run parameters.. */
uint32 timing_test_iters = 0;
/*...scalars and fixed-size arrays... */
uint32 i,j,k = 0;
/* TODO: some of these need to become 64-bit: */
uint32 dum = 0,findex = 0,ierr = 0,ilo = 0,ihi = 0,iseed,isprime,kblocks = 0,maxiter = 0,n = 0,npad = 0;
uint64 itmp64,cy, s1 = 0ull,s2 = 0ull,s3 = 0ull; // s1,2,3: Triply-redundant whole-array checksum on b,c-arrays used in the G-check
uint32 mode_flag = 0, first_sub, last_sub;
/* Exponent of number to be tested - note that for trial-factoring, we represent p
strictly in string[STR_MAX_LEN] form in this module, only converting it to numeric
form in the factoring module. For all other types of assignments uint64 should suffice: */
uint64 p = 0, i1,i2,i3, rmodb,mmodb;
uint32 nbits_in_p = 0, nfac;
/* Res64 and Selfridge-Hurwitz residues: */
uint64 Res64, Res35m1, Res36m1;
/*...Known Mersenne prime exponents. This array must be null-terminated. */
// Dec 2018: Including M51, there are (31, 19) p = 1,3 (mod 4), resp., vs (25.8, 24.5) predicted
// (for p < 10^8) by the Lenstra/Wagstaff heuristic (cf. est_num_mp_in_interval() in util.c):
const uint32 knowns[] = {2,3,5,7,13,17,19,31,61,89,107,127,521,607,1279,2203,2281,3217,4253,4423,9689,9941,11213,19937,21701 // M#1-25
,23209,44497,86243,110503,132049,216091,756839,859433,1257787,1398269,2976221,3021377,6972593,13466917,20996011 // M#26-40
,24036583,25964951,30402457,32582657,37156667,42643801,43112609,57885161,74207281,77232917,82589933,0x0}; // M#41-51
/*...What a bunch of characters... */
char *cptr = 0x0, *endp, gcd_str[STR_MAX_LEN], aid[33] = "\0"; // 32-hexit Primenet assignment id needs 33rd char for \0
/*...initialize logicals and factoring parameters... */
int restart = FALSE, use_lowmem = 0, check_interval = 0;
#if INCLUDE_TF
uint32 bit_depth_todo = 0;
uint64 factor_k_start = 0;
uint32 factor_pass_start = 0, factor_pass_hi = 0;
double log2_min_factor = 0, log2_max_factor = 0;
#endif
double tests_saved = 0.0; // v21: make this a dfloat to allow fractional-parts
uint32 pm1_done = FALSE, split_curr_assignment = FALSE, s2_continuation = FALSE, s2_partial = FALSE;
uint32 pm1_bigstep = 0, pm1_stage2_mem_multiple = 0, psmall = 0;
/*...allocatable data arrays and associated params: */
static uint64 nbytes = 0, nalloc = 0, arrtmp_alloc = 0, s1p_alloc = 0;
static double *a_ptmp = 0x0, *a = 0x0, *b = 0x0, *c = 0x0, *d = 0x0, *e = 0x0;
// uint64 scratch array and 4 pointers used to store cast-to-(uint64*) version of above b,c,d,e-pointers
static uint64 *arrtmp = 0x0, *b_uint64_ptr = 0x0, *c_uint64_ptr = 0x0, *d_uint64_ptr = 0x0, *e_uint64_ptr = 0x0;
double final_res_offset;
/*...time-related stuff. clock_t is typically an int (signed 32-bit)
and a typical value of CLOCKS_PER_SEC (e.g. as defined in <machine/machtime.h>
on Alpha TruUnix) is 1000000, so we should accumulate clock_t-stored time
differences at least roughly every half hour (preferably every second or so)
to avoid weirdness due to flipping of the sign bit or integer overflow.
*/
/*clock_t clock1, clock2; Moved these to [mers|fermat]_mod_square.c */
double tdiff,tdif2;
#define SIZE 256
time_t calendar_time;
struct tm *local_time, *gm_time;
char timebuffer[SIZE];
/*...entry point for one or more Lucas-Lehmer tests is here. */
MODULUS_TYPE = mod_type;
TEST_TYPE = test_type;
INTERACT = FALSE;
RANGE_BEG:
p = 0ull; ierr = 0;
USE_SHORT_CY_CHAIN = 0; // v19: Reset carry-chain length fiddler to default (faster/lower-accuracy) at start of each run:
ROE_ITER = 0; ROE_VAL = 0.0;
NERR_GCHECK = NERR_ROE = 0; // v20: Add counters for Gerbicz-check errors and dangerously high ROEs encountered
// during test - if a restart, will re-read actual cumulative values from checkpoint file.
// Clear out any FFT-radix or known-factor data that might remain from a just-completed run:
for(i = 0; i < 10; i++) { RADIX_VEC[i] = 0; }
nfac = 0; mi64_clear(KNOWN_FACTORS,40);
NRADICES = 0;
RESTARTFILE[0] = STATFILE[0] = '\0';
restart = FALSE;
B1 = 0; B2 = B2_start = 0ull; gcd_str[0] = '\0'; split_curr_assignment = s2_continuation = s2_partial = FALSE;
pm1_bigstep = pm1_stage2_mem_multiple = psmall = 0;
// Check for user-set value of various flags. Failure-to-find-or-parse results in isNaN(dtmp) = TRUE, print nothing in that case:
double dtmp = mlucas_getOptVal(MLUCAS_INI_FILE,"LowMem");
if(dtmp != 0) {
if(dtmp != dtmp) { // isNaN is C99, want something that also works on pre-C99 platforms
sprintf(cbuf,"User did not set LowMem in %s ... allowing all test types.\n",MLUCAS_INI_FILE);
} else if(dtmp == 1) {
sprintf(cbuf,"User set LowMem = 1 in %s ... this allows PRP-testing but excludes p-1 stage 2.\n",MLUCAS_INI_FILE); use_lowmem = 1;
} else if(dtmp == 2) {
sprintf(cbuf,"User set LowMem = 2 in %s ... this excludes both PRP-testing and p-1 stage 2.\n",MLUCAS_INI_FILE); use_lowmem = 2;
} else {
sprintf(cbuf,"User set unsupported value LowMem = %f in %s ... ignoring.\n",dtmp,MLUCAS_INI_FILE);
}
mlucas_fprint(cbuf,1);
}
dtmp = mlucas_getOptVal(MLUCAS_INI_FILE,"CheckInterval");
if(dtmp != 0) {
if(dtmp != dtmp) {
sprintf(cbuf,"User did not set CheckInterval in %s ... using default.\n",MLUCAS_INI_FILE);
} else if(dtmp < 1000 || dtmp > 1000000) {
sprintf(cbuf,"User set CheckInterval = %f in %s ... values < 10^3 or > 10^6 are not supported, ignoring.\n",dtmp,MLUCAS_INI_FILE);
} else if(DNINT(dtmp) != dtmp) {
sprintf(cbuf,"User set non-whole-number CheckInterval = %f in %s ... ignoring.\n",dtmp,MLUCAS_INI_FILE);
} else {
sprintf(cbuf,"User set CheckInterval = %d in %s.\n",(int)dtmp,MLUCAS_INI_FILE); check_interval = (int)dtmp;
}
mlucas_fprint(cbuf,1);
}
/* ...If multithreading enabled, set max. # of threads based on # of available (logical) processors,
with the default #threads = 1 and affinity set to logical core 0, unless user overrides those via -nthread or -cpu:
*/
#ifdef MULTITHREAD
#ifdef USE_OMP
// OpenMP not currently supported (attempting to build with this #define enabled barfs in
// preprocessing via #error in platform.h), this is merely placeholder for possible future use:
ASSERT(HERE, MAX_THREADS = omp_get_num_procs(), "Illegal #Cores value stored in MAX_THREADS");
#elif(defined(USE_PTHREAD))
ASSERT(HERE, MAX_THREADS = get_num_cores(), "Illegal #Cores value stored in MAX_THREADS");
#else
#error Unrecognized multithreading model!
#endif
// MAX_THREADS based on number of processing cores will most often be a power of 2, but don't assume that.
ASSERT(HERE, MAX_THREADS > 0,"MAX_THREADS must be > 0");
ASSERT(HERE, MAX_THREADS <= MAX_CORES,"MAX_THREADS exceeds the MAX_CORES setting in Mdata.h .");
if(!NTHREADS) {
NTHREADS = 1;
fprintf(stderr,"No CPU set or threadcount specified ... running single-threaded.\n");
// Use the same affinity-setting code here as for the -cpu option, but simply for cores [0:NTHREADS-1]:
sprintf(cbuf,"0:%d",NTHREADS-1);
parseAffinityString(cbuf);
} else if(NTHREADS > MAX_CORES) {
sprintf(cbuf,"ERROR: NTHREADS = %d exceeds the MAX_CORES setting in Mdata.h = %d\n", NTHREADS, MAX_CORES);
ASSERT(HERE, 0, cbuf);
} else { // In timing-test mode, allow #threads > #cores
if(NTHREADS > MAX_THREADS) {
fprintf(stderr,"WARN: NTHREADS = %d exceeds number of cores = %d\n", NTHREADS, MAX_THREADS);
}
fprintf(stderr,"NTHREADS = %d\n", NTHREADS);
}
#else
MAX_THREADS = NTHREADS = 1;
#endif // #ifdef MULTITHREAD ?
/* Make number of iterations between checkpoints dependent on #threads -
don't want excessively frequent savefile writes, at most 1 or 2 an hour is needed:
*/
// Oct 2021: Allow this to be set via CheckInterval option in mlucas.ini and read at program-start time,
// albeit for PRP-tests subject to constraints related to Gerbicz checking:
if(!check_interval) {
if(NTHREADS > 4)
ITERS_BETWEEN_CHECKPOINTS = 100000;
else
ITERS_BETWEEN_CHECKPOINTS = 10000;
} else if(check_interval < 1000) {
ASSERT(HERE,0,"User-set value of check_interval must >= 1000.");
} else
ITERS_BETWEEN_CHECKPOINTS = check_interval;
fprintf(stderr,"Setting ITERS_BETWEEN_CHECKPOINTS = %u.\n",ITERS_BETWEEN_CHECKPOINTS);
i = ITERS_BETWEEN_GCHECKS;
j = ITERS_BETWEEN_GCHECK_UPDATES;
ASSERT(HERE, i == j*j, "#iterations between Gerbicz-checksum updates must = sqrt(#iterations between residue-integrity checks)");
// v19: If PRP test, make sure Gerbicz-checkproduct interval divides checkpoint-writing one.
// If not true, merely warn here because user may be doing LL/DC/p-1 and not PRP-tests:
k = ITERS_BETWEEN_CHECKPOINTS;
if(i%k != 0 || k%j != 0)
fprintf(stderr,"WARN: G-checkproduct update interval must divide savefile-update one, which must divide the G-check interval ... this will trigger an assertion-exit if a PRP-test is attempted.");
// Alloc bitwise multiply-by-base array, needed to support P-1 factoring and PRP testing:
if(!BASE_MULTIPLIER_BITS) {
j = ((ITERS_BETWEEN_CHECKPOINTS+63) >> 6) + 1; // Add 1 pad element in case compiler does not 64-bit align
BASE_MULTIPLIER_BITS = ALLOC_UINT64(BASE_MULTIPLIER_BITS, j); if(!BASE_MULTIPLIER_BITS){ sprintf(cbuf, "ERROR: unable to allocate BASE_MULTIPLIER_BITS array in main.\n"); fprintf(stderr,"%s", cbuf); ASSERT(HERE, 0,cbuf); }
BASE_MULTIPLIER_BITS = ALIGN_UINT64(BASE_MULTIPLIER_BITS); ASSERT(HERE, ((intptr_t)BASE_MULTIPLIER_BITS & 63) == 0x0,"BASE_MULTIPLIER_BITS[] not aligned on 64-byte boundary!");
for(i = 0; i < j; i++) { BASE_MULTIPLIER_BITS[i] = 0ull; } // v20: Init = 0 here, in case we jump directly into p-1 stage 2 on restart
}
/* Look for work file... */
fp = 0x0;
if (!exponent || (exponent!=0 && fft_length!=0 && iterations==0)) {
fprintf(stderr," looking for %s file...\n",WORKFILE);
fp = mlucas_fopen(WORKFILE, "r");
}
/***********************************************************************************/
/* Automated exponent-dispatch mode - Note that in the assignment-syntax comments, */
/* <> and [] mean required and optional arguments, respectively: */
/***********************************************************************************/
if (fp) {
fprintf(stderr," %s file found...reading next assignment...\n",WORKFILE);
read_next_assignment: // Read first line of worktodo.ini file into 1K character array:
if(!fgets(in_line, STR_MAX_LEN, fp)) {
fprintf(stderr,"Hit EOF while attempting to read next line of worktodo ... quitting.\n");
exit(0);
}
fprintf(stderr," %s entry: %s\n",WORKFILE,in_line);
/* Skip any whitespace at beginning of the line: */
char_addr = in_line;
while(isspace(*char_addr)) {
++char_addr;
}
// v20.1.1: Parse all lines whose 1st non-WS char is alphabetic; print "Ignoring (copy of workfile line)" for all entries not so.
// NB: Discontinue support for numeric leading char, i.e. Mersenne-exponent-only legacy format:
if(!isalpha(*char_addr)) {
fprintf(stderr," Leading non-WS char of %s entry is not alphabetic ... skipping to next entry.\n",WORKFILE);
goto read_next_assignment;
}
// Otherwise assume Prime95-style ini file format, with a possible modulus-specific leading keyword;
// Default "Test=" means Mersenne, unless "Test" preceded by an explicit modulus-type string:
MODULUS_TYPE = MODULUS_TYPE_MERSENNE;
/* Re. the recently-added-to-Primenet PRP assignment type, On Dec 19, 2017, at 5:07 PM, George Woltman wrote:
In "PRP=[aid],1,2,75869377,-1,75,0,3,4" ([aid] stands for an optional 32-hexit assignment ID)
The first four (numeric) values are k,b,n,c as in modulus = k*b^n + c
75 is how far factored
0 is the number of PRP tests that will be saved if P-1 is done and finds a factor
3 is the PRP base (used for the first-time test)
4 is the residue type (used for the first-time test)
You should only need to support residue types 1 and 5, which are identical if there are no cofactors
(==> in JSON-result, return residue type 1).
Example: PRP=C42540C352E54E906108D48FA5D89488,1,2,80340397,-1,75,0,3,1
means a PRP test of 1.2^80340397-1 = M80340397, TFed to 2^75 (ignore), 0 saved tests (ignore), PRP base 3 and type 1.
EWM: The above are PRP-DCs ... first-time tests lack the last 2 numeric args above,
and default to base = 3, residue type 1. Further, there is also a PRPDC= format for such.
PRP-CF (Mersenne-cofactor-PRP tests) have same format as PRP, but append one or more known factors:
PRP=[aid],1,2,3215747,-1,99,0,3,1,"4457025343,185822885311153245017"
[EWM: F25 cofactor-test:
PRP=n/a,1,2,33554432,+1,99,0,3,5,"2170072644496392193"
]
Ref. for residue type definitions: https://www.mersenneforum.org/showthread.php?p=468378#post468378
George: "There are (at least) 5 PRP residue types:
EWM: Needed LR-binary-modpow modmul sequence, N = M(p) = 2^p-1 = [p binary ones],
initial seed x = a, and example final result for a = 3 and M(23) (not prime), M(31) (prime):
p = 23 p = 31 bc code:
1: Fermat PRP. N is PRP if a^(N-1) = 1 mod N p-1 (x = a.x^2), followed by one (x = x^2) 5884965 +1 m = 2^p-1; x = a = 3; for(i = 0; i < (p-2); i++) { x = a*x^2 % m; }; x = x^2 % m; x
2: SPRP variant, N is PRP if a^((N-1)/2) = +/-1 mod N p-1 (x = a.x^2) 7511964 -1 m = 2^p-1; x = a = 3; for(i = 0; i < (p-2); i++) { x = a*x^2 % m; }; x
3: Fermat PRP variant. N is PRP if a^(N+1) = a^2 mod N p mod-squarings (x = x^2) 2633043 9 m = 2^p-1; x = a = 3; for(i = 0; i < p; i++) { x = x^2 % m; }; x
4: SPRP variant. N is PRP if a^((N+1)/2) = +/-a mod N p-1 mod-squarings (x = x^2) 5758678 -3 m = 2^p-1; x = a = 3; for(i = 0; i < (p-1); i++) { x = x^2 % m; }; x
5: Fermat PRP cofactor variant, N/d is PRP if a^(N-1) = a^d mod N/d
I encourage programs to return type 1 residues as that has been the standard for prime95, PFGW, LLR for many years."
EWM: Note for PRP w/Gerbicz-check we need a mod-squaring-only sequence, so we *compute* type 3, but then do a
mod-div-by-scalar to get a^(N-1) mod N from a^2 mod N and *report* the resulting type 1 residue to the server.
Example:
Compute type 3 (Fermat-PRP) residue for N = M(23) using p = 23 mod-squarings of initial seed a = 3, yielding R3 = 2633043.
Now we need to mod-divide by a^2 = 9 to get R1, which amounts to finding R1 such that a^2*R1 == R3 (mod N). 2 options here:
1: Compute the inverse of a^2 (mod N) = 1864135, R1 = (a^-2)*R3 (mod N) = 5884965 .
2: Find a small multiplier k such that R3 + k*N is divisible by a^2 = 9. Using trial and error starting from k = 0:
k R3 + k*N (mod 9)
-- -------
0 3
1 7
2 2
3 6
4 1
5 5
6 0, so R1 = (R3 + 6*N)/9 = 5884965 .
More cleverly, we compute R3 == 3 (mod 9) and N == 4 (mod 9), so we want to solve 3 + 4*k == 0 (mod 9) for the index k.
That can be done by simply rearranging to 4*k == -3 (mod 9), then finding inverse of 4 (mod 9) and multiplying both sides
by that (mod 9) to get k. Note the difference between this mod-inverse computation and the one in option [1]: This one is
not the mod-inverse w.r.to N, it's one w.r.to a^2, which is tiny compared to N.
Calling that with args x = 4, n = 9 gives our k = -3*modinv(4,9) = -3*-2 = 6, same as the above trial-and-error approach.
*/
if((char_addr = strstr(in_line, "PRP")) != 0) // This also handles the PRPDC= format
{
TEST_TYPE = TEST_TYPE_PRP;
char_addr += 3;
// Check [k,b,n,c] portion of in_line:
cptr = check_kbnc(char_addr, &p);
ASSERT(HERE, cptr != 0x0, "[k,b,n,c] portion of in_line fails to parse correctly!");
// Next 2 entries in in_line are how-far-factored and "# of PRP tests that will be saved if P-1 is done and finds a factor":
TF_BITS = 0xffffffff; tests_saved = 0.0;
if((char_addr = strstr(cptr, ",")) != 0x0) {
cptr++;
// Only check if there's an appropriate TF_BITS entry in the input line
TF_BITS = strtoul(++char_addr, &endp, 10);
ASSERT(HERE, (char_addr = strstr(cptr, ",")) != 0x0,"Expected ',' not found after TF_BITS field in assignment-specifying line!"); cptr++;
tests_saved = strtod(++char_addr, &endp);
if(tests_saved < 0 || tests_saved > 2) {
sprintf(cbuf, "ERROR: the specified tests_saved field [%10.5f] should be in the range [0,2]!\n",tests_saved); ASSERT(HERE,0,cbuf);
}
// char_addr now points to leftmost char of tests_saved field, which we will overwrite with 0;
// endp points to to-be-appended leftover portion
}
pm1_done = (tests_saved == 0);
// If there is still factoring remaining to be done, modify the assignment type appropriately.
if(pm1_done) { // pm1_done == TRUE is more or less a no-op, translating to "proceed with primality test"
cptr = char_addr; // ...but we do need to advance cptr past the ,TF_BITS,tests_saved char-block
} else {
// Create p-1 assignment, then edit original assignment line appropriately
TEST_TYPE = TEST_TYPE_PM1;
kblocks = get_default_fft_length(p);
ASSERT(HERE, pm1_set_bounds(p, kblocks<<10, TF_BITS, tests_saved), "Failed to set p-1 bounds!");
// Format the p-1 assignment into cbuf - use cptr here, as need to preserve value of char_addr:
cptr = strstr(in_line, "="); ASSERT(HERE,cptr != 0x0,"Malformed assignment!");
cptr++; while(isspace(*cptr)) { ++cptr; } // Skip any whitespace following the equals sign
if(is_hex_string(cptr, 32)) {
strncpy(aid,cptr,32); sprintf(cbuf,"Pminus1=%s,1,2,%llu,-1,%u,%llu\n",aid,p,B1,B2); // If we get here, it's a M(p), not F(m)
} else
sprintf(cbuf,"Pminus1=1,2,%llu,-1,%u,%llu\n",p,B1,B2);
// Copy up to the final (tests_saved) char of the assignment into cstr and append tests_saved = 0;
// A properly formatted tests_saved field is 1 char wide and begins at the current value of char_addr:
i = char_addr - in_line; strncpy(cstr,in_line, i); cstr[i] = '0'; cstr[i+1] = '\0';
// Append the rest of the original assignment. If original lacked a linefeed, add one to the edited copy:
strcat(cstr,endp);
if(cstr[strlen(cstr)-1] != '\n') {
strcat(cstr,"\n");
}
split_curr_assignment = TRUE; // This will trigger the corresponding code following the goto:
goto GET_NEXT_ASSIGNMENT;
} // First-time PRP test ... !cptr check is for assignments ending with [k,b,n,c] like "PRP=1,2,93018301,-1":
if(!cptr || (char_addr = strstr(cptr, ",")) == 0x0) {
PRP_BASE = 3;
TEST_TYPE = TEST_TYPE_PRP;
} else { // PRP double-check:
// NB: Hit a gcc compiler bug (which left i = 0 for e.g. char_addr = ", 3 ,...") using -O0 here ... clang compiled correctly, as did gcc -O1:
i = (int)strtol(char_addr+1, &cptr, 10); ASSERT(HERE, i == 3,"PRP-test base must be 3!");
PRP_BASE = i;
ASSERT(HERE, (char_addr = strstr(cptr, ",")) != 0x0,"Expected ',' not found in assignment-specifying line!");
i = (int)strtol(char_addr+1, &cptr, 10); ASSERT(HERE, i == 1 || i == 5,"Only PRP-tests of type 1 (PRP-only) and type 5 (PRP and subsequent cofactor-PRP check) supported!");
// Read in known prime-factors, if any supplied - resulting factors end up in KNOWN_FACTORS[]:
if(*cptr == ',') //vv--- Pass in unused file-ptr fq here in case function emits any messages:
nfac = extract_known_factors(p,cptr+1);
// Use 0-or-not-ness of KNOWN_FACTORS[0] to differentiate between PRP-only and PRP-CF:
if(KNOWN_FACTORS[0] != 0ull) {
ASSERT(HERE, i == 5,"Only PRP-CF tests of type 5 supported!");
}
}
goto GET_EXPO;
}
else if((char_addr = strstr(in_line, "Fermat")) != 0)
{
char_addr += 6;
/* Look for comma following the modulus keyword and position next-keyword search right after it: */
if(!STREQN(char_addr,",",1))
ASSERT(HERE, 0,"Expected ',' not found in input following modulus type specifier!");
else
char_addr++;
MODULUS_TYPE = MODULUS_TYPE_FERMAT;
PRP_BASE = 2; // v20: Pépin test doesn't use this as the initial seed (that defaults to 3), but rather for the random-shift
// offsets used to prevent the shift count from modding to 0 as a result of repeated doublings (mod 2^m)
}
/* "Mersenne" is the default and hence not required, but allow it: */
else if((char_addr = strstr(in_line, "Mersenne")) != 0)
{
char_addr += 8;
/* Look for comma following the modulus keyword and position next-keyword search right after it: */
if(!STREQN(char_addr,",",1))
ASSERT(HERE, 0,"Expected ',' not found in input following modulus type specifier!");
else
char_addr++;
}
// Pépin tests are assigned via "Fermat,Test=<Fermat number index>", so catch this clause by starting new if/else()
if((char_addr = strstr(in_line, "Test")) != 0)
{
TEST_TYPE = TEST_TYPE_PRIMALITY;
char_addr += 4;
}
else if((char_addr = strstr(in_line, "DoubleCheck")) != 0)
{
TEST_TYPE = TEST_TYPE_PRIMALITY;
char_addr += 11;
}
#if INCLUDE_TF
else if((char_addr = strstr(in_line, "Factor")) != 0)
{
TEST_TYPE = TEST_TYPE_TF;
char_addr += 6;
}
#endif
/* 11/23/2020: George W. re. p-1 assignment formats:
"There is no documentation on the worktodo lines (except the source code):
Pminus1=[aid,]k,b,n,c,B1,B2[,TF_BITS][,B2_start][,known_factors] (,-separated list of known factors bookended with "")
Pfactor=[aid,]k,b,n,c,TF_BITS,ll_tests_saved_if_factor_found
(***EWM: Mlucas v20 converts PRP-with-(p-1)-needed assignments to paired Pminus1|PRP
assignments to support the fused (p-1)-S1|First-part-of-PRP-test algorithm. ***)
A tests_saved value of 0.0 will bypass any P-1 factoring. The PRP residue type is defined in primenet.h .
AFAIK, the server issues Pfactor= lines not Pminus1= lines."
*/
// Use case-insensitive analog of strstr, stristr():
else if((char_addr = stristr(in_line, "pminus1")) != 0)
{
TEST_TYPE = TEST_TYPE_PM1;
char_addr += 7;
// Check [k,b,n,c] portion of in_line:
cptr = check_kbnc(char_addr, &p);
ASSERT(HERE, cptr != 0x0, "[k,b,n,c] portion of in_line fails to parse correctly!");
ASSERT(HERE, (char_addr = strstr(cptr, ",")) != 0x0 ,"Expected ',' not found in assignment-specifying line!");
B1 = (uint32)strtoul (char_addr+1, &cptr, 10);
ASSERT(HERE, (char_addr = strstr(cptr, ",")) != 0x0 ,"Expected ',' not found in assignment-specifying line!");
/* The C11 standard re. strtoull: "On success the function returns the converted integer as unsigned long long int type
and sets endPtr to point to the first character after the input number. On failure it returns 0 and sets endPtr to
point to NULL. It handles integer overflows efficiently and return ULONG_LONG_MAX on overflow."
However, in gdb-under-Mac testing of a decimal-digit input > 2^64, I found it returned ULONG_LONG_MAX, but
also set endPtr to point to the first character after the input, which leaves some ambiguity - what if the
input was in fact == ULONG_LONG_MAX? We assume here that nobody will use a p-1 stage bound so large:
*/
B2 = (uint64)strtoull(char_addr+1, &cptr, 10); ASSERT(HERE, B2 != -1ull, "strtoull() overflow detected.");
// Remaining args optional, with the 2 numerics presumed in-order, e.g. we only look for ',B2_start' field if ',TF_BITS' was present:
if((char_addr = strstr(cptr, ",")) != 0x0) {
TF_BITS = (int)strtoul(char_addr+1, &cptr, 10); ASSERT(HERE, TF_BITS < 100 ,"TF_BITS value read from assignment is out of range.");
if((char_addr = strstr(cptr, ",")) != 0x0) {
B2_start = (uint64)strtoull(char_addr+1, &cptr, 10); ASSERT(HERE, B2_start != -1ull, "strtoull() overflow detected.");
if(B2_start > B1) // It's a stage 2 continuation run
s2_continuation = TRUE;
// Read in known prime-factors, if any supplied - resulting factors end up in KNOWN_FACTORS[]:
if(*cptr == ',') nfac = extract_known_factors(p,cptr+1);
} else if((char_addr = strstr(cptr, "\"")) != 0x0) { // Known-factors list need not be preceded by TF_BITS or B2_start
nfac = extract_known_factors(p,cptr); // cptr, not cptr+1 here, since need to preserve leading " bracketing factors-list
}
}
}
else if((char_addr = stristr(in_line, "pfactor")) != 0) // Caseless substring-match as with pminus 1
{
TEST_TYPE = TEST_TYPE_PM1;
PRP_BASE = 3;
char_addr += 7;
// Check [k,b,n,c] portion of in_line:
cptr = check_kbnc(char_addr, &p);
ASSERT(HERE, cptr != 0x0, "[k,b,n,c] portion of in_line fails to parse correctly!");
ASSERT(HERE, (char_addr = strstr(cptr, ",")) != 0x0 ,"Expected ',' not found in assignment-specifying line!");
TF_BITS = (int)strtoul(char_addr+1, &cptr, 10);
ASSERT(HERE, (char_addr = strstr(cptr, ",")) != 0x0 ,"Expected ',' not found in assignment-specifying line!");
tests_saved = strtod(++char_addr, &endp);
if(tests_saved < 0 || tests_saved > 2) {
sprintf(cbuf, "ERROR: the specified tests_saved field [%10.5f] should be in the range [0,2]!\n",tests_saved); ASSERT(HERE,0,cbuf);
}
ASSERT(HERE, pm1_set_bounds(p, get_default_fft_length(p)<<10, TF_BITS, tests_saved), "Failed to set p-1 bounds!");
}
#if INCLUDE_ECM
else if(strstr(char_addr, "ECM"))
{
TEST_TYPE = TEST_TYPE_ECM;
char_addr += 3;
}
#endif
else
{
snprintf_nowarn(cbuf,STR_MAX_LEN,"WARN: Unrecognized/Unsupported option or empty assignment line. The ini file entry was %s\n",in_line);
fprintf(stderr,"%s",cbuf);
goto read_next_assignment;
}
if(!p) { // For legacy assignment types, set p here
ASSERT(HERE, (char_addr = strstr(char_addr, "=")) != 0x0,"Expected '=' not found in assignment-specifying line!");
char_addr++;
/* Skip any whitespace following the equals sign:*/
while(isspace(*char_addr)) { ++char_addr; }
/* Check for a 32-hex-digit PrimeNet v5 assignment ID preceding the exponent: */
if(is_hex_string(char_addr, 32))
char_addr += 33;
else if(STREQN_NOCASE(char_addr,"n/a",3))
char_addr = strstr(char_addr, ",") + 1;
p = strtoull(char_addr, &cptr, 10); ASSERT(HERE, p != -1ull, "strtoull() overflow detected.");
}
GET_EXPO:
// Need to init this for savefile-naming code
ASSERT(HERE, p != 0ull, "Exponent has not been set!");
sprintf(ESTRING,"%llu",p);
// In PRP-test case, have already read the exponent from the worktodo line
/* Special case of user forcing a non-default FFT length for an exponent in the worktodo file: */
if(exponent && (p != exponent)) { // || (MODULUS_TYPE != MODULUS_TYPE_MERSENNE)) 15. Oct 2012: Need same flexibility for Fermat numbers (e.g. F27 @ 7168k) as for Mersennes, so disable modulus-type part of conditional
sprintf(cbuf,"User-supplied exponent and FFT-length for full-length test requires an exponent-matching 'Test=<exponent>' or 'DoubleCheck=<exponent>' %s entry!",WORKFILE);
ASSERT(HERE, 0,cbuf);
}
/* Check #bits in the Mersenne exponent vs. the allowed maximum: */
if(MODULUS_TYPE == MODULUS_TYPE_MERSENNE) {
nbits_in_p = 64 - leadz64(p);
}
/* If it's a Fermat number, need to check size of 2^ESTRING: */
else if(MODULUS_TYPE == MODULUS_TYPE_FERMAT) {
findex = (uint32)p;
if(findex <= MAX_PRIMALITY_TEST_BITS)
p = (uint64)1 << findex;
else
ASSERT(HERE, 0,"nbits_in_p <= MAX_PRIMALITY_TEST_BITS");
// For purposes of the bits-in-p limit, treat 2^findex as having (findex) rather than (findex+1) bits:
nbits_in_p = findex;
}
else
ASSERT(HERE, 0,"MODULUS_TYPE unknown!");
ASSERT(HERE, nbits_in_p <= MAX_EXPO_BITS,"Require nbits_in_p <= MAX_EXPO_BITS");
#if INCLUDE_TF
INIT_TF:
/* If nbits_in_p > MAX_PRIMALITY_TEST_BITS, it better be a TF run: */
if(TEST_TYPE == TEST_TYPE_TF) {
/* Currently TF only supported for Mersennes: */
if(MODULUS_TYPE != MODULUS_TYPE_MERSENNE) {
sprintf(cbuf, "ERROR: Trial-factoring Currently only supported for Mersenne numbers. The ini file entry was %s\n",in_line);
fprintf(stderr,"%s",cbuf);
goto GET_NEXT_ASSIGNMENT;
}
/* For now, always start at k = 1: */
log2_min_factor = 0.0;
log2_max_factor = get_default_factoring_depth(p);
ASSERT(HERE, log2_max_factor <= MAX_FACT_BITS, "log2_max_factor > MAX_FACT_BITS!");
/* Field following the exponent is the already-factored-to depth: if none found, use defaults. */
char_addr = strstr(char_addr, ",");
if(char_addr++) {
/* Convert the ensuing numeric digits to ulong: */
TF_BITS = strtoul(char_addr, &endp, 10);
/* Specified already-factored-to depth is larger than default factor-to depth - no more factoring to be done. */
if(TF_BITS > log2_max_factor) {
sprintf(cbuf, "INFO: the specified already-factored-to depth of %u bits exceeds the default %10.4f bits - no more factoring to be done.\n", TF_BITS, log2_max_factor);
fprintf(stderr,"%s",cbuf);
goto GET_NEXT_ASSIGNMENT;
}
}
/* Mode to allow user to override normal automated-test TF default limits:
If a second ,{#bits} argument is present in the input line, then this user-set
desired factoring depth overrides the normal default for the exponent in question.
In this mode, warn if TF-to depth greater than automated-mode default, but allow: */
if(char_addr)
char_addr = strstr(char_addr, ",");
if(char_addr++) {
bit_depth_todo = strtoul(char_addr, &endp, 10);
if(bit_depth_todo > MAX_FACT_BITS) {
sprintf(cbuf, "ERROR: factor-to bit_depth of %u > max. allowed of %u. The ini file entry was %s\n",TF_BITS,MAX_FACT_BITS,in_line);
fprintf(stderr,"%s",cbuf);
goto GET_NEXT_ASSIGNMENT;
}
else if(bit_depth_todo <= TF_BITS) {
sprintf(cbuf, "ERROR: factor-to bit_depth of %u < already-done depth of %u. The ini file entry was %s\n",TF_BITS,TF_BITS,in_line);
fprintf(stderr,"%s",cbuf);
goto GET_NEXT_ASSIGNMENT;
}
else if(bit_depth_todo > log2_max_factor) {
sprintf(cbuf, "WARN: the specified factor-to depth of %u bits exceeds the default %10.4f bits - I hope you know what you're doing.\n",TF_BITS,log2_max_factor);
log2_max_factor = bit_depth_todo;
fprintf(stderr,"%s",cbuf);
}
log2_max_factor = bit_depth_todo;
}
}
else if(nbits_in_p > MAX_PRIMALITY_TEST_BITS)
{
sprintf(cbuf, "ERROR: Inputs this large only permitted for trial-factoring. The ini file entry was %s\n",in_line);
fprintf(stderr,"%s",cbuf);
goto GET_NEXT_ASSIGNMENT;
}
#endif // INCLUDE_TF
/* If "Test..." or "DoubleCheck", check for TF_BITS and pm1_done fields following the = sign:
if present and there is still factoring remaining to be done, modify the assignment type appropriately.
Assignment format:
<Test|DoubleCheck>=[aid],p,TF_BITS,pm1_done?
*/
if(TEST_TYPE == TEST_TYPE_PRIMALITY) {
/* TF_BITS: */
TF_BITS = 0xffffffff; /* Only check if there's an appropriate entry in the input line */
char_addr = strstr(char_addr, ",");
if(char_addr++) {
/* Convert the ensuing numeric digits to ulong: */
TF_BITS = strtoul(char_addr, &endp, 10);
#if INCLUDE_TF
if(TF_BITS > MAX_FACT_BITS) {
snprintf_nowarn(cbuf,STR_MAX_LEN,"ERROR: TF_BITS of %u > max. allowed of %u. The ini file entry was %s\n", TF_BITS, MAX_FACT_BITS, in_line);
fprintf(stderr,"%s",cbuf);
goto GET_NEXT_ASSIGNMENT;
}
/* If TF_BITS less than default TF depth for this exponent
and this platform is approved for trial factoring, switch task type to TF:
*/
log2_max_factor = get_default_factoring_depth(p);
if(TF_BITS < log2_max_factor) {
TEST_TYPE = TEST_TYPE_TF;
// For now, always start at k = 1:
log2_min_factor = 0.0;
}
else if(TF_BITS > log2_max_factor) {
sprintf(cbuf, "WARN: the specified already-factored-to depth of %u bits exceeds the default %10.4f bits - no more factoring to be done.\n", TF_BITS, log2_max_factor);
fprintf(stderr,"%s",cbuf);
}
#endif
// p-1 factoring already done?
char_addr = strstr(char_addr, ",");
if(char_addr++) {
/* Convert the ensuing numeric digits to ulong: */
pm1_done = strtoul(char_addr, &endp, 10);
if(pm1_done > 1) {
sprintf(cbuf, "ERROR: the specified pm1_done field [%u] should be 0 or 1!\n",pm1_done);
ASSERT(HERE,0,cbuf);
}
if(!pm1_done) { // pm1_done == TRUE is a no-op, translating to "proceed with primality test"
// Don't actually use this in pm1_set_bounds(), due to the rise of the single-shot PRP-with-proof paradigm, but for form's sake:
tests_saved = 1;
// Create p-1 assignment, then edit original assignment line appropriately
TEST_TYPE = TEST_TYPE_PM1;
kblocks = get_default_fft_length(p);
ASSERT(HERE, pm1_set_bounds(p, kblocks<<10, TF_BITS, tests_saved), "Failed to set p-1 bounds!");
// Format the p-1 assignment into cbuf:
char_addr = strstr(in_line, "="); ASSERT(HERE,char_addr != 0x0,"Malformed assignment!");
char_addr++; while(isspace(*char_addr)) { ++char_addr; } // Skip any whitespace following the equals sign
if(is_hex_string(char_addr, 32)) {
strncpy(aid,char_addr,32); sprintf(cbuf,"Pminus1=%s,1,2,%llu,-1,%u,%llu\n",aid,p,B1,B2); // If we get here, it's a M(p), not F(m)
} else
sprintf(cbuf,"Pminus1=1,2,%llu,-1,%u,%llu\n",p,B1,B2);
// Copy all but the final (pm1_done) char of the assignment into cstr and append pm1_done = 1. If in_line ends with newline, first --j:
j = strlen(in_line) - 1; j -= (in_line[j] == '\n');
strncpy(cstr,in_line,j); cstr[j] = '\0'; strcat(cstr,"1\n");
split_curr_assignment = TRUE; // This will trigger the corresponding code following the goto:
goto GET_NEXT_ASSIGNMENT;
}