src/SyrComponent/SyrComponent_CheckOfUndefined.cxx

   1 // Copyright (C) 2007-2024  CEA, EDF, OPEN CASCADE
   2 //
   3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 // This library is free software; you can redistribute it and/or
   7 // modify it under the terms of the GNU Lesser General Public
   8 // License as published by the Free Software Foundation; either
   9 // version 2.1 of the License, or (at your option) any later version.
  10 //
  11 // This library is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 // Lesser General Public License for more details.
  15 //
  16 // You should have received a copy of the GNU Lesser General Public
  17 // License along with this library; if not, write to the Free Software
  18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22
  23 //  SuperVisionTest SyrComponent : example of component performing some mathinatical operations
  24 //  File   : SyrComponent_CheckOfUndefined.cxx
  25 //  Module : SuperVisionTest
  26 //
  27 #include <stdio.h>
  28 #ifndef WIN32
  29 #include <unistd.h>
  30 #endif
  31 #include <iostream>
  32 #include <fstream>
  33 #include <strstream>
  34 #include <string>
  35 #include <math.h>
  36 #include <cstring>
  37
  38 using namespace std;
  39
  40 // ------------------------------------------------------------------
  41 // NextPrime : Compute the first prime number greater or equal than an integer
  42 // ------------------------------------------------------------------
  43
  44 #define VALUESNBR 4
  45
  46 long NextPrime (const long me )
  47 {
  48
  49  struct svalue {int signiaib ;
  50                 int nbr ;} ;
  51
  52  struct svalue values[VALUESNBR] ;
  53  long ia ;
  54  long maxia ;
  55  long ib[4] ;
  56  int n[4] ;
  57 // int signiaib[4] = { -1 , +1 , +1 , -1 } ;
  58  int signiaib[4];
  59  signiaib[0] = -1;
  60  signiaib[1] = 1;
  61  signiaib[2] = 1;
  62  signiaib[3] = -1;
  63  long remain ;
  64
  65  int nbvalues ;
  66  int loop ;
  67  int nindd ;
  68  long minn ;
  69  long maxvn ;
  70  long premret = 0 ;
  71
  72    if (me < 0 || me >
  73 #ifdef DECOSF1
  74                       127149130704178201
  75 #else
  76                       2147483647
  77 #endif
  78                                ){
  79 //      Standard_RangeError::
  80 //       Raise("Try to apply NextPrime method with negative, null or too large value.");
  81         return 0 ;
  82    }
  83
  84    if ( me <= 7 ) {
  85      if ( me <= 1 )
  86        return 1 ;
  87      else if ( me <= 2 )
  88        return 2 ;
  89      else if ( me <= 3 )
  90        return 3 ;
  91      else if ( me <= 5 )
  92        return 5 ;
  93      else if ( me <= 7 )
  94        return 7 ;
  95    }
  96
  97    minn = ( me - 1 ) / 6 ;              // n minimum
  98    while ( 6*minn+1 < me ) {
  99         minn += 1 ;
 100       }
 101
 102    maxia = long( sqrt( (double)me ) / 6 + 1 ) ;
 103
 104    maxvn = minn + VALUESNBR ;
 105
 106    nbvalues = 0 ;
 107    for ( nindd = 0 ; nindd < VALUESNBR ; nindd++ ) {
 108       if ( 6*(nindd+minn)-1 < me ) {
 109         values[nindd].nbr = 1 ;
 110         values[nindd].signiaib = -1 ;
 111         nbvalues += 1 ;
 112       }
 113       else {
 114         values[nindd].nbr = 0 ;
 115         values[nindd].signiaib = 0 ;
 116       }
 117     }
 118
 119    for ( ia = 1 ; ia <= maxia ; ia++ ) {
 120       if ( nbvalues == VALUESNBR*2 ) {
 121         break ;
 122       }
 123       remain = -VALUESNBR ;
 124       ib[0] = ( minn + ia - remain ) / (6*ia - 1) ;
 125       n[0] = int ( 6*ia*ib[0] - ia - ib[0] - minn ) ;
 126       ib[1] = ( minn - ia - remain ) / (6*ia - 1) ;
 127       n[1] = int ( 6*ia*ib[1] + ia - ib[1] - minn ) ;
 128       ib[2] = ( minn + ia - remain ) / (6*ia + 1) ;
 129       n[2] = int ( 6*ia*ib[2] - ia + ib[2] - minn ) ;
 130       ib[3] = ( minn - ia - remain ) / (6*ia + 1) ;
 131       n[3] = int ( 6*ia*ib[3] + ia + ib[3] - minn ) ;
 132       for ( loop = 0 ; loop < 4 ; loop++ ) {
 133          if ( n[loop] >= 0 && n[loop] < VALUESNBR ) {
 134            if ( ( values[n[loop]].nbr == 0 ) ||
 135                 ( values[n[loop]].signiaib == signiaib[loop] ) ) {
 136              values[n[loop]].signiaib = -signiaib[loop] ;
 137              values[n[loop]].nbr += 1 ;
 138              if ( values[n[loop]].nbr <= 2 )
 139                nbvalues += 1 ;
 140            }
 141          }
 142        }
 143     }
 144    for ( nindd = 0 ; nindd < VALUESNBR ; nindd++ ) {
 145      if ( values[nindd].nbr == 0 ) {
 146        if ( me <= 6*(nindd+minn)-1 ) {
 147           premret = 6*(nindd+minn)-1 ;
 148           break ;
 149         }
 150        else if ( me <= 6*(nindd+minn)+1 ) {
 151          premret = 6*(nindd+minn)+1 ;
 152          break ;
 153        }
 154      }
 155      else if ( values[nindd].nbr == 1 ) {
 156         if ( values[nindd].signiaib > 0 ) {
 157           if ( me <= 6*(nindd+minn)-1 ) {
 158             premret = 6*(nindd+minn)-1 ;
 159             break ;
 160           }
 161         }
 162         else {
 163           if ( me <= 6*(nindd+minn)+1 ) {
 164             premret = 6*(nindd+minn)+1 ;
 165             break ;
 166           }
 167         }
 168       }
 169    }
 170
 171    if ( premret != 0 ) {
 172      return premret ;
 173    }
 174
 175  return NextPrime ( 6*(maxvn-1)+2) ;
 176
 177 }
 178
 179 static char * base( int b , long n ) {
 180   static char c[40] ;
 181   char ccc[2] ;
 182   int i = 39 ;
 183   c[39] = '\0' ;
 184   c[38] = '0' ;
 185   if ( n ) {
 186     while ( n ) {
 187       i-- ;
 188       sprintf( ccc ,"%d" , ( n % b ) ) ;
 189       c[i] = ccc[0] ;
 190       n = n / b ;
 191     }
 192   }
 193   else
 194     i = 38 ;
 195   return &c[i] ;
 196 }
 197
 198 static char * B2( long n ) {
 199   return base( 2 , n ) ;
 200 }
 201
 202 static char * B3( long n ) {
 203   return base( 3 , n ) ;
 204 }
 205
 206 static char * B4( long n ) {
 207   return base( 4 , n ) ;
 208 }
 209
 210 static char * B5( long n ) {
 211   return base( 5 , n ) ;
 212 }
 213
 214 static char * B7( long n ) {
 215   return base( 7 , n ) ;
 216 }
 217
 218 static char * B8( long n ) {
 219   return base( 8 , n ) ;
 220 }
 221
 222 static long S( long f , long n ) {
 223   int i ;
 224   long res = 0 ;
 225   for ( i = 0 ; i <= n ; i++ ) {
 226     res += 1 << (f*i) ;
 227   }
 228   return res ;
 229 }
 230
 231 static long S2( long n ) {
 232   return S( 2 , n ) ;
 233 }
 234
 235 static long S6( long n ) {
 236   return S( 6 , n ) ;
 237 }
 238
 239 static long rP( long odd ) {
 240   if ( (( odd - 1 ) % 6 ) == 0 )
 241     return 1 ;
 242   if ( (( odd - 3 ) % 6 ) == 0 )
 243     return 3 ;
 244   if ( (( odd - 5 ) % 6 ) == 0 )
 245     return 5 ;
 246   return 0 ;
 247 }
 248 static long nP( long odd ) {
 249   return ( odd / 6 ) ;
 250 }
 251
 252 static long rQ( long N ) {
 253   if ( ( N & 7 ) == 5 ) {
 254     return 5 ;
 255   }
 256   return ( N & 3 ) ;
 257 }
 258 static long nQ( long N ) {
 259   if ( ( N & 3 ) == 3 ) {
 260     return ( N / 4 ) ;
 261   }
 262   else {
 263     return ( N / 8 ) ;
 264   }
 265 }
 266
 267 static long rT( long N ) {
 268   return ( nP( N ) % 3 ) ;
 269 }
 270 static long nT( long N ) {
 271   return ( nP( N ) / 3 ) ;
 272 }
 273
 274 static char * fact( int n ) {
 275   static char chn[132] ;
 276   if ( n == 0 ) {
 277     strcpy( chn , "0" ) ;
 278   }
 279   else if ( n == 1 ) {
 280     strcpy( chn , "1" ) ;
 281   }
 282   else {
 283     int pow2 = 0 ;
 284     int pow3 = 0 ;
 285     while ( (n & 1 ) == 0 ) {
 286       pow2 += 1 ;
 287       n = (n >> 1 ) ;
 288     }
 289     while ( ( n % 3 ) == 0 ) {
 290       pow3 += 1 ;
 291       n = n / 3 ;
 292     }
 293     int pos = 0 ;
 294     if ( pow2 ) {
 295       if ( pow2 == 1 ) {
 296         sprintf( &chn[pos] , "%d" , 2 ) ;
 297       }
 298       else {
 299         sprintf( &chn[pos] , "2**%d" , pow2 ) ;
 300       }
 301       pos = strlen( chn ) ;
 302     }
 303     if ( pow3 ) {
 304       if ( pow2 ) {
 305         strcat( chn , "*" ) ;
 306         pos = strlen( chn ) ;
 307       }
 308       if ( pow3 == 1 ) {
 309         sprintf( &chn[pos] , "%d" , 3 ) ;
 310       }
 311       else {
 312         sprintf( &chn[pos] , "3**%d" , pow3 ) ;
 313       }
 314       pos = strlen( chn ) ;
 315     }
 316     if ( n != 1 ) {
 317       if ( pow2 || pow3 ) {
 318         sprintf( &chn[pos] , "*%d" , n ) ;
 319       }
 320       else {
 321         sprintf( &chn[pos] , "%d" , n ) ;
 322       }
 323     }
 324   }
 325   return chn ;
 326 }
 327
 328 static long Next( const long prev , long &pow2 ) {
 329   long next = ( 3*prev + 1 )/2 ;
 330   pow2 = 1 ;
 331   while ( ( next & 1 ) == 0 ) {
 332     next = ( next >> 1 ) ;
 333     pow2 += 1 ;
 334   }
 335   return next ;
 336 }
 337
 338 int T1a( int n0 , int n1 ) {
 339   return (1<<(6*n1+2))*(3*n0+0) + (1<<(6*n1+1)) + S2(3*n1+0) ;
 340 }
 341 int T1b( int n0 , int n1 ) {
 342   return (1<<(6*n1+6))*(3*n0+1) + (1<<(6*n1+5)) + S2(3*n1+2) ;
 343 }
 344 int T1c( int n0 , int n1 ) {
 345   return (1<<(6*n1+4))*(3*n0+2) + (1<<(6*n1+3)) + S2(3*n1+1) ;
 346 }
 347 int T2a( int n0 , int n1 ) {
 348   return (1<<(6*n1+7))*(3*n0+0) + S2(3*n1+2) ;
 349 }
 350 int T2b( int n0 , int n1 ) {
 351   return (1<<(6*n1+3))*(3*n0+1) + S2(3*n1+0) ;
 352 }
 353 int T2c( int n0 , int n1 ) {
 354   return (1<<(6*n1+5))*(3*n0+2) + S2(3*n1+1) ;
 355 }
 356
 357 int P1a1( int n0 , int n1 ) {
 358   return (1<<(6*n1+4))*(3*n0+0) + (1<<(6*n1+3)) + S2(3*n1+1) ;
 359 }
 360 int P1a2( int n0 , int n1 ) {
 361   return (1<<(6*n1+6))*(3*n0+0) + (1<<(6*n1+5)) + S2(3*n1+2) ;
 362 }
 363 int P1b1( int n0 , int n1 ) {
 364   return (1<<(6*n1+2))*(3*n0+1) + (1<<(6*n1+1)) + S2(3*n1+0) ;
 365 }
 366 int P1b2( int n0 , int n1 ) {
 367   return (1<<(6*n1+4))*(3*n0+1) + (1<<(6*n1+3)) + S2(3*n1+1) ;
 368 }
 369 int P1c1( int n0 , int n1 ) {
 370   return (1<<(6*n1+2))*(3*n0+2) + (1<<(6*n1+1)) + S2(3*n1+0) ;
 371 }
 372 int P1c2( int n0 , int n1 ) {
 373   return (1<<(6*n1+6))*(3*n0+2) + (1<<(6*n1+5)) + S2(3*n1+2) ;
 374 }
 375
 376 int P2a1( int n0 , int n1 ) {
 377   return (1<<(6*n1+3))*(3*n0+0) + S2(3*n1+0) ;
 378 }
 379 int P2a2( int n0 , int n1 ) {
 380   return (1<<(6*n1+5))*(3*n0+0) + S2(3*n1+1) ;
 381 }
 382 int P2b1( int n0 , int n1 ) {
 383   return (1<<(6*n1+5))*(3*n0+1) + S2(3*n1+1) ;
 384 }
 385 int P2b2( int n0 , int n1 ) {
 386   return (1<<(6*n1+7))*(3*n0+1) + S2(3*n1+2) ;
 387 }
 388 int P2c1( int n0 , int n1 ) {
 389   return (1<<(6*n1+3))*(3*n0+2) + S2(3*n1+0) ;
 390 }
 391 int P2c2( int n0 , int n1 ) {
 392   return (1<<(6*n1+7))*(3*n0+2) + S2(3*n1+2) ;
 393 }
 394
 395 typedef int (*fct)(int,int);
 396
 397 ostream * _Trace ;
 398
 399 void Compute( char * fctname , fct afct , int PT , int PP , int PPNext ) {
 400   int n0 ;
 401   int n1 ;
 402   long value = 0 ;
 403   long next = 0 ;
 404   for ( n1 = 0 ; n1 <= 3 ; n1++ ) {
 405     bool pn1 = true ;
 406     for ( n0 = 0 ; n0 <= 500 && value < 3000 ; n0++ ) {
 407       value = afct( n0 , n1 ) ;
 408       long pow2 ;
 409       next = Next( value , pow2 ) ;
 410       long kp2 ;
 411       if ( ( rP( value ) == 3 && rP( Next( next , kp2 ) ) == PPNext ) ||
 412              ( rT( value ) == PT && rP( value ) == PP &&
 413                rP( next ) == PPNext ) ) {
 414         if ( rQ( value ) == 3 ||
 415              ( rQ( value ) == 1 && ( nQ( value ) & 1 ) == 0 ) ) {
 416 //          if ( pn1 ) {
 417 //            *_Trace << "  n1 " << n1 ;
 418 //            *_Trace << " n0 " << n0 << " " ;
 419 //            pn1 = false ;
 420 //          }
 421 //          else {
 422             if ( value < 10 )
 423               *_Trace << "000" ;
 424             else if ( value < 100 )
 425               *_Trace << "00" ;
 426             else if ( value < 1000 )
 427               *_Trace << "0" ;
 428             *_Trace << value << " " ;
 429             long lownext = (next-1)/4 ;
 430             if ( 4*lownext+1 == next && ( lownext & 1 ) == 1 ) {
 431               *_Trace << "---" ;
 432             }
 433             *_Trace << fctname << " n0 " << n0 << " " ;
 434 //          }
 435           *_Trace << value << " = 6[3*" << nT( value ) << "+"
 436                   << rT( value ) << "]+" << rP( value ) << " --> "
 437                   << next << " (>>" << pow2 << ") = 6[3*"
 438                   << nT( next ) << "+" << rT( next )
 439                   << "] + " << rP( next ) << endl ;
 440         }
 441 #if 0
 442         else {
 443           if ( pn1 ) {
 444             *_Trace << "  n1 " << n1 ;
 445             *_Trace << " n0 " << n0 << " " ;
 446             pn1 = false ;
 447           }
 448           else {
 449             *_Trace << "       n0 " << n0 << " " ;
 450           }
 451           *_Trace << "---" << value << " = 6[3*" << nT( value ) << "+"
 452                   << rT( value ) << "]+" << rP( value ) << " = 4*["
 453                   << "2*" << nQ( value )/2 << "+1]+1" ;
 454           if ( Next( (value-1)/4 , kp2 ) != next )
 455             *_Trace << " " << Next( (value-1)/4 , kp2 ) << " != " << next ;
 456           *_Trace << endl ;
 457         }
 458 #endif
 459       }
 460     }
 461     value = 0 ;
 462   }
 463 }
 464
 465 int SeqSyr( int n ) {
 466   long pow2 ;
 467   int nstep = 0 ;
 468   int EQ ;
 469   int r ;
 470   int q ;
 471   int Sr ;
 472   int Sq ;
 473   int nn = n ;
 474   *_Trace << "   " << n ;
 475   while ( n != 1 ) {
 476     n = Next( n , pow2 ) ;
 477     nstep += 1 ;
 478     int EQ ;
 479     r = rQ( n ) ;
 480     q = nQ( n ) ;
 481     if ( r == 3 ) {
 482       EQ = 4 ;
 483     }
 484     else {
 485       EQ = 8 ;
 486     }
 487     if ( ( nstep % 2 ) == 0 ) {
 488       *_Trace << endl << "   " ;
 489     }
 490     Sr = rP( n ) ;
 491     Sq = nP( n ) ;
 492     *_Trace << " -> " << n << " = " << EQ << "*" << q << "+" << r
 493             << " = 6*" << Sq << "+" << Sr ;
 494     if ( r == 5 ) {
 495       break ;
 496     }
 497   }
 498   *_Trace << endl ;
 499   if ( r == 5 && q > nn/8 ) {
 500     return n ;
 501   }
 502   return 0 ;
 503 }
 504
 505 void little( int depth , int prt , int n , int r ) {
 506
 507   int n0 ;
 508   int r0 ;
 509   int n1 ;
 510   int r1 ;
 511   int r2 ;
 512   depth += 1 ;
 513
 514   if ( depth == prt )
 515     *_Trace << endl << depth << " little2 n = " << n << " r = " << r << endl ;
 516
 517   n0 = n*2 ;
 518   r0 = r*2 - 1 ;
 519   if ( ( r0 % 3 ) == 0 ) {
 520     r1 = r0/3 ;
 521     if ( depth == prt ) {
 522       *_Trace << "n->3n " << n0 << "n + " << r1 << " = 4[ "
 523               << n0/4 << "n + " << r1/4 << " ]+3" << endl ;
 524       SeqSyr( r1 ) ;
 525     }
 526   }
 527   else if ( ( ( n0 + r0 ) % 3 ) == 0 ) {
 528     r1 = ( n0 + r0 ) / 3 ;
 529     if ( depth == prt ) {
 530       *_Trace << "n->3n+1 " << n0 << "n + " << r1 << " = 4[ "
 531               << n0/4 << "n + " << r1/4 << " ]+3" << endl ;
 532       SeqSyr( r1 ) ;
 533     }
 534   }
 535   else if ( ( ( 2*n0 + r0 ) % 3 ) == 0 ) {
 536     r1 = ( 2*n0 + r0 ) / 3 ;
 537     if ( depth == prt ) {
 538       *_Trace << "n->3n+2 " << n0 << "n + " << r1 << " = 4[ "
 539               << n0/4 << "n + " << r1/4 << " ]+3" << endl ;
 540       SeqSyr( r1 ) ;
 541     }
 542   }
 543
 544   if ( depth == prt )
 545     *_Trace << endl << depth << " little4 n = " << n << " r = " << r << endl ;
 546
 547   n1 = n0*2 ;
 548   r2 = r0*2 + 1 ;
 549   if ( ( r2 % 3 ) == 0 ) {
 550     r2 = r2/3 ;
 551     if ( depth == prt ) {
 552       *_Trace << "n->3n " << n1 << "n + " << r2 << " = 8[ "
 553               << n1/8 << "n + " << r2/8 << " ]+1" << endl ;
 554       SeqSyr( r2 ) ;
 555     }
 556   }
 557   else if ( ( ( n1 + r2 ) % 3 ) == 0 ) {
 558     r2 = ( n1 + r2 ) / 3 ;
 559     if ( depth == prt ) {
 560       *_Trace << "n->3n+1 " << n1 << "n + " << r2 << " = 8[ "
 561               << n1/8 << "n + " << r2/8 << " ]+1" << endl ;
 562       SeqSyr( r2 ) ;
 563     }
 564   }
 565   else if ( ( ( 2*n1 + r2 ) % 3 ) == 0 ) {
 566     r2 = ( 2*n1 + r2 ) / 3 ;
 567     if ( depth == prt ) {
 568       *_Trace << "n->3n+2 " << n1 << "n + " << r2 << " = 8[ "
 569               << n1/8 << "n + " << r2/8 << " ]+1" << endl ;
 570       SeqSyr( r2 ) ;
 571     }
 572   }
 573
 574   if ( depth < prt ) {
 575     little( depth , prt , n0 , r1 ) ;
 576     little( depth , prt , n1 , r2 ) ;
 577   }
 578
 579 }
 580
 581 int main(int argc, char **argv) {
 582
 583   int k ;
 584   int i ;
 585   int j ;
 586   int l ;
 587   int l1 ;
 588   int n = 0 ;
 589   int p ;
 590   int P ;
 591   int Prem = 37 ;
 592
 593   _Trace = new ofstream( "/cas_01/jr/linux/SyrP.log" );
 594
 595   for ( l = 0 ; l <= 17 ; l++ ) {
 596     for ( l1 = -1 ; l1 <= 1 ; l1 += 2 ) {
 597       P = 6*l + l1 ;
 598       if ( P > 0 && NextPrime( P-1 ) == P ) {
 599         for ( k = -1 ; k <= 1 ; k += 2 ) {
 600           for ( i = 1 ; i <= 30 ; i++ ) {
 601             for ( j = 1 ; j <= 30 ; j++ ) {
 602               n = int( pow( (double)2 , (double)i ) * pow( (double)3 , (double)j )*P ) + k ;
 603               if ( n < 0 || n > pow( (double)2 , (double)30 ) ) {
 604                 break ;
 605               }
 606               if ( NextPrime( n-1 ) == n ) {
 607                 *_Trace << n << " = 2**" << i << " * 3**" << j << " * " << P << " + "
 608                         << k << endl ;
 609               }
 610 #if 0
 611               if ( NextPrime( n-1 ) != n && ( n % P ) == 0 ) {
 612                 *_Trace << n << " = 2**" << i << " * 3**" << j << " * " << P << " + "
 613                         << k << " = " ;
 614                 p = 5 ;
 615                 while ( p*p <= n ) {
 616                   while ( ( n % p ) == 0 ) {
 617                     *_Trace << p << "*" ;
 618                     n = n / p ;
 619                   }
 620                   p = NextPrime( p + 1 ) ;
 621                 }
 622                 *_Trace << n << endl ;
 623               }
 624 #endif
 625             }
 626             *_Trace << endl ;
 627           }
 628         }
 629       }
 630       *_Trace << endl ;
 631     }
 632   }
 633
 634 #if 0
 635   _Trace = new ofstream( "/cas_01/jr/linux/Syr.log" );
 636
 637   int n = 1 ;
 638   int prt ;
 639
 640 #if 0
 641   for ( prt = 1 ; prt <= 5 ; prt++ ) {
 642     *_Trace << endl << prt << ". 8n + 1 :" << endl ;
 643     little( 0 , prt , 8 , 1 ) ;
 644   }
 645
 646   for ( prt = 1 ; prt <= 5 ; prt++ ) {
 647     *_Trace << endl << prt << ". 4n + 3 :" << endl ;
 648     little( 0 , prt , 4 , 3 ) ;
 649   }
 650
 651   for ( prt = 1 ; prt <= 6 ; prt++ ) {
 652     *_Trace << endl << prt << ". 8n + 5 :" << endl ;
 653     little( 0 , prt , 8 , 5 ) ;
 654   }
 655 #endif
 656
 657   int i ;
 658   int max = 256 ;
 659   for ( i = 0 ; i <= max ; i++ ) {
 660     *_Trace << endl << endl << "8*" << i << "+5 = 4[2*" << i << "+1]+1 = " ;
 661     int ii = 2*i+1 ;
 662     while ( ((ii-1) & 7) == 4 ) {
 663       *_Trace << endl << "        8*" << ii/8 << "+5 = 4[2*" << ii/4
 664               << "+1]+1 = " ;
 665       ii = 2*(ii/4)+1 ;
 666     }
 667     *_Trace << "6*" << nP(8*i + 5) << "+" << rP(8*i + 5) << " :" << endl ;
 668     n = SeqSyr( 8*i + 5 ) ;
 669     n = n/8 ;
 670     while ( n && n > max ) {
 671       *_Trace << endl << "   ==> 8*" << n << "+5 = 4[2*" << n << "+1]+1 = " ;
 672       int ii = 2*n+1 ;
 673       while ( ((ii-1) & 7) == 4 ) {
 674         *_Trace << endl << "8*" << ii/8 << "+5 = 4[2*" << ii/4 << "+1]+1 = " ;
 675         ii = 2*(ii/4)+1 ;
 676       }
 677       *_Trace << "6*" << nP(8*n + 5) << "+" << rP(8*n + 5) << " :" << endl ;
 678       n = SeqSyr( 8*n + 5 ) ;
 679       n = n/8 ;
 680     }
 681   }
 682
 683 #if 0
 684   int x ;
 685   int rx ;
 686   int a ;
 687   int q ;
 688   int r ;
 689   int X ;
 690   int n ;
 691   int rn ;
 692   int xx ;
 693   int y ;
 694   int p ;
 695   int pow ;
 696
 697   *_Trace << "2**[x(n+1)] = (2**x-1)*S[x](n)+1" << endl << endl ;
 698
 699   *_Trace << "2**[(2x+1)(2n+rn+1)] = (2**(2x+1)-1)*S[2x+1](2n+rn)+1" << endl ;
 700   *_Trace << "                     = 2**(2x+1){(2**(2x+1)-1)*S[2x+1](2n+rn)-1}"
 701           << endl << endl ;
 702   *_Trace << "2**[2x(3n+rn+1)] = (2**2x-1)*S[2x](3n+rn)+1" << endl ;
 703   *_Trace << "                 = 2**2x{(2**2x-1)*S[2x](3n+rn-1)-1}"
 704           << endl << endl ;
 705   *_Trace << "2**(2x+1)-1 = 6S2(x-1)+1" << endl << endl ;
 706
 707   *_Trace << "2**(2(3x+1))-1 = 3(6*2*7*S6(x-1)+1) *4 :" << endl ;
 708   *_Trace << "2**(2(3x+2))-1 = 3(6*2**3*7*S6(x-1)+5) *4 :" << endl ;
 709   *_Trace << "2**(2(3x+3))-1 = 3(6*2**5*7*S6(x-1)+21)" << endl ;
 710   *_Trace << "               = 3**2(2**6*7*S6(x-1)+7)" << endl ;
 711   *_Trace << "               = 3**2*7(2**6*S6(x-1)+1)" << endl ;
 712   *_Trace << "               = 3**2*7*S6(x)" << endl ;
 713   *_Trace << "2**(2(3x))-1 = 3**2*7*S6(x-1)" << endl ;
 714   *_Trace << "2**(2(3x+4))-1 = 3(6*2*7*S6(x)+1) *4 :" << endl ;
 715   *_Trace << "2**(2(3x+5))-1 = 3(6*2**3*7*S6(x)+5)" << endl << endl ;
 716   *_Trace << "2**(2(3x+rx))-1 = 3**2*2**(2rx)*7*S6(x-1)+2**(2rx)-1" << endl
 717           << endl ;
 718
 719   *_Trace << "2**[(2x+1)(2n+rn+1)] = (6S2(x-1)+1)*S[2x+1](2n+rn)+1" << endl ;
 720   *_Trace << "2**[2(3x+rx)(3n+rn+1)] = {3**2*2**(2rx)*7*S6(x-1)+2**(2rx)-1}*"
 721           << endl << "                 S[2(3x+rx)](3n+rn)+1" << endl
 722           << endl ;
 723
 724   *_Trace << "(6x+5)y+x = 2**3*7*S6(x-1)" << endl ;
 725   *_Trace << "x = 0 ; y = 0" << endl ;
 726   *_Trace << "x = 1 ; y = 5" << endl ;
 727   *_Trace << "x = 2 ; y = 214 = 2*107 = 2*O153" << endl ;
 728   *_Trace << "x = 3 ; y = 10131 = 3*3377 = 3*O6461" << endl ;
 729   *_Trace << "x = 4 ; y = 514244 = 4*128561 = 4*O373061" << endl ;
 730   *_Trace << "x = 5 ; y = ?" << endl ;
 731   *_Trace << "x = 6 ; y = 1489853154 = 6*248308859 = 6*O1663162173" << endl
 732           << endl ;
 733
 734   for ( pow = 1 ; pow <= 31 ; pow++ ) {
 735     *_Trace << endl ;
 736     for ( x = 0 ; x <= 30 ; x++ ) {
 737       for ( n = 0 ; n <= 30 ; n++ ) {
 738         for ( rn = 0 ; rn <= 1 ; rn++ ) {
 739           if ( (2*x+1)*(2*n+rn+1) == pow && (2*x+1)*(2*n+rn+1) <= 31 ) {
 740             int val = (1 << (2*x+1)*(2*n+rn+1)) ;
 741             *_Trace << "2**(" << 2*x+1 << "*(" << 2*n+rn << "+1)) = " ;
 742             *_Trace << "2**[(2*" << x << "+1)(2*" << n << "+" << rn
 743                     << "+1)] = (6*" ;
 744             int SIG = S(2,x-1) ;
 745             int pow3 = 0 ;
 746             while ( SIG && (SIG % 3) == 0 ) {
 747               pow3 += 1 ;
 748               SIG = SIG/3 ;
 749             }
 750             *_Trace << "3**" << pow3 << "*" << SIG
 751                     << "+1) * " << fact( S(2*x+1,2*n+rn) )
 752                     << " + 1 = " << endl
 753                     << "      (6*S2(" << x-1 << ")+1)*S" << 2*x+1 << "("
 754                     << 2*n+rn << ")+1 = "
 755                     << "(6*S2(" << x << "-1)+1)*S[2*" << x << "+1](2*"
 756                     << n << "+" << rn << ")+1 = " << val << endl ;
 757             if ( val != (6*S(2,x-1)+1)*S(2*x+1,2*n+rn)+1 ) {
 758               *_Trace << val << " != " << (6*S(2,x-1)+1)*S(2*x+1,2*n+rn)+1
 759                       << endl ;
 760             }
 761           }
 762         }
 763       }
 764     }
 765     for ( x = 0 ; x <= 30 ; x++ ) {
 766       for ( rx = 0 ; rx <= 2 ; rx++ ) {
 767         for ( n = 0 ; n <= 30 ; n++ ) {
 768           for ( rn = 0 ; rn <= 2 ; rn++ ) {
 769             if ( 3*x+rx != 0 && 2*(3*x+rx)*(3*n+rn+1) == pow
 770                  && 2*(3*x+rx)*(3*n+rn+1) <= 31 ) {
 771               *_Trace << "2**[" << 2*(3*x+rx) << "*(" << 3*n+rn << "+1)] = " ;
 772               *_Trace << "2**[2(3*" << x << "+" << rx << ")(3*" << n << "+"
 773                       << rn << "+1)] = [9*" << (1 << (2*rx)) << "*7*"
 774                       << S(6,x-1) << "+" ;
 775               int DRX = (1 << (2*rx))-1 ;
 776               int pow3 = 0 ;
 777               while ( DRX && (DRX % 3) == 0 ) {
 778                 pow3 += 1 ;
 779                 DRX = DRX/3 ;
 780               }
 781               *_Trace << "3**" << pow3 << "*" << DRX << "] * " ;
 782               int SIG = S(2*(3*x+rx),3*n+rn) ;
 783               pow3 = 0 ;
 784               while ( SIG && (SIG % 3) == 0 ) {
 785                 pow3 += 1 ;
 786                 SIG = SIG/3 ;
 787               }
 788               *_Trace << "3**" << pow3 << "*" << fact( SIG )
 789                       << " + 1 = " << endl
 790                       << "   [3**2*2**" << 2*rx << "*7*S6("
 791                       << x << "-1)+" << (1 << (2*rx))-1 << "] * S"
 792                       << 2*(3*x+rx) << "(" << 3*n+rn << ")+1 = " << endl
 793                       << "      {3**2*2**(2*" << rx << ")*7*S6("
 794                       << x << "-1)+2**(2*" << rx << ")-1} * S[2(3*"
 795                       << x << "+" << rx << ")](3*" << n << "+" << rn << ")+1"
 796                       << endl ;
 797               int val = (1 << (2*(3*x+rx)*(3*n+rn+1)) ) ;
 798               if ( val != (3*3*(1<<(2*rx))*7*S(6,x-1)+(1<<(2*rx))-1)*S(2*(3*x+rx),3*n+rn)+1 ) {
 799                 *_Trace << val << " != " << (3*3*(1<<(2*rx))*7*S(6,x-1)+(1<<(2*rx))-1) * S(2*(3*x+rx),3*n+rn)+1 << endl ;
 800               }
 801             }
 802           }
 803         }
 804       }
 805     }
 806   }
 807   return 0;
 808
 809
 810   *_Trace << endl ;
 811   for ( x = 1 ; x <= 30 ; x++ ) {
 812     for ( n = 0 ; n <= 30 ; n++ ) {
 813       if ( x*(n+1) <= 30 ) {
 814         X = ((1 << x ) - 1)*S(x,n) + 1 ;
 815         if ( (1 << (x*(n+1)) ) == X ) {
 816           *_Trace << "2**[" << x << "(" << n << "+1)] = 2**[[" ;
 817           int pow2 = 0 ;
 818           int pow3 = 0 ;
 819           int p ;
 820           p = x ;
 821           while ( p && (p & 1) == 0 ) {
 822             pow2 += 1 ;
 823             p = (p >> 1 ) ;
 824           }
 825           while ( p && (p % 3) == 0 ) {
 826             pow3 += 1 ;
 827             p = p/3 ;
 828           }
 829           if ( pow2 ) {
 830             *_Trace << "2" ;
 831             if ( pow2 != 1 ) {
 832               *_Trace << "**" << pow2 ;
 833             }
 834             if ( pow3 || p != 1 ) {
 835               *_Trace << "*" ;
 836             }
 837           }
 838           if ( pow3 ) {
 839             *_Trace << "3" ;
 840             if ( pow3 != 1 ) {
 841               *_Trace << "**" << pow3 ;
 842             }
 843             if ( p != 1 ) {
 844               *_Trace << "*" ;
 845             }
 846           }
 847           if ( !((pow2 || pow3) && p == 1) ) {
 848             *_Trace << p ;
 849           }
 850           *_Trace << "][" ;
 851           p = n+1 ;
 852           pow2 = 0 ;
 853           pow3 = 0 ;
 854           while ( p && (p & 1) == 0 ) {
 855             pow2 += 1 ;
 856             p = (p >> 1 ) ;
 857           }
 858           while ( p && (p % 3) == 0 ) {
 859             pow3 += 1 ;
 860             p = p/3 ;
 861           }
 862           if ( pow2 ) {
 863             *_Trace << "2" ;
 864             if ( pow2 != 1 ) {
 865               *_Trace << "**" << pow2 ;
 866             }
 867             if ( pow3 || p != 1 ) {
 868               *_Trace << "*" ;
 869             }
 870           }
 871           if ( pow3 ) {
 872             *_Trace << "3" ;
 873             if ( pow3 != 1 ) {
 874               *_Trace << "**" << pow3 ;
 875             }
 876             if ( p != 1 ) {
 877               *_Trace << "*" ;
 878             }
 879           }
 880           if ( !((pow2 || pow3) && p == 1) ) {
 881             *_Trace << p ;
 882           }
 883           *_Trace << "]]" ;
 884           int d = (1 << x) - 1 ;
 885           pow3 = 0 ;
 886           while ( (d % 3) == 0 ) {
 887             pow3 += 1 ;
 888             d = d / 3 ;
 889           }
 890           *_Trace << " = (2**" << x << "-1)*S" ;
 891           if ( x & 1 ) {
 892             *_Trace << x << "[2*" << n/2 ;
 893             if ( n & 1 ) {
 894               *_Trace << "+1" ;
 895             }
 896             *_Trace << "]+1 = " ;
 897           }
 898           else {
 899             *_Trace << x << "[3*" << n/3 << "+" << (n % 3) << "]+1 = " ;
 900           }
 901           *_Trace << "{" ;
 902           if ( pow3 ) {
 903             *_Trace << "3" ;
 904             if ( pow3 != 1 ) {
 905               *_Trace << "**" << pow3 ;
 906             }
 907             *_Trace << "*" ;
 908           }
 909           if ( d == 1 ) {
 910           }
 911           else if ( d == 5 ) {
 912             *_Trace << "[5]" ;
 913           }
 914           else {
 915             if ( nP(d) ) {
 916               *_Trace << "[6*" << nP(d) << "+" << rP(d) << "]" ;
 917             }
 918             else {
 919               *_Trace << "[" << rP(d) << "]" ;
 920             }
 921           }
 922           *_Trace << "}*{" ;
 923           int Y = S(x,n) ;
 924           pow2 = 0 ;
 925           while ( (Y % 2) == 0 ) {
 926             pow2 += 1 ;
 927             Y = Y / 2 ;
 928           }
 929           pow3 = 0 ;
 930           while ( (Y % 3) == 0 ) {
 931             pow3 += 1 ;
 932             Y = Y / 3 ;
 933           }
 934           if ( pow2 ) {
 935             *_Trace << "2" ;
 936             if ( pow2 != 1 ) {
 937               *_Trace << "**" << pow2 ;
 938             }
 939             *_Trace << "*" ;
 940           }
 941           if ( pow3 ) {
 942             *_Trace << "3" ;
 943             if ( pow3 != 1 ) {
 944               *_Trace << "**" << pow3 ;
 945             }
 946             *_Trace << "*" ;
 947           }
 948           if ( nP(Y) && ( (nP(Y) & 1) == 0 || (nP(Y) % 3) == 0 ) ) {
 949             *_Trace << "[6*" ;
 950             int YY = nP(Y) ;
 951             pow2 = 0 ;
 952             while ( (YY % 2) == 0 ) {
 953               pow2 += 1 ;
 954               YY = YY / 2 ;
 955             }
 956             pow3 = 0 ;
 957             while ( (YY % 3) == 0 ) {
 958               pow3 += 1 ;
 959               YY = YY / 3 ;
 960             }
 961             if ( pow2 ) {
 962               *_Trace << "2" ;
 963               if ( pow2 != 1 ) {
 964                 *_Trace << "**" << pow2 ;
 965               }
 966               *_Trace << "*" ;
 967             }
 968             if ( pow3 ) {
 969               *_Trace << "3" ;
 970               if ( pow3 != 1 ) {
 971                 *_Trace << "**" << pow3 ;
 972               }
 973               *_Trace << "*" ;
 974             }
 975             *_Trace << YY << "+" << rP(Y) << "]}+1" << endl ;
 976           }
 977           else {
 978             if ( nP(Y) ) {
 979               *_Trace << "[6*" << nP(Y) << "+" << rP(Y) << "]}+1" << endl ;
 980             }
 981             else {
 982               *_Trace << "[" << rP(Y) << "]}+1" << endl ;
 983             }
 984           }
 985         }
 986         else {
 987           *_Trace << endl << (1 << (x*(n+1)) ) << " = 2**[" << x << "(" << n
 988                   << "+1)] != " << X << "(2**" << x
 989                   << "-1)*S" << x << "[" << n << "]+1" << endl ;
 990         }
 991       }
 992     }
 993   }
 994   return 0;
 995
 996   *_Trace << "S[x]( n ) = Sigma((2**x)*i) i = 0 , n" << endl ;
 997   *_Trace << "S[x]( aq + r ) = X*S[a*x]( q - 1 ) + Sx( r )" << endl ;
 998   for ( x = 2 ; x <= 10 ; x++ ) {
 999     *_Trace << endl << endl ;
1000     for ( a = 2 ; a <= 12 ; a++ ) {
1001       *_Trace << endl ;
1002       for ( q = 1 ; q <= 12 ; q++ ) {
1003         for ( r = 0 ; r < a ; r++ ) {
1004           if ( a*q + r <= 24 ) {
1005             int i = S( x , a*q+r ) - S( x , r ) ;
1006             int X = i / S( a*x , q-1 ) ;
1007             int p2 = 0 ;
1008             if ( X*S( a*x , q-1 ) == i && (i % S( a*x , q-1 ) ) == 0 ) {
1009               while ( ( X & 1 ) == 0 ) {
1010                 p2 += 1 ;
1011                 X = ( X >> 1 ) ;
1012               }
1013               int p3 = 0 ;
1014               while ( X > 0 && ( X % 3 ) == 0 ) {
1015                 p3 += 1 ;
1016                 X = ( X / 3 ) ;
1017               }
1018               *_Trace << "S( " << x << " , " << a << "*" << q << " + " << r
1019                       << " ) = 2**" << p2 << "*3**" << p3 << "*"<< X
1020                       << "*S[" << a << "*" << x << "]( " << q << "-1 ) + S( "
1021                       << x << " , " << r << " )" << endl ;
1022             }
1023           }
1024         }
1025       }
1026     }
1027   }
1028   return 0;
1029
1030   int PT ;
1031   int PP ;
1032   int PPNext ;
1033   for ( PP = 1 ; PP <= 5 ; PP+=4 ) {
1034     for ( PPNext = 1 ; PPNext <= 5 ; PPNext+=4 ) {
1035       for ( PT = 0 ; PT <= 2 ; PT++ ) {
1036         if ( PT == 0 && PP == 5 ) {
1037           *_Trace << endl << "T1a 2**(6n1+2)(3n0) + 2**(6n1+1) + S2(3n1+0) --> 6(3n0+0)+5"
1038                   << endl ;
1039           Compute( "T1a" , T1a , PT , PP , PPNext ) ;
1040         }
1041         if ( PT == 1 && PP == 1 ) {
1042           *_Trace << endl << "T2b 2**(6n1+3)(3n0+1) + S2(3n1+0) --> 6(3n0+1)+1"
1043                   << endl ;
1044           Compute( "T2b" , T2b , PT , PP , PPNext ) ;
1045         }
1046         if ( PT == 0 && PP == 1 ) {
1047           *_Trace << endl << "T2a 2**(6n1+7)(3n0) + S2(3n1+2) --> 6(3n0+0)+1"
1048                   << endl ;
1049           Compute( "T2a" , T2a , PT , PP , PPNext ) ;
1050         }
1051         if ( PT == 2 && PP == 5 ) {
1052           *_Trace << endl << "T1c 2**(6n1+4)(3n0+2) + 2**(6n1+3) + S2(3n1+1) --> 6(3n0+2)+5"
1053                   << endl ;
1054           Compute( "T1c" , T1c , PT , PP , PPNext ) ;
1055         }
1056         if ( PT == 2 && PP == 1 ) {
1057           *_Trace << endl << "T2c 2**(6n1+5)(3n0+2) + S2(3n1+1) --> 6(3n0+2)+1"
1058                   << endl ;
1059           Compute( "T2c" , T2c , PT , PP , PPNext ) ;
1060         }
1061         if ( PT == 1 && PP == 5 ) {
1062           *_Trace << endl << "T1b 2**(6n1+6)(3n0+1) + 2**(6n1+5) + S2(3n1+2) --> 6(3n0+1)+5"
1063                   << endl ;
1064           Compute( "T1b" , T1b , PT , PP , PPNext ) ;
1065         }
1066
1067 //      for ( PPNext = 1 ; PPNext <= 5 ; PPNext+=4 ) {
1068         *_Trace << endl << "P2a1 (1<<(6*n1+3))*(3*n0+0) + S2(3*n1+0) --> 6(3n0+0)+1"
1069                 << endl ;
1070         Compute( "P2a1" , P2a1 , PT , PP , PPNext ) ;
1071
1072         *_Trace << endl << "P2a2 (1<<(6*n1+5))*(3*n0+0) + S2(3*n1+1) --> 6(3n0+0)+1"
1073                 << endl ;
1074         Compute( "P2a2" , P2a2 , PT , PP , PPNext ) ;
1075
1076         *_Trace << endl << "P1b1 (1<<(6*n1+2))*(3*n0+1) + (1<<(6*n1+1)) + S2(3*n1+0) --> 6(3n0+1)+5"
1077                 << endl ;
1078         Compute( "P1b1" , P1b1 , PT , PP , PPNext ) ;
1079
1080         *_Trace << endl << "P1c1 (1<<(6*n1+2))*(3*n0+2) + (1<<(6*n1+1)) + S2(3*n1+0) --> 6(3n0+2)+5"
1081                 << endl ;
1082         Compute( "P1c1" , P1c1 , PT , PP , PPNext ) ;
1083
1084         *_Trace << endl << "P1a1 (1<<(6*n1+4))*(3*n0+0) + (1<<(6*n1+3)) + S2(3*n1+1) --> 6(3n0+0)+5"
1085                 << endl ;
1086         Compute( "P1a1" , P1a1 , PT , PP , PPNext ) ;
1087
1088         *_Trace << endl << "P2c1 (1<<(6*n1+3))*(3*n0+2) + S2(3*n1+0) --> 6(3n0+2)+1"
1089                 << endl ;
1090         Compute( "P2c1" , P2c1 , PT , PP , PPNext ) ;
1091
1092         *_Trace << endl << "P1b2 (1<<(6*n1+4))*(3*n0+1) + (1<<(6*n1+3)) + S2(3*n1+1) --> 6(3n0+1)+5"
1093                 << endl ;
1094         Compute( "P1b2" , P1b2 , PT , PP , PPNext ) ;
1095
1096         *_Trace << endl << "P2b1 (1<<(6*n1+5))*(3*n0+1) + S2(3*n1+1) --> 6(3n0+1)+1"
1097                 << endl ;
1098         Compute( "P2b1" , P2b1 , PT , PP , PPNext ) ;
1099
1100         *_Trace << endl << "P1a2 (1<<(6*n1+6))*(3*n0+0) + (1<<(6*n1+5)) + S2(3*n1+2) --> 6(3n0)+5"
1101                 << endl ;
1102         Compute( "P1a2" , P1a2 , PT , PP , PPNext ) ;
1103
1104         *_Trace << endl << "P2b2 (1<<(6*n1+7))*(3*n0+1) + S2(3*n1+2) --> 6(3n0+1)+1"
1105                 << endl ;
1106         Compute( "P2b2" , P2b2 , PT , PP , PPNext ) ;
1107
1108         *_Trace << endl << "P1c2 (1<<(6*n1+6))*(3*n0+2) + (1<<(6*n1+5)) + S2(3*n1+2) --> 6(3n0+2)+5"
1109                 << endl ;
1110         Compute( "P1c2" , P1c2 , PT , PP , PPNext ) ;
1111
1112         *_Trace << endl << "P2c2 (1<<(6*n1+7))*(3*n0+2) + S2(3*n1+2) --> 6(3n0+2)+1"
1113                 << endl ;
1114         Compute( "P2c2" , P2c2 , PT , PP , PPNext ) ;
1115       }
1116
1117     }
1118   }
1119 #endif
1120
1121 #endif
1122
1123   return 0;
1124 }
1125