Mercurial > hg > openjdk > aarch64-port > hotspot
changeset 5571:be0600ec1102
Merge
| author | kvn |
|---|---|
| date | Thu, 27 Jun 2013 11:12:19 -0700 |
| parents | f8972b867ded (current diff) 6a0ead6dc6db (diff) |
| children | 2b9380b0bf0b e0fe0c9a88da 3ea89789ba39 |
| files | |
| diffstat | 11 files changed, 800 insertions(+), 363 deletions(-) [+] |
line wrap: on
line diff
--- a/make/linux/makefiles/gcc.make Thu Jun 27 10:56:10 2013 +0200 +++ b/make/linux/makefiles/gcc.make Thu Jun 27 11:12:19 2013 -0700 @@ -350,9 +350,9 @@ ifeq ($(DEBUG_CFLAGS/$(BUILDARCH)),) ifeq ($(USE_CLANG), true) # Clang doesn't understand -gstabs - OPT_CFLAGS += -g + DEBUG_CFLAGS += -g else - OPT_CFLAGS += -gstabs + DEBUG_CFLAGS += -gstabs endif endif @@ -365,9 +365,9 @@ ifeq ($(FASTDEBUG_CFLAGS/$(BUILDARCH)),) ifeq ($(USE_CLANG), true) # Clang doesn't understand -gstabs - OPT_CFLAGS += -g + FASTDEBUG_CFLAGS += -g else - OPT_CFLAGS += -gstabs + FASTDEBUG_CFLAGS += -gstabs endif endif
--- a/src/cpu/sparc/vm/macroAssembler_sparc.cpp Thu Jun 27 10:56:10 2013 +0200 +++ b/src/cpu/sparc/vm/macroAssembler_sparc.cpp Thu Jun 27 11:12:19 2013 -0700 @@ -1161,12 +1161,6 @@ while (offset() % modulus != 0) nop(); } - -void MacroAssembler::safepoint() { - relocate(breakpoint_Relocation::spec(breakpoint_Relocation::safepoint)); -} - - void RegistersForDebugging::print(outputStream* s) { FlagSetting fs(Debugging, true); int j;
--- a/src/cpu/sparc/vm/relocInfo_sparc.cpp Thu Jun 27 10:56:10 2013 +0200 +++ b/src/cpu/sparc/vm/relocInfo_sparc.cpp Thu Jun 27 11:12:19 2013 -0700 @@ -193,36 +193,6 @@ return *(address*)addr(); } - -int Relocation::pd_breakpoint_size() { - // minimum breakpoint size, in short words - return NativeIllegalInstruction::instruction_size / sizeof(short); -} - -void Relocation::pd_swap_in_breakpoint(address x, short* instrs, int instrlen) { - Untested("pd_swap_in_breakpoint"); - // %%% probably do not need a general instrlen; just use the trap size - if (instrs != NULL) { - assert(instrlen * sizeof(short) == NativeIllegalInstruction::instruction_size, "enough instrlen in reloc. data"); - for (int i = 0; i < instrlen; i++) { - instrs[i] = ((short*)x)[i]; - } - } - NativeIllegalInstruction::insert(x); -} - - -void Relocation::pd_swap_out_breakpoint(address x, short* instrs, int instrlen) { - Untested("pd_swap_out_breakpoint"); - assert(instrlen * sizeof(short) == sizeof(int), "enough buf"); - union { int l; short s[1]; } u; - for (int i = 0; i < instrlen; i++) { - u.s[i] = instrs[i]; - } - NativeInstruction* ni = nativeInstruction_at(x); - ni->set_long_at(0, u.l); -} - void poll_Relocation::fix_relocation_after_move(const CodeBuffer* src, CodeBuffer* dest) { }
--- a/src/cpu/x86/vm/relocInfo_x86.cpp Thu Jun 27 10:56:10 2013 +0200 +++ b/src/cpu/x86/vm/relocInfo_x86.cpp Thu Jun 27 11:12:19 2013 -0700 @@ -177,30 +177,6 @@ return *pd_address_in_code(); } -int Relocation::pd_breakpoint_size() { - // minimum breakpoint size, in short words - return NativeIllegalInstruction::instruction_size / sizeof(short); -} - -void Relocation::pd_swap_in_breakpoint(address x, short* instrs, int instrlen) { - Untested("pd_swap_in_breakpoint"); - if (instrs != NULL) { - assert(instrlen * sizeof(short) == NativeIllegalInstruction::instruction_size, "enough instrlen in reloc. data"); - for (int i = 0; i < instrlen; i++) { - instrs[i] = ((short*)x)[i]; - } - } - NativeIllegalInstruction::insert(x); -} - - -void Relocation::pd_swap_out_breakpoint(address x, short* instrs, int instrlen) { - Untested("pd_swap_out_breakpoint"); - assert(NativeIllegalInstruction::instruction_size == sizeof(short), "right address unit for update"); - NativeInstruction* ni = nativeInstruction_at(x); - *(short*)ni->addr_at(0) = instrs[0]; -} - void poll_Relocation::fix_relocation_after_move(const CodeBuffer* src, CodeBuffer* dest) { #ifdef _LP64 if (!Assembler::is_polling_page_far()) {
--- a/src/cpu/zero/vm/relocInfo_zero.cpp Thu Jun 27 10:56:10 2013 +0200 +++ b/src/cpu/zero/vm/relocInfo_zero.cpp Thu Jun 27 11:12:19 2013 -0700 @@ -52,22 +52,6 @@ return (address *) addr(); } -int Relocation::pd_breakpoint_size() { - ShouldNotCallThis(); -} - -void Relocation::pd_swap_in_breakpoint(address x, - short* instrs, - int instrlen) { - ShouldNotCallThis(); -} - -void Relocation::pd_swap_out_breakpoint(address x, - short* instrs, - int instrlen) { - ShouldNotCallThis(); -} - void poll_Relocation::fix_relocation_after_move(const CodeBuffer* src, CodeBuffer* dst) { ShouldNotCallThis();
--- a/src/share/vm/code/nmethod.cpp Thu Jun 27 10:56:10 2013 +0200 +++ b/src/share/vm/code/nmethod.cpp Thu Jun 27 11:12:19 2013 -0700 @@ -1081,11 +1081,6 @@ metadata_Relocation* reloc = iter.metadata_reloc(); reloc->fix_metadata_relocation(); } - - // There must not be any interfering patches or breakpoints. - assert(!(iter.type() == relocInfo::breakpoint_type - && iter.breakpoint_reloc()->active()), - "no active breakpoint"); } }
--- a/src/share/vm/code/relocInfo.cpp Thu Jun 27 10:56:10 2013 +0200 +++ b/src/share/vm/code/relocInfo.cpp Thu Jun 27 11:12:19 2013 -0700 @@ -338,31 +338,6 @@ _limit = limit; } - -void PatchingRelocIterator:: prepass() { - // turn breakpoints off during patching - _init_state = (*this); // save cursor - while (next()) { - if (type() == relocInfo::breakpoint_type) { - breakpoint_reloc()->set_active(false); - } - } - (RelocIterator&)(*this) = _init_state; // reset cursor for client -} - - -void PatchingRelocIterator:: postpass() { - // turn breakpoints back on after patching - (RelocIterator&)(*this) = _init_state; // reset cursor again - while (next()) { - if (type() == relocInfo::breakpoint_type) { - breakpoint_Relocation* bpt = breakpoint_reloc(); - bpt->set_active(bpt->enabled()); - } - } -} - - // All the strange bit-encodings are in here. // The idea is to encode relocation data which are small integers // very efficiently (a single extra halfword). Larger chunks of @@ -704,51 +679,6 @@ _target = address_from_scaled_offset(offset, base); } - -void breakpoint_Relocation::pack_data_to(CodeSection* dest) { - short* p = (short*) dest->locs_end(); - address point = dest->locs_point(); - - *p++ = _bits; - - assert(_target != NULL, "sanity"); - - if (internal()) normalize_address(_target, dest); - - jint target_bits = - (jint)( internal() ? scaled_offset (_target, point) - : runtime_address_to_index(_target) ); - if (settable()) { - // save space for set_target later - p = add_jint(p, target_bits); - } else { - p = add_var_int(p, target_bits); - } - - for (int i = 0; i < instrlen(); i++) { - // put placeholder words until bytes can be saved - p = add_short(p, (short)0x7777); - } - - dest->set_locs_end((relocInfo*) p); -} - - -void breakpoint_Relocation::unpack_data() { - _bits = live_bits(); - - int targetlen = datalen() - 1 - instrlen(); - jint target_bits = 0; - if (targetlen == 0) target_bits = 0; - else if (targetlen == 1) target_bits = *(data()+1); - else if (targetlen == 2) target_bits = relocInfo::jint_from_data(data()+1); - else { ShouldNotReachHere(); } - - _target = internal() ? address_from_scaled_offset(target_bits, addr()) - : index_to_runtime_address (target_bits); -} - - //// miscellaneous methods oop* oop_Relocation::oop_addr() { int n = _oop_index; @@ -933,81 +863,6 @@ return target; } - -breakpoint_Relocation::breakpoint_Relocation(int kind, address target, bool internal) { - bool active = false; - bool enabled = (kind == initialization); - bool removable = (kind != safepoint); - bool settable = (target == NULL); - - int bits = kind; - if (enabled) bits |= enabled_state; - if (internal) bits |= internal_attr; - if (removable) bits |= removable_attr; - if (settable) bits |= settable_attr; - - _bits = bits | high_bit; - _target = target; - - assert(this->kind() == kind, "kind encoded"); - assert(this->enabled() == enabled, "enabled encoded"); - assert(this->active() == active, "active encoded"); - assert(this->internal() == internal, "internal encoded"); - assert(this->removable() == removable, "removable encoded"); - assert(this->settable() == settable, "settable encoded"); -} - - -address breakpoint_Relocation::target() const { - return _target; -} - - -void breakpoint_Relocation::set_target(address x) { - assert(settable(), "must be settable"); - jint target_bits = - (jint)(internal() ? scaled_offset (x, addr()) - : runtime_address_to_index(x)); - short* p = &live_bits() + 1; - p = add_jint(p, target_bits); - assert(p == instrs(), "new target must fit"); - _target = x; -} - - -void breakpoint_Relocation::set_enabled(bool b) { - if (enabled() == b) return; - - if (b) { - set_bits(bits() | enabled_state); - } else { - set_active(false); // remove the actual breakpoint insn, if any - set_bits(bits() & ~enabled_state); - } -} - - -void breakpoint_Relocation::set_active(bool b) { - assert(!b || enabled(), "cannot activate a disabled breakpoint"); - - if (active() == b) return; - - // %%% should probably seize a lock here (might not be the right lock) - //MutexLockerEx ml_patch(Patching_lock, true); - //if (active() == b) return; // recheck state after locking - - if (b) { - set_bits(bits() | active_state); - if (instrlen() == 0) - fatal("breakpoints in original code must be undoable"); - pd_swap_in_breakpoint (addr(), instrs(), instrlen()); - } else { - set_bits(bits() & ~active_state); - pd_swap_out_breakpoint(addr(), instrs(), instrlen()); - } -} - - //--------------------------------------------------------------------------------- // Non-product code
--- a/src/share/vm/code/relocInfo.hpp Thu Jun 27 10:56:10 2013 +0200 +++ b/src/share/vm/code/relocInfo.hpp Thu Jun 27 11:12:19 2013 -0700 @@ -49,9 +49,6 @@ // RelocIterator // A StackObj which iterates over the relocations associated with // a range of code addresses. Can be used to operate a copy of code. -// PatchingRelocIterator -// Specialized subtype of RelocIterator which removes breakpoints -// temporarily during iteration, then restores them. // BoundRelocation // An _internal_ type shared by packers and unpackers of relocations. // It pastes together a RelocationHolder with some pointers into @@ -204,15 +201,6 @@ // immediate field must not straddle a unit of memory coherence. // //%note reloc_3 // -// relocInfo::breakpoint_type -- a conditional breakpoint in the code -// Value: none -// Instruction types: any whatsoever -// Data: [b [T]t i...] -// The b is a bit-packed word representing the breakpoint's attributes. -// The t is a target address which the breakpoint calls (when it is enabled). -// The i... is a place to store one or two instruction words overwritten -// by a trap, so that the breakpoint may be subsequently removed. -// // relocInfo::static_stub_type -- an extra stub for each static_call_type // Value: none // Instruction types: a virtual call: { set_oop; jump; } @@ -271,8 +259,8 @@ section_word_type = 9, // internal, but a cross-section reference poll_type = 10, // polling instruction for safepoints poll_return_type = 11, // polling instruction for safepoints at return - breakpoint_type = 12, // an initialization barrier or safepoint - metadata_type = 13, // metadata that used to be oops + metadata_type = 12, // metadata that used to be oops + yet_unused_type_1 = 13, // Still unused yet_unused_type_2 = 14, // Still unused data_prefix_tag = 15, // tag for a prefix (carries data arguments) type_mask = 15 // A mask which selects only the above values @@ -312,7 +300,6 @@ visitor(internal_word) \ visitor(poll) \ visitor(poll_return) \ - visitor(breakpoint) \ visitor(section_word) \ @@ -454,7 +441,7 @@ public: enum { // Conservatively large estimate of maximum length (in shorts) - // of any relocation record (probably breakpoints are largest). + // of any relocation record. // Extended format is length prefix, data words, and tag/offset suffix. length_limit = 1 + 1 + (3*BytesPerWord/BytesPerShort) + 1, have_format = format_width > 0 @@ -571,8 +558,6 @@ void initialize(nmethod* nm, address begin, address limit); - friend class PatchingRelocIterator; - // make an uninitialized one, for PatchingRelocIterator: RelocIterator() { initialize_misc(); } public: @@ -779,9 +764,6 @@ void pd_verify_data_value (address x, intptr_t off) { pd_set_data_value(x, off, true); } address pd_call_destination (address orig_addr = NULL); void pd_set_call_destination (address x); - void pd_swap_in_breakpoint (address x, short* instrs, int instrlen); - void pd_swap_out_breakpoint (address x, short* instrs, int instrlen); - static int pd_breakpoint_size (); // this extracts the address of an address in the code stream instead of the reloc data address* pd_address_in_code (); @@ -1302,87 +1284,6 @@ void fix_relocation_after_move(const CodeBuffer* src, CodeBuffer* dest); }; - -class breakpoint_Relocation : public Relocation { - relocInfo::relocType type() { return relocInfo::breakpoint_type; } - - enum { - // attributes which affect the interpretation of the data: - removable_attr = 0x0010, // buffer [i...] allows for undoing the trap - internal_attr = 0x0020, // the target is an internal addr (local stub) - settable_attr = 0x0040, // the target is settable - - // states which can change over time: - enabled_state = 0x0100, // breakpoint must be active in running code - active_state = 0x0200, // breakpoint instruction actually in code - - kind_mask = 0x000F, // mask for extracting kind - high_bit = 0x4000 // extra bit which is always set - }; - - public: - enum { - // kinds: - initialization = 1, - safepoint = 2 - }; - - // If target is NULL, 32 bits are reserved for a later set_target(). - static RelocationHolder spec(int kind, address target = NULL, bool internal_target = false) { - RelocationHolder rh = newHolder(); - new(rh) breakpoint_Relocation(kind, target, internal_target); - return rh; - } - - private: - // We require every bits value to NOT to fit into relocInfo::datalen_width, - // because we are going to actually store state in the reloc, and so - // cannot allow it to be compressed (and hence copied by the iterator). - - short _bits; // bit-encoded kind, attrs, & state - address _target; - - breakpoint_Relocation(int kind, address target, bool internal_target); - - friend class RelocIterator; - breakpoint_Relocation() { } - - short bits() const { return _bits; } - short& live_bits() const { return data()[0]; } - short* instrs() const { return data() + datalen() - instrlen(); } - int instrlen() const { return removable() ? pd_breakpoint_size() : 0; } - - void set_bits(short x) { - assert(live_bits() == _bits, "must be the only mutator of reloc info"); - live_bits() = _bits = x; - } - - public: - address target() const; - void set_target(address x); - - int kind() const { return bits() & kind_mask; } - bool enabled() const { return (bits() & enabled_state) != 0; } - bool active() const { return (bits() & active_state) != 0; } - bool internal() const { return (bits() & internal_attr) != 0; } - bool removable() const { return (bits() & removable_attr) != 0; } - bool settable() const { return (bits() & settable_attr) != 0; } - - void set_enabled(bool b); // to activate, you must also say set_active - void set_active(bool b); // actually inserts bpt (must be enabled 1st) - - // data is packed as 16 bits, followed by the target (1 or 2 words), followed - // if necessary by empty storage for saving away original instruction bytes. - void pack_data_to(CodeSection* dest); - void unpack_data(); - - // during certain operations, breakpoints must be out of the way: - void fix_relocation_after_move(const CodeBuffer* src, CodeBuffer* dest) { - assert(!active(), "cannot perform relocation on enabled breakpoints"); - } -}; - - // We know all the xxx_Relocation classes, so now we can define these: #define EACH_CASE(name) \ inline name##_Relocation* RelocIterator::name##_reloc() { \ @@ -1401,25 +1302,4 @@ initialize(nm, begin, limit); } -// if you are going to patch code, you should use this subclass of -// RelocIterator -class PatchingRelocIterator : public RelocIterator { - private: - RelocIterator _init_state; - - void prepass(); // deactivates all breakpoints - void postpass(); // reactivates all enabled breakpoints - - // do not copy these puppies; it would have unpredictable side effects - // these are private and have no bodies defined because they should not be called - PatchingRelocIterator(const RelocIterator&); - void operator=(const RelocIterator&); - - public: - PatchingRelocIterator(nmethod* nm, address begin = NULL, address limit = NULL) - : RelocIterator(nm, begin, limit) { prepass(); } - - ~PatchingRelocIterator() { postpass(); } -}; - #endif // SHARE_VM_CODE_RELOCINFO_HPP
--- a/src/share/vm/interpreter/bytecodeInterpreter.cpp Thu Jun 27 10:56:10 2013 +0200 +++ b/src/share/vm/interpreter/bytecodeInterpreter.cpp Thu Jun 27 11:12:19 2013 -0700 @@ -481,9 +481,9 @@ // So we have a second version of the assertion which handles the case where EnableInvokeDynamic was // switched off because of the wrong classes. if (EnableInvokeDynamic || FLAG_IS_CMDLINE(EnableInvokeDynamic)) { - assert(abs(istate->_stack_base - istate->_stack_limit) == (istate->_method->max_stack() + 1), "bad stack limit"); + assert(labs(istate->_stack_base - istate->_stack_limit) == (istate->_method->max_stack() + 1), "bad stack limit"); } else { - const int extra_stack_entries = Method::extra_stack_entries_for_indy; + const int extra_stack_entries = Method::extra_stack_entries_for_jsr292; assert(labs(istate->_stack_base - istate->_stack_limit) == (istate->_method->max_stack() + extra_stack_entries + 1), "bad stack limit"); } @@ -2233,7 +2233,7 @@ } Method* method = cache->f1_as_method(); - VERIFY_OOP(method); + if (VerifyOops) method->verify(); if (cache->has_appendix()) { ConstantPool* constants = METHOD->constants(); @@ -2265,8 +2265,7 @@ } Method* method = cache->f1_as_method(); - - VERIFY_OOP(method); + if (VerifyOops) method->verify(); if (cache->has_appendix()) { ConstantPool* constants = METHOD->constants();
--- a/src/share/vm/opto/memnode.cpp Thu Jun 27 10:56:10 2013 +0200 +++ b/src/share/vm/opto/memnode.cpp Thu Jun 27 11:12:19 2013 -0700 @@ -2943,11 +2943,19 @@ Node* my_mem = in(MemBarNode::Precedent); // The MembarAquire may keep an unused LoadNode alive through the Precedent edge if ((my_mem != NULL) && (opc == Op_MemBarAcquire) && (my_mem->outcnt() == 1)) { - assert(my_mem->unique_out() == this, "sanity"); - phase->hash_delete(this); - del_req(Precedent); - phase->is_IterGVN()->_worklist.push(my_mem); // remove dead node later - my_mem = NULL; + // if the Precedent is a decodeN and its input (a Load) is used at more than one place, + // replace this Precedent (decodeN) with the Load instead. + if ((my_mem->Opcode() == Op_DecodeN) && (my_mem->in(1)->outcnt() > 1)) { + Node* load_node = my_mem->in(1); + set_req(MemBarNode::Precedent, load_node); + phase->is_IterGVN()->_worklist.push(my_mem); + my_mem = load_node; + } else { + assert(my_mem->unique_out() == this, "sanity"); + del_req(Precedent); + phase->is_IterGVN()->_worklist.push(my_mem); // remove dead node later + my_mem = NULL; + } } if (my_mem != NULL && my_mem->is_Mem()) { const TypeOopPtr* t_oop = my_mem->in(MemNode::Address)->bottom_type()->isa_oopptr();
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test/compiler/8005956/PolynomialRoot.java Thu Jun 27 11:12:19 2013 -0700 @@ -0,0 +1,776 @@ +//package com.polytechnik.utils; +/* + * (C) Vladislav Malyshkin 2010 + * This file is under GPL version 3. + * + */ + +/** Polynomial root. + * @version $Id: PolynomialRoot.java,v 1.105 2012/08/18 00:00:05 mal Exp $ + * @author Vladislav Malyshkin mal@gromco.com + */ + +/** +* @test +* @bug 8005956 +* @summary C2: assert(!def_outside->member(r)) failed: Use of external LRG overlaps the same LRG defined in this block +* +* @run main PolynomialRoot +*/ + +public class PolynomialRoot { + + +public static int findPolynomialRoots(final int n, + final double [] p, + final double [] re_root, + final double [] im_root) +{ + if(n==4) + { + return root4(p,re_root,im_root); + } + else if(n==3) + { + return root3(p,re_root,im_root); + } + else if(n==2) + { + return root2(p,re_root,im_root); + } + else if(n==1) + { + return root1(p,re_root,im_root); + } + else + { + throw new RuntimeException("n="+n+" is not supported yet"); + } +} + + + +static final double SQRT3=Math.sqrt(3.0),SQRT2=Math.sqrt(2.0); + + +private static final boolean PRINT_DEBUG=false; + +public static int root4(final double [] p,final double [] re_root,final double [] im_root) +{ + if(PRINT_DEBUG) System.err.println("=====================root4:p="+java.util.Arrays.toString(p)); + final double vs=p[4]; + if(PRINT_DEBUG) System.err.println("p[4]="+p[4]); + if(!(Math.abs(vs)>EPS)) + { + re_root[0]=re_root[1]=re_root[2]=re_root[3]= + im_root[0]=im_root[1]=im_root[2]=im_root[3]=Double.NaN; + return -1; + } + +/* zsolve_quartic.c - finds the complex roots of + * x^4 + a x^3 + b x^2 + c x + d = 0 + */ + final double a=p[3]/vs,b=p[2]/vs,c=p[1]/vs,d=p[0]/vs; + if(PRINT_DEBUG) System.err.println("input a="+a+" b="+b+" c="+c+" d="+d); + + + final double r4 = 1.0 / 4.0; + final double q2 = 1.0 / 2.0, q4 = 1.0 / 4.0, q8 = 1.0 / 8.0; + final double q1 = 3.0 / 8.0, q3 = 3.0 / 16.0; + final int mt; + + /* Deal easily with the cases where the quartic is degenerate. The + * ordering of solutions is done explicitly. */ + if (0 == b && 0 == c) + { + if (0 == d) + { + re_root[0]=-a; + im_root[0]=im_root[1]=im_root[2]=im_root[3]=0; + re_root[1]=re_root[2]=re_root[3]=0; + return 4; + } + else if (0 == a) + { + if (d > 0) + { + final double sq4 = Math.sqrt(Math.sqrt(d)); + re_root[0]=sq4*SQRT2/2; + im_root[0]=re_root[0]; + re_root[1]=-re_root[0]; + im_root[1]=re_root[0]; + re_root[2]=-re_root[0]; + im_root[2]=-re_root[0]; + re_root[3]=re_root[0]; + im_root[3]=-re_root[0]; + if(PRINT_DEBUG) System.err.println("Path a=0 d>0"); + } + else + { + final double sq4 = Math.sqrt(Math.sqrt(-d)); + re_root[0]=sq4; + im_root[0]=0; + re_root[1]=0; + im_root[1]=sq4; + re_root[2]=0; + im_root[2]=-sq4; + re_root[3]=-sq4; + im_root[3]=0; + if(PRINT_DEBUG) System.err.println("Path a=0 d<0"); + } + return 4; + } + } + + if (0.0 == c && 0.0 == d) + { + root2(new double []{p[2],p[3],p[4]},re_root,im_root); + re_root[2]=im_root[2]=re_root[3]=im_root[3]=0; + return 4; + } + + if(PRINT_DEBUG) System.err.println("G Path c="+c+" d="+d); + final double [] u=new double[3]; + + if(PRINT_DEBUG) System.err.println("Generic Path"); + /* For non-degenerate solutions, proceed by constructing and + * solving the resolvent cubic */ + final double aa = a * a; + final double pp = b - q1 * aa; + final double qq = c - q2 * a * (b - q4 * aa); + final double rr = d - q4 * a * (c - q4 * a * (b - q3 * aa)); + final double rc = q2 * pp , rc3 = rc / 3; + final double sc = q4 * (q4 * pp * pp - rr); + final double tc = -(q8 * qq * q8 * qq); + if(PRINT_DEBUG) System.err.println("aa="+aa+" pp="+pp+" qq="+qq+" rr="+rr+" rc="+rc+" sc="+sc+" tc="+tc); + final boolean flag_realroots; + + /* This code solves the resolvent cubic in a convenient fashion + * for this implementation of the quartic. If there are three real + * roots, then they are placed directly into u[]. If two are + * complex, then the real root is put into u[0] and the real + * and imaginary part of the complex roots are placed into + * u[1] and u[2], respectively. */ + { + final double qcub = (rc * rc - 3 * sc); + final double rcub = (rc*(2 * rc * rc - 9 * sc) + 27 * tc); + + final double Q = qcub / 9; + final double R = rcub / 54; + + final double Q3 = Q * Q * Q; + final double R2 = R * R; + + final double CR2 = 729 * rcub * rcub; + final double CQ3 = 2916 * qcub * qcub * qcub; + + if(PRINT_DEBUG) System.err.println("CR2="+CR2+" CQ3="+CQ3+" R="+R+" Q="+Q); + + if (0 == R && 0 == Q) + { + flag_realroots=true; + u[0] = -rc3; + u[1] = -rc3; + u[2] = -rc3; + } + else if (CR2 == CQ3) + { + flag_realroots=true; + final double sqrtQ = Math.sqrt (Q); + if (R > 0) + { + u[0] = -2 * sqrtQ - rc3; + u[1] = sqrtQ - rc3; + u[2] = sqrtQ - rc3; + } + else + { + u[0] = -sqrtQ - rc3; + u[1] = -sqrtQ - rc3; + u[2] = 2 * sqrtQ - rc3; + } + } + else if (R2 < Q3) + { + flag_realroots=true; + final double ratio = (R >= 0?1:-1) * Math.sqrt (R2 / Q3); + final double theta = Math.acos (ratio); + final double norm = -2 * Math.sqrt (Q); + + u[0] = norm * Math.cos (theta / 3) - rc3; + u[1] = norm * Math.cos ((theta + 2.0 * Math.PI) / 3) - rc3; + u[2] = norm * Math.cos ((theta - 2.0 * Math.PI) / 3) - rc3; + } + else + { + flag_realroots=false; + final double A = -(R >= 0?1:-1)*Math.pow(Math.abs(R)+Math.sqrt(R2-Q3),1.0/3.0); + final double B = Q / A; + + u[0] = A + B - rc3; + u[1] = -0.5 * (A + B) - rc3; + u[2] = -(SQRT3*0.5) * Math.abs (A - B); + } + if(PRINT_DEBUG) System.err.println("u[0]="+u[0]+" u[1]="+u[1]+" u[2]="+u[2]+" qq="+qq+" disc="+((CR2 - CQ3) / 2125764.0)); + } + /* End of solution to resolvent cubic */ + + /* Combine the square roots of the roots of the cubic + * resolvent appropriately. Also, calculate 'mt' which + * designates the nature of the roots: + * mt=1 : 4 real roots + * mt=2 : 0 real roots + * mt=3 : 2 real roots + */ + + + final double w1_re,w1_im,w2_re,w2_im,w3_re,w3_im,mod_w1w2,mod_w1w2_squared; + if (flag_realroots) + { + mod_w1w2=-1; + mt = 2; + int jmin=0; + double vmin=Math.abs(u[jmin]); + for(int j=1;j<3;j++) + { + final double vx=Math.abs(u[j]); + if(vx<vmin) + { + vmin=vx; + jmin=j; + } + } + final double u1=u[(jmin+1)%3],u2=u[(jmin+2)%3]; + mod_w1w2_squared=Math.abs(u1*u2); + if(u1>=0) + { + w1_re=Math.sqrt(u1); + w1_im=0; + } + else + { + w1_re=0; + w1_im=Math.sqrt(-u1); + } + if(u2>=0) + { + w2_re=Math.sqrt(u2); + w2_im=0; + } + else + { + w2_re=0; + w2_im=Math.sqrt(-u2); + } + if(PRINT_DEBUG) System.err.println("u1="+u1+" u2="+u2+" jmin="+jmin); + } + else + { + mt = 3; + final double w_mod2_sq=u[1]*u[1]+u[2]*u[2],w_mod2=Math.sqrt(w_mod2_sq),w_mod=Math.sqrt(w_mod2); + if(w_mod2_sq<=0) + { + w1_re=w1_im=0; + } + else + { + // calculate square root of a complex number (u[1],u[2]) + // the result is in the (w1_re,w1_im) + final double absu1=Math.abs(u[1]),absu2=Math.abs(u[2]),w; + if(absu1>=absu2) + { + final double t=absu2/absu1; + w=Math.sqrt(absu1*0.5 * (1.0 + Math.sqrt(1.0 + t * t))); + if(PRINT_DEBUG) System.err.println(" Path1 "); + } + else + { + final double t=absu1/absu2; + w=Math.sqrt(absu2*0.5 * (t + Math.sqrt(1.0 + t * t))); + if(PRINT_DEBUG) System.err.println(" Path1a "); + } + if(u[1]>=0) + { + w1_re=w; + w1_im=u[2]/(2*w); + if(PRINT_DEBUG) System.err.println(" Path2 "); + } + else + { + final double vi = (u[2] >= 0) ? w : -w; + w1_re=u[2]/(2*vi); + w1_im=vi; + if(PRINT_DEBUG) System.err.println(" Path2a "); + } + } + final double absu0=Math.abs(u[0]); + if(w_mod2>=absu0) + { + mod_w1w2=w_mod2; + mod_w1w2_squared=w_mod2_sq; + w2_re=w1_re; + w2_im=-w1_im; + } + else + { + mod_w1w2=-1; + mod_w1w2_squared=w_mod2*absu0; + if(u[0]>=0) + { + w2_re=Math.sqrt(absu0); + w2_im=0; + } + else + { + w2_re=0; + w2_im=Math.sqrt(absu0); + } + } + if(PRINT_DEBUG) System.err.println("u[0]="+u[0]+"u[1]="+u[1]+" u[2]="+u[2]+" absu0="+absu0+" w_mod="+w_mod+" w_mod2="+w_mod2); + } + + /* Solve the quadratic in order to obtain the roots + * to the quartic */ + if(mod_w1w2>0) + { + // a shorcut to reduce rounding error + w3_re=qq/(-8)/mod_w1w2; + w3_im=0; + } + else if(mod_w1w2_squared>0) + { + // regular path + final double mqq8n=qq/(-8)/mod_w1w2_squared; + w3_re=mqq8n*(w1_re*w2_re-w1_im*w2_im); + w3_im=-mqq8n*(w1_re*w2_im+w2_re*w1_im); + } + else + { + // typically occur when qq==0 + w3_re=w3_im=0; + } + + final double h = r4 * a; + if(PRINT_DEBUG) System.err.println("w1_re="+w1_re+" w1_im="+w1_im+" w2_re="+w2_re+" w2_im="+w2_im+" w3_re="+w3_re+" w3_im="+w3_im+" h="+h); + + re_root[0]=w1_re+w2_re+w3_re-h; + im_root[0]=w1_im+w2_im+w3_im; + re_root[1]=-(w1_re+w2_re)+w3_re-h; + im_root[1]=-(w1_im+w2_im)+w3_im; + re_root[2]=w2_re-w1_re-w3_re-h; + im_root[2]=w2_im-w1_im-w3_im; + re_root[3]=w1_re-w2_re-w3_re-h; + im_root[3]=w1_im-w2_im-w3_im; + + return 4; +} + + + + static void setRandomP(final double [] p,final int n,java.util.Random r) + { + if(r.nextDouble()<0.1) + { + // integer coefficiens + for(int j=0;j<p.length;j++) + { + if(j<=n) + { + p[j]=(r.nextInt(2)<=0?-1:1)*r.nextInt(10); + } + else + { + p[j]=0; + } + } + } + else + { + // real coefficiens + for(int j=0;j<p.length;j++) + { + if(j<=n) + { + p[j]=-1+2*r.nextDouble(); + } + else + { + p[j]=0; + } + } + } + if(Math.abs(p[n])<1e-2) + { + p[n]=(r.nextInt(2)<=0?-1:1)*(0.1+r.nextDouble()); + } + } + + + static void checkValues(final double [] p, + final int n, + final double rex, + final double imx, + final double eps, + final String txt) + { + double res=0,ims=0,sabs=0; + final double xabs=Math.abs(rex)+Math.abs(imx); + for(int k=n;k>=0;k--) + { + final double res1=(res*rex-ims*imx)+p[k]; + final double ims1=(ims*rex+res*imx); + res=res1; + ims=ims1; + sabs+=xabs*sabs+p[k]; + } + sabs=Math.abs(sabs); + if(false && sabs>1/eps? + (!(Math.abs(res/sabs)<=eps)||!(Math.abs(ims/sabs)<=eps)) + : + (!(Math.abs(res)<=eps)||!(Math.abs(ims)<=eps))) + { + throw new RuntimeException( + getPolinomTXT(p)+"\n"+ + "\t x.r="+rex+" x.i="+imx+"\n"+ + "res/sabs="+(res/sabs)+" ims/sabs="+(ims/sabs)+ + " sabs="+sabs+ + "\nres="+res+" ims="+ims+" n="+n+" eps="+eps+" "+ + " sabs>1/eps="+(sabs>1/eps)+ + " f1="+(!(Math.abs(res/sabs)<=eps)||!(Math.abs(ims/sabs)<=eps))+ + " f2="+(!(Math.abs(res)<=eps)||!(Math.abs(ims)<=eps))+ + " "+txt); + } + } + + static String getPolinomTXT(final double [] p) + { + final StringBuilder buf=new StringBuilder(); + buf.append("order="+(p.length-1)+"\t"); + for(int k=0;k<p.length;k++) + { + buf.append("p["+k+"]="+p[k]+";"); + } + return buf.toString(); + } + + static String getRootsTXT(int nr,final double [] re,final double [] im) + { + final StringBuilder buf=new StringBuilder(); + for(int k=0;k<nr;k++) + { + buf.append("x."+k+"("+re[k]+","+im[k]+")\n"); + } + return buf.toString(); + } + + static void testRoots(final int n, + final int n_tests, + final java.util.Random rn, + final double eps) + { + final double [] p=new double [n+1]; + final double [] rex=new double [n],imx=new double [n]; + for(int i=0;i<n_tests;i++) + { + for(int dg=n;dg-->-1;) + { + for(int dr=3;dr-->0;) + { + setRandomP(p,n,rn); + for(int j=0;j<=dg;j++) + { + p[j]=0; + } + if(dr==0) + { + p[0]=-1+2.0*rn.nextDouble(); + } + else if(dr==1) + { + p[0]=p[1]=0; + } + + findPolynomialRoots(n,p,rex,imx); + + for(int j=0;j<n;j++) + { + //System.err.println("j="+j); + checkValues(p,n,rex[j],imx[j],eps," t="+i); + } + } + } + } + System.err.println("testRoots(): n_tests="+n_tests+" OK, dim="+n); + } + + + + + static final double EPS=0; + + public static int root1(final double [] p,final double [] re_root,final double [] im_root) + { + if(!(Math.abs(p[1])>EPS)) + { + re_root[0]=im_root[0]=Double.NaN; + return -1; + } + re_root[0]=-p[0]/p[1]; + im_root[0]=0; + return 1; + } + + public static int root2(final double [] p,final double [] re_root,final double [] im_root) + { + if(!(Math.abs(p[2])>EPS)) + { + re_root[0]=re_root[1]=im_root[0]=im_root[1]=Double.NaN; + return -1; + } + final double b2=0.5*(p[1]/p[2]),c=p[0]/p[2],d=b2*b2-c; + if(d>=0) + { + final double sq=Math.sqrt(d); + if(b2<0) + { + re_root[1]=-b2+sq; + re_root[0]=c/re_root[1]; + } + else if(b2>0) + { + re_root[0]=-b2-sq; + re_root[1]=c/re_root[0]; + } + else + { + re_root[0]=-b2-sq; + re_root[1]=-b2+sq; + } + im_root[0]=im_root[1]=0; + } + else + { + final double sq=Math.sqrt(-d); + re_root[0]=re_root[1]=-b2; + im_root[0]=sq; + im_root[1]=-sq; + } + return 2; + } + + public static int root3(final double [] p,final double [] re_root,final double [] im_root) + { + final double vs=p[3]; + if(!(Math.abs(vs)>EPS)) + { + re_root[0]=re_root[1]=re_root[2]= + im_root[0]=im_root[1]=im_root[2]=Double.NaN; + return -1; + } + final double a=p[2]/vs,b=p[1]/vs,c=p[0]/vs; + /* zsolve_cubic.c - finds the complex roots of x^3 + a x^2 + b x + c = 0 + */ + final double q = (a * a - 3 * b); + final double r = (a*(2 * a * a - 9 * b) + 27 * c); + + final double Q = q / 9; + final double R = r / 54; + + final double Q3 = Q * Q * Q; + final double R2 = R * R; + + final double CR2 = 729 * r * r; + final double CQ3 = 2916 * q * q * q; + final double a3=a/3; + + if (R == 0 && Q == 0) + { + re_root[0]=re_root[1]=re_root[2]=-a3; + im_root[0]=im_root[1]=im_root[2]=0; + return 3; + } + else if (CR2 == CQ3) + { + /* this test is actually R2 == Q3, written in a form suitable + for exact computation with integers */ + + /* Due to finite precision some double roots may be missed, and + will be considered to be a pair of complex roots z = x +/- + epsilon i close to the real axis. */ + + final double sqrtQ = Math.sqrt (Q); + + if (R > 0) + { + re_root[0] = -2 * sqrtQ - a3; + re_root[1]=re_root[2]=sqrtQ - a3; + im_root[0]=im_root[1]=im_root[2]=0; + } + else + { + re_root[0]=re_root[1] = -sqrtQ - a3; + re_root[2]=2 * sqrtQ - a3; + im_root[0]=im_root[1]=im_root[2]=0; + } + return 3; + } + else if (R2 < Q3) + { + final double sgnR = (R >= 0 ? 1 : -1); + final double ratio = sgnR * Math.sqrt (R2 / Q3); + final double theta = Math.acos (ratio); + final double norm = -2 * Math.sqrt (Q); + final double r0 = norm * Math.cos (theta/3) - a3; + final double r1 = norm * Math.cos ((theta + 2.0 * Math.PI) / 3) - a3; + final double r2 = norm * Math.cos ((theta - 2.0 * Math.PI) / 3) - a3; + + re_root[0]=r0; + re_root[1]=r1; + re_root[2]=r2; + im_root[0]=im_root[1]=im_root[2]=0; + return 3; + } + else + { + final double sgnR = (R >= 0 ? 1 : -1); + final double A = -sgnR * Math.pow (Math.abs (R) + Math.sqrt (R2 - Q3), 1.0 / 3.0); + final double B = Q / A; + + re_root[0]=A + B - a3; + im_root[0]=0; + re_root[1]=-0.5 * (A + B) - a3; + im_root[1]=-(SQRT3*0.5) * Math.abs(A - B); + re_root[2]=re_root[1]; + im_root[2]=-im_root[1]; + return 3; + } + + } + + + static void root3a(final double [] p,final double [] re_root,final double [] im_root) + { + if(Math.abs(p[3])>EPS) + { + final double v=p[3], + a=p[2]/v,b=p[1]/v,c=p[0]/v, + a3=a/3,a3a=a3*a, + pd3=(b-a3a)/3, + qd2=a3*(a3a/3-0.5*b)+0.5*c, + Q=pd3*pd3*pd3+qd2*qd2; + if(Q<0) + { + // three real roots + final double SQ=Math.sqrt(-Q); + final double th=Math.atan2(SQ,-qd2); + im_root[0]=im_root[1]=im_root[2]=0; + final double f=2*Math.sqrt(-pd3); + re_root[0]=f*Math.cos(th/3)-a3; + re_root[1]=f*Math.cos((th+2*Math.PI)/3)-a3; + re_root[2]=f*Math.cos((th+4*Math.PI)/3)-a3; + //System.err.println("3r"); + } + else + { + // one real & two complex roots + final double SQ=Math.sqrt(Q); + final double r1=-qd2+SQ,r2=-qd2-SQ; + final double v1=Math.signum(r1)*Math.pow(Math.abs(r1),1.0/3), + v2=Math.signum(r2)*Math.pow(Math.abs(r2),1.0/3), + sv=v1+v2; + // real root + re_root[0]=sv-a3; + im_root[0]=0; + // complex roots + re_root[1]=re_root[2]=-0.5*sv-a3; + im_root[1]=(v1-v2)*(SQRT3*0.5); + im_root[2]=-im_root[1]; + //System.err.println("1r2c"); + } + } + else + { + re_root[0]=re_root[1]=re_root[2]=im_root[0]=im_root[1]=im_root[2]=Double.NaN; + } + } + + + static void printSpecialValues() + { + for(int st=0;st<6;st++) + { + //final double [] p=new double []{8,1,3,3.6,1}; + final double [] re_root=new double [4],im_root=new double [4]; + final double [] p; + final int n; + if(st<=3) + { + if(st<=0) + { + p=new double []{2,-4,6,-4,1}; + //p=new double []{-6,6,-6,8,-2}; + } + else if(st==1) + { + p=new double []{0,-4,8,3,-9}; + } + else if(st==2) + { + p=new double []{-1,0,2,0,-1}; + } + else + { + p=new double []{-5,2,8,-2,-3}; + } + root4(p,re_root,im_root); + n=4; + } + else + { + p=new double []{0,2,0,1}; + if(st==4) + { + p[1]=-p[1]; + } + root3(p,re_root,im_root); + n=3; + } + System.err.println("======== n="+n); + for(int i=0;i<=n;i++) + { + if(i<n) + { + System.err.println(String.valueOf(i)+"\t"+ + p[i]+"\t"+ + re_root[i]+"\t"+ + im_root[i]); + } + else + { + System.err.println(String.valueOf(i)+"\t"+p[i]+"\t"); + } + } + } + } + + + + public static void main(final String [] args) + { + final long t0=System.currentTimeMillis(); + final double eps=1e-6; + //checkRoots(); + final java.util.Random r=new java.util.Random(-1381923); + printSpecialValues(); + + final int n_tests=10000000; + //testRoots(2,n_tests,r,eps); + //testRoots(3,n_tests,r,eps); + testRoots(4,n_tests,r,eps); + final long t1=System.currentTimeMillis(); + System.err.println("PolynomialRoot.main: "+n_tests+" tests OK done in "+(t1-t0)+" milliseconds. ver=$Id: PolynomialRoot.java,v 1.105 2012/08/18 00:00:05 mal Exp $"); + } + + + +}