/*
References:
- saffah_cc_v41_260924, https://duck.ac/submission/118929: directly reused FFT,
carry formatter, parser and startup. Public source has no separate license.
- https://duck.ac/problem/1004: both operands have exactly one million digits.
Idea: Replace the first 1024 runtime sin/cos pairs with precomputed hexadecimal
double constants generated from the identical arguments using Python math.
Keep both sine and cosine arrays independently to preserve rounding; keep
runtime generation for the other root tables. Also skip guaranteed-length
full-input boundary scans. This removes transcendental calls, not FFT work.
Experiment purpose: Measure exact-root precomputation on the official workload.
*/
#define DIAGTILE NC
#define TIL DIAGTILE
#pragma GCC target("arch=skylake","tune=skylake")
#ifndef FUSE16
#define FUSE16 0
#endif
#ifndef NTPARSE
#define NTPARSE 1
#endif
#ifndef NOMEMSET
#define NOMEMSET 1
#endif
#ifndef PFV
#define PFV 4
#endif
#ifndef PFW
#define PFW 2
#endif
#ifndef PFV
#define PFV 0
#endif
#ifndef PFW
#define PFW 0
#endif
#ifndef ROWPF
#define ROWPF 0
#endif
#ifndef NTPARSE
#define NTPARSE 0
#endif
#ifndef NTCARRY
#define NTCARRY 0
#endif
#ifndef FUSE16
#define FUSE16 0
#endif
#ifndef NOMEMSET
#define NOMEMSET 0
#endif
#pragma GCC optimize("O3","unroll-loops","inline-functions","ipa-cp-clone","ipa-sra")
#include <cstdint>
#include <cstring>
#include <cmath>
#include <immintrin.h>
#ifndef CBW
#define CBW 16
#endif
#ifndef ZSKIP
#define ZSKIP 1
#endif
#ifndef CBW2
#define CBW2 16
#endif
#ifndef PF
#define PF 1
#endif
#define LOGN 19
#define NN (1<<LOGN)
#define NR (1<<10)
#define NC (1<<9)
#define STRIDE (NC+16)
#define ARRSZ ((size_t)NR*STRIDE)
typedef uint32_t u32; typedef uint64_t u64; typedef int64_t i64;
typedef double f64;
struct DI {
unsigned long abi;
const char *s; unsigned long sn;
char *o; unsigned long ol; unsigned long os;
char *e; unsigned long el; unsigned long es;
const char *IB; unsigned long IBl;
char *OB; unsigned long OBl;
unsigned long tsc;
} __attribute__((packed));
#define TGT
#define FT
alignas(4096) static f64 RE_[2*ARRSZ + 256 + 2048];
static f64 *const RE = RE_;
static f64 *const IM = RE_ + ARRSZ + 128 + 2048;
static f64 TWR[2][NC], TWI[2][NC];
alignas(64) static f64 RDr[NC], RDi[NC], RIr[NC], RIi[NC];
static f64 CDr[NR], CDi[NR], CIr[NR], CIi[NR];
static f64 DQr[NR], DQi[NR], DTr[NC], DTi[NC];
static u32 BREV[NR], BREVC[NC], BREV9[512], CJ0[NC], CJ1[NC];
static f64 DQ18[512], DQ18i[512], DT18[NC], DT18i[NC];
#define NN2 ((1<<18))
static u32 RVIX[NR], CVIX0[NC], CVIX1[NC];
static f64 LIN[(size_t)NN + 8];
static char OUTB[(size_t)2*1000000 + 2*1000000/8 + 4096];
static u32 D8[(1<<19)];
static char DIG2T[201] =
"00010203040506070809101112131415161718192021222324252627282930313233343536373839"
"40414243444546474849505152535455565758596061626364656667686970717273747576777879"
"8081828384858687888990919293949596979899";
FT static void build_r4(void);
FT static void build_all(void) {
static const double cr[1024]={0x1.0000000000000p+0,0x1.fffd8858e8a92p-1,0x1.fff62169b92dbp-1,0x1.ffe9cb44b51a1p-1,0x1.ffd886084cd0dp-1,0x1.ffc251df1d3f8p-1,0x1.ffa72effef75dp-1,0x1.ff871dadb81dfp-1,0x1.ff621e3796d7ep-1,0x1.ff3830f8d575cp-1,0x1.ff095658e71adp-1,0x1.fed58ecb673c4p-1,0x1.fe9cdad01883ap-1,0x1.fe5f3af2e3940p-1,0x1.fe1cafcbd5b09p-1,0x1.fdd539ff1f456p-1,0x1.fd88da3d12526p-1,0x1.fd37914220b84p-1,0x1.fce15fd6da67bp-1,0x1.fc8646cfeb721p-1,0x1.fc26470e19fd3p-1,0x1.fbc1617e44186p-1,0x1.fb5797195d741p-1,0x1.fae8e8e46cfbbp-1,0x1.fa7557f08a517p-1,0x1.f9fce55adb2c8p-1,0x1.f97f924c9099bp-1,0x1.f8fd5ffae41dbp-1,0x1.f8764fa714ba9p-1,0x1.f7ea629e63d6ep-1,0x1.f7599a3a12077p-1,0x1.f6c3f7df5bbb7p-1,0x1.f6297cff75cb0p-1,0x1.f58a2b1789e84p-1,0x1.f4e603b0b2f2dp-1,0x1.f43d085ff92ddp-1,0x1.f38f3ac64e589p-1,0x1.f2dc9c9089a9dp-1,0x1.f2252f7763adap-1,0x1.f168f53f7205dp-1,0x1.f0a7efb9230d7p-1,0x1.efe220c0b95edp-1,0x1.ef178a3e473c2p-1,0x1.ee482e25a9dbcp-1,0x1.ed740e7684963p-1,0x1.ec9b2d3c3bf84p-1,0x1.ebbd8c8df0b74p-1,0x1.eadb2e8e7a88ep-1,0x1.e9f4156c62ddap-1,0x1.e9084361df7f3p-1,0x1.e817bab4cd10dp-1,0x1.e7227db6a9744p-1,0x1.e6288ec48e112p-1,0x1.e529f04729ffcp-1,0x1.e426a4b2bc17ep-1,0x1.e31eae870ce25p-1,0x1.e212104f686e5p-1,0x1.e100cca2980acp-1,0x1.dfeae622dbe2bp-1,0x1.ded05f7de47dap-1,0x1.ddb13b6ccc23dp-1,0x1.dc8d7cb410260p-1,0x1.db6526238a09bp-1,0x1.da383a9668988p-1,0x1.d906bcf328d46p-1,0x1.d7d0b02b8ecfap-1,0x1.d696173c9e68bp-1,0x1.d556f52e93eb1p-1,0x1.d4134d14dc93ap-1,0x1.d2cb220e0ef9fp-1,0x1.d17e7743e35dcp-1,0x1.d02d4feb2bd92p-1,0x1.ced7af43cc773p-1,0x1.cd7d9898b32f6p-1,0x1.cc1f0f3fcfc5cp-1,0x1.cabc169a0b901p-1,0x1.c954b213411f5p-1,0x1.c7e8e52233cf3p-1,0x1.c678b3488739bp-1,0x1.c5042012b6907p-1,0x1.c38b2f180bdb1p-1,0x1.c20de3fa971b0p-1,0x1.c08c426725549p-1,0x1.bf064e15377ddp-1,0x1.bd7c0ac6f952ap-1,0x1.bbed7c49380eap-1,0x1.ba5aa673590d3p-1,0x1.b8c38d27504e9p-1,0x1.b728345196e3ep-1,0x1.b5889fe921405p-1,0x1.b3e4d3ef55712p-1,0x1.b23cd470013b4p-1,0x1.b090a58150200p-1,0x1.aee04b43c1474p-1,0x1.ad2bc9e21d511p-1,0x1.ab7325916c0d4p-1,0x1.a9b66290ea1a3p-1,0x1.a7f58529fe69dp-1,0x1.a63091b02fae2p-1,0x1.a4678c8119ac8p-1,0x1.a29a7a0462782p-1,0x1.a0c95eabaf937p-1,0x1.9ef43ef29af94p-1,0x1.9d1b1f5ea80d6p-1,0x1.9b3e047f38741p-1,0x1.995cf2ed80d22p-1,0x1.9777ef4c7d742p-1,0x1.958efe48e6dd7p-1,0x1.93a22499263fcp-1,0x1.91b166fd49da2p-1,0x1.8fbcca3ef940dp-1,0x1.8dc45331698ccp-1,0x1.8bc806b151741p-1,0x1.89c7e9a4dd4abp-1,0x1.87c400fba2ebfp-1,0x1.85bc51ae958ccp-1,0x1.83b0e0bff976ep-1,0x1.81a1b33b57accp-1,0x1.7f8ece3571771p-1,0x1.7d7836cc33db3p-1,0x1.7b5df226aafafp-1,0x1.79400574f55e5p-1,0x1.771e75f037261p-1,0x1.74f948da8d28dp-1,0x1.72d0837efff97p-1,0x1.70a42b3176d7ap-1,0x1.6e74454eaa8aep-1,0x1.6c40d73c18275p-1,0x1.6a09e667f3bcdp-1,0x1.67cf78491af10p-1,0x1.6591925f0783ep-1,0x1.63503a31c1be9p-1,0x1.610b7551d2cdfp-1,0x1.5ec3495837074p-1,0x1.5c77bbe65018dp-1,0x1.5a28d2a5d7251p-1,0x1.57d69348cec9fp-1,0x1.5581038975137p-1,0x1.5328292a35596p-1,0x1.50cc09f59a09bp-1,0x1.4e6cabbe3e5e9p-1,0x1.4c0a145ec0005p-1,0x1.49a449b9b0939p-1,0x1.473b51b987347p-1,0x1.44cf325091dd6p-1,0x1.425ff178e6bb2p-1,0x1.3fed9534556d5p-1,0x1.3d78238c58344p-1,0x1.3affa292050b9p-1,0x1.3884185dfeb22p-1,0x1.36058b10659f3p-1,0x1.338400d0c8e57p-1,0x1.30ff7fce17036p-1,0x1.2e780e3e8ea17p-1,0x1.2bedb25faf3eap-1,0x1.2960727629ca8p-1,0x1.26d054cdd12dfp-1,0x1.243d5fb98ac20p-1,0x1.21a799933eb59p-1,0x1.1f0f08bbc861bp-1,0x1.1c73b39ae68c9p-1,0x1.19d5a09f2b9b8p-1,0x1.1734d63dedb49p-1,0x1.14915af336cecp-1,0x1.11eb3541b4b23p-1,0x1.0f426bb2a8e7fp-1,0x1.0c9704d5d898ep-1,0x1.09e907417c5e0p-1,0x1.073879922ffedp-1,0x1.0485626ae221ap-1,0x1.01cfc874c3eb7p-1,0x1.fe2f64be71210p-2,0x1.f8ba4dbf89abbp-2,0x1.f3405963fd069p-2,0x1.edc1952ef78d7p-2,0x1.e83e0eaf85116p-2,0x1.e2b5d3806f63ep-2,0x1.dd28f1481cc57p-2,0x1.d79775b86e389p-2,0x1.d2016e8e9db5bp-2,0x1.cc66e9931c45ep-2,0x1.c6c7f4997000bp-2,0x1.c1249d8011ee8p-2,0x1.bb7cf2304bd02p-2,0x1.b5d1009e15cc2p-2,0x1.b020d6c7f400bp-2,0x1.aa6c82b6d3fccp-2,0x1.a4b4127dea1e4p-2,0x1.9ef7943a8ed8ap-2,0x1.993716141bdfep-2,0x1.9372a63bc93d7p-2,0x1.8daa52ec8a4b0p-2,0x1.87de2a6aea964p-2,0x1.820e3b04eaac5p-2,0x1.7c3a9311dcce8p-2,0x1.766340f2418f8p-2,0x1.7088530fa45a1p-2,0x1.6aa9d7dc77e19p-2,0x1.64c7ddd3f27c5p-2,0x1.5ee27379ea693p-2,0x1.58f9a75ab1fddp-2,0x1.530d880af3c24p-2,0x1.4d1e24278e76bp-2,0x1.472b8a5571055p-2,0x1.4135c94176603p-2,0x1.3b3cefa0414b9p-2,0x1.35410c2e18154p-2,0x1.2f422daec0389p-2,0x1.294062ed59f05p-2,0x1.233bbabc3bb71p-2,0x1.1d3443f4cdb3dp-2,0x1.172a0d7765177p-2,0x1.111d262b1f678p-2,0x1.0b0d9cfdbdb91p-2,0x1.04fb80e37fdafp-2,0x1.fdcdc1adfedfcp-3,0x1.f19f97b215f1ep-3,0x1.e56ca1e101a20p-3,0x1.d934fe5454317p-3,0x1.ccf8cb312b284p-3,0x1.c0b826a7e4f62p-3,0x1.b4732ef3d6722p-3,0x1.a82a025b00451p-3,0x1.9bdcbf2dc4368p-3,0x1.8f8b83c69a60dp-3,0x1.83366e89c64c8p-3,0x1.76dd9de50bf35p-3,0x1.6a81304f64ab6p-3,0x1.5e214448b3fcbp-3,0x1.51bdf8597c5f8p-3,0x1.45576b1293e58p-3,0x1.38edbb0cd8d13p-3,0x1.2c8106e8e613ap-3,0x1.20116d4ec7bcfp-3,0x1.139f0cedaf578p-3,0x1.072a047ba831fp-3,0x1.f564e56a97314p-4,0x1.dc70ecbae9fd1p-4,0x1.c3785c79ec2dep-4,0x1.aa7b724495c0ep-4,0x1.917a6bc29b438p-4,0x1.787586a5d5b1fp-4,0x1.5f6d00a9aa418p-4,0x1.4661179272096p-4,0x1.2d52092ce19f8p-4,0x1.1440134d709b6p-4,0x1.f656e79f820eap-5,0x1.c428d12c0d7f0p-5,0x1.91f65f10dd824p-5,0x1.5fc00d290cd57p-5,0x1.2d865759455e4p-5,0x1.f693731d1cef4p-6,0x1.92155f7a36677p-6,0x1.2d936bbe30efdp-6,0x1.921d1fcdec78fp-7,0x1.921f0fe67009fp-8,0x1.1a62633145c07p-54,-0x1.921f0fe670012p-8,-0x1.921d1fcdec749p-7,-0x1.2d936bbe30ed9p-6,-0x1.92155f7a36654p-6,-0x1.f693731d1ced1p-6,-0x1.2d865759455d2p-5,-0x1.5fc00d290cd45p-5,-0x1.91f65f10dd813p-5,-0x1.c428d12c0d7dfp-5,-0x1.f656e79f820d9p-5,-0x1.1440134d709adp-4,-0x1.2d52092ce19efp-4,-0x1.466117927208ep-4,-0x1.5f6d00a9aa40fp-4,-0x1.787586a5d5b16p-4,-0x1.917a6bc29b42fp-4,-0x1.aa7b724495c05p-4,-0x1.c3785c79ec2d5p-4,-0x1.dc70ecbae9fc8p-4,-0x1.f564e56a9730bp-4,-0x1.072a047ba831bp-3,-0x1.139f0cedaf574p-3,-0x1.20116d4ec7bcbp-3,-0x1.2c8106e8e6136p-3,-0x1.38edbb0cd8d0fp-3,-0x1.45576b1293e54p-3,-0x1.51bdf8597c5f3p-3,-0x1.5e214448b3fc7p-3,-0x1.6a81304f64ab2p-3,-0x1.76dd9de50bf30p-3,-0x1.83366e89c64c4p-3,-0x1.8f8b83c69a608p-3,-0x1.9bdcbf2dc4363p-3,-0x1.a82a025b0044dp-3,-0x1.b4732ef3d671ep-3,-0x1.c0b826a7e4f5ep-3,-0x1.ccf8cb312b280p-3,-0x1.d934fe5454312p-3,-0x1.e56ca1e101a1cp-3,-0x1.f19f97b215f1ap-3,-0x1.fdcdc1adfedf7p-3,-0x1.04fb80e37fdadp-2,-0x1.0b0d9cfdbdb8fp-2,-0x1.111d262b1f676p-2,-0x1.172a0d7765175p-2,-0x1.1d3443f4cdb3bp-2,-0x1.233bbabc3bb6fp-2,-0x1.294062ed59f02p-2,-0x1.2f422daec0387p-2,-0x1.35410c2e18152p-2,-0x1.3b3cefa0414b7p-2,-0x1.4135c94176600p-2,-0x1.472b8a5571053p-2,-0x1.4d1e24278e769p-2,-0x1.530d880af3c22p-2,-0x1.58f9a75ab1fdbp-2,-0x1.5ee27379ea691p-2,-0x1.64c7ddd3f27c3p-2,-0x1.6aa9d7dc77e17p-2,-0x1.7088530fa459fp-2,-0x1.766340f2418f6p-2,-0x1.7c3a9311dcce6p-2,-0x1.820e3b04eaac3p-2,-0x1.87de2a6aea962p-2,-0x1.8daa52ec8a4aep-2,-0x1.9372a63bc93d5p-2,-0x1.993716141bdfcp-2,-0x1.9ef7943a8ed88p-2,-0x1.a4b4127dea1e2p-2,-0x1.aa6c82b6d3fc6p-2,-0x1.b020d6c7f4009p-2,-0x1.b5d1009e15cbcp-2,-0x1.bb7cf2304bd00p-2,-0x1.c1249d8011ee2p-2,-0x1.c6c7f49970009p-2,-0x1.cc66e9931c460p-2,-0x1.d2016e8e9db59p-2,-0x1.d79775b86e38ap-2,-0x1.dd28f1481cc55p-2,-0x1.e2b5d3806f63cp-2,-0x1.e83e0eaf85110p-2,-0x1.edc1952ef78d5p-2,-0x1.f3405963fd063p-2,-0x1.f8ba4dbf89ab9p-2,-0x1.fe2f64be7120bp-2,-0x1.01cfc874c3eb6p-1,-0x1.0485626ae221bp-1,-0x1.073879922ffedp-1,-0x1.09e907417c5e1p-1,-0x1.0c9704d5d898dp-1,-0x1.0f426bb2a8e7ep-1,-0x1.11eb3541b4b21p-1,-0x1.14915af336cebp-1,-0x1.1734d63dedb47p-1,-0x1.19d5a09f2b9b7p-1,-0x1.1c73b39ae68c6p-1,-0x1.1f0f08bbc861ap-1,-0x1.21a799933eb59p-1,-0x1.243d5fb98ac1ep-1,-0x1.26d054cdd12dfp-1,-0x1.2960727629ca7p-1,-0x1.2bedb25faf3eap-1,-0x1.2e780e3e8ea15p-1,-0x1.30ff7fce17035p-1,-0x1.338400d0c8e55p-1,-0x1.36058b10659f2p-1,-0x1.3884185dfeb23p-1,-0x1.3affa292050b8p-1,-0x1.3d78238c58344p-1,-0x1.3fed9534556d3p-1,-0x1.425ff178e6bb2p-1,-0x1.44cf325091dd5p-1,-0x1.473b51b987347p-1,-0x1.49a449b9b0937p-1,-0x1.4c0a145ec0004p-1,-0x1.4e6cabbe3e5e7p-1,-0x1.50cc09f59a09bp-1,-0x1.5328292a35597p-1,-0x1.5581038975136p-1,-0x1.57d69348ceca0p-1,-0x1.5a28d2a5d724fp-1,-0x1.5c77bbe65018cp-1,-0x1.5ec3495837073p-1,-0x1.610b7551d2cdfp-1,-0x1.63503a31c1be7p-1,-0x1.6591925f0783dp-1,-0x1.67cf78491af0ep-1,-0x1.6a09e667f3bccp-1,-0x1.6c40d73c18276p-1,-0x1.6e74454eaa8aep-1,-0x1.70a42b3176d7ap-1,-0x1.72d0837efff95p-1,-0x1.74f948da8d28dp-1,-0x1.771e75f037260p-1,-0x1.79400574f55e5p-1,-0x1.7b5df226aafadp-1,-0x1.7d7836cc33db2p-1,-0x1.7f8ece357176fp-1,-0x1.81a1b33b57acbp-1,-0x1.83b0e0bff976ep-1,-0x1.85bc51ae958cbp-1,-0x1.87c400fba2ebfp-1,-0x1.89c7e9a4dd4a9p-1,-0x1.8bc806b151741p-1,-0x1.8dc45331698cbp-1,-0x1.8fbcca3ef940dp-1,-0x1.91b166fd49da0p-1,-0x1.93a22499263fbp-1,-0x1.958efe48e6dd5p-1,-0x1.9777ef4c7d741p-1,-0x1.995cf2ed80d23p-1,-0x1.9b3e047f38740p-1,-0x1.9d1b1f5ea80d6p-1,-0x1.9ef43ef29af93p-1,-0x1.a0c95eabaf937p-1,-0x1.a29a7a0462781p-1,-0x1.a4678c8119ac8p-1,-0x1.a63091b02fae0p-1,-0x1.a7f58529fe69cp-1,-0x1.a9b66290ea1a4p-1,-0x1.ab7325916c0d4p-1,-0x1.ad2bc9e21d511p-1,-0x1.aee04b43c1473p-1,-0x1.b090a58150200p-1,-0x1.b23cd470013b3p-1,-0x1.b3e4d3ef55712p-1,-0x1.b5889fe921404p-1,-0x1.b728345196e3dp-1,-0x1.b8c38d27504e7p-1,-0x1.ba5aa673590d2p-1,-0x1.bbed7c49380ebp-1,-0x1.bd7c0ac6f9529p-1,-0x1.bf064e15377ddp-1,-0x1.c08c426725548p-1,-0x1.c20de3fa971b0p-1,-0x1.c38b2f180bdb0p-1,-0x1.c5042012b6907p-1,-0x1.c678b3488739ap-1,-0x1.c7e8e52233cf3p-1,-0x1.c954b213411f4p-1,-0x1.cabc169a0b900p-1,-0x1.cc1f0f3fcfc5dp-1,-0x1.cd7d9898b32f5p-1,-0x1.ced7af43cc773p-1,-0x1.d02d4feb2bd92p-1,-0x1.d17e7743e35dcp-1,-0x1.d2cb220e0ef9ep-1,-0x1.d4134d14dc93ap-1,-0x1.d556f52e93eb0p-1,-0x1.d696173c9e68bp-1,-0x1.d7d0b02b8ecf8p-1,-0x1.d906bcf328d46p-1,-0x1.da383a9668988p-1,-0x1.db6526238a09ap-1,-0x1.dc8d7cb410260p-1,-0x1.ddb13b6ccc23cp-1,-0x1.ded05f7de47dap-1,-0x1.dfeae622dbe2ap-1,-0x1.e100cca2980acp-1,-0x1.e212104f686e4p-1,-0x1.e31eae870ce25p-1,-0x1.e426a4b2bc17dp-1,-0x1.e529f04729ffcp-1,-0x1.e6288ec48e112p-1,-0x1.e7227db6a9744p-1,-0x1.e817bab4cd10dp-1,-0x1.e9084361df7f2p-1,-0x1.e9f4156c62ddap-1,-0x1.eadb2e8e7a88dp-1,-0x1.ebbd8c8df0b74p-1,-0x1.ec9b2d3c3bf83p-1,-0x1.ed740e7684963p-1,-0x1.ee482e25a9dbbp-1,-0x1.ef178a3e473c2p-1,-0x1.efe220c0b95edp-1,-0x1.f0a7efb9230d7p-1,-0x1.f168f53f7205dp-1,-0x1.f2252f7763ad9p-1,-0x1.f2dc9c9089a9dp-1,-0x1.f38f3ac64e588p-1,-0x1.f43d085ff92ddp-1,-0x1.f4e603b0b2f2cp-1,-0x1.f58a2b1789e84p-1,-0x1.f6297cff75cb0p-1,-0x1.f6c3f7df5bbb7p-1,-0x1.f7599a3a12077p-1,-0x1.f7ea629e63d6ep-1,-0x1.f8764fa714ba9p-1,-0x1.f8fd5ffae41dap-1,-0x1.f97f924c9099bp-1,-0x1.f9fce55adb2c8p-1,-0x1.fa7557f08a517p-1,-0x1.fae8e8e46cfbap-1,-0x1.fb5797195d741p-1,-0x1.fbc1617e44186p-1,-0x1.fc26470e19fd3p-1,-0x1.fc8646cfeb721p-1,-0x1.fce15fd6da67bp-1,-0x1.fd37914220b84p-1,-0x1.fd88da3d12525p-1,-0x1.fdd539ff1f456p-1,-0x1.fe1cafcbd5b09p-1,-0x1.fe5f3af2e3940p-1,-0x1.fe9cdad01883ap-1,-0x1.fed58ecb673c4p-1,-0x1.ff095658e71adp-1,-0x1.ff3830f8d575cp-1,-0x1.ff621e3796d7ep-1,-0x1.ff871dadb81dfp-1,-0x1.ffa72effef75dp-1,-0x1.ffc251df1d3f8p-1,-0x1.ffd886084cd0dp-1,-0x1.ffe9cb44b51a1p-1,-0x1.fff62169b92dbp-1,-0x1.fffd8858e8a92p-1,-0x1.0000000000000p+0,-0x1.fffd8858e8a92p-1,-0x1.fff62169b92dbp-1,-0x1.ffe9cb44b51a1p-1,-0x1.ffd886084cd0dp-1,-0x1.ffc251df1d3f8p-1,-0x1.ffa72effef75dp-1,-0x1.ff871dadb81dfp-1,-0x1.ff621e3796d7ep-1,-0x1.ff3830f8d575cp-1,-0x1.ff095658e71adp-1,-0x1.fed58ecb673c4p-1,-0x1.fe9cdad01883ap-1,-0x1.fe5f3af2e3941p-1,-0x1.fe1cafcbd5b09p-1,-0x1.fdd539ff1f456p-1,-0x1.fd88da3d12526p-1,-0x1.fd37914220b85p-1,-0x1.fce15fd6da67bp-1,-0x1.fc8646cfeb722p-1,-0x1.fc26470e19fd3p-1,-0x1.fbc1617e44186p-1,-0x1.fb5797195d741p-1,-0x1.fae8e8e46cfbbp-1,-0x1.fa7557f08a517p-1,-0x1.f9fce55adb2c8p-1,-0x1.f97f924c9099bp-1,-0x1.f8fd5ffae41dbp-1,-0x1.f8764fa714baap-1,-0x1.f7ea629e63d6fp-1,-0x1.f7599a3a12078p-1,-0x1.f6c3f7df5bbb7p-1,-0x1.f6297cff75cb0p-1,-0x1.f58a2b1789e84p-1,-0x1.f4e603b0b2f2dp-1,-0x1.f43d085ff92ddp-1,-0x1.f38f3ac64e589p-1,-0x1.f2dc9c9089a9dp-1,-0x1.f2252f7763adap-1,-0x1.f168f53f7205ep-1,-0x1.f0a7efb9230d7p-1,-0x1.efe220c0b95edp-1,-0x1.ef178a3e473c3p-1,-0x1.ee482e25a9dbcp-1,-0x1.ed740e7684964p-1,-0x1.ec9b2d3c3bf84p-1,-0x1.ebbd8c8df0b75p-1,-0x1.eadb2e8e7a88ep-1,-0x1.e9f4156c62ddbp-1,-0x1.e9084361df7f3p-1,-0x1.e817bab4cd10ep-1,-0x1.e7227db6a9745p-1,-0x1.e6288ec48e113p-1,-0x1.e529f04729ffdp-1,-0x1.e426a4b2bc17ep-1,-0x1.e31eae870ce25p-1,-0x1.e212104f686e5p-1,-0x1.e100cca2980adp-1,-0x1.dfeae622dbe2bp-1,-0x1.ded05f7de47dbp-1,-0x1.ddb13b6ccc23dp-1,-0x1.dc8d7cb410261p-1,-0x1.db6526238a09bp-1,-0x1.da383a9668989p-1,-0x1.d906bcf328d47p-1,-0x1.d7d0b02b8ecf9p-1,-0x1.d696173c9e68bp-1,-0x1.d556f52e93eb1p-1,-0x1.d4134d14dc93ap-1,-0x1.d2cb220e0ef9fp-1,-0x1.d17e7743e35ddp-1,-0x1.d02d4feb2bd93p-1,-0x1.ced7af43cc774p-1,-0x1.cd7d9898b32f6p-1,-0x1.cc1f0f3fcfc5ep-1,-0x1.cabc169a0b901p-1,-0x1.c954b213411f5p-1,-0x1.c7e8e52233cf4p-1,-0x1.c678b3488739bp-1,-0x1.c5042012b6908p-1,-0x1.c38b2f180bdb1p-1,-0x1.c20de3fa971b1p-1,-0x1.c08c426725549p-1,-0x1.bf064e15377dep-1,-0x1.bd7c0ac6f952ap-1,-0x1.bbed7c49380ecp-1,-0x1.ba5aa673590d3p-1,-0x1.b8c38d27504e9p-1,-0x1.b728345196e3ep-1,-0x1.b5889fe921405p-1,-0x1.b3e4d3ef55713p-1,-0x1.b23cd470013b4p-1,-0x1.b090a58150201p-1,-0x1.aee04b43c1474p-1,-0x1.ad2bc9e21d512p-1,-0x1.ab7325916c0d5p-1,-0x1.a9b66290ea1a5p-1,-0x1.a7f58529fe69dp-1,-0x1.a63091b02fae2p-1,-0x1.a4678c8119ac9p-1,-0x1.a29a7a0462782p-1,-0x1.a0c95eabaf938p-1,-0x1.9ef43ef29af94p-1,-0x1.9d1b1f5ea80d7p-1,-0x1.9b3e047f38741p-1,-0x1.995cf2ed80d24p-1,-0x1.9777ef4c7d742p-1,-0x1.958efe48e6dd6p-1,-0x1.93a22499263fcp-1,-0x1.91b166fd49da2p-1,-0x1.8fbcca3ef940ep-1,-0x1.8dc45331698ccp-1,-0x1.8bc806b151742p-1,-0x1.89c7e9a4dd4abp-1,-0x1.87c400fba2ec0p-1,-0x1.85bc51ae958cdp-1,-0x1.83b0e0bff9770p-1,-0x1.81a1b33b57acdp-1,-0x1.7f8ece3571770p-1,-0x1.7d7836cc33db3p-1,-0x1.7b5df226aafafp-1,-0x1.79400574f55e6p-1,-0x1.771e75f037261p-1,-0x1.74f948da8d28fp-1,-0x1.72d0837efff97p-1,-0x1.70a42b3176d7cp-1,-0x1.6e74454eaa8afp-1,-0x1.6c40d73c18277p-1,-0x1.6a09e667f3bcep-1,-0x1.67cf78491af10p-1,-0x1.6591925f0783fp-1,-0x1.63503a31c1be9p-1,-0x1.610b7551d2ce0p-1,-0x1.5ec3495837074p-1,-0x1.5c77bbe65018ep-1,-0x1.5a28d2a5d7251p-1,-0x1.57d69348ceca2p-1,-0x1.5581038975138p-1,-0x1.5328292a35598p-1,-0x1.50cc09f59a09cp-1,-0x1.4e6cabbe3e5ecp-1,-0x1.4c0a145ec0002p-1,-0x1.49a449b9b0939p-1,-0x1.473b51b987349p-1,-0x1.44cf325091ddap-1,-0x1.425ff178e6bb0p-1,-0x1.3fed9534556d5p-1,-0x1.3d78238c58346p-1,-0x1.3affa292050bdp-1,-0x1.3884185dfeb21p-1,-0x1.36058b10659f4p-1,-0x1.338400d0c8e5ap-1,-0x1.30ff7fce17033p-1,-0x1.2e780e3e8ea17p-1,-0x1.2bedb25faf3ecp-1,-0x1.2960727629cacp-1,-0x1.26d054cdd12dep-1,-0x1.243d5fb98ac20p-1,-0x1.21a799933eb5bp-1,-0x1.1f0f08bbc861fp-1,-0x1.1c73b39ae68c8p-1,-0x1.19d5a09f2b9b9p-1,-0x1.1734d63dedb4cp-1,-0x1.14915af336ceap-1,-0x1.11eb3541b4b23p-1,-0x1.0f426bb2a8e80p-1,-0x1.0c9704d5d8993p-1,-0x1.09e907417c5dfp-1,-0x1.073879922ffeep-1,-0x1.0485626ae221dp-1,-0x1.01cfc874c3ebbp-1,-0x1.fe2f64be7120fp-2,-0x1.f8ba4dbf89abdp-2,-0x1.f3405963fd06ep-2,-0x1.edc1952ef78d2p-2,-0x1.e83e0eaf85114p-2,-0x1.e2b5d3806f63fp-2,-0x1.dd28f1481cc60p-2,-0x1.d79775b86e387p-2,-0x1.d2016e8e9db5dp-2,-0x1.cc66e9931c464p-2,-0x1.c6c7f49970014p-2,-0x1.c1249d8011ee6p-2,-0x1.bb7cf2304bd04p-2,-0x1.b5d1009e15cc7p-2,-0x1.b020d6c7f4006p-2,-0x1.aa6c82b6d3fcap-2,-0x1.a4b4127dea1e9p-2,-0x1.9ef7943a8ed93p-2,-0x1.993716141bdfdp-2,-0x1.9372a63bc93d9p-2,-0x1.8daa52ec8a4b6p-2,-0x1.87de2a6aea96dp-2,-0x1.820e3b04eaac3p-2,-0x1.7c3a9311dccebp-2,-0x1.766340f2418fep-2,-0x1.7088530fa459bp-2,-0x1.6aa9d7dc77e18p-2,-0x1.64c7ddd3f27cbp-2,-0x1.5ee27379ea69cp-2,-0x1.58f9a75ab1fdbp-2,-0x1.530d880af3c26p-2,-0x1.4d1e24278e771p-2,-0x1.472b8a557104fp-2,-0x1.4135c94176601p-2,-0x1.3b3cefa0414bbp-2,-0x1.35410c2e1815ap-2,-0x1.2f422daec0384p-2,-0x1.294062ed59f07p-2,-0x1.233bbabc3bb77p-2,-0x1.1d3443f4cdb47p-2,-0x1.172a0d7765175p-2,-0x1.111d262b1f67ap-2,-0x1.0b0d9cfdbdb97p-2,-0x1.04fb80e37fdaap-2,-0x1.fdcdc1adfedf8p-3,-0x1.f19f97b215f23p-3,-0x1.e56ca1e101a2cp-3,-0x1.d934fe545430bp-3,-0x1.ccf8cb312b289p-3,-0x1.c0b826a7e4f6ep-3,-0x1.b4732ef3d6736p-3,-0x1.a82a025b0044ep-3,-0x1.9bdcbf2dc436cp-3,-0x1.8f8b83c69a619p-3,-0x1.83366e89c64bdp-3,-0x1.76dd9de50bf31p-3,-0x1.6a81304f64abbp-3,-0x1.5e214448b3fd7p-3,-0x1.51bdf8597c5ecp-3,-0x1.45576b1293e5dp-3,-0x1.38edbb0cd8d1fp-3,-0x1.2c8106e8e614ep-3,-0x1.20116d4ec7bccp-3,-0x1.139f0cedaf57dp-3,-0x1.072a047ba832cp-3,-0x1.f564e56a972fdp-4,-0x1.dc70ecbae9fc9p-4,-0x1.c3785c79ec2e7p-4,-0x1.aa7b724495c27p-4,-0x1.917a6bc29b421p-4,-0x1.787586a5d5b28p-4,-0x1.5f6d00a9aa431p-4,-0x1.46611792720bfp-4,-0x1.2d52092ce19f1p-4,-0x1.1440134d709bfp-4,-0x1.f656e79f8211cp-5,-0x1.c428d12c0d7c2p-5,-0x1.91f65f10dd816p-5,-0x1.5fc00d290cd69p-5,-0x1.2d86575945615p-5,-0x1.f693731d1ced8p-6,-0x1.92155f7a3669bp-6,-0x1.2d936bbe30f60p-6,-0x1.921d1fcdec8d6p-7,-0x1.921f0fe67002dp-8,-0x1.a79394c9e8a0ap-53,0x1.921f0fe66fe85p-8,0x1.921d1fcdec802p-7,0x1.2d936bbe30ef6p-6,0x1.92155f7a36631p-6,0x1.f693731d1ce6ep-6,0x1.2d865759455e0p-5,0x1.5fc00d290cd34p-5,0x1.91f65f10dd7e1p-5,0x1.c428d12c0d78dp-5,0x1.f656e79f820e7p-5,0x1.1440134d709a5p-4,0x1.2d52092ce19d6p-4,0x1.46611792720a5p-4,0x1.5f6d00a9aa416p-4,0x1.787586a5d5b0dp-4,0x1.917a6bc29b407p-4,0x1.aa7b724495c0cp-4,0x1.c3785c79ec2ccp-4,0x1.dc70ecbae9fafp-4,0x1.f564e56a972e3p-4,0x1.072a047ba831fp-3,0x1.139f0cedaf56fp-3,0x1.20116d4ec7bbfp-3,0x1.2c8106e8e6141p-3,0x1.38edbb0cd8d12p-3,0x1.45576b1293e50p-3,0x1.51bdf8597c5dfp-3,0x1.5e214448b3fcap-3,0x1.6a81304f64aaep-3,0x1.76dd9de50bf24p-3,0x1.83366e89c64b0p-3,0x1.8f8b83c69a60cp-3,0x1.9bdcbf2dc435fp-3,0x1.a82a025b00441p-3,0x1.b4732ef3d6729p-3,0x1.c0b826a7e4f61p-3,0x1.ccf8cb312b27cp-3,0x1.d934fe54542fep-3,0x1.e56ca1e101a1fp-3,0x1.f19f97b215f16p-3,0x1.fdcdc1adfedebp-3,0x1.04fb80e37fda3p-2,0x1.0b0d9cfdbdb90p-2,0x1.111d262b1f674p-2,0x1.172a0d776516fp-2,0x1.1d3443f4cdb41p-2,0x1.233bbabc3bb70p-2,0x1.294062ed59f00p-2,0x1.2f422daec037dp-2,0x1.35410c2e18153p-2,0x1.3b3cefa0414b5p-2,0x1.4135c941765fbp-2,0x1.472b8a5571049p-2,0x1.4d1e24278e76ap-2,0x1.530d880af3c20p-2,0x1.58f9a75ab1fd5p-2,0x1.5ee27379ea696p-2,0x1.64c7ddd3f27c5p-2,0x1.6aa9d7dc77e11p-2,0x1.7088530fa4595p-2,0x1.766340f2418f8p-2,0x1.7c3a9311dcce4p-2,0x1.820e3b04eaabdp-2,0x1.87de2a6aea967p-2,0x1.8daa52ec8a4b0p-2,0x1.9372a63bc93d3p-2,0x1.993716141bdf7p-2,0x1.9ef7943a8ed8dp-2,0x1.a4b4127dea1e3p-2,0x1.aa6c82b6d3fc4p-2,0x1.b020d6c7f4000p-2,0x1.b5d1009e15cc1p-2,0x1.bb7cf2304bcfep-2,0x1.c1249d8011ee0p-2,0x1.c6c7f4997000ep-2,0x1.cc66e9931c45ep-2,0x1.d2016e8e9db57p-2,0x1.d79775b86e381p-2,0x1.dd28f1481cc5ap-2,0x1.e2b5d3806f63ap-2,0x1.e83e0eaf8510fp-2,0x1.edc1952ef78ccp-2,0x1.f3405963fd068p-2,0x1.f8ba4dbf89ab7p-2,0x1.fe2f64be71209p-2,0x1.01cfc874c3eb9p-1,0x1.0485626ae221ap-1,0x1.073879922ffecp-1,0x1.09e907417c5ddp-1,0x1.0c9704d5d8990p-1,0x1.0f426bb2a8e7dp-1,0x1.11eb3541b4b20p-1,0x1.14915af336ce7p-1,0x1.1734d63dedb49p-1,0x1.19d5a09f2b9b7p-1,0x1.1c73b39ae68c5p-1,0x1.1f0f08bbc861dp-1,0x1.21a799933eb58p-1,0x1.243d5fb98ac1dp-1,0x1.26d054cdd12dbp-1,0x1.2960727629ca9p-1,0x1.2bedb25faf3e9p-1,0x1.2e780e3e8ea14p-1,0x1.30ff7fce17031p-1,0x1.338400d0c8e57p-1,0x1.36058b10659f1p-1,0x1.3884185dfeb1fp-1,0x1.3affa292050bbp-1,0x1.3d78238c58343p-1,0x1.3fed9534556d2p-1,0x1.425ff178e6baep-1,0x1.44cf325091dd7p-1,0x1.473b51b987346p-1,0x1.49a449b9b0936p-1,0x1.4c0a145ec0000p-1,0x1.4e6cabbe3e5e9p-1,0x1.50cc09f59a09ap-1,0x1.5328292a35593p-1,0x1.5581038975139p-1,0x1.57d69348cec9fp-1,0x1.5a28d2a5d724ep-1,0x1.5c77bbe650189p-1,0x1.5ec3495837075p-1,0x1.610b7551d2cdep-1,0x1.63503a31c1be7p-1,0x1.6591925f07839p-1,0x1.67cf78491af10p-1,0x1.6a09e667f3bcbp-1,0x1.6c40d73c18272p-1,0x1.6e74454eaa8b0p-1,0x1.70a42b3176d79p-1,0x1.72d0837efff94p-1,0x1.74f948da8d28ap-1,0x1.771e75f037262p-1,0x1.79400574f55e4p-1,0x1.7b5df226aafadp-1,0x1.7d7836cc33daep-1,0x1.7f8ece3571771p-1,0x1.81a1b33b57acap-1,0x1.83b0e0bff976bp-1,0x1.85bc51ae958cdp-1,0x1.87c400fba2ebep-1,0x1.89c7e9a4dd4a9p-1,0x1.8bc806b15173ep-1,0x1.8dc45331698cdp-1,0x1.8fbcca3ef940cp-1,0x1.91b166fd49da0p-1,0x1.93a22499263f8p-1,0x1.958efe48e6dd7p-1,0x1.9777ef4c7d740p-1,0x1.995cf2ed80d20p-1,0x1.9b3e047f38742p-1,0x1.9d1b1f5ea80d5p-1,0x1.9ef43ef29af93p-1,0x1.a0c95eabaf934p-1,0x1.a29a7a0462782p-1,0x1.a4678c8119ac7p-1,0x1.a63091b02fae0p-1,0x1.a7f58529fe699p-1,0x1.a9b66290ea1a3p-1,0x1.ab7325916c0d3p-1,0x1.ad2bc9e21d50ep-1,0x1.aee04b43c1474p-1,0x1.b090a581501ffp-1,0x1.b23cd470013b3p-1,0x1.b3e4d3ef5570fp-1,0x1.b5889fe921405p-1,0x1.b728345196e3dp-1,0x1.b8c38d27504e7p-1,0x1.ba5aa673590cfp-1,0x1.bbed7c49380eap-1,0x1.bd7c0ac6f9529p-1,0x1.bf064e15377dbp-1,0x1.c08c42672554ap-1,0x1.c20de3fa971afp-1,0x1.c38b2f180bdafp-1,0x1.c5042012b6904p-1,0x1.c678b3488739bp-1,0x1.c7e8e52233cf2p-1,0x1.c954b213411f3p-1,0x1.cabc169a0b8fep-1,0x1.cc1f0f3fcfc5cp-1,0x1.cd7d9898b32f5p-1,0x1.ced7af43cc771p-1,0x1.d02d4feb2bd93p-1,0x1.d17e7743e35dcp-1,0x1.d2cb220e0ef9dp-1,0x1.d4134d14dc938p-1,0x1.d556f52e93eb1p-1,0x1.d696173c9e68ap-1,0x1.d7d0b02b8ecf8p-1,0x1.d906bcf328d44p-1,0x1.da383a9668988p-1,0x1.db6526238a09ap-1,0x1.dc8d7cb41025ep-1,0x1.ddb13b6ccc23dp-1,0x1.ded05f7de47d9p-1,0x1.dfeae622dbe2ap-1,0x1.e100cca2980aap-1,0x1.e212104f686e5p-1,0x1.e31eae870ce24p-1,0x1.e426a4b2bc17dp-1,0x1.e529f04729ffdp-1,0x1.e6288ec48e112p-1,0x1.e7227db6a9744p-1,0x1.e817bab4cd10bp-1,0x1.e9084361df7f3p-1,0x1.e9f4156c62ddap-1,0x1.eadb2e8e7a88dp-1,0x1.ebbd8c8df0b73p-1,0x1.ec9b2d3c3bf84p-1,0x1.ed740e7684963p-1,0x1.ee482e25a9dbbp-1,0x1.ef178a3e473c3p-1,0x1.efe220c0b95ecp-1,0x1.f0a7efb9230d7p-1,0x1.f168f53f7205cp-1,0x1.f2252f7763adap-1,0x1.f2dc9c9089a9dp-1,0x1.f38f3ac64e588p-1,0x1.f43d085ff92dbp-1,0x1.f4e603b0b2f2dp-1,0x1.f58a2b1789e84p-1,0x1.f6297cff75cafp-1,0x1.f6c3f7df5bbb8p-1,0x1.f7599a3a12077p-1,0x1.f7ea629e63d6ep-1,0x1.f8764fa714ba8p-1,0x1.f8fd5ffae41dbp-1,0x1.f97f924c9099ap-1,0x1.f9fce55adb2c8p-1,0x1.fa7557f08a516p-1,0x1.fae8e8e46cfbbp-1,0x1.fb5797195d741p-1,0x1.fbc1617e44185p-1,0x1.fc26470e19fd3p-1,0x1.fc8646cfeb721p-1,0x1.fce15fd6da67bp-1,0x1.fd37914220b84p-1,0x1.fd88da3d12526p-1,0x1.fdd539ff1f456p-1,0x1.fe1cafcbd5b09p-1,0x1.fe5f3af2e3940p-1,0x1.fe9cdad01883ap-1,0x1.fed58ecb673c4p-1,0x1.ff095658e71adp-1,0x1.ff3830f8d575cp-1,0x1.ff621e3796d7ep-1,0x1.ff871dadb81dfp-1,0x1.ffa72effef75cp-1,0x1.ffc251df1d3f8p-1,0x1.ffd886084cd0dp-1,0x1.ffe9cb44b51a1p-1,0x1.fff62169b92dbp-1,0x1.fffd8858e8a92p-1};
static const double ci[1024]={-0x0.0p+0,-0x1.921f0fe670071p-8,-0x1.921d1fcdec784p-7,-0x1.2d936bbe30efdp-6,-0x1.92155f7a3667ep-6,-0x1.f693731d1cf01p-6,-0x1.2d865759455cdp-5,-0x1.5fc00d290cd43p-5,-0x1.91f65f10dd814p-5,-0x1.c428d12c0d7e2p-5,-0x1.f656e79f820e0p-5,-0x1.1440134d709b2p-4,-0x1.2d52092ce19f6p-4,-0x1.4661179272096p-4,-0x1.5f6d00a9aa419p-4,-0x1.787586a5d5b21p-4,-0x1.917a6bc29b42cp-4,-0x1.aa7b724495c04p-4,-0x1.c3785c79ec2d5p-4,-0x1.dc70ecbae9fc8p-4,-0x1.f564e56a9730ep-4,-0x1.072a047ba831dp-3,-0x1.139f0cedaf576p-3,-0x1.20116d4ec7bcep-3,-0x1.2c8106e8e613ap-3,-0x1.38edbb0cd8d14p-3,-0x1.45576b1293e5ap-3,-0x1.51bdf8597c5f2p-3,-0x1.5e214448b3fc6p-3,-0x1.6a81304f64ab2p-3,-0x1.76dd9de50bf31p-3,-0x1.83366e89c64c5p-3,-0x1.8f8b83c69a60ap-3,-0x1.9bdcbf2dc4366p-3,-0x1.a82a025b00451p-3,-0x1.b4732ef3d6722p-3,-0x1.c0b826a7e4f63p-3,-0x1.ccf8cb312b286p-3,-0x1.d934fe5454311p-3,-0x1.e56ca1e101a1bp-3,-0x1.f19f97b215f1ap-3,-0x1.fdcdc1adfedf8p-3,-0x1.04fb80e37fdaep-2,-0x1.0b0d9cfdbdb90p-2,-0x1.111d262b1f677p-2,-0x1.172a0d7765177p-2,-0x1.1d3443f4cdb3dp-2,-0x1.233bbabc3bb72p-2,-0x1.294062ed59f05p-2,-0x1.2f422daec0386p-2,-0x1.35410c2e18152p-2,-0x1.3b3cefa0414b7p-2,-0x1.4135c94176602p-2,-0x1.472b8a5571054p-2,-0x1.4d1e24278e76ap-2,-0x1.530d880af3c24p-2,-0x1.58f9a75ab1fddp-2,-0x1.5ee27379ea693p-2,-0x1.64c7ddd3f27c6p-2,-0x1.6aa9d7dc77e16p-2,-0x1.7088530fa459ep-2,-0x1.766340f2418f6p-2,-0x1.7c3a9311dcce7p-2,-0x1.820e3b04eaac4p-2,-0x1.87de2a6aea963p-2,-0x1.8daa52ec8a4afp-2,-0x1.9372a63bc93d7p-2,-0x1.993716141bdfep-2,-0x1.9ef7943a8ed8ap-2,-0x1.a4b4127dea1e4p-2,-0x1.aa6c82b6d3fc9p-2,-0x1.b020d6c7f4009p-2,-0x1.b5d1009e15cc0p-2,-0x1.bb7cf2304bd01p-2,-0x1.c1249d8011ee7p-2,-0x1.c6c7f4997000ap-2,-0x1.cc66e9931c45dp-2,-0x1.d2016e8e9db5bp-2,-0x1.d79775b86e389p-2,-0x1.dd28f1481cc58p-2,-0x1.e2b5d3806f63bp-2,-0x1.e83e0eaf85113p-2,-0x1.edc1952ef78d5p-2,-0x1.f3405963fd068p-2,-0x1.f8ba4dbf89abap-2,-0x1.fe2f64be71210p-2,-0x1.01cfc874c3eb7p-1,-0x1.0485626ae221ap-1,-0x1.073879922ffedp-1,-0x1.09e907417c5e1p-1,-0x1.0c9704d5d898fp-1,-0x1.0f426bb2a8e7dp-1,-0x1.11eb3541b4b22p-1,-0x1.14915af336cebp-1,-0x1.1734d63dedb49p-1,-0x1.19d5a09f2b9b8p-1,-0x1.1c73b39ae68c8p-1,-0x1.1f0f08bbc861bp-1,-0x1.21a799933eb58p-1,-0x1.243d5fb98ac1fp-1,-0x1.26d054cdd12dfp-1,-0x1.2960727629ca8p-1,-0x1.2bedb25faf3eap-1,-0x1.2e780e3e8ea16p-1,-0x1.30ff7fce17035p-1,-0x1.338400d0c8e57p-1,-0x1.36058b10659f3p-1,-0x1.3884185dfeb22p-1,-0x1.3affa292050b9p-1,-0x1.3d78238c58343p-1,-0x1.3fed9534556d4p-1,-0x1.425ff178e6bb1p-1,-0x1.44cf325091dd6p-1,-0x1.473b51b987347p-1,-0x1.49a449b9b0938p-1,-0x1.4c0a145ec0004p-1,-0x1.4e6cabbe3e5e9p-1,-0x1.50cc09f59a09bp-1,-0x1.5328292a35596p-1,-0x1.5581038975137p-1,-0x1.57d69348cec9fp-1,-0x1.5a28d2a5d7250p-1,-0x1.5c77bbe65018cp-1,-0x1.5ec3495837074p-1,-0x1.610b7551d2cdep-1,-0x1.63503a31c1be9p-1,-0x1.6591925f0783ep-1,-0x1.67cf78491af10p-1,-0x1.6a09e667f3bccp-1,-0x1.6c40d73c18275p-1,-0x1.6e74454eaa8aep-1,-0x1.70a42b3176d7ap-1,-0x1.72d0837efff96p-1,-0x1.74f948da8d28dp-1,-0x1.771e75f037261p-1,-0x1.79400574f55e4p-1,-0x1.7b5df226aafafp-1,-0x1.7d7836cc33db2p-1,-0x1.7f8ece3571770p-1,-0x1.81a1b33b57acbp-1,-0x1.83b0e0bff976dp-1,-0x1.85bc51ae958ccp-1,-0x1.87c400fba2ebfp-1,-0x1.89c7e9a4dd4aap-1,-0x1.8bc806b151741p-1,-0x1.8dc45331698ccp-1,-0x1.8fbcca3ef940cp-1,-0x1.91b166fd49da2p-1,-0x1.93a22499263fbp-1,-0x1.958efe48e6dd7p-1,-0x1.9777ef4c7d741p-1,-0x1.995cf2ed80d22p-1,-0x1.9b3e047f38740p-1,-0x1.9d1b1f5ea80d5p-1,-0x1.9ef43ef29af94p-1,-0x1.a0c95eabaf936p-1,-0x1.a29a7a0462782p-1,-0x1.a4678c8119ac8p-1,-0x1.a63091b02fae2p-1,-0x1.a7f58529fe69dp-1,-0x1.a9b66290ea1a3p-1,-0x1.ab7325916c0d4p-1,-0x1.ad2bc9e21d510p-1,-0x1.aee04b43c1473p-1,-0x1.b090a581501ffp-1,-0x1.b23cd470013b3p-1,-0x1.b3e4d3ef55712p-1,-0x1.b5889fe921405p-1,-0x1.b728345196e3ep-1,-0x1.b8c38d27504e9p-1,-0x1.ba5aa673590d2p-1,-0x1.bbed7c49380eap-1,-0x1.bd7c0ac6f9529p-1,-0x1.bf064e15377ddp-1,-0x1.c08c426725549p-1,-0x1.c20de3fa971afp-1,-0x1.c38b2f180bdb0p-1,-0x1.c5042012b6907p-1,-0x1.c678b3488739bp-1,-0x1.c7e8e52233cf3p-1,-0x1.c954b213411f5p-1,-0x1.cabc169a0b900p-1,-0x1.cc1f0f3fcfc5cp-1,-0x1.cd7d9898b32f6p-1,-0x1.ced7af43cc773p-1,-0x1.d02d4feb2bd92p-1,-0x1.d17e7743e35dbp-1,-0x1.d2cb220e0ef9fp-1,-0x1.d4134d14dc93ap-1,-0x1.d556f52e93eb1p-1,-0x1.d696173c9e68bp-1,-0x1.d7d0b02b8ecf9p-1,-0x1.d906bcf328d46p-1,-0x1.da383a9668987p-1,-0x1.db6526238a09ap-1,-0x1.dc8d7cb410260p-1,-0x1.ddb13b6ccc23cp-1,-0x1.ded05f7de47d9p-1,-0x1.dfeae622dbe2bp-1,-0x1.e100cca2980acp-1,-0x1.e212104f686e5p-1,-0x1.e31eae870ce25p-1,-0x1.e426a4b2bc17ep-1,-0x1.e529f04729ffcp-1,-0x1.e6288ec48e112p-1,-0x1.e7227db6a9744p-1,-0x1.e817bab4cd10cp-1,-0x1.e9084361df7f2p-1,-0x1.e9f4156c62ddbp-1,-0x1.eadb2e8e7a88ep-1,-0x1.ebbd8c8df0b74p-1,-0x1.ec9b2d3c3bf84p-1,-0x1.ed740e7684963p-1,-0x1.ee482e25a9dbcp-1,-0x1.ef178a3e473c2p-1,-0x1.efe220c0b95ecp-1,-0x1.f0a7efb9230d7p-1,-0x1.f168f53f7205dp-1,-0x1.f2252f7763ad9p-1,-0x1.f2dc9c9089a9dp-1,-0x1.f38f3ac64e589p-1,-0x1.f43d085ff92ddp-1,-0x1.f4e603b0b2f2dp-1,-0x1.f58a2b1789e84p-1,-0x1.f6297cff75cb0p-1,-0x1.f6c3f7df5bbb7p-1,-0x1.f7599a3a12077p-1,-0x1.f7ea629e63d6ep-1,-0x1.f8764fa714ba9p-1,-0x1.f8fd5ffae41dbp-1,-0x1.f97f924c9099bp-1,-0x1.f9fce55adb2c8p-1,-0x1.fa7557f08a517p-1,-0x1.fae8e8e46cfbbp-1,-0x1.fb5797195d741p-1,-0x1.fbc1617e44186p-1,-0x1.fc26470e19fd3p-1,-0x1.fc8646cfeb721p-1,-0x1.fce15fd6da67bp-1,-0x1.fd37914220b84p-1,-0x1.fd88da3d12525p-1,-0x1.fdd539ff1f456p-1,-0x1.fe1cafcbd5b09p-1,-0x1.fe5f3af2e3940p-1,-0x1.fe9cdad01883ap-1,-0x1.fed58ecb673c4p-1,-0x1.ff095658e71adp-1,-0x1.ff3830f8d575cp-1,-0x1.ff621e3796d7ep-1,-0x1.ff871dadb81dfp-1,-0x1.ffa72effef75dp-1,-0x1.ffc251df1d3f8p-1,-0x1.ffd886084cd0dp-1,-0x1.ffe9cb44b51a1p-1,-0x1.fff62169b92dbp-1,-0x1.fffd8858e8a92p-1,-0x1.0000000000000p+0,-0x1.fffd8858e8a92p-1,-0x1.fff62169b92dbp-1,-0x1.ffe9cb44b51a1p-1,-0x1.ffd886084cd0dp-1,-0x1.ffc251df1d3f8p-1,-0x1.ffa72effef75dp-1,-0x1.ff871dadb81dfp-1,-0x1.ff621e3796d7ep-1,-0x1.ff3830f8d575cp-1,-0x1.ff095658e71adp-1,-0x1.fed58ecb673c4p-1,-0x1.fe9cdad01883ap-1,-0x1.fe5f3af2e3941p-1,-0x1.fe1cafcbd5b09p-1,-0x1.fdd539ff1f456p-1,-0x1.fd88da3d12526p-1,-0x1.fd37914220b84p-1,-0x1.fce15fd6da67bp-1,-0x1.fc8646cfeb721p-1,-0x1.fc26470e19fd3p-1,-0x1.fbc1617e44186p-1,-0x1.fb5797195d741p-1,-0x1.fae8e8e46cfbbp-1,-0x1.fa7557f08a517p-1,-0x1.f9fce55adb2c8p-1,-0x1.f97f924c9099bp-1,-0x1.f8fd5ffae41dbp-1,-0x1.f8764fa714ba9p-1,-0x1.f7ea629e63d6ep-1,-0x1.f7599a3a12077p-1,-0x1.f6c3f7df5bbb7p-1,-0x1.f6297cff75cb0p-1,-0x1.f58a2b1789e84p-1,-0x1.f4e603b0b2f2dp-1,-0x1.f43d085ff92ddp-1,-0x1.f38f3ac64e589p-1,-0x1.f2dc9c9089a9dp-1,-0x1.f2252f7763adap-1,-0x1.f168f53f7205dp-1,-0x1.f0a7efb9230d7p-1,-0x1.efe220c0b95edp-1,-0x1.ef178a3e473c2p-1,-0x1.ee482e25a9dbcp-1,-0x1.ed740e7684963p-1,-0x1.ec9b2d3c3bf84p-1,-0x1.ebbd8c8df0b75p-1,-0x1.eadb2e8e7a88ep-1,-0x1.e9f4156c62ddbp-1,-0x1.e9084361df7f2p-1,-0x1.e817bab4cd10dp-1,-0x1.e7227db6a9744p-1,-0x1.e6288ec48e112p-1,-0x1.e529f04729ffdp-1,-0x1.e426a4b2bc17fp-1,-0x1.e31eae870ce25p-1,-0x1.e212104f686e5p-1,-0x1.e100cca2980acp-1,-0x1.dfeae622dbe2bp-1,-0x1.ded05f7de47dap-1,-0x1.ddb13b6ccc23cp-1,-0x1.dc8d7cb410260p-1,-0x1.db6526238a09bp-1,-0x1.da383a9668988p-1,-0x1.d906bcf328d46p-1,-0x1.d7d0b02b8ecfap-1,-0x1.d696173c9e68bp-1,-0x1.d556f52e93eb1p-1,-0x1.d4134d14dc93ap-1,-0x1.d2cb220e0ef9fp-1,-0x1.d17e7743e35ddp-1,-0x1.d02d4feb2bd92p-1,-0x1.ced7af43cc774p-1,-0x1.cd7d9898b32f6p-1,-0x1.cc1f0f3fcfc5dp-1,-0x1.cabc169a0b901p-1,-0x1.c954b213411f4p-1,-0x1.c7e8e52233cf4p-1,-0x1.c678b3488739bp-1,-0x1.c5042012b6908p-1,-0x1.c38b2f180bdb1p-1,-0x1.c20de3fa971b0p-1,-0x1.c08c426725549p-1,-0x1.bf064e15377dep-1,-0x1.bd7c0ac6f952ap-1,-0x1.bbed7c49380ebp-1,-0x1.ba5aa673590d3p-1,-0x1.b8c38d27504e8p-1,-0x1.b728345196e3ep-1,-0x1.b5889fe921405p-1,-0x1.b3e4d3ef55712p-1,-0x1.b23cd470013b4p-1,-0x1.b090a58150201p-1,-0x1.aee04b43c1474p-1,-0x1.ad2bc9e21d512p-1,-0x1.ab7325916c0d5p-1,-0x1.a9b66290ea1a5p-1,-0x1.a7f58529fe69dp-1,-0x1.a63091b02fae1p-1,-0x1.a4678c8119ac9p-1,-0x1.a29a7a0462782p-1,-0x1.a0c95eabaf938p-1,-0x1.9ef43ef29af94p-1,-0x1.9d1b1f5ea80d7p-1,-0x1.9b3e047f38741p-1,-0x1.995cf2ed80d24p-1,-0x1.9777ef4c7d742p-1,-0x1.958efe48e6dd6p-1,-0x1.93a22499263fcp-1,-0x1.91b166fd49da1p-1,-0x1.8fbcca3ef940ep-1,-0x1.8dc45331698ccp-1,-0x1.8bc806b151742p-1,-0x1.89c7e9a4dd4aap-1,-0x1.87c400fba2ec0p-1,-0x1.85bc51ae958cdp-1,-0x1.83b0e0bff976fp-1,-0x1.81a1b33b57accp-1,-0x1.7f8ece3571770p-1,-0x1.7d7836cc33db3p-1,-0x1.7b5df226aafaep-1,-0x1.79400574f55e6p-1,-0x1.771e75f037261p-1,-0x1.74f948da8d28ep-1,-0x1.72d0837efff96p-1,-0x1.70a42b3176d7bp-1,-0x1.6e74454eaa8afp-1,-0x1.6c40d73c18277p-1,-0x1.6a09e667f3bcdp-1,-0x1.67cf78491af0fp-1,-0x1.6591925f0783ep-1,-0x1.63503a31c1be9p-1,-0x1.610b7551d2ce0p-1,-0x1.5ec3495837074p-1,-0x1.5c77bbe65018ep-1,-0x1.5a28d2a5d7250p-1,-0x1.57d69348ceca1p-1,-0x1.5581038975138p-1,-0x1.5328292a35598p-1,-0x1.50cc09f59a09cp-1,-0x1.4e6cabbe3e5e8p-1,-0x1.4c0a145ec0005p-1,-0x1.49a449b9b0938p-1,-0x1.473b51b987348p-1,-0x1.44cf325091dd6p-1,-0x1.425ff178e6bb3p-1,-0x1.3fed9534556d4p-1,-0x1.3d78238c58346p-1,-0x1.3affa292050bap-1,-0x1.3884185dfeb24p-1,-0x1.36058b10659f4p-1,-0x1.338400d0c8e56p-1,-0x1.30ff7fce17036p-1,-0x1.2e780e3e8ea16p-1,-0x1.2bedb25faf3ebp-1,-0x1.2960727629ca8p-1,-0x1.26d054cdd12e0p-1,-0x1.243d5fb98ac1fp-1,-0x1.21a799933eb5bp-1,-0x1.1f0f08bbc861cp-1,-0x1.1c73b39ae68c8p-1,-0x1.19d5a09f2b9b9p-1,-0x1.1734d63dedb48p-1,-0x1.14915af336cecp-1,-0x1.11eb3541b4b22p-1,-0x1.0f426bb2a8e7fp-1,-0x1.0c9704d5d898fp-1,-0x1.09e907417c5e2p-1,-0x1.073879922ffeep-1,-0x1.0485626ae221dp-1,-0x1.01cfc874c3eb8p-1,-0x1.fe2f64be7120ep-2,-0x1.f8ba4dbf89abcp-2,-0x1.f3405963fd066p-2,-0x1.edc1952ef78d8p-2,-0x1.e83e0eaf85113p-2,-0x1.e2b5d3806f63fp-2,-0x1.dd28f1481cc58p-2,-0x1.d79775b86e38dp-2,-0x1.d2016e8e9db5cp-2,-0x1.cc66e9931c463p-2,-0x1.c6c7f4997000cp-2,-0x1.c1249d8011ee5p-2,-0x1.bb7cf2304bd03p-2,-0x1.b5d1009e15cbfp-2,-0x1.b020d6c7f400cp-2,-0x1.aa6c82b6d3fc9p-2,-0x1.a4b4127dea1e8p-2,-0x1.9ef7943a8ed8bp-2,-0x1.993716141be03p-2,-0x1.9372a63bc93d8p-2,-0x1.8daa52ec8a4b5p-2,-0x1.87de2a6aea965p-2,-0x1.820e3b04eaac2p-2,-0x1.7c3a9311dcceap-2,-0x1.766340f2418f5p-2,-0x1.7088530fa45a2p-2,-0x1.6aa9d7dc77e17p-2,-0x1.64c7ddd3f27cap-2,-0x1.5ee27379ea694p-2,-0x1.58f9a75ab1fe2p-2,-0x1.530d880af3c25p-2,-0x1.4d1e24278e770p-2,-0x1.472b8a5571056p-2,-0x1.4135c94176600p-2,-0x1.3b3cefa0414bap-2,-0x1.35410c2e18151p-2,-0x1.2f422daec038ap-2,-0x1.294062ed59f06p-2,-0x1.233bbabc3bb76p-2,-0x1.1d3443f4cdb3fp-2,-0x1.172a0d776517cp-2,-0x1.111d262b1f679p-2,-0x1.0b0d9cfdbdb96p-2,-0x1.04fb80e37fdb0p-2,-0x1.fdcdc1adfedf6p-3,-0x1.f19f97b215f21p-3,-0x1.e56ca1e101a1ap-3,-0x1.d934fe5454319p-3,-0x1.ccf8cb312b286p-3,-0x1.c0b826a7e4f6cp-3,-0x1.b4732ef3d6724p-3,-0x1.a82a025b0045bp-3,-0x1.9bdcbf2dc436ap-3,-0x1.8f8b83c69a617p-3,-0x1.83366e89c64cbp-3,-0x1.76dd9de50bf2fp-3,-0x1.6a81304f64ab9p-3,-0x1.5e214448b3fc6p-3,-0x1.51bdf8597c5fap-3,-0x1.45576b1293e5bp-3,-0x1.38edbb0cd8d1dp-3,-0x1.2c8106e8e613cp-3,-0x1.20116d4ec7bdap-3,-0x1.139f0cedaf57ap-3,-0x1.072a047ba831ap-3,-0x1.f564e56a97319p-4,-0x1.dc70ecbae9fc5p-4,-0x1.c3785c79ec2e2p-4,-0x1.aa7b724495c03p-4,-0x1.917a6bc29b43cp-4,-0x1.787586a5d5b23p-4,-0x1.5f6d00a9aa42cp-4,-0x1.466117927209bp-4,-0x1.2d52092ce1a0cp-4,-0x1.1440134d709bbp-4,-0x1.f656e79f820d3p-5,-0x1.c428d12c0d7f9p-5,-0x1.91f65f10dd80dp-5,-0x1.5fc00d290cd60p-5,-0x1.2d865759455cdp-5,-0x1.f693731d1cf46p-6,-0x1.92155f7a36689p-6,-0x1.2d936bbe30f4ep-6,-0x1.921d1fcdec7b3p-7,-0x1.921f0fe6701e6p-8,-0x1.1a62633145c07p-53,0x1.921f0fe6700ccp-8,0x1.921d1fcdec726p-7,0x1.2d936bbe30f08p-6,0x1.92155f7a36642p-6,0x1.f693731d1ceffp-6,0x1.2d865759455a9p-5,0x1.5fc00d290cd3dp-5,0x1.91f65f10dd7eap-5,0x1.c428d12c0d7d6p-5,0x1.f656e79f820b0p-5,0x1.1440134d709a9p-4,0x1.2d52092ce19fbp-4,0x1.4661179272089p-4,0x1.5f6d00a9aa41bp-4,0x1.787586a5d5b12p-4,0x1.917a6bc29b42bp-4,0x1.aa7b724495bf1p-4,0x1.c3785c79ec2d1p-4,0x1.dc70ecbae9fb4p-4,0x1.f564e56a97307p-4,0x1.072a047ba8311p-3,0x1.139f0cedaf572p-3,0x1.20116d4ec7bd1p-3,0x1.2c8106e8e6134p-3,0x1.38edbb0cd8d14p-3,0x1.45576b1293e52p-3,0x1.51bdf8597c5f1p-3,0x1.5e214448b3fbdp-3,0x1.6a81304f64ab0p-3,0x1.76dd9de50bf26p-3,0x1.83366e89c64c2p-3,0x1.8f8b83c69a60ep-3,0x1.9bdcbf2dc4361p-3,0x1.a82a025b00452p-3,0x1.b4732ef3d671cp-3,0x1.c0b826a7e4f63p-3,0x1.ccf8cb312b27ep-3,0x1.d934fe5454310p-3,0x1.e56ca1e101a12p-3,0x1.f19f97b215f18p-3,0x1.fdcdc1adfedeep-3,0x1.04fb80e37fdacp-2,0x1.0b0d9cfdbdb92p-2,0x1.111d262b1f675p-2,0x1.172a0d7765178p-2,0x1.1d3443f4cdb3ap-2,0x1.233bbabc3bb71p-2,0x1.294062ed59f01p-2,0x1.2f422daec0386p-2,0x1.35410c2e1814dp-2,0x1.3b3cefa0414b6p-2,0x1.4135c941765fcp-2,0x1.472b8a5571052p-2,0x1.4d1e24278e76bp-2,0x1.530d880af3c21p-2,0x1.58f9a75ab1fdep-2,0x1.5ee27379ea690p-2,0x1.64c7ddd3f27c6p-2,0x1.6aa9d7dc77e12p-2,0x1.7088530fa459ep-2,0x1.766340f2418f1p-2,0x1.7c3a9311dcce5p-2,0x1.820e3b04eaabep-2,0x1.87de2a6aea961p-2,0x1.8daa52ec8a4b1p-2,0x1.9372a63bc93d4p-2,0x1.993716141bdffp-2,0x1.9ef7943a8ed87p-2,0x1.a4b4127dea1e4p-2,0x1.aa6c82b6d3fc5p-2,0x1.b020d6c7f4008p-2,0x1.b5d1009e15cbbp-2,0x1.bb7cf2304bcffp-2,0x1.c1249d8011ee1p-2,0x1.c6c7f49970008p-2,0x1.cc66e9931c45fp-2,0x1.d2016e8e9db58p-2,0x1.d79775b86e38ap-2,0x1.dd28f1481cc54p-2,0x1.e2b5d3806f63bp-2,0x1.e83e0eaf8510fp-2,0x1.edc1952ef78d4p-2,0x1.f3405963fd062p-2,0x1.f8ba4dbf89ab8p-2,0x1.fe2f64be7120ap-2,0x1.01cfc874c3eb6p-1,0x1.0485626ae221bp-1,0x1.073879922ffecp-1,0x1.09e907417c5e1p-1,0x1.0c9704d5d898dp-1,0x1.0f426bb2a8e7dp-1,0x1.11eb3541b4b20p-1,0x1.14915af336cebp-1,0x1.1734d63dedb46p-1,0x1.19d5a09f2b9b7p-1,0x1.1c73b39ae68c6p-1,0x1.1f0f08bbc861ap-1,0x1.21a799933eb59p-1,0x1.243d5fb98ac1ep-1,0x1.26d054cdd12dfp-1,0x1.2960727629ca6p-1,0x1.2bedb25faf3eap-1,0x1.2e780e3e8ea14p-1,0x1.30ff7fce17034p-1,0x1.338400d0c8e54p-1,0x1.36058b10659f2p-1,0x1.3884185dfeb22p-1,0x1.3affa292050b8p-1,0x1.3d78238c58344p-1,0x1.3fed9534556d3p-1,0x1.425ff178e6bb1p-1,0x1.44cf325091dd4p-1,0x1.473b51b987347p-1,0x1.49a449b9b0937p-1,0x1.4c0a145ec0003p-1,0x1.4e6cabbe3e5e7p-1,0x1.50cc09f59a09ap-1,0x1.5328292a35596p-1,0x1.5581038975136p-1,0x1.57d69348ceca0p-1,0x1.5a28d2a5d724fp-1,0x1.5c77bbe65018cp-1,0x1.5ec3495837072p-1,0x1.610b7551d2cdep-1,0x1.63503a31c1be7p-1,0x1.6591925f0783dp-1,0x1.67cf78491af0ep-1,0x1.6a09e667f3bccp-1,0x1.6c40d73c18275p-1,0x1.6e74454eaa8adp-1,0x1.70a42b3176d7ap-1,0x1.72d0837efff95p-1,0x1.74f948da8d28dp-1,0x1.771e75f03725fp-1,0x1.79400574f55e4p-1,0x1.7b5df226aafadp-1,0x1.7d7836cc33db1p-1,0x1.7f8ece357176ep-1,0x1.81a1b33b57acbp-1,0x1.83b0e0bff976bp-1,0x1.85bc51ae958cep-1,0x1.87c400fba2ebfp-1,0x1.89c7e9a4dd4a9p-1,0x1.8bc806b15173ep-1,0x1.8dc45331698cdp-1,0x1.8fbcca3ef940cp-1,0x1.91b166fd49da0p-1,0x1.93a22499263f8p-1,0x1.958efe48e6dd7p-1,0x1.9777ef4c7d741p-1,0x1.995cf2ed80d20p-1,0x1.9b3e047f38742p-1,0x1.9d1b1f5ea80d5p-1,0x1.9ef43ef29af93p-1,0x1.a0c95eabaf934p-1,0x1.a29a7a0462783p-1,0x1.a4678c8119ac8p-1,0x1.a63091b02fae0p-1,0x1.a7f58529fe69ap-1,0x1.a9b66290ea1a3p-1,0x1.ab7325916c0d4p-1,0x1.ad2bc9e21d50fp-1,0x1.aee04b43c1475p-1,0x1.b090a581501ffp-1,0x1.b23cd470013b3p-1,0x1.b3e4d3ef5570fp-1,0x1.b5889fe921406p-1,0x1.b728345196e3dp-1,0x1.b8c38d27504e7p-1,0x1.ba5aa673590d0p-1,0x1.bbed7c49380eap-1,0x1.bd7c0ac6f9529p-1,0x1.bf064e15377dbp-1,0x1.c08c42672554ap-1,0x1.c20de3fa971afp-1,0x1.c38b2f180bdb0p-1,0x1.c5042012b6905p-1,0x1.c678b3488739cp-1,0x1.c7e8e52233cf3p-1,0x1.c954b213411f4p-1,0x1.cabc169a0b8fep-1,0x1.cc1f0f3fcfc5cp-1,0x1.cd7d9898b32f5p-1,0x1.ced7af43cc771p-1,0x1.d02d4feb2bd93p-1,0x1.d17e7743e35dcp-1,0x1.d2cb220e0ef9ep-1,0x1.d4134d14dc938p-1,0x1.d556f52e93eb1p-1,0x1.d696173c9e68ap-1,0x1.d7d0b02b8ecf8p-1,0x1.d906bcf328d44p-1,0x1.da383a9668988p-1,0x1.db6526238a09ap-1,0x1.dc8d7cb41025fp-1,0x1.ddb13b6ccc23dp-1,0x1.ded05f7de47dap-1,0x1.dfeae622dbe2ap-1,0x1.e100cca2980aap-1,0x1.e212104f686e5p-1,0x1.e31eae870ce24p-1,0x1.e426a4b2bc17dp-1,0x1.e529f04729ffdp-1,0x1.e6288ec48e112p-1,0x1.e7227db6a9744p-1,0x1.e817bab4cd10cp-1,0x1.e9084361df7f3p-1,0x1.e9f4156c62ddap-1,0x1.eadb2e8e7a88dp-1,0x1.ebbd8c8df0b73p-1,0x1.ec9b2d3c3bf84p-1,0x1.ed740e7684963p-1,0x1.ee482e25a9dbbp-1,0x1.ef178a3e473c3p-1,0x1.efe220c0b95ecp-1,0x1.f0a7efb9230d7p-1,0x1.f168f53f7205cp-1,0x1.f2252f7763adap-1,0x1.f2dc9c9089a9dp-1,0x1.f38f3ac64e588p-1,0x1.f43d085ff92dcp-1,0x1.f4e603b0b2f2dp-1,0x1.f58a2b1789e84p-1,0x1.f6297cff75cafp-1,0x1.f6c3f7df5bbb8p-1,0x1.f7599a3a12077p-1,0x1.f7ea629e63d6ep-1,0x1.f8764fa714ba8p-1,0x1.f8fd5ffae41dbp-1,0x1.f97f924c9099ap-1,0x1.f9fce55adb2c8p-1,0x1.fa7557f08a516p-1,0x1.fae8e8e46cfbbp-1,0x1.fb5797195d741p-1,0x1.fbc1617e44186p-1,0x1.fc26470e19fd4p-1,0x1.fc8646cfeb721p-1,0x1.fce15fd6da67bp-1,0x1.fd37914220b84p-1,0x1.fd88da3d12526p-1,0x1.fdd539ff1f456p-1,0x1.fe1cafcbd5b09p-1,0x1.fe5f3af2e3940p-1,0x1.fe9cdad01883ap-1,0x1.fed58ecb673c4p-1,0x1.ff095658e71adp-1,0x1.ff3830f8d575cp-1,0x1.ff621e3796d7ep-1,0x1.ff871dadb81dfp-1,0x1.ffa72effef75cp-1,0x1.ffc251df1d3f8p-1,0x1.ffd886084cd0dp-1,0x1.ffe9cb44b51a1p-1,0x1.fff62169b92dbp-1,0x1.fffd8858e8a92p-1,0x1.0000000000000p+0,0x1.fffd8858e8a92p-1,0x1.fff62169b92dbp-1,0x1.ffe9cb44b51a1p-1,0x1.ffd886084cd0dp-1,0x1.ffc251df1d3f8p-1,0x1.ffa72effef75dp-1,0x1.ff871dadb81dfp-1,0x1.ff621e3796d7ep-1,0x1.ff3830f8d575cp-1,0x1.ff095658e71adp-1,0x1.fed58ecb673c4p-1,0x1.fe9cdad01883ap-1,0x1.fe5f3af2e3940p-1,0x1.fe1cafcbd5b09p-1,0x1.fdd539ff1f456p-1,0x1.fd88da3d12526p-1,0x1.fd37914220b84p-1,0x1.fce15fd6da67bp-1,0x1.fc8646cfeb722p-1,0x1.fc26470e19fd4p-1,0x1.fbc1617e44186p-1,0x1.fb5797195d741p-1,0x1.fae8e8e46cfbbp-1,0x1.fa7557f08a516p-1,0x1.f9fce55adb2c8p-1,0x1.f97f924c9099bp-1,0x1.f8fd5ffae41dcp-1,0x1.f8764fa714ba9p-1,0x1.f7ea629e63d6fp-1,0x1.f7599a3a12078p-1,0x1.f6c3f7df5bbb8p-1,0x1.f6297cff75cb0p-1,0x1.f58a2b1789e85p-1,0x1.f4e603b0b2f2ep-1,0x1.f43d085ff92dcp-1,0x1.f38f3ac64e589p-1,0x1.f2dc9c9089a9ep-1,0x1.f2252f7763adbp-1,0x1.f168f53f7205dp-1,0x1.f0a7efb9230d7p-1,0x1.efe220c0b95edp-1,0x1.ef178a3e473c4p-1,0x1.ee482e25a9dbcp-1,0x1.ed740e7684964p-1,0x1.ec9b2d3c3bf85p-1,0x1.ebbd8c8df0b74p-1,0x1.eadb2e8e7a88ep-1,0x1.e9f4156c62ddbp-1,0x1.e9084361df7f4p-1,0x1.e817bab4cd10dp-1,0x1.e7227db6a9745p-1,0x1.e6288ec48e113p-1,0x1.e529f04729ffep-1,0x1.e426a4b2bc17ep-1,0x1.e31eae870ce25p-1,0x1.e212104f686e6p-1,0x1.e100cca2980abp-1,0x1.dfeae622dbe2bp-1,0x1.ded05f7de47dbp-1,0x1.ddb13b6ccc23ep-1,0x1.dc8d7cb410260p-1,0x1.db6526238a09bp-1,0x1.da383a9668989p-1,0x1.d906bcf328d45p-1,0x1.d7d0b02b8ecf9p-1,0x1.d696173c9e68cp-1,0x1.d556f52e93eb2p-1,0x1.d4134d14dc939p-1,0x1.d2cb220e0ef9fp-1,0x1.d17e7743e35ddp-1,0x1.d02d4feb2bd94p-1,0x1.ced7af43cc773p-1,0x1.cd7d9898b32f7p-1,0x1.cc1f0f3fcfc5ep-1,0x1.cabc169a0b8ffp-1,0x1.c954b213411f5p-1,0x1.c7e8e52233cf4p-1,0x1.c678b3488739dp-1,0x1.c5042012b6906p-1,0x1.c38b2f180bdb1p-1,0x1.c20de3fa971b1p-1,0x1.c08c42672554cp-1,0x1.bf064e15377ddp-1,0x1.bd7c0ac6f952ap-1,0x1.bbed7c49380ecp-1,0x1.ba5aa673590d1p-1,0x1.b8c38d27504e9p-1,0x1.b728345196e3fp-1,0x1.b5889fe921407p-1,0x1.b3e4d3ef55711p-1,0x1.b23cd470013b5p-1,0x1.b090a58150201p-1,0x1.aee04b43c1477p-1,0x1.ad2bc9e21d510p-1,0x1.ab7325916c0d5p-1,0x1.a9b66290ea1a5p-1,0x1.a7f58529fe69cp-1,0x1.a63091b02fae2p-1,0x1.a4678c8119acap-1,0x1.a29a7a0462785p-1,0x1.a0c95eabaf936p-1,0x1.9ef43ef29af95p-1,0x1.9d1b1f5ea80d7p-1,0x1.9b3e047f38744p-1,0x1.995cf2ed80d22p-1,0x1.9777ef4c7d743p-1,0x1.958efe48e6dd9p-1,0x1.93a22499263fap-1,0x1.91b166fd49da2p-1,0x1.8fbcca3ef940ep-1,0x1.8dc45331698cfp-1,0x1.8bc806b151740p-1,0x1.89c7e9a4dd4abp-1,0x1.87c400fba2ec1p-1,0x1.85bc51ae958d0p-1,0x1.83b0e0bff976ep-1,0x1.81a1b33b57acdp-1,0x1.7f8ece3571773p-1,0x1.7d7836cc33db1p-1,0x1.7b5df226aafafp-1,0x1.79400574f55e6p-1,0x1.771e75f037264p-1,0x1.74f948da8d28cp-1,0x1.72d0837efff97p-1,0x1.70a42b3176d7cp-1,0x1.6e74454eaa8b2p-1,0x1.6c40d73c18275p-1,0x1.6a09e667f3bcep-1,0x1.67cf78491af13p-1,0x1.6591925f0783cp-1,0x1.63503a31c1be9p-1,0x1.610b7551d2ce1p-1,0x1.5ec3495837078p-1,0x1.5c77bbe65018cp-1,0x1.5a28d2a5d7251p-1,0x1.57d69348ceca2p-1,0x1.558103897513cp-1,0x1.5328292a35596p-1,0x1.50cc09f59a09dp-1,0x1.4e6cabbe3e5ecp-1,0x1.4c0a145ec0003p-1,0x1.49a449b9b0939p-1,0x1.473b51b987349p-1,0x1.44cf325091ddap-1,0x1.425ff178e6bb1p-1,0x1.3fed9534556d5p-1,0x1.3d78238c58347p-1,0x1.3affa292050bep-1,0x1.3884185dfeb22p-1,0x1.36058b10659f4p-1,0x1.338400d0c8e5ap-1,0x1.30ff7fce17034p-1,0x1.2e780e3e8ea17p-1,0x1.2bedb25faf3ecp-1,0x1.2960727629cacp-1,0x1.26d054cdd12dep-1,0x1.243d5fb98ac20p-1,0x1.21a799933eb5cp-1,0x1.1f0f08bbc8620p-1,0x1.1c73b39ae68c8p-1,0x1.19d5a09f2b9bap-1,0x1.1734d63dedb4cp-1,0x1.14915af336ceap-1,0x1.11eb3541b4b23p-1,0x1.0f426bb2a8e80p-1,0x1.0c9704d5d8993p-1,0x1.09e907417c5e0p-1,0x1.073879922ffefp-1,0x1.0485626ae221ep-1,0x1.01cfc874c3ebcp-1,0x1.fe2f64be71210p-2,0x1.f8ba4dbf89abep-2,0x1.f3405963fd06fp-2,0x1.edc1952ef78d3p-2,0x1.e83e0eaf85115p-2,0x1.e2b5d3806f640p-2,0x1.dd28f1481cc61p-2,0x1.d79775b86e388p-2,0x1.d2016e8e9db5ep-2,0x1.cc66e9931c465p-2,0x1.c6c7f49970015p-2,0x1.c1249d8011ee7p-2,0x1.bb7cf2304bd05p-2,0x1.b5d1009e15cc8p-2,0x1.b020d6c7f4007p-2,0x1.aa6c82b6d3fcbp-2,0x1.a4b4127dea1eap-2,0x1.9ef7943a8ed94p-2,0x1.993716141bdfep-2,0x1.9372a63bc93dap-2,0x1.8daa52ec8a4b7p-2,0x1.87de2a6aea96ep-2,0x1.820e3b04eaac4p-2,0x1.7c3a9311dccecp-2,0x1.766340f2418ffp-2,0x1.7088530fa459cp-2,0x1.6aa9d7dc77e19p-2,0x1.64c7ddd3f27ccp-2,0x1.5ee27379ea69dp-2,0x1.58f9a75ab1fdcp-2,0x1.530d880af3c27p-2,0x1.4d1e24278e772p-2,0x1.472b8a5571050p-2,0x1.4135c94176602p-2,0x1.3b3cefa0414bcp-2,0x1.35410c2e1815bp-2,0x1.2f422daec0385p-2,0x1.294062ed59f08p-2,0x1.233bbabc3bb78p-2,0x1.1d3443f4cdb48p-2,0x1.172a0d7765177p-2,0x1.111d262b1f67bp-2,0x1.0b0d9cfdbdb98p-2,0x1.04fb80e37fdabp-2,0x1.fdcdc1adfedfap-3,0x1.f19f97b215f25p-3,0x1.e56ca1e101a2ep-3,0x1.d934fe545430dp-3,0x1.ccf8cb312b28bp-3,0x1.c0b826a7e4f70p-3,0x1.b4732ef3d6738p-3,0x1.a82a025b00450p-3,0x1.9bdcbf2dc436ep-3,0x1.8f8b83c69a61bp-3,0x1.83366e89c64bfp-3,0x1.76dd9de50bf33p-3,0x1.6a81304f64abdp-3,0x1.5e214448b3fdap-3,0x1.51bdf8597c5eep-3,0x1.45576b1293e5fp-3,0x1.38edbb0cd8d21p-3,0x1.2c8106e8e6151p-3,0x1.20116d4ec7bcep-3,0x1.139f0cedaf57fp-3,0x1.072a047ba832ep-3,0x1.f564e56a97302p-4,0x1.dc70ecbae9fcep-4,0x1.c3785c79ec2ebp-4,0x1.aa7b724495c2bp-4,0x1.917a6bc29b425p-4,0x1.787586a5d5b2cp-4,0x1.5f6d00a9aa435p-4,0x1.46611792720c3p-4,0x1.2d52092ce19f5p-4,0x1.1440134d709c3p-4,0x1.f656e79f82125p-5,0x1.c428d12c0d7cbp-5,0x1.91f65f10dd81fp-5,0x1.5fc00d290cd71p-5,0x1.2d8657594561ep-5,0x1.f693731d1cee9p-6,0x1.92155f7a366acp-6,0x1.2d936bbe30f72p-6,0x1.921d1fcdec8f9p-7,0x1.921f0fe670073p-8};
for(int q=0;q<NR;q++){DQr[q]=cr[q];DQi[q]=ci[q];}
for (int t=0;t<NC;t++){ f64 a=-2.0*M_PI*(f64)t/(f64)NN; DTr[t]=cos(a); DTi[t]=sin(a); }
for (int t=0;t<NC;t++){ f64 a=-2.0*M_PI*(f64)t/(f64)NN2; DT18[t]=cos(a); DT18i[t]=sin(a); }
build_r4();
for (int i=0;i<NR;i++){ u32 r=0; for(int b=0;b<10;b++) if(i&(1u<<b)) r|=1u<<(9-b); BREV[i]=r; }
for (int i=0;i<512;i++){ u32 r=0; for(int b=0;b<9;b++) if(i&(1u<<b)) r|=1u<<(8-b); BREV9[i]=r; }
for (int q=0;q<512;q++){ DQ18[q]=DQr[2*q]; DQ18i[q]=DQi[2*q]; }
for (int j=0;j<NC;j++){ u32 r=0; for(int b=0;b<9;b++) if(j&(1u<<b)) r|=1u<<(8-b); BREVC[j]=r; }
for (int i=0;i<NR;i++){ u32 k1=BREV[i]; u32 k1n=(NR-k1)%NR; RVIX[i]=BREV[k1n]; }
for (int j=0;j<NC;j++){
u32 k2=BREVC[j];
CVIX0[j]=BREVC[(NC-k2)%NC];
CVIX1[j]=BREVC[(NC-k2+NC-1)%NC];
}
for (int j=0;j<NC;j+=2){
u32 k2 = BREVC[j];
for (int par=0; par<2; par++) {
u32 k2p = par + 2*k2;
u32 cc = BREVC[k2p & 511];
if (par) CJ1[j]=cc; else CJ0[j]=cc;
}
}
{ int n=NC, k=0;
for (int h=n>>1;h>=1;h>>=1){ for(int j=0;j<h;j++){ int q=512*j/h; RDr[k+j]=DQr[q]; RDi[k+j]=DQi[q];} k+=h; }
k=0;
for(int h=1;h<n;h<<=1){ for(int j=0;j<h;j++){ RIr[k+j]=RDr[k+j]; RIi[k+j]=-RDi[k+j]; } k+=h; } }
for (int kq=0;kq<2;kq++)
for (int c=0;c<(NC>>1);c++){ u32 eh = kq + 2u*BREVC[2*c]; TWR[kq][c]=DQr[eh]; TWI[kq][c]=DQi[eh]; }
{ int n=NR, k=0;
for (int h=n>>1;h>=1;h>>=1){ for(int j=0;j<h;j++){ int q=512*j/h; CDr[k+j]=DQr[q]; CDi[k+j]=DQi[q];} k+=h; }
k=0;
for(int h=1;h<n;h<<=1){ for(int j=0;j<h;j++){ CIr[k+j]=CDr[k+j]; CIi[k+j]=-CDi[k+j]; } k+=h; } }
}
FT static inline void low3_dif(f64*pr, f64*pi) {
__m256d r0=_mm256_load_pd(pr), r1=_mm256_load_pd(pr+4);
__m256d i0=_mm256_load_pd(pi), i1=_mm256_load_pd(pi+4);
{ const f64 twr[4]={1.0, 0.7071067811865475244, 0.0, -0.7071067811865475244};
const f64 twi[4]={0.0,-0.7071067811865475244,-1.0, -0.7071067811865475244};
__m256d wr=_mm256_load_pd(twr), wi=_mm256_load_pd(twi);
__m256d sr=_mm256_add_pd(r0,r1), si=_mm256_add_pd(i0,i1);
__m256d dr=_mm256_sub_pd(r0,r1), di=_mm256_sub_pd(i0,i1);
__m256d t2=_mm256_mul_pd(di,wi), drr=_mm256_fmsub_pd(dr,wr,t2);
__m256d t3=_mm256_mul_pd(dr,wi), dii=_mm256_fmadd_pd(di,wr,t3);
r0=sr;i0=si;r1=drr;i1=dii; }
#define H2D(V,W) do{ \
__m256d Av=_mm256_permute2f128_pd(V,V,1), Aw=_mm256_permute2f128_pd(W,W,1); \
__m256d P=_mm256_add_pd(V,Av), Q=_mm256_sub_pd(V,Av), Q2=_mm256_sub_pd(Av,V); \
__m256d R2=_mm256_add_pd(W,Aw), S3=_mm256_sub_pd(Aw,W); \
__m256d X=_mm256_blend_pd(P,Q2,0xC), Y=_mm256_blend_pd(R2,S3,0xC); \
V=_mm256_blend_pd(X,Y,0x8); \
W=_mm256_blend_pd(Y,Q,0x8); \
}while(0)
H2D(r0,i0); H2D(r1,i1);
#undef H2D
#define H1D(V,W) do{ \
V=_mm256_fmsubadd_pd(_mm256_permute_pd(V,5),_mm256_set1_pd(1.0),V); \
W=_mm256_fmsubadd_pd(_mm256_permute_pd(W,5),_mm256_set1_pd(1.0),W); \
}while(0)
H1D(r0,i0); H1D(r1,i1);
#undef H1D
_mm256_store_pd(pr,r0); _mm256_store_pd(pi,i0);
_mm256_store_pd(pr+4,r1); _mm256_store_pd(pi+4,i1);
}
FT static inline void low3_dit(f64*pr, f64*pi) {
__m256d r0=_mm256_load_pd(pr), r1=_mm256_load_pd(pr+4);
__m256d i0=_mm256_load_pd(pi), i1=_mm256_load_pd(pi+4);
#define H1I(V,W) do{ \
V=_mm256_fmsubadd_pd(_mm256_permute_pd(V,5),_mm256_set1_pd(1.0),V); \
W=_mm256_fmsubadd_pd(_mm256_permute_pd(W,5),_mm256_set1_pd(1.0),W); \
}while(0)
H1I(r0,i0); H1I(r1,i1);
#undef H1I
#define H2I(V,W) do{ \
__m256d nW=_mm256_sub_pd(_mm256_setzero_pd(),W); \
__m256d P=_mm256_blend_pd(V,nW,0x8), R=_mm256_blend_pd(W,V,0x8); \
__m256d Ps=_mm256_permute2f128_pd(P,P,1), Rs=_mm256_permute2f128_pd(R,R,1); \
V=_mm256_blend_pd(_mm256_add_pd(P,Ps),_mm256_sub_pd(Ps,P),0xC); \
W=_mm256_blend_pd(_mm256_add_pd(R,Rs),_mm256_sub_pd(Rs,R),0xC); \
}while(0)
H2I(r0,i0); H2I(r1,i1);
#undef H2I
{ const f64 twr[4]={1.0, 0.7071067811865475244, 0.0, -0.7071067811865475244};
const f64 twi[4]={0.0, 0.7071067811865475244, 1.0, 0.7071067811865475244};
__m256d wr=_mm256_load_pd(twr), wi=_mm256_load_pd(twi);
__m256d t2=_mm256_mul_pd(i1,wi), prr=_mm256_fmsub_pd(r1,wr,t2);
__m256d t3=_mm256_mul_pd(r1,wi), pii=_mm256_fmadd_pd(i1,wr,t3);
__m256d s=_mm256_add_pd(r0,prr), d=_mm256_sub_pd(r0,prr);
__m256d s2=_mm256_add_pd(i0,pii), d2=_mm256_sub_pd(i0,pii);
r0=s;i0=s2;r1=d;i1=d2; }
_mm256_store_pd(pr,r0); _mm256_store_pd(pi,i0);
_mm256_store_pd(pr+4,r1); _mm256_store_pd(pi+4,i1);
}
#define R4BODY(R0,R1,R2,R3,I0,I1,I2,I3,W1R,W1I,W2R,W2I,W3R,W3I) do{ \
__m256d x0r=_mm256_load_pd(R0+c),x0i=_mm256_load_pd(I0+c); \
__m256d x1r=_mm256_load_pd(R1+c),x1i=_mm256_load_pd(I1+c); \
__m256d x2r=_mm256_load_pd(R2+c),x2i=_mm256_load_pd(I2+c); \
__m256d x3r=_mm256_load_pd(R3+c),x3i=_mm256_load_pd(I3+c); \
__m256d ar=_mm256_add_pd(x0r,x2r), ai=_mm256_add_pd(x0i,x2i); \
__m256d br=_mm256_sub_pd(x0r,x2r), bi=_mm256_sub_pd(x0i,x2i); \
__m256d cr=_mm256_add_pd(x1r,x3r), ci=_mm256_add_pd(x1i,x3i); \
__m256d dr=_mm256_sub_pd(x1r,x3r), di=_mm256_sub_pd(x1i,x3i); \
_mm256_store_pd(R0+c,_mm256_add_pd(ar,cr)); \
_mm256_store_pd(I0+c,_mm256_add_pd(ai,ci)); \
{ __m256d y2r=_mm256_sub_pd(ar,cr), y2i=_mm256_sub_pd(ai,ci); \
__m256d u=_mm256_mul_pd(y2i,W2I), z2r=_mm256_fmsub_pd(y2r,W2R,u); \
__m256d v=_mm256_mul_pd(y2r,W2I), z2i=_mm256_fmadd_pd(y2i,W2R,v); \
_mm256_store_pd(R1+c,z2r); _mm256_store_pd(I1+c,z2i); } \
{ __m256d u1r=_mm256_add_pd(br,di), u1i=_mm256_sub_pd(bi,dr); \
__m256d u=_mm256_mul_pd(u1i,W1I), y1r=_mm256_fmsub_pd(u1r,W1R,u); \
__m256d v=_mm256_mul_pd(u1r,W1I), y1i=_mm256_fmadd_pd(u1i,W1R,v); \
_mm256_store_pd(R2+c,y1r); _mm256_store_pd(I2+c,y1i); } \
{ __m256d u3r=_mm256_sub_pd(br,di), u3i=_mm256_add_pd(bi,dr); \
__m256d u=_mm256_mul_pd(u3i,W3I), y3r=_mm256_fmsub_pd(u3r,W3R,u); \
__m256d v=_mm256_mul_pd(u3r,W3I), y3i=_mm256_fmadd_pd(u3i,W3R,v); \
_mm256_store_pd(R3+c,y3r); _mm256_store_pd(I3+c,y3i); } }while(0)
#define R4BODY0(R0,R1,R2,R3,I0,I1,I2,I3,W1R,W1I,W2R,W2I,W3R,W3I) do{ \
__m256d x0r=_mm256_load_pd(R0+c),x0i=_mm256_load_pd(I0+c); \
__m256d x1r=_mm256_load_pd(R1+c),x1i=_mm256_load_pd(I1+c); \
__m256d ar=x0r, ai=x0i, br=x0r, bi=x0i; \
__m256d cr=x1r, ci=x1i, dr=x1r, di=x1i; \
_mm256_store_pd(R0+c,_mm256_add_pd(ar,cr)); \
_mm256_store_pd(I0+c,_mm256_add_pd(ai,ci)); \
{ __m256d y2r=_mm256_sub_pd(ar,cr), y2i=_mm256_sub_pd(ai,ci); \
__m256d u=_mm256_mul_pd(y2i,W2I), z2r=_mm256_fmsub_pd(y2r,W2R,u); \
__m256d v=_mm256_mul_pd(y2r,W2I), z2i=_mm256_fmadd_pd(y2i,W2R,v); \
_mm256_store_pd(R1+c,z2r); _mm256_store_pd(I1+c,z2i); } \
{ __m256d u1r=_mm256_add_pd(br,di), u1i=_mm256_sub_pd(bi,dr); \
__m256d u=_mm256_mul_pd(u1i,W1I), y1r=_mm256_fmsub_pd(u1r,W1R,u); \
__m256d v=_mm256_mul_pd(u1r,W1I), y1i=_mm256_fmadd_pd(u1i,W1R,v); \
_mm256_store_pd(R2+c,y1r); _mm256_store_pd(I2+c,y1i); } \
{ __m256d u3r=_mm256_sub_pd(br,di), u3i=_mm256_add_pd(bi,dr); \
__m256d u=_mm256_mul_pd(u3i,W3I), y3r=_mm256_fmsub_pd(u3r,W3R,u); \
__m256d v=_mm256_mul_pd(u3r,W3I), y3i=_mm256_fmadd_pd(u3i,W3R,v); \
_mm256_store_pd(R3+c,y3r); _mm256_store_pd(I3+c,y3i); } }while(0)
#define R4DIT1(R0,R1,R2,R3,I0,I1,I2,I3) do{ \
__m256d x0r=_mm256_load_pd(R0+c),x0i=_mm256_load_pd(I0+c); \
__m256d x1r=_mm256_load_pd(R1+c),x1i=_mm256_load_pd(I1+c); \
__m256d x2r=_mm256_load_pd(R2+c),x2i=_mm256_load_pd(I2+c); \
__m256d x3r=_mm256_load_pd(R3+c),x3i=_mm256_load_pd(I3+c); \
__m256d Ar=_mm256_add_pd(x0r,x1r), Ai=_mm256_add_pd(x0i,x1i); \
__m256d Br=_mm256_sub_pd(x0r,x1r), Bi=_mm256_sub_pd(x0i,x1i); \
__m256d Cr=_mm256_add_pd(x2r,x3r), Ci=_mm256_add_pd(x2i,x3i); \
__m256d Dr=_mm256_sub_pd(x2r,x3r), Di=_mm256_sub_pd(x2i,x3i); \
_mm256_store_pd(R0+c,_mm256_add_pd(Ar,Cr)); _mm256_store_pd(I0+c,_mm256_add_pd(Ai,Ci)); \
_mm256_store_pd(R2+c,_mm256_sub_pd(Ar,Cr)); _mm256_store_pd(I2+c,_mm256_sub_pd(Ai,Ci)); \
_mm256_store_pd(R1+c,_mm256_sub_pd(Br,Di)); _mm256_store_pd(I1+c,_mm256_add_pd(Bi,Dr)); \
_mm256_store_pd(R3+c,_mm256_add_pd(Br,Di)); _mm256_store_pd(I3+c,_mm256_sub_pd(Bi,Dr)); }while(0)
#define R4DIT(R0,R1,R2,R3,I0,I1,I2,I3,W1R,W1I,W2R,W2I) do{ \
__m256d x0r=_mm256_load_pd(R0+c),x0i=_mm256_load_pd(I0+c); \
__m256d x1r=_mm256_load_pd(R1+c),x1i=_mm256_load_pd(I1+c); \
__m256d x2r=_mm256_load_pd(R2+c),x2i=_mm256_load_pd(I2+c); \
__m256d x3r=_mm256_load_pd(R3+c),x3i=_mm256_load_pd(I3+c); \
__m256d a1r=_mm256_fmsub_pd(x1r,W1R,_mm256_mul_pd(x1i,W1I)); \
__m256d a1i=_mm256_fmadd_pd(x1r,W1I,_mm256_mul_pd(x1i,W1R)); \
__m256d a3r=_mm256_fmsub_pd(x3r,W1R,_mm256_mul_pd(x3i,W1I)); \
__m256d a3i=_mm256_fmadd_pd(x3r,W1I,_mm256_mul_pd(x3i,W1R)); \
__m256d Ar=_mm256_add_pd(x0r,a1r), Ai=_mm256_add_pd(x0i,a1i); \
__m256d Br=_mm256_sub_pd(x0r,a1r), Bi=_mm256_sub_pd(x0i,a1i); \
__m256d Cr=_mm256_add_pd(x2r,a3r), Ci=_mm256_add_pd(x2i,a3i); \
__m256d Dr=_mm256_sub_pd(x2r,a3r), Di=_mm256_sub_pd(x2i,a3i); \
__m256d c2r=_mm256_fmsub_pd(Cr,W2R,_mm256_mul_pd(Ci,W2I)); \
__m256d c2i=_mm256_fmadd_pd(Cr,W2I,_mm256_mul_pd(Ci,W2R)); \
__m256d Pd=_mm256_fmsub_pd(Dr,W2R,_mm256_mul_pd(Di,W2I)); \
__m256d Qd=_mm256_fmadd_pd(Dr,W2I,_mm256_mul_pd(Di,W2R)); \
_mm256_store_pd(R0+c,_mm256_add_pd(Ar,c2r)); \
_mm256_store_pd(I0+c,_mm256_add_pd(Ai,c2i)); \
_mm256_store_pd(R2+c,_mm256_sub_pd(Ar,c2r)); \
_mm256_store_pd(I2+c,_mm256_sub_pd(Ai,c2i)); \
_mm256_store_pd(R1+c,_mm256_sub_pd(Br,Qd)); \
_mm256_store_pd(I1+c,_mm256_add_pd(Bi,Pd)); \
_mm256_store_pd(R3+c,_mm256_add_pd(Br,Qd)); \
_mm256_store_pd(I3+c,_mm256_sub_pd(Bi,Pd)); }while(0)
#define R4BODYT(R0,R1,R2,R3,I0,I1,I2,I3) do{ \
__m256d x0r=_mm256_load_pd(R0+c),x0i=_mm256_load_pd(I0+c); \
__m256d x1r=_mm256_load_pd(R1+c),x1i=_mm256_load_pd(I1+c); \
__m256d x2r=_mm256_load_pd(R2+c),x2i=_mm256_load_pd(I2+c); \
__m256d x3r=_mm256_load_pd(R3+c),x3i=_mm256_load_pd(I3+c); \
__m256d ar=_mm256_add_pd(x0r,x2r), ai=_mm256_add_pd(x0i,x2i); \
__m256d br=_mm256_sub_pd(x0r,x2r), bi=_mm256_sub_pd(x0i,x2i); \
__m256d cr=_mm256_add_pd(x1r,x3r), ci=_mm256_add_pd(x1i,x3i); \
__m256d dr=_mm256_sub_pd(x1r,x3r), di=_mm256_sub_pd(x1i,x3i); \
_mm256_store_pd(R0+c,_mm256_add_pd(ar,cr)); \
_mm256_store_pd(I0+c,_mm256_add_pd(ai,ci)); \
_mm256_store_pd(R1+c,_mm256_sub_pd(ar,cr)); \
_mm256_store_pd(I1+c,_mm256_sub_pd(ai,ci)); \
_mm256_store_pd(R2+c,_mm256_add_pd(br,di)); \
_mm256_store_pd(I2+c,_mm256_sub_pd(bi,dr)); \
_mm256_store_pd(R3+c,_mm256_sub_pd(br,di)); \
_mm256_store_pd(I3+c,_mm256_add_pd(bi,dr)); }while(0)
#ifndef NST
#define NST 5
#endif
alignas(64) static f64 RT1r[176],RT1i[176],RT2r[176],RT2i[176],RT3r[176],RT3i[176];
alignas(64) static f64 RDT1r[168],RDT1i[168],RDT2r[168],RDT2i[168],RDT3r[168],RDT3i[168];
alignas(64) static f64 CDT1r[85],CDT1i[85],CDT2r[85],CDT2i[85],CDT3r[85],CDT3i[85];
__attribute__((noinline)) static void row_dif(f64 *pr, f64 *pi) {
const int n=NC;
{ const int h=128;
for (int j=0;j<h;j+=4) {
__m256d w1r=_mm256_load_pd(RT1r+j), w1i=_mm256_load_pd(RT1i+j);
__m256d w2r=_mm256_load_pd(RT2r+j), w2i=_mm256_load_pd(RT2i+j);
__m256d w3r=_mm256_load_pd(RT3r+j), w3i=_mm256_load_pd(RT3i+j);
for (int s=0;s<n;s+=4*h) {
f64 *a=pr+s+j, *b=pi+s+j;
const int c=0;
R4BODY(a,a+h,a+2*h,a+3*h, b,b+h,b+2*h,b+3*h, w1r,w1i,w2r,w2i,w3r,w3i);
}
} }
{ const int h=32;
for (int j=0;j<h;j+=4) {
__m256d w1r=_mm256_load_pd(RT1r+128+j), w1i=_mm256_load_pd(RT1i+128+j);
__m256d w2r=_mm256_load_pd(RT2r+128+j), w2i=_mm256_load_pd(RT2i+128+j);
__m256d w3r=_mm256_load_pd(RT3r+128+j), w3i=_mm256_load_pd(RT3i+128+j);
for (int s=0;s<n;s+=4*h) {
f64 *a=pr+s+j, *b=pi+s+j;
const int c=0;
R4BODY(a,a+h,a+2*h,a+3*h, b,b+h,b+2*h,b+3*h, w1r,w1i,w2r,w2i,w3r,w3i);
}
} }
{ const int h=8;
for (int j=0;j<h;j+=4) {
__m256d w1r=_mm256_load_pd(RT1r+160+j), w1i=_mm256_load_pd(RT1i+160+j);
__m256d w2r=_mm256_load_pd(RT2r+160+j), w2i=_mm256_load_pd(RT2i+160+j);
__m256d w3r=_mm256_load_pd(RT3r+160+j), w3i=_mm256_load_pd(RT3i+160+j);
for (int s=0;s<n;s+=4*h) {
f64 *a=pr+s+j, *b=pi+s+j;
const int c=0;
R4BODY(a,a+h,a+2*h,a+3*h, b,b+h,b+2*h,b+3*h, w1r,w1i,w2r,w2i,w3r,w3i);
}
} }
for (int s=0;s<n;s+=16) { low3_dif(pr+s,pi+s); low3_dif(pr+s+8,pi+s+8); }
}
static f64 T1r[400],T1i[400],T2r[400],T2i[400],T3r[400],T3i[400];
FT static void build_r4(void){
int k=0;
for (int h=NR>>2; h>=1; h>>=2){
for(int j=0;j<h;j++){
int q1=256*j/h, q2=512*j/h, q3=768*j/h;
T1r[k+j]=DQr[q1]; T1i[k+j]=DQi[q1];
T2r[k+j]=DQr[q2]; T2i[k+j]=DQi[q2];
T3r[k+j]=DQr[q3]; T3i[k+j]=DQi[q3];
} k+=h; }
k=0;
for (int h=NC>>2; h>=8; h>>=2){
for(int j=0;j<h;j++){
int q1=256*j/h, q2=512*j/h, q3=768*j/h;
RT1r[k+j]=DQr[q1]; RT1i[k+j]=DQi[q1];
RT2r[k+j]=DQr[q2]; RT2i[k+j]=DQi[q2];
RT3r[k+j]=DQr[q3]; RT3i[k+j]=DQi[q3];
} k+=h; }
{ int m=0;
for (int h=8; h<NC; h<<=2){
for(int j=0;j<h;j++){ int q1=512*j/h, q2=256*j/h;
RDT1r[m+j]=DQr[q1]; RDT1i[m+j]=-DQi[q1];
RDT2r[m+j]=DQr[q2]; RDT2i[m+j]=-DQi[q2];
RDT3r[m+j]=DQi[q2]; RDT3i[m+j]=DQr[q2];
} m+=h; } }
{ int m=0;
for (int h=1; h*4<=512; h<<=2){
for(int j=0;j<h;j++){ int q1=512*j/h, q2=256*j/h;
CDT1r[m+j]=DQr[q1]; CDT1i[m+j]=-DQi[q1];
CDT2r[m+j]=DQr[q2]; CDT2i[m+j]=-DQi[q2];
CDT3r[m+j]=DQi[q2]; CDT3i[m+j]=DQr[q2];
} m+=h; } }
}
__attribute__((noinline)) static void col_dif(int j0){ const int n=NR;
{
const int h=n>>2; const int K=0;
for(int j=0;j<h;j++){
__m256d w1r=_mm256_set1_pd(T1r[K+j]),w1i=_mm256_set1_pd(T1i[K+j]);
__m256d w2r=_mm256_set1_pd(T2r[K+j]),w2i=_mm256_set1_pd(T2i[K+j]);
__m256d w3r=_mm256_set1_pd(T3r[K+j]),w3i=_mm256_set1_pd(T3i[K+j]);
const size_t st1=(size_t)h*STRIDE, st4=(size_t)4*h*STRIDE;
f64 *r0=RE+(size_t)j*STRIDE+j0,*i0=IM+(size_t)j*STRIDE+j0;
#if 1
if (j+4<h){
const char *p0=(const char*)(RE+(size_t)(j+4)*STRIDE+j0);
const char *q0=(const char*)(IM+(size_t)(j+4)*STRIDE+j0);
const size_t stp=(size_t)h*STRIDE*8;
for (int o=0;o<CBW*8;o+=64){
_mm_prefetch(p0+o,_MM_HINT_T0); _mm_prefetch(p0+stp+o,_MM_HINT_T0);
_mm_prefetch(p0+2*stp+o,_MM_HINT_T0); _mm_prefetch(p0+3*stp+o,_MM_HINT_T0);
_mm_prefetch(q0+o,_MM_HINT_T0); _mm_prefetch(q0+stp+o,_MM_HINT_T0);
_mm_prefetch(q0+2*stp+o,_MM_HINT_T0); _mm_prefetch(q0+3*stp+o,_MM_HINT_T0); } }
#endif
for(int s=0;s<n;s+=4*h){
f64 *r1=r0+st1,*i1=i0+st1, *r2=r1+st1,*i2=i1+st1, *r3=r2+st1,*i3=i2+st1;
#if ZSKIP
for(int c=0;c<CBW;c+=4) R4BODY0(r0,r1,r2,r3,i0,i1,i2,i3,w1r,w1i,w2r,w2i,w3r,w3i);
#else
for(int c=0;c<CBW;c+=4) R4BODY(r0,r1,r2,r3,i0,i1,i2,i3,w1r,w1i,w2r,w2i,w3r,w3i);
#endif
r0+=st4; i0+=st4; }
}
}
const int GW=n>>2, K64=n>>2, K16=K64+(n>>4), K4=K16+(n>>6);
for (int g=0; g<n; g+=GW) {
{ const int h=GW>>2;
const size_t st1=(size_t)h*STRIDE, st4=(size_t)4*h*STRIDE;
for(int j=0;j<h;j++){
__m256d w1r=_mm256_set1_pd(T1r[K64+j]),w1i=_mm256_set1_pd(T1i[K64+j]);
__m256d w2r=_mm256_set1_pd(T2r[K64+j]),w2i=_mm256_set1_pd(T2i[K64+j]);
__m256d w3r=_mm256_set1_pd(T3r[K64+j]),w3i=_mm256_set1_pd(T3i[K64+j]);
f64 *r0=RE+(size_t)(g+j)*STRIDE+j0,*i0=IM+(size_t)(g+j)*STRIDE+j0;
#if 1
if (j+1<h){
const size_t stp=(size_t)h*STRIDE*8;
const char*p0=(const char*)(RE+(size_t)(g+j+1)*STRIDE+j0);
const char*q0=(const char*)(IM+(size_t)(g+j+1)*STRIDE+j0);
for (int o=0;o<CBW*8;o+=64){
_mm_prefetch(p0+o,_MM_HINT_T0); _mm_prefetch(p0+stp+o,_MM_HINT_T0);
_mm_prefetch(p0+2*stp+o,_MM_HINT_T0); _mm_prefetch(p0+3*stp+o,_MM_HINT_T0);
_mm_prefetch(q0+o,_MM_HINT_T0); _mm_prefetch(q0+stp+o,_MM_HINT_T0);
_mm_prefetch(q0+2*stp+o,_MM_HINT_T0); _mm_prefetch(q0+3*stp+o,_MM_HINT_T0); }
}
#endif
for(int s=g;s<g+GW;s+=4*h){
f64 *r1=r0+st1,*i1=i0+st1, *r2=r1+st1,*i2=i1+st1, *r3=r2+st1,*i3=i2+st1;
for(int c=0;c<CBW;c+=4) R4BODY(r0,r1,r2,r3,i0,i1,i2,i3,w1r,w1i,w2r,w2i,w3r,w3i);
r0+=st4; i0+=st4; }
}
}
for (int g2=g; g2<g+GW; g2+=64) {
{ const int h=16;
const size_t st1=(size_t)h*STRIDE;
for(int j=0;j<h;j++){
__m256d w1r=_mm256_set1_pd(T1r[K16+j]),w1i=_mm256_set1_pd(T1i[K16+j]);
__m256d w2r=_mm256_set1_pd(T2r[K16+j]),w2i=_mm256_set1_pd(T2i[K16+j]);
__m256d w3r=_mm256_set1_pd(T3r[K16+j]),w3i=_mm256_set1_pd(T3i[K16+j]);
f64 *r0=RE+(size_t)(g2+j)*STRIDE+j0,*i0=IM+(size_t)(g2+j)*STRIDE+j0;
f64 *r1=r0+st1,*i1=i0+st1, *r2=r1+st1,*i2=i1+st1, *r3=r2+st1,*i3=i2+st1;
#if 1
if (j+1<h){
const size_t stp=(size_t)h*STRIDE*8;
const char*p0=(const char*)(RE+(size_t)(g2+j+1)*STRIDE+j0);
const char*q0=(const char*)(IM+(size_t)(g2+j+1)*STRIDE+j0);
for (int o=0;o<CBW*8;o+=64){
_mm_prefetch(p0+o,_MM_HINT_T0); _mm_prefetch(p0+stp+o,_MM_HINT_T0);
_mm_prefetch(p0+2*stp+o,_MM_HINT_T0); _mm_prefetch(p0+3*stp+o,_MM_HINT_T0);
_mm_prefetch(q0+o,_MM_HINT_T0); _mm_prefetch(q0+stp+o,_MM_HINT_T0);
_mm_prefetch(q0+2*stp+o,_MM_HINT_T0); _mm_prefetch(q0+3*stp+o,_MM_HINT_T0); }
}
#endif
for(int c=0;c<CBW;c+=4) R4BODY(r0,r1,r2,r3,i0,i1,i2,i3,w1r,w1i,w2r,w2i,w3r,w3i);
}
}
for (int g3=g2; g3<g2+64; g3+=16) {
{ const int h=4;
const size_t st1=(size_t)h*STRIDE;
for(int j=0;j<h;j++){
__m256d w1r=_mm256_set1_pd(T1r[K4+j]),w1i=_mm256_set1_pd(T1i[K4+j]);
__m256d w2r=_mm256_set1_pd(T2r[K4+j]),w2i=_mm256_set1_pd(T2i[K4+j]);
__m256d w3r=_mm256_set1_pd(T3r[K4+j]),w3i=_mm256_set1_pd(T3i[K4+j]);
f64 *r0=RE+(size_t)(g3+j)*STRIDE+j0,*i0=IM+(size_t)(g3+j)*STRIDE+j0;
f64 *r1=r0+st1,*i1=i0+st1, *r2=r1+st1,*i2=i1+st1, *r3=r2+st1,*i3=i2+st1;
for(int c=0;c<CBW;c+=4) R4BODY(r0,r1,r2,r3,i0,i1,i2,i3,w1r,w1i,w2r,w2i,w3r,w3i);
}
}
for (int g4=g3; g4<g3+16; g4+=4) {
f64 *r0=RE+(size_t)g4*STRIDE+j0,*i0=IM+(size_t)g4*STRIDE+j0;
f64 *r1=r0+STRIDE,*i1=i0+STRIDE, *r2=r1+STRIDE,*i2=i1+STRIDE, *r3=r2+STRIDE,*i3=i2+STRIDE;
for(int c=0;c<CBW;c+=4) R4BODYT(r0,r1,r2,r3,i0,i1,i2,i3);
}
}
}
}
}
__attribute__((noinline)) static void row_dit(f64 *pr, f64 *pi) {
const int n=NC;
for (int s=0;s<n;s+=16) { low3_dit(pr+s,pi+s); low3_dit(pr+s+8,pi+s+8); }
{ const int h=8; const int k=0;
for (int j=0;j<h;j+=4) {
__m256d w1r=_mm256_load_pd(RDT1r+k+j), w1i=_mm256_load_pd(RDT1i+k+j);
__m256d w2r=_mm256_load_pd(RDT2r+k+j), w2i=_mm256_load_pd(RDT2i+k+j);
for (int s=0;s<n;s+=4*h) {
f64 *a=pr+s+j, *b=pi+s+j;
const int c=0;
R4DIT(a,a+h,a+2*h,a+3*h, b,b+h,b+2*h,b+3*h, w1r,w1i,w2r,w2i);
}
} }
{ const int h=32; const int k=8;
for (int j=0;j<h;j+=4) {
__m256d w1r=_mm256_load_pd(RDT1r+k+j), w1i=_mm256_load_pd(RDT1i+k+j);
__m256d w2r=_mm256_load_pd(RDT2r+k+j), w2i=_mm256_load_pd(RDT2i+k+j);
for (int s=0;s<n;s+=4*h) {
f64 *a=pr+s+j, *b=pi+s+j;
const int c=0;
R4DIT(a,a+h,a+2*h,a+3*h, b,b+h,b+2*h,b+3*h, w1r,w1i,w2r,w2i);
}
} }
{ const int h=128; const int k=40;
for (int j=0;j<h;j+=4) {
__m256d w1r=_mm256_load_pd(RDT1r+k+j), w1i=_mm256_load_pd(RDT1i+k+j);
__m256d w2r=_mm256_load_pd(RDT2r+k+j), w2i=_mm256_load_pd(RDT2i+k+j);
for (int s=0;s<n;s+=4*h) {
f64 *a=pr+s+j, *b=pi+s+j;
const int c=0;
R4DIT(a,a+h,a+2*h,a+3*h, b,b+h,b+2*h,b+3*h, w1r,w1i,w2r,w2i);
}
} }
}
FT static void col_dit(int j0) {
const int n=NR; int k=n-1;
for (int h=1;h<n;h<<=1) { k-=h;
for (int s=0;s<n;s+=2*h)
for (int j=0;j<h;j++) {
f64 *ar=RE+(size_t)(s+j)*STRIDE+j0, *ai=IM+(size_t)(s+j)*STRIDE+j0;
f64 *br=ar+(size_t)h*STRIDE, *bi=ai+(size_t)h*STRIDE;
__m256d wr=_mm256_set1_pd(CIr[k+j]), wi=_mm256_set1_pd(CIi[k+j]);
for (int c=0;c<128;c+=4) {
__m256d xr=_mm256_load_pd(ar+c), xi=_mm256_load_pd(ai+c);
__m256d yr=_mm256_load_pd(br+c), yi=_mm256_load_pd(bi+c);
__m256d t2=_mm256_mul_pd(yi,wi);
__m256d rr=_mm256_fmsub_pd(yr,wr,t2);
__m256d t3=_mm256_mul_pd(yr,wi);
__m256d ii=_mm256_fmadd_pd(yi,wr,t3);
_mm256_store_pd(ar+c,_mm256_add_pd(xr,rr));
_mm256_store_pd(ai+c,_mm256_add_pd(xi,ii));
_mm256_store_pd(br+c,_mm256_sub_pd(xr,rr));
_mm256_store_pd(bi+c,_mm256_sub_pd(xi,ii));
}
}
}
}
static inline void cw(u32 e, f64 &r, f64 &i, int sign) {
f64 rr = DQr[e>>9]*DTr[e&(NC-1)] - DQi[e>>9]*DTi[e&(NC-1)];
f64 ii = DQr[e>>9]*DTi[e&(NC-1)] + DQi[e>>9]*DTr[e&(NC-1)];
if (sign<0) ii = -ii;
r=rr; i=ii;
}
FT static inline void diagonal_row(int i, int sign, f64 sc) {
{
u32 k1 = BREV[i];
f64 *pr=RE+(size_t)i*STRIDE, *pi=IM+(size_t)i*STRIDE;
f64 rr,ri; cw(k1, rr, ri, sign);
f64 r2r=rr*rr-ri*ri, r2i=2.0*rr*ri;
f64 r3r=r2r*rr-r2i*ri, r3i=r2r*ri+r2i*rr;
f64 r4r=r2r*r2r-r2i*r2i, r4i=2.0*r2r*r2i;
__m256d Cr=_mm256_set_pd(r3r, r2r, rr, 1.0);
__m256d Ci=_mm256_set_pd(r3i, r2i, ri, 0.0);
f64 r8r=r4r*r4r-r4i*r4i, r8i=2.0*r4r*r4i;
f64 r16r=r8r*r8r-r8i*r8i, r16i=2.0*r8r*r8i;
for (int j=0;j<NC;j+=DIAGTILE) {
u32 be = (u32)(((u64)k1*(u32)j) & (u32)(NN-1));
f64 base_r, base_i; cw(be, base_r, base_i, sign);
base_r *= sc; base_i *= sc;
f64 b1r=base_r*r4r-base_i*r4i, b1i=base_r*r4i+base_i*r4r;
f64 b2r=b1r*r4r-b1i*r4i, b2i=b1r*r4i+b1i*r4r;
f64 b3r=b2r*r4r-b2i*r4i, b3i=b2r*r4i+b2i*r4r;
#if PFDIAG>0
if (j+64<NC) { _mm_prefetch((const char*)(pr+j+64),_MM_HINT_T0); _mm_prefetch((const char*)(pi+j+64),_MM_HINT_T0);
_mm_prefetch((const char*)(pr+j+96),_MM_HINT_T0); _mm_prefetch((const char*)(pi+j+96),_MM_HINT_T0); }
#endif
__m256d br4=_mm256_set_pd(b3r,b2r,b1r,base_r);
__m256d bi4=_mm256_set_pd(b3i,b2i,b1i,base_i);
const __m256d R16r=_mm256_set1_pd(r16r), R16i=_mm256_set1_pd(r16i);
for (int t=0;t<DIAGTILE;t+=16) {
#define DCHAIN(BR4,BI4,C,OFF) do{ \
__m256d br=_mm256_permute4x64_pd(BR4,_MM_SHUFFLE(C,C,C,C)); \
__m256d bi=_mm256_permute4x64_pd(BI4,_MM_SHUFFLE(C,C,C,C)); \
__m256d wr=_mm256_fmsub_pd(br,Cr,_mm256_mul_pd(bi,Ci)); \
__m256d wi=_mm256_fmadd_pd(br,Ci,_mm256_mul_pd(bi,Cr)); \
__m256d xr=_mm256_load_pd(pr+j+t+OFF), xi=_mm256_load_pd(pi+j+t+OFF); \
__m256d q2=_mm256_mul_pd(xi,wi); \
__m256d orr=_mm256_fmsub_pd(xr,wr,q2); \
__m256d q3=_mm256_mul_pd(xr,wi); \
__m256d oii=_mm256_fmadd_pd(xi,wr,q3); \
_mm256_store_pd(pr+j+t+OFF,orr); \
_mm256_store_pd(pi+j+t+OFF,oii); }while(0)
DCHAIN(br4,bi4,0,0); DCHAIN(br4,bi4,1,4);
DCHAIN(br4,bi4,2,8); DCHAIN(br4,bi4,3,12);
#undef DCHAIN
__m256d nre=_mm256_fmsub_pd(br4,R16r,_mm256_mul_pd(bi4,R16i));
__m256d nim=_mm256_fmadd_pd(br4,R16i,_mm256_mul_pd(bi4,R16r));
br4=nre; bi4=nim;
}
}
}
}
FT static void diagonal(int sign, f64 sc) {
for (int i=0;i<NR;i++) diagonal_row(i, sign, sc);
}
FT static inline void cw18(u32 e, f64 &r, f64 &i, int sign) {
r = DQ18[e>>9]*DT18[e&(NC-1)] - DQ18i[e>>9]*DT18i[e&(NC-1)];
i = DQ18[e>>9]*DT18i[e&(NC-1)] + DQ18i[e>>9]*DT18[e&(NC-1)];
if (sign<0) i = -i;
}
FT static inline void diagonal18_row(int i, int sign, f64 sc) {
{
u32 k1 = BREV9[i];
f64 *pr=RE+(size_t)i*STRIDE, *pi=IM+(size_t)i*STRIDE;
f64 rr,ri; cw18(k1, rr, ri, sign);
f64 r2r=rr*rr-ri*ri, r2i=2.0*rr*ri;
f64 r3r=r2r*rr-r2i*ri, r3i=r2r*ri+r2i*rr;
f64 r4r=r2r*r2r-r2i*r2i, r4i=2.0*r2r*r2i;
__m256d Cr=_mm256_set_pd(r3r, r2r, rr, 1.0);
__m256d Ci=_mm256_set_pd(r3i, r2i, ri, 0.0);
f64 r8r=r4r*r4r-r4i*r4i, r8i=2.0*r4r*r4i;
f64 r16r=r8r*r8r-r8i*r8i, r16i=2.0*r8r*r8i;
for (int j=0;j<NC;j+=DIAGTILE) {
u32 be = (u32)(((u64)k1*(u32)j) & (u32)(NN2-1));
f64 base_r, base_i; cw18(be, base_r, base_i, sign);
base_r *= sc; base_i *= sc;
f64 b1r=base_r*r4r-base_i*r4i, b1i=base_r*r4i+base_i*r4r;
f64 b2r=b1r*r4r-b1i*r4i, b2i=b1r*r4i+b1i*r4r;
f64 b3r=b2r*r4r-b2i*r4i, b3i=b2r*r4i+b2i*r4r;
__m256d br4=_mm256_set_pd(b3r,b2r,b1r,base_r);
__m256d bi4=_mm256_set_pd(b3i,b2i,b1i,base_i);
const __m256d R16r=_mm256_set1_pd(r16r), R16i=_mm256_set1_pd(r16i);
for (int t=0;t<DIAGTILE;t+=16) {
#define ECHAIN(BR4,BI4,C,OFF) do{ \
__m256d br=_mm256_permute4x64_pd(BR4,_MM_SHUFFLE(C,C,C,C)); \
__m256d bi=_mm256_permute4x64_pd(BI4,_MM_SHUFFLE(C,C,C,C)); \
__m256d wr=_mm256_fmsub_pd(br,Cr,_mm256_mul_pd(bi,Ci)); \
__m256d wi=_mm256_fmadd_pd(br,Ci,_mm256_mul_pd(bi,Cr)); \
__m256d xr=_mm256_load_pd(pr+j+t+OFF), xi=_mm256_load_pd(pi+j+t+OFF); \
__m256d q2=_mm256_mul_pd(xi,wi); \
__m256d orr=_mm256_fmsub_pd(xr,wr,q2); \
__m256d q3=_mm256_mul_pd(xr,wi); \
__m256d oii=_mm256_fmadd_pd(xi,wr,q3); \
_mm256_store_pd(pr+j+t+OFF,orr); \
_mm256_store_pd(pi+j+t+OFF,oii); }while(0)
ECHAIN(br4,bi4,0,0); ECHAIN(br4,bi4,1,4);
ECHAIN(br4,bi4,2,8); ECHAIN(br4,bi4,3,12);
#undef ECHAIN
__m256d nre=_mm256_fmsub_pd(br4,R16r,_mm256_mul_pd(bi4,R16i));
__m256d nim=_mm256_fmadd_pd(br4,R16i,_mm256_mul_pd(bi4,R16r));
br4=nre; bi4=nim;
}
}
}
}
FT static void diagonal18(int sign, f64 sc) {
for (int i=0;i<512;i++) diagonal18_row(i, sign, sc);
}
__attribute__((noinline)) static void col512_dit(int j0) {
const int n=512;
int m=1;
{
const int m1=1, m16=1+4;
const size_t s4=(size_t)4*STRIDE, s16=(size_t)16*STRIDE;
for (int W=0;W<n;W+=64) {
for (int k=0;k<4;k++) {
const int w=W+16*k;
for (int q=0;q<4;q++) {
f64 *ar=RE+(size_t)(w+4*q)*STRIDE+j0, *ai=IM+(size_t)(w+4*q)*STRIDE+j0;
f64 *br=ar+(size_t)STRIDE, *bi=ai+(size_t)STRIDE;
f64 *cr=br+(size_t)STRIDE, *ci=bi+(size_t)STRIDE;
f64 *dr=cr+(size_t)STRIDE, *di=ci+(size_t)STRIDE;
for (int c=0;c<CBW2;c+=4) R4DIT1(ar,br,cr,dr,ai,bi,ci,di);
}
for (int j4=0;j4<4;j4++) {
__m256d w1r=_mm256_set1_pd(CDT1r[m1+j4]), w1i=_mm256_set1_pd(CDT1i[m1+j4]);
__m256d w2r=_mm256_set1_pd(CDT2r[m1+j4]), w2i=_mm256_set1_pd(CDT2i[m1+j4]);
f64 *ar=RE+(size_t)(w+j4)*STRIDE+j0, *ai=IM+(size_t)(w+j4)*STRIDE+j0;
f64 *br=ar+s4, *bi=ai+s4, *cr=br+s4, *ci=bi+s4, *dr=cr+s4, *di=ci+s4;
for (int c=0;c<CBW2;c+=4) R4DIT(ar,br,cr,dr,ai,bi,ci,di,w1r,w1i,w2r,w2i);
}
}
for (int j16=0;j16<16;j16++) {
__m256d w1r=_mm256_set1_pd(CDT1r[m16+j16]), w1i=_mm256_set1_pd(CDT1i[m16+j16]);
__m256d w2r=_mm256_set1_pd(CDT2r[m16+j16]), w2i=_mm256_set1_pd(CDT2i[m16+j16]);
f64 *ar=RE+(size_t)(W+j16)*STRIDE+j0, *ai=IM+(size_t)(W+j16)*STRIDE+j0;
f64 *br=ar+s16, *bi=ai+s16, *cr=br+s16, *ci=bi+s16, *dr=cr+s16, *di=ci+s16;
for (int c=0;c<CBW2;c+=4) R4DIT(ar,br,cr,dr,ai,bi,ci,di,w1r,w1i,w2r,w2i);
}
}
}
m=21;
for (int h=16; h*4<=n; h<<=2) { if (h==16) { continue; }
for (int j=0;j<h;j++) {
__m256d w1r=_mm256_set1_pd(CDT1r[m+j]), w1i=_mm256_set1_pd(CDT1i[m+j]);
__m256d w2r=_mm256_set1_pd(CDT2r[m+j]), w2i=_mm256_set1_pd(CDT2i[m+j]);
if (j+5<h) {
const size_t stp=(size_t)h*STRIDE*8;
const char*p0=(const char*)(RE+(size_t)(j+5)*STRIDE+j0);
const char*q0=(const char*)(IM+(size_t)(j+5)*STRIDE+j0);
for (int o=0;o<CBW2*8;o+=64){
_mm_prefetch(p0+o,_MM_HINT_T0); _mm_prefetch(p0+1*stp+o,_MM_HINT_T0); _mm_prefetch(p0+2*stp+o,_MM_HINT_T0); _mm_prefetch(p0+3*stp+o,_MM_HINT_T0); _mm_prefetch(q0+o,_MM_HINT_T0); _mm_prefetch(q0+1*stp+o,_MM_HINT_T0); _mm_prefetch(q0+2*stp+o,_MM_HINT_T0); _mm_prefetch(q0+3*stp+o,_MM_HINT_T0); }
}
for (int s=0;s<n;s+=4*h) {
f64 *ar=RE+(size_t)(s+j)*STRIDE+j0, *ai=IM+(size_t)(s+j)*STRIDE+j0;
f64 *br=ar+(size_t)h*STRIDE, *bi=ai+(size_t)h*STRIDE;
f64 *cr=br+(size_t)h*STRIDE, *ci=bi+(size_t)h*STRIDE;
f64 *dr=cr+(size_t)h*STRIDE, *di=ci+(size_t)h*STRIDE;
for (int c=0;c<CBW2;c+=4)
R4DIT(ar,br,cr,dr,ai,bi,ci,di,w1r,w1i,w2r,w2i);
}
if (h==64) {
for (int k=0;k<4;k++) {
const int jj=j+(h*k);
f64 *ar=RE+(size_t)jj*STRIDE+j0, *ai=IM+(size_t)jj*STRIDE+j0;
f64 *br=ar+(size_t)256*STRIDE, *bi=ai+(size_t)256*STRIDE;
__m256d wr=_mm256_set1_pd(RIr[jj]), wi=_mm256_set1_pd(RIi[jj]);
for (int c=0;c<CBW2;c+=4) {
__m256d xr=_mm256_load_pd(ar+c), xi=_mm256_load_pd(ai+c);
__m256d yr=_mm256_load_pd(br+c), yi=_mm256_load_pd(bi+c);
__m256d t2=_mm256_mul_pd(yi,wi);
__m256d rr=_mm256_fmsub_pd(yr,wr,t2);
__m256d t3=_mm256_mul_pd(yr,wi);
__m256d ii=_mm256_fmadd_pd(yi,wr,t3);
_mm256_store_pd(ar+c,_mm256_add_pd(xr,rr));
_mm256_store_pd(ai+c,_mm256_add_pd(xi,ii));
_mm256_store_pd(br+c,_mm256_sub_pd(xr,rr));
_mm256_store_pd(bi+c,_mm256_sub_pd(xi,ii));
}
}
}
}
m+=h;
}
}
alignas(64) static f64 SVr[2][NC], SVi[2][NC];
#ifndef R18VEC
#define R18VEC 1
#endif
#ifndef PFR18
#define PFR18 1
#endif
#ifndef PFDIAG
#define PFDIAG 1
#endif
#ifndef PFPAIR
#define PFPAIR 0
#endif
#if R18VEC>0
FT static inline void reformat18_row(int i) {
const f64 *pr = RE+(size_t)i*STRIDE, *pi = IM+(size_t)i*STRIDE;
int r = i>>1, par = i&1;
f64 *wr = RE+(size_t)r*STRIDE + (par?256:0), *wi = IM+(size_t)r*STRIDE + (par?256:0);
u32 k1 = BREV[i]; u32 kq = k1>>9;
f64 dr = DTr[k1 & (NC-1)], di = DTi[k1 & (NC-1)];
if (i<2) { pr = SVr[i]; pi = SVi[i]; }
const f64 *TR = TWR[kq], *TI = TWI[kq];
const __m256d sg=_mm256_set1_pd(-0.0);
const __m256d Vr=_mm256_set1_pd(dr), Vi=_mm256_set1_pd(di);
for (int c=0;c<(NC>>1);c+=4) {
#if PFR18>0
if (c+32 < (NC>>1)) { _mm_prefetch((const char*)(pr+2*c+64),_MM_HINT_T0); _mm_prefetch((const char*)(pi+2*c+64),_MM_HINT_T0);
_mm_prefetch((const char*)(pr+2*c+96),_MM_HINT_T0); _mm_prefetch((const char*)(pi+2*c+96),_MM_HINT_T0); }
#endif
__m256d A=_mm256_load_pd(pr+2*c), B=_mm256_load_pd(pr+2*c+4);
__m256d P=_mm256_permute4x64_pd(A,0xD8), Q=_mm256_permute4x64_pd(B,0xD8);
__m256d z1r=_mm256_permute2f128_pd(P,Q,0x20);
__m256d z2r=_mm256_permute2f128_pd(P,Q,0x31);
A=_mm256_load_pd(pi+2*c); B=_mm256_load_pd(pi+2*c+4);
P=_mm256_permute4x64_pd(A,0xD8); Q=_mm256_permute4x64_pd(B,0xD8);
__m256d z1i=_mm256_permute2f128_pd(P,Q,0x20);
__m256d z2i=_mm256_permute2f128_pd(P,Q,0x31);
__m256d qr0=_mm256_load_pd(TR+c), qi0=_mm256_load_pd(TI+c);
__m256d er=_mm256_add_pd(z1r,z2r), dr2=_mm256_sub_pd(z1r,z2r);
__m256d ei=_mm256_add_pd(z1i,z2i), di2=_mm256_sub_pd(z1i,z2i);
__m256d t0=_mm256_mul_pd(qr0,Vr), tr=_mm256_fnmadd_pd(qi0,Vi,t0);
__m256d t1=_mm256_mul_pd(qr0,Vi), ti=_mm256_fnmadd_pd(qi0,Vr,_mm256_xor_pd(t1,sg));
__m256d t2=_mm256_mul_pd(tr,dr2), orr=_mm256_fnmadd_pd(ti,di2,t2);
__m256d t3=_mm256_mul_pd(tr,di2), oii=_mm256_fmadd_pd(ti,dr2,t3);
_mm256_store_pd(wr+c,_mm256_sub_pd(er,oii));
_mm256_store_pd(wi+c,_mm256_add_pd(ei,orr));
}
}
#else
FT static inline void reformat18_row_scalar(int i);
FT static inline void reformat18_row(int i) { reformat18_row_scalar(i); }
#endif
FT static inline void reformat18_row_scalar(int i) {
{
const f64 *pr = RE+(size_t)i*STRIDE, *pi = IM+(size_t)i*STRIDE;
int r = i>>1, par = i&1;
const u32 *CJ = par ? CJ1 : CJ0;
f64 *wr = RE+(size_t)r*STRIDE, *wi = IM+(size_t)r*STRIDE;
u32 k1 = BREV[i];
u32 kq = k1>>9;
f64 dr = DTr[k1 & (NC-1)], di = DTi[k1 & (NC-1)];
if (i<2) { pr = SVr[i]; pi = SVi[i]; }
for (int j=0;j<NC;j+=2) {
f64 z1r,z1i,z2r,z2i;
{ z1r=pr[j]; z1i=pi[j]; z2r=pr[j+1]; z2i=pi[j+1]; }
u32 eh = kq + 2u*BREVC[j];
f64 qr0=DQr[eh], qi0=DQi[eh];
f64 tr = qr0*dr - qi0*di;
f64 ti = -(qr0*di + qi0*dr);
f64 er=z1r+z2r, ei=z1i+z2i;
f64 dr=z1r-z2r, di=z1i-z2i;
f64 orr=tr*dr-ti*di, oii=tr*di+ti*dr;
int ccol = (int)CJ[j];
wr[ccol] = er - oii;
wi[ccol] = ei + orr;
}
}
}
FT static void reformat18(void) {
memcpy(SVr[0], RE, sizeof(f64)*NC); memcpy(SVi[0], IM, sizeof(f64)*NC);
memcpy(SVr[1], RE+STRIDE, sizeof(f64)*NC); memcpy(SVi[1], IM+STRIDE, sizeof(f64)*NC);
for (int i=0;i<NR;i++) reformat18_row(i);
}
FT static int parse_digits2(const char *s, long len, f64 *dst, int stride) {
int m = (int)((len + 3) / 4);
long lo = len & 3;
if (lo) { u32 v=0; for(long j=0;j<lo;j++) v=v*10+(u32)(s[j]-'0'); long kk=m-1;
dst[(size_t)(kk>>9)*stride + (kk & 511)] = (f64)v; }
const __m256i c10 = _mm256_setr_epi8(10,1,10,1,10,1,10,1,10,1,10,1,10,1,10,1,
10,1,10,1,10,1,10,1,10,1,10,1,10,1,10,1);
const __m256i c100 = _mm256_setr_epi16(100,1,100,1,100,1,100,1,100,1,100,1,100,1,100,1);
const __m256i adj = _mm256_set1_epi32(53328);
const __m256i rv = _mm256_setr_epi32(7,6,5,4,3,2,1,0);
long i = len;
f64 *wp = dst;
int col = 0;
while (i - 64 >= lo) {
__m256i d0 = _mm256_loadu_si256((const __m256i*)(s+i-32));
__m256i d1 = _mm256_loadu_si256((const __m256i*)(s+i-64));
i -= 64;
{ long pa = i-320; if (pa < 0) pa = 0; _mm_prefetch((const char*)(s+pa), _MM_HINT_T0);
long pb = i-448; if (pb < 0) pb = 0; _mm_prefetch((const char*)(s+pb), _MM_HINT_T0); }
__m256i v0 = _mm256_sub_epi32(_mm256_madd_epi16(_mm256_maddubs_epi16(d0,c10),c100), adj);
__m256i v1 = _mm256_sub_epi32(_mm256_madd_epi16(_mm256_maddubs_epi16(d1,c10),c100), adj);
v0 = _mm256_permutevar8x32_epi32(v0, rv);
v1 = _mm256_permutevar8x32_epi32(v1, rv);
__m256d q0=_mm256_cvtepi32_pd(_mm256_castsi256_si128(v0));
__m256d q1=_mm256_cvtepi32_pd(_mm256_extracti128_si256(v0,1));
__m256d q2=_mm256_cvtepi32_pd(_mm256_castsi256_si128(v1));
__m256d q3=_mm256_cvtepi32_pd(_mm256_extracti128_si256(v1,1));
#if NTPARSE>0
_mm256_stream_pd(wp,q0); _mm256_stream_pd(wp+4,q1);
wp += 8; col += 8; if (col == 512) { col = 0; wp += stride - 512; }
_mm256_stream_pd(wp,q2); _mm256_stream_pd(wp+4,q3);
#else
_mm256_store_pd(wp,q0); _mm256_store_pd(wp+4,q1);
wp += 8; col += 8; if (col == 512) { col = 0; wp += stride - 512; }
_mm256_store_pd(wp,q2); _mm256_store_pd(wp+4,q3);
#endif
wp += 8; col += 8; if (col == 512) { col = 0; wp += stride - 512; }
}
while (i - 32 >= lo) {
i -= 32;
__m256i v = _mm256_sub_epi32(_mm256_madd_epi16(
_mm256_maddubs_epi16(_mm256_loadu_si256((const __m256i*)(s+i)),c10),c100), adj);
v = _mm256_permutevar8x32_epi32(v, rv);
__m256d q0=_mm256_cvtepi32_pd(_mm256_castsi256_si128(v));
__m256d q1=_mm256_cvtepi32_pd(_mm256_extracti128_si256(v,1));
#if NTPARSE>0
_mm256_stream_pd(wp,q0); _mm256_stream_pd(wp+4,q1);
#else
_mm256_store_pd(wp,q0); _mm256_store_pd(wp+4,q1);
#endif
wp += 8; col += 8;
if (col == 512) { col = 0; wp += stride - 512; }
}
#if NTPARSE>0
_mm_sfence();
#endif
int k = (int)((i - lo) / 4);
while (i - 4 >= lo) {
i -= 4;
long kk = --k;
dst[(size_t)(kk>>9)*stride + (kk & 511)] =
(f64)((((s[i]-'0')*10 + (s[i+1]-'0'))*10 + (s[i+2]-'0'))*10 + (s[i+3]-'0'));
}
return m;
}
TGT static inline long kb_scan(const char *s, long i, long sn, int want_digit) {
const __m256i lo = _mm256_set1_epi8('0'-1), hi = _mm256_set1_epi8('9'+1);
if (i < sn) {
while ((((uintptr_t)(s + i)) & 31) != 0) {
if (i >= sn) break;
if (((s[i] >= '0' && s[i] <= '9') ? 1 : 0) == want_digit) return i;
i++;
}
long n = sn - 32;
if (want_digit) {
for (; i <= n; i += 32) {
__m256i v = _mm256_load_si256((const __m256i *)(s + i));
unsigned m = (unsigned)_mm256_movemask_epi8(
_mm256_and_si256(_mm256_cmpgt_epi8(v, lo), _mm256_cmpgt_epi8(hi, v)));
if (m) return i + __builtin_ctz(m);
}
while (i < sn && (s[i] < '0' || s[i] > '9')) i++;
} else {
for (; i <= n; i += 32) {
__m256i v = _mm256_load_si256((const __m256i *)(s + i));
unsigned m = (unsigned)_mm256_movemask_epi8(
_mm256_and_si256(_mm256_cmpgt_epi8(v, lo), _mm256_cmpgt_epi8(hi, v)));
if (~m) return i + __builtin_ctz(~m);
}
while (i < sn && s[i] >= '0' && s[i] <= '9') i++;
}
}
return i;
}
TGT static inline int split_lines(const char *s, long sn, const char **b0, long *n0,
const char **b1, long *n1) {
long i = kb_scan(s, 0, sn, 1);
long j = kb_scan(s, i, sn, 0);
if (j == i) return 0;
*b0 = s + i; *n0 = j - i;
long p = kb_scan(s, j, sn, 1);
long k = kb_scan(s, p, sn, 0);
if (k == p) return 0;
*b1 = s + p; *n1 = k - p;
return 1;
}
#define COEF(k) RE[(size_t)((k)>>9)*STRIDE + ((k) & 511)]
static u32 TAB4[10000];
FT static void build_tab4(void) {
for (u32 d = 0; d < 10000; d++)
TAB4[d] = (u32)('0'+d/1000) | ((u32)('0'+(d/100)%10) << 8)
| ((u32)('0'+(d/10)%10) << 16) | ((u32)('0'+d%10) << 24);
}
FT static long carry_format_m(int tot, char *o, char *end) {
const int nm = (tot + 1) >> 1;
int nfull = ((tot - 1) + 1) >> 1;
if (nfull > nm) nfull = nm;
if (nfull < 0) nfull = 0;
u64 carry = 0;
char *p = end;
const f64 *pr = RE, *pi = IM;
int m = 0;
while (m < nfull) {
int cnt = nfull - m; if (cnt > NC) cnt = NC;
const f64 *qr = pr, *qi = pi;
for (int k = 0; k < cnt; k++) {
u64 C = (u64)_mm_cvtsd_si64(_mm_load_sd((const f64*)(qr + k)))
+ (u64)_mm_cvtsd_si64(_mm_load_sd((const f64*)(qi + k))) * 10000ull;
u64 total = C + carry;
u64 t4 = total / 10000ull;
u64 ex = total / 100000000ull;
carry = ex;
u32 lo = (u32)(total - t4 * 10000ull);
u32 hi = (u32)t4 - (u32)ex * 10000u;
p -= 8;
*(u32*)p = TAB4[hi];
*(u32*)(p + 4) = TAB4[lo];
}
pr += STRIDE; pi += STRIDE; m += cnt;
}
pr = RE + (size_t)(nfull/NC)*STRIDE + (nfull%NC);
pi = IM + (size_t)(nfull/NC)*STRIDE + (nfull%NC);
for (; m < nm; m++) {
u64 C = (u64)_mm_cvtsd_si64(_mm_load_sd(pr));
u64 total = C + carry;
u64 t4 = total / 10000ull;
u64 ex = total / 100000000ull;
carry = ex;
u32 lo = (u32)(total - t4 * 10000ull);
u32 hi = (u32)t4 - (u32)ex * 10000u;
p -= 8;
*(u32*)p = TAB4[hi];
*(u32*)(p + 4) = TAB4[lo];
pr++; pi++;
}
while (carry) { u32 d = (u32)(carry % 100000000ull); carry /= 100000000ull;
u32 hi = d / 10000u, lo = d % 10000u;
p -= 8; *(u32*)p = TAB4[hi]; *(u32*)(p+4) = TAB4[lo]; }
while (p < end && *p == '0') p++;
if (p == end) *(--p) = '0';
(void)o;
return (long)(end - p);
}
FT static inline void pointwise_row(int i);
static unsigned char FWD_DONE[NR];
#if FUSE16>0
FT static void fft_forward2(int sign, f64 sc);
#endif
int g_kstop=3;
FT static void fft_forward(void) {
#ifdef REVCD
for (int j0=NC-CBW;j0>=0;j0-=CBW) col_dif(j0);
#else
for (int j0=0;j0<NC;j0+=CBW) col_dif(j0);
#endif
memset(FWD_DONE, 0, sizeof FWD_DONE);
for (int r=0;r<(NR>>1);r++) {
const int i0 = 2*r;
for (int t=0;t<2;t++) {
int x = i0+t;
if (FWD_DONE[x]) continue;
int m = (int)RVIX[x];
#if ROWPF>0
{ int n1=x+1; if(n1<NR && !FWD_DONE[n1]){ const char*q0=(const char*)(RE+(size_t)n1*STRIDE); const char*w0=(const char*)(IM+(size_t)n1*STRIDE);
for(int o=0;o<NC*8;o+=64){ _mm_prefetch(q0+o,_MM_HINT_T0); _mm_prefetch(w0+o,_MM_HINT_T0); } }
int n2=(int)RVIX[x<NR-1?x+1:x]; if(n2!=x+1 && !FWD_DONE[n2]){ const char*q1=(const char*)(RE+(size_t)n2*STRIDE); const char*w1=(const char*)(IM+(size_t)n2*STRIDE);
for(int o=0;o<NC*8;o+=64){ _mm_prefetch(q1+o,_MM_HINT_T0); _mm_prefetch(w1+o,_MM_HINT_T0); } } }
#endif
diagonal_row(x, +1, 1.0);
row_dif(RE+(size_t)x*STRIDE, IM+(size_t)x*STRIDE);
if (m != x) {
diagonal_row(m, +1, 1.0);
row_dif(RE+(size_t)m*STRIDE, IM+(size_t)m*STRIDE);
FWD_DONE[m] = 1;
}
FWD_DONE[x] = 1;
pointwise_row(x);
if (m != x) pointwise_row(m);
}
if (r == 0) {
{ f64 zr=RE[1], zi=IM[1]; RE[1]=4.0*zr*zi; IM[1]=0.0; }
memcpy(SVr[0], RE, sizeof(f64)*NC); memcpy(SVi[0], IM, sizeof(f64)*NC);
memcpy(SVr[1], RE+STRIDE, sizeof(f64)*NC); memcpy(SVi[1], IM+STRIDE, sizeof(f64)*NC);
}
reformat18_row(i0);
reformat18_row(i0+1);
row_dit(RE+(size_t)r*STRIDE, IM+(size_t)r*STRIDE);
diagonal18_row(r, -1, 1.0/(4.0*(f64)NN));
}
for (int j0=0;j0<NC;j0+=CBW2) col512_dit(j0);
}
FT static void fft_inverse(void) {
for (int i=0;i<NR;i++) row_dit(RE+(size_t)i*STRIDE, IM+(size_t)i*STRIDE);
diagonal(-1, 1.0/(f64)NN);
for (int j0=0; !g_kstop && j0<NC;j0+=128) col_dit(j0);
}
#if FUSE16>0
FT static void fused_row_loop(int sign, f64 sc) {
for (int i=0;i<NR;i++) {
if (!FWD_DONE[i]) {
int m = (int)RVIX[i];
diagonal_row(i, +1, 1.0);
row_dif(RE+(size_t)i*STRIDE, IM+(size_t)i*STRIDE);
if (m != i) { diagonal_row(m, +1, 1.0); row_dif(RE+(size_t)m*STRIDE, IM+(size_t)m*STRIDE); }
FWD_DONE[i] = 1; FWD_DONE[m] = 1;
pointwise_row(i);
if (m != i) pointwise_row(m);
}
if (i & 1) {
int r = i>>1;
if (i == 1) { f64 zr=RE[1], zi=IM[1]; RE[1]=4.0*zr*zi; IM[1]=0.0; }
reformat18_row(i-1);
reformat18_row(i);
row_dit(RE+(size_t)r*STRIDE, IM+(size_t)r*STRIDE);
diagonal18_row(r, sign, sc);
}
}
}
#endif
__attribute__((unused)) FT static void fft_inverse18(int sign, f64 sc) {
memcpy(SVr[0], RE, sizeof(f64)*NC); memcpy(SVi[0], IM, sizeof(f64)*NC);
memcpy(SVr[1], RE+STRIDE, sizeof(f64)*NC); memcpy(SVi[1], IM+STRIDE, sizeof(f64)*NC);
for (int i=0;i<NR;i++) {
reformat18_row(i);
if (i & 1) { int r = i>>1;
row_dit(RE+(size_t)r*STRIDE, IM+(size_t)r*STRIDE);
diagonal18_row(r, sign, sc); }
}
for (int j0=0;j0<NC;j0+=CBW2) col512_dit(j0);
}
#ifndef PWVEC
#define PWVEC 1
#endif
#if PWVEC>0
#define BLD(X,Y,M) _mm256_blend_pd(X,Y,M)
#define PWC(z1r,z1i,z2r,z2i,cr,ci) do{ \
__m256d a1r=_mm256_add_pd(z1r,z2r), a1i=_mm256_sub_pd(z1i,z2i); \
__m256d b1r=_mm256_add_pd(z1i,z2i), b1i=_mm256_sub_pd(z2r,z1r); \
\
__m256d u=_mm256_mul_pd(a1i,b1i); cr=_mm256_fmsub_pd(a1r,b1r,u); \
__m256d v=_mm256_mul_pd(a1i,b1r); ci=_mm256_fmadd_pd(a1r,b1i,v); }while(0)
FT static void pointwise_vec(int i) {
unsigned ri = RVIX[i];
f64 *pr=RE+(size_t)i*STRIDE, *pi=IM+(size_t)i*STRIDE;
f64 *qr=RE+(size_t)ri*STRIDE, *qi=IM+(size_t)ri*STRIDE;
const __m256d sg=_mm256_set1_pd(-0.0);
for(int A=0;A<256;A+=8){
int B=504-A;
__m256d Ur0=_mm256_load_pd(pr+A), Ur1=_mm256_load_pd(pr+A+4);
__m256d Ui0=_mm256_load_pd(pi+A), Ui1=_mm256_load_pd(pi+A+4);
__m256d Vr0=_mm256_load_pd(qr+B), Vr1=_mm256_load_pd(qr+B+4);
__m256d Vi0=_mm256_load_pd(qi+B), Vi1=_mm256_load_pd(qi+B+4);
__m256d z2r0=_mm256_permute4x64_pd(Vr1,0x1B),z2r1=_mm256_permute4x64_pd(Vr0,0x1B);
__m256d z2i0=_mm256_permute4x64_pd(Vi1,0x1B),z2i1=_mm256_permute4x64_pd(Vi0,0x1B);
__m256d cr0,ci0,cr1,ci1;
PWC(Ur0,Ui0,z2r0,z2i0,cr0,ci0);
PWC(Ur1,Ui1,z2r1,z2i1,cr1,ci1);
_mm256_store_pd(pr+A, cr0); _mm256_store_pd(pr+A+4,cr1);
_mm256_store_pd(pi+A, ci0); _mm256_store_pd(pi+A+4,ci1);
_mm256_store_pd(qr+B, _mm256_permute4x64_pd(cr1,0x1B));
_mm256_store_pd(qr+B+4, _mm256_permute4x64_pd(cr0,0x1B));
_mm256_store_pd(qi+B, _mm256_xor_pd(_mm256_permute4x64_pd(ci1,0x1B),sg));
_mm256_store_pd(qi+B+4, _mm256_xor_pd(_mm256_permute4x64_pd(ci0,0x1B),sg));
}
}
#endif
FT static inline void pointwise_row_scalar(int i);
FT static inline void pointwise_row(int i) {
#if PWVEC>0
unsigned ri = RVIX[i];
if (BREV[i] && ri != (unsigned)i) { pointwise_vec(i); return; }
#endif
pointwise_row_scalar(i);
}
FT static inline void pointwise_row_scalar(int i) {
{
unsigned ri = RVIX[i];
const u32 *CX = BREV[i] ? CVIX1 : CVIX0;
f64 *pr=RE+(size_t)i*STRIDE, *pi=IM+(size_t)i*STRIDE;
f64 *qr=RE+(size_t)ri*STRIDE, *qi=IM+(size_t)ri*STRIDE;
for (int j=0;j<NC;j+=2) {
int jr = (int)CX[j];
f64 z1r=pr[j], z1i=pi[j], z2r=qr[jr], z2i=qi[jr];
f64 a1r=(z1r+z2r), a1i=(z1i-z2i);
f64 b1r=(z1i+z2i), b1i=-(z1r-z2r);
f64 cr=a1r*b1r - a1i*b1i, ci=a1r*b1i + a1i*b1r;
pr[j]=cr; pi[j]=ci;
qr[jr]=cr; qi[jr]=-ci;
}
}
}
FT static void pointwise(void) {
for (int i=0;i<NR;i++) pointwise_row(i);
{ f64 zr=RE[1], zi=IM[1]; RE[1]=4.0*zr*zi; IM[1]=0.0; }
}
FT static void run_job(DI *d) {
const char *s; long sn;
s = d->s; sn = (long)d->sn;
const char *b0,*b1; long n0,n1;
const long first=kb_scan(s,0,sn,1);
const long second=kb_scan(s,first+1000000,sn,1);
if(second+1000000>sn)return;
b0=s+first; b1=s+second; n0=n1=1000000;
build_all();
const int m0 = (int)((n0 + 3) / 4), m1 = (int)((n1 + 3) / 4);
parse_digits2(b0, n0, RE, STRIDE);
parse_digits2(b1, n1, IM, STRIDE);
int tot = m0 + m1 - 1;
#if FUSE16>0
fft_forward();
memcpy(SVr[0], RE, sizeof(f64)*NC); memcpy(SVi[0], IM, sizeof(f64)*NC);
memcpy(SVr[1], RE+STRIDE, sizeof(f64)*NC); memcpy(SVi[1], IM+STRIDE, sizeof(f64)*NC);
fused_row_loop(-1, 1.0/(4.0*(f64)NN));
for (int j0=0;j0<NC;j0+=CBW2) col512_dit(j0);
#else
fft_forward();
#endif
long olim = (long)d->ol - 1; if (olim > 2000000) olim = 2000000;
char *cap = (char *)d->o + olim;
build_tab4();
long L = carry_format_m(tot, cap, cap);
if (cap - L != (char *)d->o) memmove(d->o, cap - L, (size_t)L);
((char *)d->o)[L] = '\n';
d->os = L + 1;
}
#ifndef NO_LOCAL_TEST
extern "C" void __libc_start_main(void *m, int argc, char **argv) {
(void)m;
unsigned long *p = (unsigned long *)(argv + argc + 1);
while (*p) p++;
p++;
DI *d = 0;
for (int i = 0; i < 32 && p[0]; i++, p += 2)
if (p[0] == 0x6b637564UL) { d = (DI *)p[1]; break; }
if (d) run_job(d);
__asm__ volatile("syscall" ::"a"(60), "D"(0) : "rcx", "r11", "memory");
for (;;);
}
int main(){ return 0; }
#endif
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 5.258 ms | 10 MB + 368 KB | Accepted | Score: 100 | 显示更多 |