src/SyrComponent/SyrComponent_CheckOfUndefined.cxx

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